text:
## real double shift code

##################################################
## Bulge things

## Absorb Ut
## We have the bulge is generated by
## Ut * A * U in the singleshift case and
## Ut * Vt * A * V * U in the DoubleShift Case
## This combines the Ut with A (or Ut, Vt, A into W,A)

function absorb_Ut(QF::QFactorization{T,S, VV}, RF, storage, ctr) where {T,S <: Real, VV}

    U, V = storage.VU[2], storage.VU[1]
    Ut, Vt = U', V'

    i = idx(U)
    j  = i + 1

    W = QF.Q[i]

    p = i == 1 ? one(S) : get_parity(QF,i-1)

    W = dflip(W, p)
    W, Ut, Vt = turnover(Ut, Vt, W)
    fuse!(Vt, QF)
    Ut  = dflip(Ut, p)

    QF.Q[i] = Ut
    storage.limb[1] = W

    return nothing
end

## For double shift we have V then U
function passthrough_triu(QF::QFactorization{T, S, Vt}, RF, storage, ctr,  dir::Val{:right}) where {T, S <: Real, Vt}

    U, V = storage.VU[2], storage.VU[1]
    j = idx(V)

    flag = false

    if  j <= ctr.tr
        # this leverages Ct*B*U = I*U when Ct, B not touched
        # so V, U are not touched, but Ct, B are
        # flag is true if done, false if not
        flag = simple_passthrough!(RF, U, V)
    end

    if j > ctr.tr || !flag

        V = passthrough!(RF, V)
        U = passthrough!(RF, U)
        storage.VU[1], storage.VU[2]  = V, U

    end

    ctr.tr -= 2

    return nothing

end

## When Ct and B are identical, we can update just one and leave U,V alone
function simple_passthrough!(RF::RFactorizationRankOne{T, S}, U, V) where {T, S <: Real}

    j = idx(V)
    N = length(RF)

    _ = passthrough!(RF.B, V)
    _ = passthrough!(RF.B, U)

    for k in -1:1
        a,b = vals(RF.B[j+k])
        jj = N+1-(j+k)
        RF.Ct[jj] = Rotator(a, -b, j+k)  ## adjoint
    end

    true
end

## passthrough Q
##
## we have Q U or (W,Q) V U to pass through
## This depends on Twisted but not pencil
##
## If we update indices and use a unitary transform, return false (not absorbed)
## else return true (was absorved)
function passthrough_Q(QF::QFactorization{T, S, Vt}, RF, storage, ctr, dir::Val{:right}) where {T, S <: Real, Vt}

    U, V, W = storage.VU[2], storage.VU[1],  storage.limb[1]
    i = idx(U); j = i + 1

    if j < ctr.stop_index
        V = passthrough!(QF, V)
        U = passthrough!(QF, U)
        V, U, W = turnover(W, V, U)
        storage.VU[2], storage.VU[1],  storage.limb[1] = U, V, W

        return false

    else

        p = get_parity(QF,j+1) # 1 or -1 possibly ## XXX
        V = dflip(V, p)
        fuse!(QF, V)

        # now turnover U, merge with W, unitary over, pass through, and fuse
        U = passthrough!(QF.Q, U)
        i += 1 # after turnover, U moves down
        U, Di = fuse(W, U)

        # pass U through triangle then fuse
        U = passthrough!(RF, U)

        p = get_parity(QF,i+1)
        U = dflip(U, p)
        fuse!(QF, U)

        return true

    end
end

##################################################
function deflate(QF::AbstractQFactorization{T, S}, k, ctr) where {T,S <: Real}

    c,s = vals(QF.Q[k])
    i = idx(QF.Q[k])
    QF.Q[k] = Rotator(sign(c), zero(T), i) # ± 1, not just 1

end

## The zero_index and stop_index+1 point at diagonal rotators [1 0; 0 1] or [-1 0; 0 -1]
## this recoves 1 or -1
function get_parity(QF::QFactorization{T, S}, k) where {T, S <: Real}
    if k > length(QF)
        return one(S)
    else
        c::S, s::T = vals(QF.Q[k])
        ## @assert iszero(s)
        sign(c)
    end
end

##################################################
## Transformations

## Let a D matrix be one of [1 0; 0 1] or [-1 0; 0 1] (D^2 = I). Then we have this move
## D    --->   D  (we update the rotator)
##   [       [
##
## this is `dflip`
## @vars c s
## D=[-1 0 0; 0 -1 0; 0 0 1]
## U = [1 0 0; 0 c s; 0 -s c]
## U1 = [1 0 0; 0 c -s; 0 s c]
## D*U - U1*D
function dflip(a::Rotator{T,S}, d=one(T)) where {T, S<:Real}
    c::S, s::T = vals(a)
    p = sign(d)
    i = idx(a)
    return Rotator(c, p*s, i)
end


text:
"""
    emojify(path::String)

Go recursively over all the files contained in path and replace
all possible occurence of functions with random emoji aliases
"""
function emojify(path::String; overwrite=👌)
    if isdir(path)
        for subpath in readdir(path)
            emojify(joinpath(path, subpath); overwrite=overwrite)
        end
    elseif isfile(path) && endswith(path, ".jl")
        return emojify_file(path; overwrite=overwrite)
    end
    return ⬛
end

function emojify_file(filepath::String; overwrite=👍)
    str = String(read(filepath))
    str = emojify_string(str)
    if overwrite
        write(filepath, str)
        return ⬛
    else
        return str
    end
end

function emojify_string(str::String)
    for func in 🔑(func_to_emojis)
        str = replace(str, Regex("\\b" * string(func) * "\\b") => RandString(to_string.(func_to_emojis[func])))
    end
    return str
end

function to_string(emoji::Symbol)
    string(emoji)
end

function to_string(emoji::Tuple)
    return string(emoji[1])
end

## Allow to get a random string every time it's printed
struct RandString{TS}
    strings::TS
end

function 🖨(io::IO, rs::RandString)
    🖨(io, 🎲(rs.strings))
end

ext: jl


text:
@testset "c\"#include ...\"" begin
	@eval c`-I$(@__DIR__)`
	ctx = C.CONTEXT_CACHE[@__MODULE__]

	@test length(ctx.hdrs) == 0
	@eval c"""
	#include <stddef.h>
	"""
	@test length(ctx.hdrs) == 1
	@eval c"""
	#include "test.h"
	"""i
	@test length(ctx.hdrs) == 1
	@eval c"""
	#include <inttypes.h>
	"""s
	@test length(ctx.hdrs) == 0

	str = String(take!(ctx.src))
	@test occursin("<stddef.h>", str)
	@test occursin("\"test.h\"", str)
	@test occursin("<inttypes.h>", str)
end


text:
include("./imports.jl")

@testset ExtendedTestSet "Normalize" begin
    @testset ExtendedTestSet "normalize" begin
        image = rand(Float32, 10, 10, 3)
        means = SVector{3}(rand(Float32, 1, 1, 3))
        stds = SVector{3}(rand(Float32, 1, 1, 3))
        @test denormalize(normalize(copy(image), means, stds), means, stds) ≈ image
    end

    @testset ExtendedTestSet "Normalize" begin
        item = ArrayItem(zeros(10, 10, 3))
        means =  SVector{3}(ones(1, 1, 3))
        stds =  SVector{3}(ones(1, 1, 3))
        tfm = Normalize(means, stds)
        @test itemdata(apply(tfm, item)) ≈ -1 .* ones(10, 10, 3)
        testapply(tfm, item)
        item1 = ArrayItem(rand(10, 10, 3))
        item2 = ArrayItem(rand(10, 10, 3))
        testapply!(tfm, item1, item2)
    end
end

@testset ExtendedTestSet "NormalizeIntensity" begin
    item = ArrayItem(Float32.(collect(1:5)))
    ground_truth = [-1.2649, -0.6325, 0, 0.6325, 1.2649]
    tfm = NormalizeIntensity()
    @test itemdata(apply(tfm, item)) ≈ ground_truth
end

@testset ExtendedTestSet "ToEltype" begin
    data = rand(RGB, 10, 10)
    imagergb = Image(data)
    imagegray = Image((Gray.(data)))
    @test itemdata(apply(ToEltype(Gray), imagegray)) === itemdata(imagegray)
    @test itemdata(apply(ToEltype(Gray), imagergb)) ≈ itemdata(imagegray)

    tfm = ToEltype(Gray)
    testapply(tfm, Image)
    testapply!(tfm, Image)
end

@testset ExtendedTestSet "ImageToTensor" begin
    @testset ExtendedTestSet "imagetotensor,tensortoimage" begin
        data = rand(RGB, 10, 10)
        @test size(imagetotensor(data)) == (10, 10, 3)
        @test data |> imagetotensor |> tensortoimage ≈ data

    end

    @testset ExtendedTestSet "ImageToTensor" begin
        image = Image(rand(RGB, 10, 10))
        tfm = ImageToTensor()
        @test_nowarn apply(tfm, image)
        a = itemdata(apply(tfm, image))
        @test size(a) == (10, 10, 3)
        testapply(tfm, Image)
        testapply!(tfm, Image)
    end
end

@testset ExtendedTestSet "OneHot" begin
    tfm = OneHot()
    mask = rand(1:4, 10, 10)
    item1 = MaskMulti(mask, 1:4)
    item2 = MaskMulti(rand(1:3, 10, 10), 1:4)
    @test_nowarn apply(tfm, item1)
    aitem = apply(tfm, item1)
    @test size(itemdata(aitem)) == (10, 10, 4)

    testapply(tfm, item1)
    testapply!(tfm, item1, item2)

end

@testset ExtendedTestSet "Image pipeline" begin
    item1 = Image(rand(RGB, 64, 64))
    item2 = Image(rand(RGB, 64, 64))

    tfm = compose(
        ToEltype(RGB),
        RandomResizeCrop((32, 48)),
        ImageToTensor(),
        Normalize([0, 0, 0], [1, 1, 1])
    )

    res = apply(tfm, item1)
    a = itemdata(res)
    @test size(a) == (32, 48, 3)
    @test eltype(a) == Float32

    testapply(tfm, item1)
    testapply!(tfm, item1, item2)
end


text:
using Libdl

cd(dirname(@__FILE__))

if (!ispath("lib"))
  run(`mkdir lib`)
end

cd("src")
run(`make`)
run(`mv libgiac_c.$(dlext) ../lib`)


text:
@testset "cplx2real - Exponent" begin
    z1 = MPC.Variable("z1", Complex)
    z2 = MPC.Variable("z2", Complex)
    x = MPC.Variable("x", Real)
    b = MPC.Variable("b", Bool)

    expo_cplx = Exponent(SortedDict(z1=>Degree(1,0),
                                    x=>Degree(4,0),
                                    b=>Degree(1,0)))

    z1_re = MPC.Variable("z1_Re", Real)
    z1_im = MPC.Variable("z1_Im", Real)

    expo_real = Exponent(SortedDict(z1_re=>Degree(1,0),
                                    x=>Degree(4,0),
                                    b=>Degree(1,0)))
    expo_imag = Exponent(SortedDict(z1_im=>Degree(1,0),
                                    x=>Degree(4,0),
                                    b=>Degree(1,0)))

    preal, pimag = cplx2real(expo_cplx, 1.)
    @test preal == Polynomial(SortedDict{Exponent, Number}(expo_real=>1.0))
    @test pimag == Polynomial(SortedDict{Exponent, Number}(expo_imag=>1.0))

    expo_cplx = Exponent(SortedDict(z1=>Degree(0,1),
                                    x=>Degree(4,0),
                                    b=>Degree(1,0)))

    preal, pimag = cplx2real(expo_cplx, 1.)
    @test preal == Polynomial(SortedDict{Exponent, Number}(expo_real=>1.0))
    @test pimag == Polynomial(SortedDict{Exponent, Number}(expo_imag=>-1.0))
end

@testset "cplx2real - Polynomial" begin
    z1 = MPC.Variable("z1", Complex)
    z2 = MPC.Variable("z2", Complex)
    x = MPC.Variable("x", Real)
    b = MPC.Variable("b", Bool)

    expo_cplx = Exponent(SortedDict(z1=>Degree(1,0),
                                    x=>Degree(4,0),
                                    b=>Degree(1,0)))
    p = Polynomial(SortedDict{Exponent, Number}(expo_cplx=>3+5im,
                                                Exponent()=>5im))

    # Corresponding real quantities
    z1_re = MPC.Variable("z1_Re", Real)
    z1_im = MPC.Variable("z1_Im", Real)
    expo_real = Exponent(SortedDict(z1_re=>Degree(1,0),
                                    x=>Degree(4,0),
                                    b=>Degree(1,0)))
    expo_imag = Exponent(SortedDict(z1_im=>Degree(1,0),
                                    x=>Degree(4,0),
                                    b=>Degree(1,0)))

    preal, pimag = cplx2real(p)
    @test preal == Polynomial(SortedDict{Exponent,Number}(expo_real=>3.0,
                                                          expo_imag=>-5.0))
    @test pimag == Polynomial(SortedDict{Exponent,Number}(Exponent()=>5,
                                                          expo_imag=>3.0,
                                                          expo_real=>5.0))

end

@testset "cplx2real - Point" begin
    z1 = MPC.Variable("z1", Complex)
    z2 = MPC.Variable("z2", Complex)
    x = MPC.Variable("x", Real)
    b = MPC.Variable("b", Bool)

    pt_c = Point(SortedDict(z1=>1+2im,
                            z2=>3im,
                            x=>3.5,
                            b=>1))

    z1_Re = MPC.Variable("z1_Re", Real)
    z1_Im = MPC.Variable("z1_Im", Real)
    z2_Im = MPC.Variable("z2_Im", Real)

    pt_re = Point(SortedDict(z1_Re=>1,
                             z1_Im=>2,
                             z2_Im=>3,
                             x=>3.5,
                             b=>1))

    @test cplx2real(pt_c) == pt_re

    @test real2cplx(pt_re) == pt_c
end


text:
include("../utils.jl")
include("group_utils.jl")

using ManifoldsBase: VeeOrthogonalBasis

Random.seed!(10)

@testset "Special Euclidean group" begin
    @testset "SpecialEuclidean($n)" for n in (2, 3, 4)
        G = SpecialEuclidean(n)
        @test isa(G, SpecialEuclidean{n})
        @test repr(G) == "SpecialEuclidean($n)"
        M = base_manifold(G)
        @test M === TranslationGroup(n) × SpecialOrthogonal(n)
        @test submanifold(G, 1) === TranslationGroup(n)
        @test submanifold(G, 2) === SpecialOrthogonal(n)
        Rn = Rotations(n)
        p = Matrix(I, n, n)

        if n == 2
            t = Vector{Float64}.([1:2, 2:3, 3:4])
            ω = [[1.0], [2.0], [1.0]]
            tuple_pts = [(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)]
            tuple_X = [([-1.0, 2.0], hat(Rn, p, [1.0])), ([1.0, -2.0], hat(Rn, p, [0.5]))]
        elseif n == 3
            t = Vector{Float64}.([1:3, 2:4, 4:6])
            ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
            tuple_pts = [(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)]
            tuple_X = [
                ([-1.0, 2.0, 1.0], hat(Rn, p, [1.0, 0.5, -0.5])),
                ([-2.0, 1.0, 0.5], hat(Rn, p, [-1.0, -0.5, 1.1])),
            ]
        else # n == 4
            t = Vector{Float64}.([1:4, 2:5, 3:6])
            ω = [
                [1.0, 2.0, 3.0, 1.0, 2.0, 3.0],
                [3.0, 2.0, 1.0, 1.0, 2.0, 3.0],
                [1.0, 3.0, 2.0, 1.0, 2.0, 3.0],
            ]
            tuple_pts = [(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)]
            tuple_X = [
                ([-1.0, 2.0, 1.0, 3.0], hat(Rn, p, [1.0, 0.5, -0.5, 0.0, 2.0, 1.0])),
                ([-2.0, 1.5, -1.0, 2.0], hat(Rn, p, [1.0, -0.5, 0.5, 1.0, 0.0, 1.0])),
            ]
        end

        @testset "product repr" begin
            pts = [ProductRepr(tp...) for tp in tuple_pts]
            X_pts = [ProductRepr(tX...) for tX in tuple_X]

            g1, g2 = pts[1:2]
            t1, R1 = g1.parts
            t2, R2 = g2.parts
            g1g2 = ProductRepr(R1 * t2 + t1, R1 * R2)
            @test isapprox(G, compose(G, g1, g2), g1g2)
            g1g2mat = affine_matrix(G, g1g2)
            @test g1g2mat ≈ affine_matrix(G, g1) * affine_matrix(G, g2)
            @test affine_matrix(G, g1g2mat) === g1g2mat
            @test affine_matrix(G, Identity(G)) isa SDiagonal{n,Float64}
            @test affine_matrix(G, Identity(G)) == SDiagonal{n,Float64}(I)

            w = translate_diff(G, pts[1], Identity(G), X_pts[1])
            w2 = allocate(w)
            w2.parts[1] .= w.parts[1]
            w2.parts[2] .= pts[1].parts[2] * w.parts[2]
            w2mat = screw_matrix(G, w2)
            @test w2mat ≈ affine_matrix(G, pts[1]) * screw_matrix(G, X_pts[1])
            @test screw_matrix(G, w2mat) === w2mat

            test_group(
                G,
                pts,
                X_pts,
                X_pts;
                test_diff=true,
                test_lie_bracket=true,
                test_adjoint_action=true,
                diff_convs=[(), (LeftAction(),), (RightAction(),)],
            )
        end

        @testset "affine matrix" begin
            pts = [affine_matrix(G, ProductRepr(tp...)) for tp in tuple_pts]
            X_pts = [screw_matrix(G, ProductRepr(tX...)) for tX in tuple_X]
            test_group(
                G,
                pts,
                X_pts,
                X_pts;
                test_diff=true,
                test_lie_bracket=true,
                diff_convs=[(), (LeftAction(),), (RightAction(),)],
                atol=1e-9,
            )
            # specific affine tests
            p = copy(G, pts[1])
            X = copy(G, p, X_pts[1])
            X[n + 1, n + 1] = 0.1
            @test_throws DomainError is_vector(G, p, X, true)
            X2 = zeros(n + 2, n + 2)
            # nearly correct just too large (and the error from before)
            X2[1:n, 1:n] .= X[1:n, 1:n]
            X2[1:n, end] .= X[1:n, end]
            X2[end, end] = X[end, end]
            @test_throws DomainError is_vector(G, p, X2, true)
            p[n + 1, n + 1] = 0.1
            @test_throws DomainError is_point(G, p, true)
            p2 = zeros(n + 2, n + 2)
            # nearly correct just too large (and the error from before)
            p2[1:n, 1:n] .= p[1:n, 1:n]
            p2[1:n, end] .= p[1:n, end]
            p2[end, end] = p[end, end]
            @test_throws DomainError is_point(G, p2, true)
            # exp/log_lie for ProductGroup on arrays
            X = copy(G, p, X_pts[1])
            p3 = exp_lie(G, X)
            X3 = log_lie(G, p3)
            isapprox(G, Identity(G), X, X3)
        end

        @testset "hat/vee" begin
            p = ProductRepr(tuple_pts[1]...)
            X = ProductRepr(tuple_X[1]...)
            Xexp = [X.parts[1]; vee(Rn, p.parts[2], X.parts[2])]
            Xc = vee(G, p, X)
            @test Xc ≈ Xexp
            @test isapprox(G, p, hat(G, p, Xc), X)

            Xc = vee(G, affine_matrix(G, p), screw_matrix(G, X))
            @test Xc ≈ Xexp
            @test hat(G, affine_matrix(G, p), Xc) ≈ screw_matrix(G, X)

            e = Identity(G)
            Xe = log_lie(G, p)
            Xc = vee(G, e, Xe)
            @test_throws MethodError vee(M, e, Xe)
            w = similar(Xc)
            vee!(G, w, e, Xe)
            @test isapprox(Xc, w)
            @test_throws MethodError vee!(M, w, e, Xe)

            w = similar(Xc)
            vee!(G, w, identity_element(G), Xe)
            @test isapprox(Xc, w)

            Ye = hat(G, e, Xc)
            @test_throws MethodError hat(M, e, Xc)
            isapprox(G, e, Xe, Ye)
            Ye2 = copy(G, p, X)
            hat!(G, Ye2, e, Xc)
            @test_throws MethodError hat!(M, Ye, e, Xc)
            @test isapprox(G, e, Ye, Ye2)

            Ye2 = copy(G, p, X)
            hat!(G, Ye2, identity_element(G), Xc)
            @test isapprox(G, e, Ye, Ye2)
        end
    end

    G = SpecialEuclidean(11)
    @test affine_matrix(G, Identity(G)) isa Diagonal{Float64,Vector{Float64}}
    @test affine_matrix(G, Identity(G)) == Diagonal(ones(11))

    @testset "Explicit embedding in GL(n+1)" begin
        G = SpecialEuclidean(3)
        t = Vector{Float64}.([1:3, 2:4, 4:6])
        ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
        p = Matrix(I, 3, 3)
        Rn = Rotations(3)
        pts = [ProductRepr(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)]
        X = ProductRepr([-1.0, 2.0, 1.0], hat(Rn, p, [1.0, 0.5, -0.5]))
        q = ProductRepr([0.0, 0.0, 0.0], p)

        GL = GeneralLinear(4)
        SEGL = EmbeddedManifold(G, GL)
        @test Manifolds.SpecialEuclideanInGeneralLinear(3) === SEGL
        pts_gl = [embed(SEGL, pp) for pp in pts]
        q_gl = embed(SEGL, q)
        X_gl = embed(SEGL, pts_gl[1], X)

        q_gl2 = allocate(q_gl)
        embed!(SEGL, q_gl2, q)
        @test isapprox(SEGL, q_gl2, q_gl)

        q2 = allocate(q)
        project!(SEGL, q2, q_gl)
        @test isapprox(G, q, q2)

        @test isapprox(G, pts[1], project(SEGL, pts_gl[1]))
        @test isapprox(G, pts[1], X, project(SEGL, pts_gl[1], X_gl))

        X_gl2 = allocate(X_gl)
        embed!(SEGL, X_gl2, pts_gl[1], X)
        @test isapprox(SEGL, pts_gl[1], X_gl2, X_gl)

        X2 = allocate(X)
        project!(SEGL, X2, pts_gl[1], X_gl)
        @test isapprox(G, pts[1], X, X2)

        for conv in [LeftAction(), RightAction()]
            tpgl = translate(GL, pts_gl[2], pts_gl[1], conv)
            tXgl = translate_diff(GL, pts_gl[2], pts_gl[1], X_gl, conv)
            tpse = translate(G, pts[2], pts[1], conv)
            tXse = translate_diff(G, pts[2], pts[1], X, conv)
            @test isapprox(G, tpse, project(SEGL, tpgl))
            @test isapprox(G, tpse, tXse, project(SEGL, tpgl, tXgl))

            @test isapprox(
                G,
                pts_gl[1],
                X_gl,
                translate_diff(G, Identity(G), pts_gl[1], X_gl, conv),
            )
        end
    end

    @testset "Adjoint action on 𝔰𝔢(3)" begin
        G = SpecialEuclidean(3)
        t = Vector{Float64}.([1:3, 2:4, 4:6])
        ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
        p = Matrix(I, 3, 3)
        Rn = Rotations(3)
        pts = [ProductRepr(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)]
        X = ProductRepr([-1.0, 2.0, 1.0], hat(Rn, p, [1.0, 0.5, -0.5]))
        q = ProductRepr([0.0, 0.0, 0.0], p)

        # adjoint action of SE(3)
        fX = TFVector(vee(G, q, X), VeeOrthogonalBasis())
        fXp = adjoint_action(G, pts[1], fX)
        fXp2 = adjoint_action(G, pts[1], X)
        @test isapprox(G, pts[1], hat(G, pts[1], fXp.data), fXp2)
    end
end

max_stars_repo_name: JuliaNLSolvers/ManifoldMuseum

max_issues_repo_name: JuliaNLSolvers/ManifoldMuseum

max_forks_repo_name: JuliaNLSolvers/ManifoldMuseum


text:
# Low-level Interfaces
# ====================

function iserror(code::Csize_t)
    return LibZstd.ZSTD_isError(code) != 0
end

function zstderror(stream, code::Csize_t)
    ptr = LibZstd.ZSTD_getErrorName(code)
    error("zstd error: ", unsafe_string(ptr))
end

function max_clevel()
    return LibZstd.ZSTD_maxCLevel()
end

const MAX_CLEVEL = max_clevel()

const InBuffer = LibZstd.ZSTD_inBuffer
InBuffer() = InBuffer(C_NULL, 0, 0)
const OutBuffer = LibZstd.ZSTD_outBuffer
OutBuffer() = OutBuffer(C_NULL, 0, 0)

# ZSTD_CStream
mutable struct CStream
    ptr::Ptr{Cvoid}
    ibuffer::InBuffer
    obuffer::OutBuffer

    function CStream()
        ptr = LibZstd.ZSTD_createCStream()
        if ptr == C_NULL
            throw(OutOfMemoryError())
        end
        return new(ptr, InBuffer(), OutBuffer())
    end
end

function initialize!(cstream::CStream, level::Integer)
    return LibZstd.ZSTD_initCStream(cstream.ptr, level)
end

function reset!(cstream::CStream, srcsize::Integer)
    # ZSTD_resetCStream is deprecated
    # https://github.com/facebook/zstd/blob/9d2a45a705e22ad4817b41442949cd0f78597154/lib/zstd.h#L2253-L2272
    res = LibZstd.ZSTD_CCtx_reset(cstream.ptr, LibZstd.ZSTD_reset_session_only)
    if iserror(res)
        return res
    end
    if srcsize == 0
        # From zstd.h:
        # Note: ZSTD_resetCStream() interprets pledgedSrcSize == 0 as ZSTD_CONTENTSIZE_UNKNOWN, but
        # ZSTD_CCtx_setPledgedSrcSize() does not do the same, so ZSTD_CONTENTSIZE_UNKNOWN must be
        # explicitly specified.
        srcsize = ZSTD_CONTENTSIZE_UNKNOWN
    end
    return LibZstd.ZSTD_CCtx_setPledgedSrcSize(cstream.ptr, srcsize)
    #return ccall((:ZSTD_resetCStream, libzstd), Csize_t, (Ptr{Cvoid}, Culonglong), cstream.ptr, srcsize)

end

function compress!(cstream::CStream)
    return LibZstd.ZSTD_compressStream(cstream.ptr, pointer_from_objref(cstream.obuffer), pointer_from_objref(cstream.ibuffer))
end

function finish!(cstream::CStream)
    return LibZstd.ZSTD_endStream(cstream.ptr, pointer_from_objref(cstream.obuffer))
end

function free!(cstream::CStream)
    return LibZstd.ZSTD_freeCStream(cstream.ptr)
end

# ZSTD_DStream
mutable struct DStream
    ptr::Ptr{Cvoid}
    ibuffer::InBuffer
    obuffer::OutBuffer

    function DStream()
        ptr = LibZstd.ZSTD_createDStream()
        if ptr == C_NULL
            throw(OutOfMemoryError())
        end
        return new(ptr, InBuffer(), OutBuffer())
    end
end

function initialize!(dstream::DStream)
    return LibZstd.ZSTD_initDStream(dstream.ptr)
end

function reset!(dstream::DStream)
    # LibZstd.ZSTD_resetDStream is deprecated
    # https://github.com/facebook/zstd/blob/9d2a45a705e22ad4817b41442949cd0f78597154/lib/zstd.h#L2332-L2339
    return LibZstd.ZSTD_DCtx_reset(dstream.ptr, LibZstd.ZSTD_reset_session_only)
end

function decompress!(dstream::DStream)
    return LibZstd.ZSTD_decompressStream(dstream.ptr, pointer_from_objref(dstream.obuffer), pointer_from_objref(dstream.ibuffer))
end

function free!(dstream::DStream)
    return LibZstd.ZSTD_freeDStream(dstream.ptr)
end

# Misc. functions
# ---------------

const ZSTD_CONTENTSIZE_UNKNOWN = Culonglong(0) - 1
const ZSTD_CONTENTSIZE_ERROR   = Culonglong(0) - 2

function find_decompressed_size(src::Ptr, size::Integer)
    return LibZstd.ZSTD_findDecompressedSize(src, size)
end


text:
addprocs(3) # Remove addprocs if you want to use a single process only
using IncrementalInference
using KernelDensityEstimate, Gadfly # for vstack
using Cairo, Fontconfig # for drawing PNG/PDF

gt = Dict{ASCIIString, Array{Float64,2}}()
# HMM computed ground truth for first 3 poses only
gt["x1"]=[[-100.0; 1.96]';[0.0; 1.96]']'
gt["x2"]=[[-50.0; 3.1]';[50.0; 3.1]']'
gt["x3"]=[[100.0; 3.05]';[0.0; 3.05]']'

fg = emptyFactorGraph()

N=100

doors = [-100.0;0.0;100.0;300.0]'
cov = [3.0]
# pd = kde!(doors,cov)
# pd = resample(pd,N);
# bws = getBW(pd)[:,1]
# doors2 = getPoints(pd);

v1 = addNode!(fg,"x1",doors,N=N)
f1  = addFactor!(fg,[v1], Obsv2(doors, cov', [1.0]))

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,"x2", tem, N=N)
addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

# monocular sighting would look something like
#addFactor!(fg, Mono, ["x3","l1"], [14.0], [1.0], [1.0])
#addFactor!(fg, Mono, ["x4","l1"], [11.0], [1.0], [1.0])

v3=addNode!(fg,"x3",4.0*randn(1,N)+getVal(v2)+50.0, N=N)
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))
f2 = addFactor!(fg,[v3], Obsv2(doors, cov', [1.0]))

if true

    v4=addNode!(fg,"x4",2.0*randn(1,N)+getVal(v3)+50.0, N=N)
    addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

    if true
        l1=addNode!(fg, "l1", 0.5*randn(1,N)+getVal(v3)+64.0, N=N)
        addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))
        addFactor!(fg, [v4,l1], Ranged([16.0],[0.5],[1.0]))
    end

    v5=addNode!(fg,"x5",2.0*randn(1,N)+getVal(v4)+50.0, N=N)
    addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

    v6=addNode!(fg,"x6",1.25*randn(1,N)+getVal(v5)+40.0, N=N)
    addFactor!(fg,[v5;v6],Odo([40.0]',[1.25]',[1.0]))

    v7=addNode!(fg,"x7",2.0*randn(1,N)+getVal(v6) +60.0, N=N)
    addFactor!(fg,[v6;v7],Odo([60.0]',[2.0]',[1.0]))

    f3 = addFactor!(fg,[v7], Obsv2(doors, cov', [1.0]))

    # HMM computed ground truth, extended for 7 poses with landmark
    gt = Dict{ASCIIString, Array{Float64,2}}()
    gt["x1"]=([0.0;1.97304 ]')' # -0.0342366
    gt["x2"]=([50.0; 2.83153 ]')' # 49.8797
    gt["x3"]=([100.0; 1.65557 ]')' # 99.8351
    gt["x4"]=([150.0; 1.64945 ]')' # 148.637
    gt["x5"]=([200.0; 1.77992 ]')' # 198.62
    gt["x6"]=([240.0; 2.20466 ]')' # 238.492
    gt["x7"]=([300.0; 2.14353 ]')' # 298.467
    gt["l1"]=([165.0; 1.17284 ]')' # 164.102
end

writeGraphPdf(fg);
tree = prepBatchTree!(fg,drawpdf=true);

# spyCliqMat(tree.cliques[1])

# list vertices in fg
@show xx,ll = ls(fg)

# from = fg.v[fg.IDs["l1"]]
# to = fg.v[fg.IDs["x7"]]
# fgs = subgraphShortestPath(fg, from=from, to=to)

# using Graphs
# el = shortest_path(fg.g, ones(19),from, to)

# do belief propagation inference over tree once
[inferOverTree!(fg, tree, N=100) for i in 1:1];

# draw all beliefs
DOYTICKS = false
xx,ll = ls(fg)
msgPlots = drawHorBeliefsList(fg, xx, gt=gt,nhor=2);
evalstr = ""
for i in 1:length(msgPlots)
    evalstr = string(evalstr, ",msgPlots[$(i)]")
end
pl = eval(parse(string("vstack(",evalstr[2:end],")")));
Gadfly.draw(PNG("4doors.png",15cm,20cm),pl) # can also do PNG
# Gadfly.set_default_plot_size(15cm, 20cm)
# pl

# if false
#
#   # draw upward messages
#   msgPlots=drawTreeUpwardMsgs(fg, tree, N=500); # to init memory for eval(parse(string))
#   println("Upward messages for all cliques except root")
#   # vvMsgs = vstackedDensities(fg, tree, msgPlots)
#   Gadfly.set_default_plot_size(17cm, 15cm)
#   Gadfly.draw(PNG("results/testMsgs.png",17cm,15cm),vvMsgs)
#   # vvMsgs
# end

#
# using ProtoBuf
#
# iob = PipeBuffer()
# dat = IncrementalInference.dlapi.getvertex(fg,1).attributes["data"]
# pd = convert(PackedVariableNodeData, dat) #fg.v[1].attributes["data"]
# writeproto(iob, pd)
# dd = readproto(iob, PackedVariableNodeData());
# unpckd = convert(VariableNodeData, dd)
#


text:
using MultiFilePatterns
using Base.Test

@test get("A000a0001", "A") == 0
@test get("A000a0001", "a") == 1
@test get("A123B456", "A") == 123
@test get("A123B456", "B") == 456
@test get("A123B456", "A1") == 23
@test get("A123B456", "3B45") == 6
@test get("A123-456", "A") == 123
@test get("A123-456", "-") == 456
@test get("A123 456", "A") == 123
@test get("A123 456", " ") == 456
@test get("A123 456", " 4") == 56
@test get("A123 456", "3 ") == 456

@test get("A123.456", "A") == 123
#@test get("A123.456", ".") == 456  # fails

@test get("A00B001C245D76", "A", "B", "C", "D") == [0, 1, 245, 76]

@test_throws Exception get("A123 456", "3")

@test get("A123A456", "A") == 123  # problematic

@test set("A0B002", "A", 1) == "A1B002"
@test set("A0B002", "A", 12) == "A12B002"
@test set("A0B002", "A", 03) == "A3B002"
@test set("A0B002", "B", 0) == "A0B000"
@test set("A0B002", "B", 12) == "A0B012"
@test set("A0B002", "B", 123) == "A0B123"
@test set("A0B002", "B", 1234) == "A0B1234"
@test set("A0B002", "B00", 1234) == "A0B001234"

@test_throws Exception get("A0B002", "B002", 6)
@test_throws Exception get("A0B002", "A0", 6)

@test set("A0B002", "A", 12, "B", 6) == "A12B006"

@test set("A0B002A0", "A", 12, "B", 6) == "A12B006A12"


text:
import NLopt

struct NLoptArgmaxAlg <: OptimAlgorithm end
"NLopt based optimization"
const NLoptArgmax = NLoptArgmaxAlg()
defΩ(::NLoptArgmaxAlg) = LinearΩ{Vector{Int}, UnitRange{Int}, Vector{Float64}}

"NLOpt style loss function (accepting vectors for input/grad) from RandVar `y`"
function nllossfunc(y, ω; usegrad = true)
  i = 1
  function innernllossfunc(ωvec::Vector, grad::Vector)
    ω = unlinearize(ωvec, ω)  # FIXME: Unlinearization requires a structure (inefficient?)
    loss = y(ω)   #
    if usegrad
      # FIXME: Set the gradient
    end
    lens(Loop, (loss = loss, i = i))
    i += 1
    loss
  end
end

function nlopt(n, lossf; alg = :LN_COBYLA, maximize = true)
  opt = NLopt.Opt(alg, n)
  NLopt.lower_bounds!(opt, 0.0)
  NLopt.upper_bounds!(opt, 1.0)
  if maximize
    NLopt.max_objective!(opt, lossf)
  else
    NLopt.min_objective!(opt, lossf)
  end
  opt
end

"ω which maximizes x(ω) using NLopt"
function Base.argmax(x::RandVar,
                     alg::NLoptArgmaxAlg,
                     ΩT::Type{OT};
                     nloptargs = ()) where {OT <: Ω}
  # initialize at random point
  ωinit = ΩT()
  x(ωinit)
  ωinitvec = linearize(ωinit)
  lossf = nllossfunc(x, ωinit)

  #  reate optimization
  n = length(ωinitvec)
  opt = nlopt(n, lossf; nloptargs...)

  # Optimize
  (minf, ωvecoptim, ret) = NLopt.optimize(opt, ωinitvec)
  unlinearize(ωvecoptim, ωinit)
end


text:
# ---
# title: 32. Longest Valid Parentheses
# id: problem32
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/longest-valid-parentheses/description/>
# hidden: true
# ---
#
# Given a string containing just the characters `'('` and `')'`, find the length
# of the longest valid (well-formed) parentheses substring.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "(()"
#     Output: 2
#     Explanation: The longest valid parentheses substring is "()".
#
#
# **Example 2:**
#
#
#
#     Input: s = ")()())"
#     Output: 4
#     Explanation: The longest valid parentheses substring is "()()".
#
#
# **Example 3:**
#
#
#
#     Input: s = ""
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 3 * 104`
#   * `s[i]` is `'('`, or `')'`.
#
#
## @lc code=start
using LeetCode

function longest_valid_parentheses(s::String)::Int
    dp = zeros(Int, length(s))
    n = 0
    for i in 2:length(s)
        if s[i] == ')'
            if s[i - 1] == '('
                ## case 1: ()()
                dp[i] = i == 2 ? 2 : dp[i - 2] + 2
            else
                ## case 2: (())
                i₍ = i - dp[i - 1] - 1
                if i₍ >= 1 && s[i₍] == '('
                    if dp[i - 1] > 0
                        dp[i] = dp[i - 1] + 2 + (i₍ == 1 ? 0 : dp[i₍ - 1])
                    end
                end
            end
            n = max(n, dp[i])
        end
    end
    return n
end
## @lc code=end


text:
struct Match
    rule::AbstractRule
    # the rhs pattern to instantiate
    pat_to_inst::Union{Nothing,Pattern}
    # the substitution
    sub::Sub
    # the id the matched the lhs
    id::EClassId
end

const MatchesBuf = Vector{Match}

function cached_ids(g::EGraph, p::Pattern)# ::Vector{Int64}
    collect(keys(g.classes))
end

# FIXME
function cached_ids(g::EGraph, p::PatTerm)
    # println("pattern $p, $(p.head)")
    # println("all ids")
    # keys(g.classes) |> println
    # println("cached symbols")
    # cached = get(g.symcache, p.head, Set{Int64}())
    # println("symbols where $(p.head) appears")
    # appears = Set{Int64}()
    # for (id, class) ∈ g.classes
    #     for n ∈ class
    #         if n.head == p.head
    #             push!(appears, id)
    #         end
    #     end
    # end
    # # println(appears)
    # if !(cached == appears)
    #     @show cached
    #     @show appears
    # end

    collect(keys(g.classes))
    # cached
    # get(g.symcache, p.head, [])
end

# function cached_ids(g::EGraph, p::PatLiteral)
#     get(g.symcache, p.val, [])
# end

function (r::SymbolicRule)(g::EGraph, id)
    ematch(g, r.ematch_program, id) .|> sub -> Match(r, r.right, sub, id)
end

function (r::DynamicRule)(g::EGraph, id)
    ematch(g, r.ematch_program, id) .|> sub -> Match(r, nothing, sub, id)
end

function (r::BidirRule)(g::EGraph, id)
    vcat(ematch(g, r.ematch_program_l, id) .|> sub -> Match(r, r.right, sub, id),
        ematch(g, r.ematch_program_r, id) .|> sub -> Match(r, r.left, sub, id))
end

function eqsat_search_threaded!(egraph::EGraph, theory::Vector{<:AbstractRule},
        scheduler::AbstractScheduler)::MatchesBuf
    matches = MatchesBuf()
    mlock = ReentrantLock()

    inequalities = filter(Base.Fix2(isa, UnequalRule), theory)
    # never skip contradiction checks
    Threads.@threads for rule ∈ inequalities
        ids = cached_ids(egraph, rule.left)

        lock(mlock) do
            for id in ids
                append!(matches, rule(egraph, id))
            end
        end
    end

    other_rules = filter(theory) do rule
        !(rule isa UnequalRule)
    end
    Threads.@threads for rule ∈ other_rules
        # don't apply banned rules
        if !cansearch(scheduler, rule)
            # println("skipping banned rule $rule")
            continue
        end

        rule_matches = []
        ids = cached_ids(egraph, rule.left)
        for id in ids
            append!(rule_matches, rule(egraph, id))
        end

        can_yield = inform!(scheduler, rule, rule_matches)
        if can_yield
            lock(mlock) do
                append!(matches, rule_matches)
            end
        end
    end
    return matches
end

function eqsat_search!(egraph::EGraph, theory::Vector{<:AbstractRule},
    scheduler::AbstractScheduler)::MatchesBuf
    matches = MatchesBuf()

    inequalities = filter(Base.Fix2(isa, UnequalRule), theory)
    # never skip contradiction checks
    for rule ∈ inequalities
        ids = cached_ids(egraph, rule.left)

        for id in ids
            append!(matches, rule(egraph, id))
        end
    end

    other_rules = filter(theory) do rule
        !(rule isa UnequalRule)
    end
    for rule ∈ other_rules
    # don't apply banned rules
        if !cansearch(scheduler, rule)
        # println("skipping banned rule $rule")
            continue
        end

        rule_matches = Match[]
        ids = cached_ids(egraph, rule.left)
        for id in ids
            append!(rule_matches, rule(egraph, id))
        end
        can_yield = inform!(scheduler, rule, rule_matches)
        if can_yield
            append!(matches, rule_matches)
        end
    end
    return matches
end


text:
"""
    FeMat{T,S}

Term with an explicit, constant matrix representation

# Fields
* `x`: matrix
* `wtx`: weighted matrix
* `piv`: pivot `Vector{Int}`` for pivoted Cholesky factorization of `wtx'wtx`
* `rank`: computational rank of `x`
* `cnames`: vector of column names
"""
mutable struct FeMat{T,S<:AbstractMatrix}
    x::S
    wtx::S
    piv::Vector{Int}
    rank::Int
    cnames::Vector{String}
end

function FeMat(X::AbstractMatrix, cnms)
    T = eltype(X)
    if size(X,2) > 0
        ch = statscholesky(Symmetric(X'X))
        pivot = ch.piv
        rank = ch.rank
        Xp = all(pivot .== 1:size(X, 2)) ? X : X[:, ch.piv]
        if rank < length(pivot)
            @warn "fixed-effects matrix is rank deficient"
        end
    else
        # although it doesn't take long for an empty matrix,
        # we can still skip the Cholesky step, which gets the rank
        # wrong anyway
        pivot = Int[]
        Xp = X
        rank = 0
    end
    # ch.rank is wrong for empty FE
    FeMat{T,typeof(X)}(Xp, Xp, pivot, rank, cnms[pivot])
end

function reweight!(A::FeMat{T}, sqrtwts::Vector{T}) where {T}
    if !isempty(sqrtwts)
        if A.x === A.wtx
            A.wtx = similar(A.x)
        end
        mul!(A.wtx, Diagonal(sqrtwts), A.x)
    end
    A
end

Base.adjoint(A::FeMat) = Adjoint(A)

Base.eltype(A::FeMat{T}) where {T} = T

fullrankwtx(A::FeMat) = rank(A) == size(A, 2) ? A.wtx : A.wtx[:, 1:rank(A)]

Base.length(A::FeMat) = length(A.wtx)

LinearAlgebra.rank(A::FeMat) = A.rank

Base.size(A::FeMat) = size(A.wtx)

Base.size(A::Adjoint{T,<:FeMat{T}}) where {T} = reverse(size(A.parent))

Base.size(A::FeMat, i) = size(A.wtx, i)

Base.copyto!(A::FeMat{T}, src::AbstractVecOrMat{T}) where {T} = copyto!(A.x, src)

*(adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T} =
    fullrankwtx(adjA.parent)' * fullrankwtx(B)

LinearAlgebra.mul!(R::StridedVecOrMat{T}, A::FeMat{T}, B::StridedVecOrMat{T}) where {T} =
    mul!(R, A.x, B)

LinearAlgebra.mul!(C, adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T} =
    mul!(C, fullrankwtx(adjA.parent)', fullrankwtx(B))

"""
    isfullrank(A::FeMat)

Does `A` have full column rank?
"""
isfullrank(A::FeMat) = A.rank == length(A.piv)

nθ(A::FeMat) = 0


text:
using .Abstract: returntype
using IRTools: Variable, returnvalue, blocks, isexpr
using IRTools.Inner: iscall

struct Trivial end

function infer(f, ::Trivial, tr = trace(typeof(f)))
  r = returntype(tr)
  r isa Abstract.Const && return Singleton(r.value)
  any(((v, st),) -> iscall(st.expr, observe), tr) && return
  r = returnvalue(blocks(tr)[end])
  r isa Variable || return
  ex = tr[r].expr
  iscall(ex, rand) && ex.args[2] isa Distribution && return ex.args[2]
  return
end


text:
@inline context() = nothing
@inline contextualized(::Type{T}) where {T} = false

macro context(Ctx)
    esc(quote
        @inline $TinyCassette.execute(ctx::$Ctx, ::typeof($context)) = ctx
        @inline $TinyCassette.execute(ctx::$Ctx, ::typeof($contextualized), ::Type{T}) where {T} = $Ctx <: T
    end)
end

macro withctx(Ctx, Fn::Expr)
    @assert isexpr(Fn, :call) "second argument must be a function call"

    fnName = Fn.args[1]
    fnArgs = Fn.args[2:end]

    esc(Expr(:call, :($TinyCassette.execute), Ctx, fnName, fnArgs...))
end


text:
#
# Autogenerated by Thrift Compiler (1.0.0-dev)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

type _enum_PartialKey
  ROW::Int32
  ROW_COLFAM::Int32
  ROW_COLFAM_COLQUAL::Int32
  ROW_COLFAM_COLQUAL_COLVIS::Int32
  ROW_COLFAM_COLQUAL_COLVIS_TIME::Int32
  ROW_COLFAM_COLQUAL_COLVIS_TIME_DEL::Int32
end
const PartialKey = _enum_PartialKey(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5))

type _enum_TablePermission
  READ::Int32
  WRITE::Int32
  BULK_IMPORT::Int32
  ALTER_TABLE::Int32
  GRANT::Int32
  DROP_TABLE::Int32
end
const TablePermission = _enum_TablePermission(Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

type _enum_SystemPermission
  GRANT::Int32
  CREATE_TABLE::Int32
  DROP_TABLE::Int32
  ALTER_TABLE::Int32
  CREATE_USER::Int32
  DROP_USER::Int32
  ALTER_USER::Int32
  SYSTEM::Int32
end
const SystemPermission = _enum_SystemPermission(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

type _enum_ScanType
  SINGLE::Int32
  BATCH::Int32
end
const ScanType = _enum_ScanType(Int32(0), Int32(1))

type _enum_ScanState
  IDLE::Int32
  RUNNING::Int32
  QUEUED::Int32
end
const ScanState = _enum_ScanState(Int32(0), Int32(1), Int32(2))

type _enum_ConditionalStatus
  ACCEPTED::Int32
  REJECTED::Int32
  VIOLATED::Int32
  UNKNOWN::Int32
  INVISIBLE_VISIBILITY::Int32
end
const ConditionalStatus = _enum_ConditionalStatus(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4))

type _enum_Durability
  DEFAULT::Int32
  NONE::Int32
  LOG::Int32
  FLUSH::Int32
  SYNC::Int32
end
const Durability = _enum_Durability(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4))

type _enum_CompactionType
  MINOR::Int32
  MERGE::Int32
  MAJOR::Int32
  FULL::Int32
end
const CompactionType = _enum_CompactionType(Int32(0), Int32(1), Int32(2), Int32(3))

type _enum_CompactionReason
  USER::Int32
  SYSTEM::Int32
  CHOP::Int32
  IDLE::Int32
  CLOSE::Int32
end
const CompactionReason = _enum_CompactionReason(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4))

type _enum_IteratorScope
  MINC::Int32
  MAJC::Int32
  SCAN::Int32
end
const IteratorScope = _enum_IteratorScope(Int32(0), Int32(1), Int32(2))

type _enum_TimeType
  LOGICAL::Int32
  MILLIS::Int32
end
const TimeType = _enum_TimeType(Int32(0), Int32(1))

type Key
  row::Vector{UInt8}
  colFamily::Vector{UInt8}
  colQualifier::Vector{UInt8}
  colVisibility::Vector{UInt8}
  timestamp::Int64
  Key() = (o=new(); fillunset(o); o)
end # type Key
meta(t::Type{Key}) = meta(t, Symbol[:timestamp], Int[], Dict{Symbol,Any}(:timestamp => Int64(9223372036854775807)))

type ColumnUpdate
  colFamily::Vector{UInt8}
  colQualifier::Vector{UInt8}
  colVisibility::Vector{UInt8}
  timestamp::Int64
  value::Vector{UInt8}
  deleteCell::Bool
  ColumnUpdate() = (o=new(); fillunset(o); o)
end # type ColumnUpdate
meta(t::Type{ColumnUpdate}) = meta(t, Symbol[:colVisibility,:timestamp,:value,:deleteCell], Int[], Dict{Symbol,Any}())

type DiskUsage
  tables::Vector{UTF8String}
  usage::Int64
  DiskUsage() = (o=new(); fillunset(o); o)
end # type DiskUsage

type KeyValue
  key::Key
  value::Vector{UInt8}
  KeyValue() = (o=new(); fillunset(o); o)
end # type KeyValue

type ScanResult
  results::Vector{KeyValue}
  more::Bool
  ScanResult() = (o=new(); fillunset(o); o)
end # type ScanResult

type Range
  start::Key
  startInclusive::Bool
  stop::Key
  stopInclusive::Bool
  Range() = (o=new(); fillunset(o); o)
end # type Range

type ScanColumn
  colFamily::Vector{UInt8}
  colQualifier::Vector{UInt8}
  ScanColumn() = (o=new(); fillunset(o); o)
end # type ScanColumn
meta(t::Type{ScanColumn}) = meta(t, Symbol[:colQualifier], Int[], Dict{Symbol,Any}())

type IteratorSetting
  priority::Int32
  name::UTF8String
  iteratorClass::UTF8String
  properties::Dict{UTF8String,UTF8String}
  IteratorSetting() = (o=new(); fillunset(o); o)
end # type IteratorSetting

type ScanOptions
  authorizations::Set{Vector{UInt8}}
  range::Range
  columns::Vector{ScanColumn}
  iterators::Vector{IteratorSetting}
  bufferSize::Int32
  ScanOptions() = (o=new(); fillunset(o); o)
end # type ScanOptions
meta(t::Type{ScanOptions}) = meta(t, Symbol[:authorizations,:range,:columns,:iterators,:bufferSize], Int[], Dict{Symbol,Any}())

type BatchScanOptions
  authorizations::Set{Vector{UInt8}}
  ranges::Vector{Range}
  columns::Vector{ScanColumn}
  iterators::Vector{IteratorSetting}
  threads::Int32
  BatchScanOptions() = (o=new(); fillunset(o); o)
end # type BatchScanOptions
meta(t::Type{BatchScanOptions}) = meta(t, Symbol[:authorizations,:ranges,:columns,:iterators,:threads], Int[], Dict{Symbol,Any}())

type KeyValueAndPeek
  keyValue::KeyValue
  hasNext::Bool
  KeyValueAndPeek() = (o=new(); fillunset(o); o)
end # type KeyValueAndPeek

type KeyExtent
  tableId::UTF8String
  endRow::Vector{UInt8}
  prevEndRow::Vector{UInt8}
  KeyExtent() = (o=new(); fillunset(o); o)
end # type KeyExtent

type Column
  colFamily::Vector{UInt8}
  colQualifier::Vector{UInt8}
  colVisibility::Vector{UInt8}
  Column() = (o=new(); fillunset(o); o)
end # type Column

type Condition
  column::Column
  timestamp::Int64
  value::Vector{UInt8}
  iterators::Vector{IteratorSetting}
  Condition() = (o=new(); fillunset(o); o)
end # type Condition
meta(t::Type{Condition}) = meta(t, Symbol[:timestamp,:value,:iterators], Int[], Dict{Symbol,Any}())

type ConditionalUpdates
  conditions::Vector{Condition}
  updates::Vector{ColumnUpdate}
  ConditionalUpdates() = (o=new(); fillunset(o); o)
end # type ConditionalUpdates
meta(t::Type{ConditionalUpdates}) = meta(t, Symbol[], Int[2,3], Dict{Symbol,Any}())

type ConditionalWriterOptions
  maxMemory::Int64
  timeoutMs::Int64
  threads::Int32
  authorizations::Set{Vector{UInt8}}
  durability::Int32
  ConditionalWriterOptions() = (o=new(); fillunset(o); o)
end # type ConditionalWriterOptions
meta(t::Type{ConditionalWriterOptions}) = meta(t, Symbol[:maxMemory,:timeoutMs,:threads,:authorizations,:durability], Int[], Dict{Symbol,Any}())

type ActiveScan
  client::UTF8String
  user::UTF8String
  table::UTF8String
  age::Int64
  idleTime::Int64
  _type::Int32
  state::Int32
  extent::KeyExtent
  columns::Vector{Column}
  iterators::Vector{IteratorSetting}
  authorizations::Vector{Vector{UInt8}}
  ActiveScan() = (o=new(); fillunset(o); o)
end # type ActiveScan

type ActiveCompaction
  extent::KeyExtent
  age::Int64
  inputFiles::Vector{UTF8String}
  outputFile::UTF8String
  _type::Int32
  reason::Int32
  localityGroup::UTF8String
  entriesRead::Int64
  entriesWritten::Int64
  iterators::Vector{IteratorSetting}
  ActiveCompaction() = (o=new(); fillunset(o); o)
end # type ActiveCompaction

type WriterOptions
  maxMemory::Int64
  latencyMs::Int64
  timeoutMs::Int64
  threads::Int32
  durability::Int32
  WriterOptions() = (o=new(); fillunset(o); o)
end # type WriterOptions
meta(t::Type{WriterOptions}) = meta(t, Symbol[:durability], Int[], Dict{Symbol,Any}())

type CompactionStrategyConfig
  className::UTF8String
  options::Dict{UTF8String,UTF8String}
  CompactionStrategyConfig() = (o=new(); fillunset(o); o)
end # type CompactionStrategyConfig

type UnknownScanner <: Exception
  msg::UTF8String
  UnknownScanner() = (o=new(); fillunset(o); o)
end # type UnknownScanner

type UnknownWriter <: Exception
  msg::UTF8String
  UnknownWriter() = (o=new(); fillunset(o); o)
end # type UnknownWriter

type NoMoreEntriesException <: Exception
  msg::UTF8String
  NoMoreEntriesException() = (o=new(); fillunset(o); o)
end # type NoMoreEntriesException

type AccumuloException <: Exception
  msg::UTF8String
  AccumuloException() = (o=new(); fillunset(o); o)
end # type AccumuloException

type AccumuloSecurityException <: Exception
  msg::UTF8String
  AccumuloSecurityException() = (o=new(); fillunset(o); o)
end # type AccumuloSecurityException

type TableNotFoundException <: Exception
  msg::UTF8String
  TableNotFoundException() = (o=new(); fillunset(o); o)
end # type TableNotFoundException

type TableExistsException <: Exception
  msg::UTF8String
  TableExistsException() = (o=new(); fillunset(o); o)
end # type TableExistsException

type MutationsRejectedException <: Exception
  msg::UTF8String
  MutationsRejectedException() = (o=new(); fillunset(o); o)
end # type MutationsRejectedException

abstract AccumuloProxyClientBase


text:
# centroid is the geometric centre

centroid(x) = x

centroid(d::Circle) = d.center
centroid(d::Interval) = mean((first(d),last(d)))

# dist and dist2 are the nearest distance and distance-squared between x and y

dist(x,y) = sqrt(dist2(x,y))
dist2(x,y) = dist(x,y)^2

dist2(x::Number,y::Number) = abs2(x-y)
dist2(d1::Domain,d2::Domain) = minimum(extrema2(d1,d2))

# diam is the diameter of x

diam(x) = sqrt(diam2(x))
diam2(x) = diam(x)^2

diam2(x::Number) = 1
diam(d::Circle) = 2d.radius
diam2(d::Interval) = abs2(first(d)-last(d))
diam2(d::UnionDomain) = maximum(extrema2(d))

# Extremal distances between domains

function Base.extrema(d1::Domain,d2::Domain)
    ext2 = extrema2(d1,d2)
    sqrt(ext2[1]),sqrt(ext2[2])
end
function extrema2(d1::Domain,d2::Domain)
    ext = extrema(d1,d2)
    ext[1]^2,ext[2]^2
end

for op in (:(Base.extrema),:extrema2)
    @eval begin
        function $op(d::UnionDomain)
            ext = [$op(d1,d2) for d1 in d,d2 in d]
            minimum(minimum(ext)),maximum(maximum(ext))
        end
        function $op(d1::UnionDomain,d2::UnionDomain)
            ext = [$op(d11,d22) for d11 in d1,d22 in d2]
            minimum(minimum(ext)),maximum(maximum(ext))
        end
        function $op(d1::Domain,d2::UnionDomain)
            ext = [$op(d1,d22) for d22 in d2]
            minimum(minimum(ext)),maximum(maximum(ext))
        end
        $op(d1::UnionDomain,d2::Domain) = $op(d2,d1)
    end
end

#
# Geometry for non-intersecting Intervals
#
function dist2(c::Number,d::Interval)
    if in(c,d)
        zero(real(c))
    else
        a,b = d.a,d.b
        x1,y1 = reim(a)
        x2,y2 = reim(b)
        x3,y3 = reim(c)
        px,py = x2-x1,y2-y1
        u = ((x3-x1)px+(y3-y1)py)/(px^2+py^2)
        u = u > 1 ? 1 : u ≥ 0 ? u : 0
        dx,dy = x1+u*px-x3,y1+u*py-y3
        dx^2+dy^2
    end
end

function extrema2(d1::Interval,d2::Interval)
    a,b = d1.a,d1.b
    c,d = d2.a,d2.b
    extrema((dist2(a,d2),dist2(b,d2),dist2(c,d1),dist2(d,d1),abs2(a-c),abs2(a-d),abs2(b-c),abs2(b-d)))
end

#
# Geometry for non-intersecting Circles
#
function dist2(z::Number,d::Circle)
    if in(z,d)
        zero(real(z))
    else
        abs2(z-d.center-d.radius*cis(angle(z-d.center)))
    end
end

function extrema2(d1::Circle,d2::Circle)
    r1,r2=d1.radius,d2.radius
    c1,c2=d1.center,d2.center
    if  d1 == d2 # Circles are equal
        (0,(2d1.radius)^2)
    elseif r1>r2&&abs(c1-c2)<r1 || r1<r2&&abs(c1-c2)<r2 # Circles are interior/exterior
        extrema((abs2(c1-c2-(r1-r2)*cis(angle(c1-c2))),
                 abs2(c1-c2+(r1+r2)*cis(angle(c1-c2))),
                 abs2(c1-c2-(r1+r2)*cis(angle(c1-c2))),
                 abs2(c1-c2+(r1-r2)*cis(angle(c1-c2)))))
    else # Circles are disjoint
        extrema((abs2(c1-c2-(r1+r2)*cis(angle(c1-c2))),
                 abs2(c1-c2+(r1+r2)*cis(angle(c1-c2)))))
    end
end

#
# Geometry between non-intersecting Intervals and Circles
#
# Consider the minimal distances from a & b to the center, adding or subtracting
# a radius in either direction. Add to these combinations the minimal distance from the center
# to the Interval, adding or subtracting a radius in either direction.
# The extremal distances are a subset.
#
function Base.extrema(d1::Interval,d2::Circle)
    a,b = d1.a,d1.b
    c,r=d2.center,d2.radius
    extrema(( dist(a,d2),dist(b,d2),dist(a,d2)+2r,dist(b,d2)+2r,dist(c,d1)+r,dist(c,d1)-r ))
end
Base.extrema(d1::Circle,d2::Interval) = extrema(d2,d1)


text:
macro wandbfinish(run)
    quote
        if !("JULIA_NO_WANDB" in keys(ENV)) || ENV["JULIA_NO_WANDB"] == "0" || ENV["JULIA_NO_WANDB"] == "false" || ENV["JULIA_NO_WANDB"] == false
            $(esc(run)).finish()
        end
    end
end

macro wandbfinish()
    quote
        if !("JULIA_NO_WANDB" in keys(ENV)) || ENV["JULIA_NO_WANDB"] == "" || ENV["JULIA_NO_WANDB"] == "false" || ENV["JULIA_NO_WANDB"] == false
            wandb.finish()
        end
    end
end


text:
using .BifurcationKit
using .Setfield

# ------------------------------------------------------------------------------------------------ #
#                           Generate a bifurcation diagram from a Process                          #
# ------------------------------------------------------------------------------------------------ #
"""Ok, so this one plots a bifurcation diagram using the equations defining a Process.
You can give the single parameter you want to vary, param, and the range over which you want to vary it, prange. If you don't supply prange, I'll take the parameter range as the extrema of the parameter function in the Process. Note that this will only work for Processes defined using ODEproblems from DifferentialEquations. We'll discard the jacobian and compute it automatically.
"""
function bifurcation_diagram(P::Process; param=1, prange=extrema(parameterseries(P, p=param)), kwargs...)
    f = constantine(P; param)
    ode = ODEProblem(f, P.X0, (P.transient_t0, P.tmax); atol = 1e-10, rtol = 1e-9)
    opts = ContinuationPar(pMin = min(prange...)|>Float64,
                            pMax = max(prange...)|>Float64,
                            ds = 0.01,
                            dsmax = 0.05,
	                        # nInversion = 16,
                            # detectBifurcation = 3,
                            # maxBisectionSteps = 50,
                            nev = 3,
                            kwargs...) # parameterseries(P, p=param)[1]
    br, = continuation((x, p) -> Array(f(x, p, 0.0)), P.X0, [min(prange...)],
                            (@lens _[1]), opts;
                            plot = true,
                            verbosity = 0)
end
export bifurcation_diagram

# jet = BifurcationKit.get3Jet(f)
# opts = ContinuationPar(pMin = min(prange...)|>Float64,
#                         pMax = max(prange...)|>Float64,
#                         ds = 0.01,
#                         dsmax = 0.05,
#                         # nInversion = 16,
#                         # detectBifurcation = 3,
#                         # maxBisectionSteps = 50,
#                         nev = 3)
# bifurcationdiagram(jet...,
#     P.X0, 10.0, (@lens _[1]),
# 	# important argument: this is the maximal
# 	# recursion level
# 	5,
#     opts)

function constantine(P::Process; param)
    # Want to turn our special non-stationary EOM to a normal, constant parameter one
    eom = P.process
    p_profiles = forcevec(deepcopy(P.parameter_profile))
    p_parameters= forcevec(deepcopy(P.parameter_profile_parameters))
    p_profiles′ = (p_profiles[1:param-1]..., constantParameter, p_profiles[param+1:end]...)
    function f(X::AbstractArray, 𝑝, 𝑡)
        p_parameters′ = (p_parameters[1:param-1]..., 𝑝[1], p_parameters[param+1:end]...)
        if p_profiles′ isa Union{Tuple, Vector} && length(p_profiles′) == 1
            p_profiles′ = p_profiles′[1]
        end
        pf = tuplef2ftuple(p_profiles′, p_parameters′)
        eom(X, pf, 𝑡)
    end
    return f
end

# optnew = NewtonPar(verbose = true, tol = 1e-12)
# opts = ContinuationPar(dsmin = 0.0001, dsmax = 0.01, ds = 0.01, pMax = 1.,
# 	newtonOptions = setproperties(optnew; maxIter = 30, tol = 1e-8),
# 	maxSteps = 300, plotEveryStep = 40,
# 	detectBifurcation = 3, nInversion = 4, tolBisectionEigenvalue = 1e-17, dsminBisection = 1e-7)

#     function optrec(x, p, l; opt = opts)
#         return setproperties(opt; maxSteps = 300, detectBifurcation = 3, detectLoop = false)
#     end

# g = (x, p) -> collect(f(x, p, 0.0))
# J = (x, p) -> x0 -> BifurcationKit.ForwardDiff.jacobian(y->g(y, p), x0)
# jet = BifurcationKit.get3Jet(g, J)
# diagram = @time bifurcationdiagram(jet...,
# 	[0.1], 1.0, (@lens _[1]),
# 	# here we specify a maximum branching level of 4
# 	4, optrec)


text:
# Line "drawing" iterator
struct LineIterator{N}
    start :: CartesianIndex{N}
    ϕ     :: Float64
    θ     :: Float64
end

RandomLineIterator(start) =
    LineIterator(start, 2π*rand(Float64), π*rand(Float64) - π/2)

Base.IteratorSize(::LineIterator) = Base.IsInfinite()
Base.iterate(iter :: LineIterator) = iter.start, 0.1

function Base.iterate(iter :: LineIterator{2}, r :: Float64)
    r    = r + sqrt(2)
    x, y = iter.start[1], iter.start[2]
    ϕ    = iter.ϕ

    xn = x + r*cos(ϕ) |> floor |> Int
    yn = y + r*sin(ϕ) |> floor |> Int
    return CartesianIndex(xn, yn), r
end

function Base.iterate(iter :: LineIterator{3}, r :: Float64)
    r       = r + sqrt(3)
    x, y, z = iter.start[1], iter.start[2], iter.start[3]
    ϕ       = iter.ϕ
    θ       = iter.θ

    xn = x + r*cos(θ)*cos(ϕ) |> floor |> Int
    yn = y + r*cos(θ)*sin(ϕ) |> floor |> Int
    zn = z + r*sin(θ)        |> floor |> Int
    return CartesianIndex(xn, yn, zn), r
end

# Theoretical two-point functions for the void phase for overlapping
# spheres with centers generated by Poisson process

function s2_theory(r, R, λ, :: Val{2})
    tmp = r/(2R)
    tmp2 = (r > 2R) ? 2 : 2/π*(π + tmp*sqrt(1 - tmp^2) - acos(tmp))
    η = λ * π * R^2
    return exp(-η*abs(tmp2))
end

function s2_theory(r, R, λ, :: Val{3})
    tmp = r/R
    tmp2 = (r > 2R) ? 2 : 1 + 3/4*tmp - 1/16*tmp^3
    η = λ * 4/3 * π * R^3
    return exp(-η*tmp2)
end

# Search for the best parameters for random spheres
function spheres_parameters(target :: AbstractArray{T, N},
                            R0     :: Float64,
                            λ0     :: Float64) where {T, N}
    @. model(x, p) = s2_theory(x, p[1], p[2], Val(N))
    s2 = Directional.s2(target, 0) |> mean
    r = range(0.0; length = length(s2))
    fit = curve_fit(model, r, s2, [R0, λ0])
    return coef(fit)
end

function draw_sphere!(array  :: AbstractArray{T, N},
                      center :: CartesianIndex{N},
                      R      :: Float64) where {T, N}
    Ri = R |> round |> Int

    sphere_indices = CartesianIndices(Tuple(-Ri:Ri for i in 1:N))
    foreach(sphere_indices) do point
        dist = mapreduce(x -> x^2, +, Tuple(point))
        coord = point + center
        if dist < R^2 && checkbounds(Bool, array, coord)
            @inbounds array[coord] = 1
        end
    end
end

function generate_spheres(shape :: NTuple{N, Int},
                          T     :: Type,
                          R     :: Float64,
                          λ     :: Float64) where N
    nspheres = pois_rand(λ * prod(shape))
    spheres  = zeros(T, shape)
    indices  = CartesianIndices(spheres)

    for n in 1:nspheres
        draw_sphere!(spheres, rand(indices), R)
    end

    return spheres
end


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля (без внутренних перегородок)

#РЕЗУЛЬТАТ: Робот - в исходном положении, в клетке с роботом стоит маркер, и все остальные клетки поля промаркированы в шахматном порядке

module ShessMark
    using HorizonSideRobots
    import Main.moves!, Main.inverse
#=
    После выполнеия include("roblib.jl") имена функций moves!, inverse оказались в пространстве имен Main
    и теперь их от туда необходимо импортировать в наш модуль
=#
    export mark_chess

    FLAG_MARK = nothing

    function mark_chess(r)
        global FLAG_MARK

        num_hor = moves!(r,West)
        num_vert = move!(r,Sud)
        #УТВ: Робот - в юго-западном углу
        FLAG_MARK = odd(num_hor+num_vert) ? true : false
        # значение FLAG_MARK определяет, нужно ли в юго-западном углу ставить маркер

        side = Ost
        mark_chess(r,side)
        while isborder(Nord)==false
            move!(r,Nord); FLAG_MARK = !FLAG_MARK
            side = inverse(side)
            mark_chess(r,side)
        end
        #Робот - у северной границы поля И маркеры расставлены в шахматном порядке

        for side in (West,Sud) moves!(r,side) end
        #Робот - в юго-западном углу

        move!(r,Ost,num_hor)
        move!(r,Nord,num_vert)
        #Робот - в исходном положении
    end
    function mark_chess(r,side)
        global FLAG_MARK
        while isborder(r,side)==false
            if FLAG_MARK == true
                putmarker(r)
            end
            move!(r,side); FLAG_MARK = !FLAG_MARK
        end
    end
end


text:
## VLLE Solver
function Obj_VLLE_pressure(model::EoSModel, F, T, v_l, v_ll, v_v, x, xx, y,ts,ps)
    x   = FractionVector(x)
    y   = FractionVector(y)
    xx  = FractionVector(xx)#julia magic, check misc.jl
    μ_l   = VT_chemical_potential(model,v_l,T,x)
    μ_ll  = VT_chemical_potential(model,v_ll,T,xx)
    μ_v   = VT_chemical_potential(model,v_v,T,y)
    p_l   = pressure(model,v_l,T,x)
    p_ll  = pressure(model,v_ll,T,xx)
    p_v   = pressure(model,v_v,T,y)
    n_c   = length(model.components)
    for i in 1:n_c
        F[i] = (μ_l[i]-μ_v[i])/(R̄*ts[i])
        F[i+n_c] = (μ_ll[i]-μ_v[i])/(R̄*ts[i])
    end
    F[end-1] = (p_l-p_v)/ps
    F[end]   = (p_ll-p_v)/ps
    return F
end

function VLLE_pressure(model::EoSModel, T; v0 =nothing)
#     TYPE = promote_type(eltype(T),eltype(x))
#     lb_v = lb_volume(model,x)
    if v0 === nothing
        v0 = x0_VLLE(model,T)
    end
    ts = T_scales(model,[v0[4],1-v0[4]])
    pmix = p_scale(model,[v0[4],1-v0[4]])
    len = length(v0[1:end])
    #xcache = zeros(eltype(x0),len)
    Fcache = zeros(eltype(v0[1:end]),len)
    f! = (F,z) -> Obj_VLLE_pressure(model, F, T, exp10(z[1]), exp10(z[2]), exp10(z[3]), z[4], z[5], z[6],ts,pmix)
    r  = Solvers.nlsolve(f!,v0[1:end],LineSearch(Newton()))
    sol = Solvers.x_sol(r)
    v_l = exp10(sol[1])
    v_ll = exp10(sol[2])
    v_v = exp10(sol[3])
    x   = sol[4]
    xx = sol[5]
    y  = sol[6]
    x = FractionVector(x)
    xx = FractionVector(xx)
    y = FractionVector(y)
    P_sat = pressure(model,v_l,T,x)
    return (P_sat, v_l, v_ll, v_v, x, xx, y)
end

function x0_VLLE(model::EoSModel, T)
    pure = split_model(model)
    sat  = saturation_pressure.(pure,T)
    x0    = [0.75,0.25]
    xx0   = [0.25,0.75]
    y0    = [0.5,0.5]
    v_l0  = sat[1][2]*x0[1]+sat[2][2]*x0[2]
    v_ll0 = sat[1][2]*xx0[1]+sat[2][2]*xx0[2]
    v_v0  = sat[1][3]*y0[1]+sat[2][3]*y0[2]
    return [log10(v_l0),log10(v_ll0),log10(v_v0),x0[1],xx0[1],y0[1]]
end

function VLLE_temperature(model,p)
    f(z) = Obj_VLLE_temperature(model,z,p)
    Ti   = T_scales(model,model.icomponents).*1.7
    T = Roots.find_zero(f,sum(Ti)/length(Ti))
    p,v_l,v_v,y = VLLE_pressure(model,T)
    return T,v_l,v_v,y
end

function Obj_VLLE_temperature(model,T,p)
    p̃,v_l,v_ll,xx = VLLE_pressure(model,T)
    return p̃-p
end


text:
group = addgroup!(SUITE, "array")

# generate some arrays
cpu_mat = rand(Float32, 512, 1000)
gpu_mat = CuArray{Float32}(undef, size(cpu_mat))
gpu_vec = reshape(gpu_mat, length(gpu_mat))
gpu_mat_ints = CuArray(rand(Int, 512, 1000))
gpu_vec_ints = reshape(gpu_mat_ints, length(gpu_mat_ints))
gpu_mat_bools = CuArray(rand(Bool, 512, 1000))
gpu_vec_bools = reshape(gpu_mat_bools, length(gpu_mat_bools))

group["construct"] = @benchmarkable CuArray{Int}(undef, 1)

group["copy"] = @async_benchmarkable copy($gpu_mat)

let group = addgroup!(group, "copyto!")
    gpu_mat2 = copy(gpu_mat)
    group["cpu_to_gpu"] = @async_benchmarkable copyto!($gpu_mat, $cpu_mat)
    group["gpu_to_cpu"] = @async_benchmarkable copyto!($cpu_mat, $gpu_mat)
    group["gpu_to_gpu"] = @async_benchmarkable copyto!($gpu_mat2, $gpu_mat)
end

let group = addgroup!(group, "iteration")
    group["scalar"] = @benchmarkable CUDA.@allowscalar [$gpu_vec[i] for i in 1:10]

    group["logical"] = @benchmarkable $gpu_vec[$gpu_vec_bools]

    group["findall"] = @benchmarkable findall($gpu_vec_bools)
    group["findall"] = @benchmarkable findall(isodd, $gpu_vec_ints)

    group["findfirst"] = @benchmarkable findfirst($gpu_vec_bools)
    group["findfirst"] = @benchmarkable findfirst(isodd, $gpu_vec_ints)

    let group = addgroup!(group, "findmin") # findmax
        group["1d"] = @async_benchmarkable findmin($gpu_vec)
        group["2d"] = @async_benchmarkable findmin($gpu_mat; dims=1)
    end
end

let group = addgroup!(group, "reverse")
    group["1d"] = @async_benchmarkable reverse($gpu_vec)
    group["2d"] = @async_benchmarkable reverse($gpu_mat; dims=1)
    group["1d_inplace"] = @async_benchmarkable reverse!($gpu_vec)
    group["2d_inplace"] = @async_benchmarkable reverse!($gpu_mat; dims=1)
end

group["broadcast"] = @async_benchmarkable $gpu_mat .= 0f0

# no need to test inplace version, which performs the same operation (but with an alloc)
let group = addgroup!(group, "accumulate")
    group["1d"] = @async_benchmarkable accumulate(+, $gpu_vec)
    group["2d"] = @async_benchmarkable accumulate(+, $gpu_mat; dims=1)
end

let group = addgroup!(group, "reductions")
    let group = addgroup!(group, "reduce")
        group["1d"] = @async_benchmarkable reduce(+, $gpu_vec)
        group["2d"] = @async_benchmarkable reduce(+, $gpu_mat; dims=1)
    end

    let group = addgroup!(group, "mapreduce")
        group["1d"] = @async_benchmarkable mapreduce(x->x+1, +, $gpu_vec)
        group["2d"] = @async_benchmarkable mapreduce(x->x+1, +, $gpu_mat; dims=1)
    end

    # used by sum, prod, minimum, maximum, all, any, count
end


text:
#############################################################################
#############################################################################
#
# This file implements polynomial multiplication
#
#############################################################################
#############################################################################

"""
Multiply two polynomials.
"""
function *(p1::Polynomial, p2::Polynomial)::Polynomial
    p_out = Polynomial()
    for t in p1
        p_out = p_out + (t * p2)
    end
    return p_out
end

"""
Power of a polynomial.
"""
function ^(p::Polynomial, n::Int)
    n < 0 && error("No negative power")
    out = one(p)
    for _ in 1:n
        out *= p
    end
    return out
end


text:
using Test
using ScottishTaxBenefitModel
using .STBParameters: MinimumWage, get_minimum_wage
using .ExampleHelpers

## FIXME don't need both
lmt = LegacyMeansTestedBenefitSystem{Float64}()
sys = get_system( scotland=true )

mwsys = MinimumWage()

@testset "Basic MW Retrieve 2019/20 values" begin
    @test get_minimum_wage(mwsys,18) == 6.15
    @test get_minimum_wage(mwsys,16) == 4.35
    @test get_minimum_wage(mwsys,1) == 0
    @test get_minimum_wage(mwsys,49) == 8.21
    @test get_minimum_wage(mwsys,99) == 8.21
    #= 20/1
    @test get_minimum_wage(mwsys,18) == 6.45
    @test get_minimum_wage(mwsys,16) == 4.55
    @test get_minimum_wage(mwsys,1) == 0
    @test get_minimum_wage(mwsys,49) == 8.72
    @test get_minimum_wage(mwsys,99) == 8.72
    =#
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
include("common.jl")
@testset "Unique times" begin
    prob = prob_dde_1delay_long

    @testset for constrained in (false, true)
        alg = MethodOfSteps(Tsit5(), constrained=constrained)
        sol = solve(prob, alg)

        @test allunique(sol.t)
    end
end


text:
using AugmentorAddons
using Test
using ImageIO
using ImageFiltering
using Augmentor

@testset "AugmentorAddons.jl" begin
    ref = testpattern() |> size
    lo = (50, 50)
    hi = (10, 10)

    sz = augment(RCutOut(100, 100)) |> size
    @test (@. sz == ref) |> all

    sz = augment(PadBorders(RGB(0, 0, 0), lo)) |> size
    @test (@. sz == (ref + 2 * lo)) |> all

    sz = augment(PadBorders(:reflect, lo...)) |> size
    @test (@. sz == (ref + 2 * lo)) |> all

    sz = augment(PadBorders(RGB(0, 0, 0), lo, hi)) |> size
    @test (@. sz == (ref + lo + hi)) |> all

    x = rand(32, 32, 3, 4)
    ref = size(x)
    lo = (50, 50, 0, 0)
    hi = (10, 10, 0, 0)
    sz = augment(x, PadBorders(0, lo, hi)) |> size
    @test (@. sz == (ref + lo + hi)) |> all
end


text:
import .MCMCChains: Chains

# A basic chains constructor that works with the Transition struct we defined.
function bundle_samples(
    rng::AbstractRNG,
    model::DensityModel,
    s::Metropolis,
    N::Integer,
    ts::Vector{T},
    chain_type::Type{Chains};
    param_names=missing,
    kwargs...
) where {ModelType<:AbstractModel, T<:AbstractTransition}
    # Turn all the transitions into a vector-of-vectors.
    vals = copy(reduce(hcat,[vcat(t.params, t.lp) for t in ts])')

    # Check if we received any parameter names.
    if ismissing(param_names)
        param_names = ["Parameter $i" for i in 1:length(s.init_params)]
    else
        # Deepcopy to be thread safe.
        param_names = deepcopy(param_names)
    end

    # Add the log density field to the parameter names.
    push!(param_names, "lp")

    # Bundle everything up and return a Chains struct.
    return Chains(vals, param_names, (internals=["lp"],))
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    Cylinder{T,N} <: ParametricSurface{T,N}

Cylinder of infinite height centered at the origin, oriented along the z-axis.
`visheight` is used for visualization purposes only, **note that this does not fully represent the surface**.

```julia
Cylinder(radius::T, visheight::T = 2.0; interface::NullOrFresnel{T} = nullinterface(T))
```
"""
struct Cylinder{T, N} <: ParametricSurface{T, N}
    radius::T
    visheight::T # Only used for visualization purposes.
    interface::NullOrFresnel{T}

    function Cylinder(
        radius::T,
        visheight::T = T(2.0);
        interface::NullOrFresnel{T} = NullInterface(T),
    ) where {T <: Real}
        @assert radius > zero(T) && visheight > zero(T)
        @assert !isnan(radius)
        return new{T, 3}(radius, visheight, interface)
    end
end

function Base.show(io::IO, a::Cylinder{T}) where {T <: Real}
    return print(io, "Cylinder{$T}($(a.radius), $(interface(a)))")
end

interface(a::Cylinder{T}) where {T <: Real} = a.interface
radius(a::Cylinder{T}) where {T <: Real} = a.radius

function uvrange(::Type{Cylinder{T, N}}) where {T <: Real, N}
    return ((-T(π), T(π)), (-one(T), one(T)))
end

function onsurface(cyl::Cylinder{T, 3}, x::T, y::T, ::T) where {T <: Real}
    return samepoint(x^2 + y^2, radius(cyl)^2)
end

function inside(cyl::Cylinder{T, 3}, x::T, y::T, ::T) where {T <: Real}
    return x^2 + y^2 - radius(cyl)^2 < zero(T)
end

function point(cyl::Cylinder{T, 3}, u::T, v::T) where {T <: Real}
    return SVector{3, T}(
        radius(cyl) * cos(u),
        radius(cyl) * sin(u),
        v * cyl.visheight / 2,
    )
end

function normal(::Cylinder{T, 3}, u::T, ::T) where {T <: Real}
    return SVector{3, T}(cos(u), sin(u), 0.0)
end

function partials(::Cylinder{T, 3}, u::T, ::T) where {T <: Real}
    return SVector{3, T}(-sin(u), cos(u), 0.0), SVector{3, T}(0.0, 0.0, 1.0)
end
function uv(cyl::Cylinder{T, 3}, x::T, y::T, z::T) where {T <: Real}
    return SVector{2, T}(atan(y, x), 2 * z / cyl.visheight)
end

# Assumes the ray has been transformed into the canonical cylinder coordinate frame which has the cylinder axis passing through (0,0,0) and aligned with the z axis.
function surfaceintersection(
    cyl::Cylinder{T, N},
    r::AbstractRay{T, N},
) where {T <: Real, N}
    # in the cylinder coordinate frame the cylinder axis is (0,0,1). project ray into plane perpendicular to cylinder axis by dropping z coordinate
    ox, oy, oz = origin(r)
    dx, dy, dz = direction(r)
    rad = radius(cyl)
    c = (ox^2 + oy^2) - rad^2

    if samepoint(one(T), abs(dz)) # ray is parallel to the cylinder axis
        if c > zero(T)
            # if strictly outside the cylinder then no intersection
            return EmptyInterval(T)
        else
            # ray is contained entirely in the cylinder
            return rayorigininterval(Infinity(T))
        end
    end

    a = dx^2 + dy^2
    b = 2 * (ox * dx + oy * dy)

    temp = quadraticroots(a, b, c)

    if temp === nothing
        return EmptyInterval(T) # no intersection with cylinder and ray not contained entirely in cylinder
    end

    t1, t2 = temp

    if isapprox(t1, t2, rtol = 1e-12, atol = 2 * eps(T))
        return EmptyInterval(T) # single root which indicates a tangent ray cylinder intersection
    end

    if t1 > zero(T)
        pt1 = point(r, t1)
    else
        pt1 = nothing
    end

    if t2 > zero(T)
        pt2 = point(r, t2)
    else
        pt2 = nothing
    end

    let int1 = nothing, int2 = nothing
        if pt1 !== nothing
            u, v = uv(cyl, pt1)
            int1 = Intersection(
                t1,
                pt1,
                SVector{3, T}(pt1[1], pt1[2], 0.0),
                u,
                v,
                interface(cyl),
            )
        end

        if pt2 !== nothing
            u, v = uv(cyl, pt2)
            int2 = Intersection(
                t2,
                pt2,
                SVector{3, T}(pt2[1], pt2[2], 0.0),
                u,
                v,
                interface(cyl),
            )
        end

        if int1 !== nothing && int2 !== nothing
            if t1 <= t2
                return Interval(int1, int2)
            else
                return Interval(int2, int1)
            end
        elseif int1 !== nothing
            return rayorigininterval(int1)
        elseif int2 !== nothing
            return rayorigininterval(int2)
        else
            return EmptyInterval(T)
        end
    end
end

function BoundingBox(cyl::Cylinder{T, 3}) where {T <: Real}
    return BoundingBox(
        -radius(cyl),
        radius(cyl),
        -radius(cyl),
        radius(cyl),
        typemin(T),
        typemax(T),
    )
end


text:
module TestHess
    using Calculus2
    using Base.Test

    # multivariate test function
    function f(x::Vector)
        return sin(x[1]) + cos(x[2])
    end
    x = [0.0, 0.0]
    n = length(x)
    H = Array(eltype(x), n, n)

    # Test in-place multivariate grad!
    hess!(H, f, x)
    @test Calculus2.maxdiff(
        H,
        [-sin(0.0) 0.0; 0.0 -cos(0.0)]
    ) < 10e-8

    # Test pure multivariate grad
    @test Calculus2.maxdiff(
        hess(f, x),
        [-sin(0.0) 0.0; 0.0 -cos(0.0)]
    ) < 10e-8

    # Test in-place multivariate gradfunc!
    fill!(H, 0.0)
    h! = hessfunc!(f, x)
    h!(H, x)
    @test Calculus2.maxdiff(
        H,
        [-sin(0.0) 0.0; 0.0 -cos(0.0)]
    ) < 10e-8

    # Test pure multivariate gradfunc
    h = hessfunc(f, x)
    @test Calculus2.maxdiff(
        h(x),
        [-sin(0.0) 0.0; 0.0 -cos(0.0)]
    ) < 10e-8

    # Test pure univariate grad
    @test isapprox(hess(sin, 1.0), -sin(1.0))

    # Test pure univariate gradfunc
    @test isapprox(hessfunc(sin, 0.0)(1.0), -sin(1.0))

    # TODO: Test automatic differentation
    # TODO: Test symbolic differentiation
    # TODO: Test direction options
end


text:
using Test
using Comonicon
using Comonicon.Tools: prompt

@testset "prompt(msg) input=$yes" for yes in ['Y', 'y', '\n', '\r']
    print(stdin.buffer, yes)
    @test prompt("input something") == true
    readavailable(stdin.buffer)
end

@testset "prompt(msg) input=$no" for no in ['N', 'n', 'a', 'b']
    print(stdin.buffer, no)
    @test prompt("input something") == false
    readavailable(stdin.buffer)
end

@testset "prompt(msg;yes) input=$input" for input in ['Y', 'y', '\n', '\r', 'N', 'n', 'a', 'b']
    print(stdin.buffer, input)
    @test prompt("input something"; yes = true) == true
    readavailable(stdin.buffer)
end

@testset "prompt(msg; require=true)" begin
    print(stdin.buffer, '\n')
    print(stdin.buffer, '\n')
    print(stdin.buffer, 'Y')
    @test prompt("input something"; require = true) == true
    readavailable(stdin.buffer)

    print(stdin.buffer, '\n')
    print(stdin.buffer, '\n')
    print(stdin.buffer, 'n')
    @test prompt("input something"; require = true) == false
    readavailable(stdin.buffer)

    print(stdin.buffer, '\n')
    print(stdin.buffer, '\n')
    print(stdin.buffer, '\n')
    @test_throws Comonicon.CommandError prompt("input something"; require = true)
    readavailable(stdin.buffer)
end


text:
export SVG, CSS, RCP, svg_css, svg_class, set_href, svg_document, SvgContent, GradientContent, radialGradient_template, linearGradient_template, svg_doc_recipe, svg_canvas_recipe

# Get default elements and attributes
const SVG = get_json(string(PKG_ROOT_DIR,"/assets/svg/","svg.json"))
const CSS = get_json(string(PKG_ROOT_DIR,"/assets/svg/","css.json"))
const RCP = get_json(string(PKG_ROOT_DIR,"/assets/svg/","recipes.json"))

"""
================================================================================
# SVG specific types and functions
================================================================================

A unique differenciator can be used for elements.
Example: "svg_rect"

SVG specific definitions:

A generic SVG document:

- svg element
- style element
- defs element
- main area

"""

# Run main function to create functions based on the SVG JSON file.
# A prefix is optional but recommended because of potential conflicts with Julia functions or types with same names.
# Example: svg filter element
PREF = "svg_"

generate(copy(SVG),PREF)

# function wrapper for template functions

function wrap(n,a=Dict(),c=string())
    f = Symbol(n)
    eval(quote $f($a,$c) end )
end

# Composite type: Minimal structure of a svg document
mutable struct SvgContent
    style
    defs
    main
    function SvgContent(style=wrap(string(PREF,"style")),defs=wrap(string(PREF,"defs")),main=wrap(string(PREF,"g")))
        new(style,defs,main)
    end
end

# Creates the actual svg document
function svg_document(a::Dict=SVG["svg"],c::SvgContent=SvgContent())
    return wrap(string(PREF,"svg"),a,string(c.style,c.defs,c.main))

end

# Creates a svg css class string
function svg_class(type::String,values::Dict,id::String=Random.randstring(8))
    c = string()
    for (k,v) in values c = string(c,k,":",v,";") end
    return string(type,".",id,"{",c,"} ")
end

function svg_css(children::Array{String,1})
    return string("<![CDATA[ ",join(children)," ]]>")
end

# Set href value string for styles and symbols. If true (default), the href value is being formatted for css, otherwise for element references.
function set_href(value,c::Bool=true)
      return c == true ? string("url(#",value,")") : string("#",value)
end

# Templates: Gradients
mutable struct GradientContent
    stop1
    stop2
    function GradientContent(stop1=wrap(string(PREF,"stop")),stop2=wrap(string(PREF,"stop")))
        new(stop1,stop2)
    end
end

function radialGradient_template(a::Dict=SVG["radialGradient"],c::GradientContent=GradientContent())
    return wrap(string(PREF,"radialGradient"),a,string(c.stop1,c.stop2))
end

function linearGradient_template(a::Dict=SVG["linearGradient"],c::GradientContent=GradientContent())
    return wrap(string(PREF,"linearGradient"),a,string(c.stop1,c.stop2))
end

"""

================================================================================
# SVG Recipes
================================================================================

Recipes are parametric document templates with a corresponding json file

"""

# A simple svg document. Just a rectangle with a radial gradient background
function svg_doc_recipe(p::Dict=get_attributes(RCP,@funcName))

    # Id's
    radialGradient_id = set_random_id()
    rect_style_id = set_random_id()

    # Gradients
    radialGradient = get_attributes(SVG,"radialGradient")
    radialGradient["id"] = radialGradient_id

    grad_content = GradientContent()
    grad_content.stop1 = svg_stop(Dict(
       "stop-color"=>p["color1"],
       "offset"=>p["offset1"]
    ))
    grad_content.stop2 = svg_stop(Dict(
       "stop-color"=>p["color2"],
       "offset"=>p["offset2"]
    ))

    # CSS
    rect_style = get_attributes(CSS,"rect")
    rect_style["fill"] = set_href(radialGradient_id)
    rect_style_class = svg_class("rect",rect_style,rect_style_id)

    # Canvas
    rect = get_attributes(SVG,"rect")
    delete!(rect, "style")
    rect["class"] = rect_style_id
    rect["x"] = p["x"]
    rect["y"] = p["y"]
    rect["width"] = p["w"]
    rect["height"] = p["h"]

    # Assemble the document

    doc = SvgContent()

    doc.style = svg_style(Dict(),svg_css([rect_style_class]))
    doc.defs = svg_defs(Dict("id"=>"defs"),radialGradient_template(radialGradient,grad_content))

    content = svg_rect(rect)
    doc.main = svg_g(Dict("id"=>"main"),content)

    svg =  get_attributes(SVG,"svg")
    svg["viewBox"] = join_str([p["x"],p["y"],p["w"],p["h"]])

    return svg_document(svg,doc)

end

# See abstract_rectangle in Abstract.jl
function svg_canvas_recipe(p::Dict=get_attributes(RCP,@funcName))

   # Build Gradients

   defs = string()

   if p["rect_class"] == "yes"

     rect_radialGradient_id = set_random_id()
     rect_radialGradient = get_attributes(SVG,"radialGradient")
     rect_radialGradient["id"] = rect_radialGradient_id

     rect_grad_content = GradientContent()
     rect_grad_content.stop1 = svg_stop(Dict(
       "stop-color"=>p["rect_color1"],
       "offset"=>p["rect_offset1"]
     ))
     rect_grad_content.stop2 = svg_stop(Dict(
       "stop-color"=>p["rect_color2"],
       "offset"=>p["rect_offset2"]
     ))

     defs = string(defs,radialGradient_template(rect_radialGradient,rect_grad_content))
   end

   if p["circle_class"] == "yes"
      circle_radialGradient_id = set_random_id()
      circle_radialGradient = get_attributes(SVG,"radialGradient")
      circle_radialGradient["id"] = circle_radialGradient_id

      circle_grad_content = GradientContent()
      circle_grad_content.stop1 = svg_stop(Dict(
        "stop-color"=>p["circle_color1"],
        "offset"=>p["circle_offset1"]
      ))
      circle_grad_content.stop2 = svg_stop(Dict(
        "stop-color"=>p["circle_color2"],
        "offset"=>p["circle_offset2"]
      ))

      defs = string(defs,radialGradient_template(circle_radialGradient,circle_grad_content))
   end
   # Add all gradients to dev element
   # defs = string(
   #    radialGradient_template(rect_radialGradient,rect_grad_content),
   #    radialGradient_template(circle_radialGradient,circle_grad_content)
   #    )

   # Build CSS Classes

   class_list = Array{String,1}()

   if p["rect_class"] == "yes"
      rect_style_id = set_random_id()
      rect_style = get_attributes(CSS,"rect")
      rect_style["fill"] = set_href(rect_radialGradient_id)
      rect_style["stroke-width"] = p["rect_stroke-width"]
      rect_style["stroke"] = p["rect_stroke"]
      push!(class_list,svg_class("rect",rect_style,rect_style_id))
   end
   if p["circle_class"] == "yes"
      circle_style_id = set_random_id()
      circle_style = get_attributes(CSS,"rect")
      circle_style["fill"] = set_href(circle_radialGradient_id)
      circle_style["stroke-width"] = p["circle_stroke-width"]
      circle_style["stroke"] = p["circle_stroke"]
      push!(class_list,svg_class("circle",circle_style,circle_style_id))
   end

   if p["line_class"] == "yes"
      line_style_id = set_random_id()
      line_style = get_attributes(CSS,"line")
      line_style["stroke-width"] = p["line_stroke-width"]
      line_style["stroke"] = p["line_stroke"]
      push!(class_list,svg_class("line",line_style,line_style_id))
   end

   # Build Elements
   el = string()
   ae = abstract_rectangle(p)
   r = string((parse(Float64,ae.r)/parse(Float64,p["factor"])))

   # Create the rectangle

   ra = get_attributes(SVG,"rect")
   if p["rect_class"] == "yes"
     delete!(ra, "style")
     ra["class"] = rect_style_id
   end
   ra["x"] = ae.points["tl"].x
   ra["y"] = ae.points["tl"].y
   ra["width"] = ae.w
   ra["height"] = ae.h
   el = string(el, svg_rect(ra))

   # Create the circles
   for (k,v) in ae.points
         ar = get_attributes(SVG,"circle")
         if p["circle_class"] == "yes"
           delete!(ar, "style")
           ar["class"] = circle_style_id
         end
         ar["cx"] = v.x
         ar["cy"] = v.y
         ar["r"] = r
         el = string(el, svg_circle(ar))
   end

   # create cross lines
   lv = get_attributes(SVG,"line")
   lh = get_attributes(SVG,"line")

   lv["x1"] = ae.points["mt"].x
   lv["y1"] = ae.points["mt"].y
   lv["x2"] = ae.points["mb"].x
   lv["y2"] = ae.points["mb"].y

   lh["x1"] = ae.points["mr"].x
   lh["y1"] = ae.points["mr"].y
   lh["x2"] = ae.points["ml"].x
   lh["y2"] = ae.points["ml"].y

   if p["line_class"] == "yes"
     delete!(lv, "style")
     lv["class"] = line_style_id
     delete!(lh, "style")
     lh["class"] = line_style_id
   end

   el = string(el,svg_line(lv),svg_line(lh))

   ga = get_attributes(SVG,"g")
   ga["type"] = @funcName

   content = svg_g(ga,el)

   # Assemble the document

   doc = SvgContent()

   if !isempty(class_list) doc.style = svg_style(Dict(),svg_css(class_list)) end
   if length(defs) > 0 doc.defs = svg_defs(Dict("id"=>"defs"),defs) end

   # doc.style = svg_style(Dict(),svg_css(class_list))
   # doc.defs = svg_defs(Dict("id"=>"defs"),defs)

   doc.main = svg_g(Dict("id"=>"main"),content)

   svg =  get_attributes(SVG,"svg")
   svg["viewBox"] = join_str([p["x"],p["y"],p["w"],p["h"]])

   return svg_document(svg,doc)

end