text:
using DashBootstrapComponents

progress = html_div([
    dbc_progress(value = 50, style = Dict("height" => "1px"), className = "mb-3"),
    dbc_progress(value = 50, style = Dict("height" => "30px")),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
function get_mesh_2d(i)
    coordinates_min = (-1, -1)
    coordinates_max = (1, 1)
    mesh1_2d = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=5,
                  n_cells_max=10_000
                  )

    coordinates_min = (-1, -1)
    coordinates_max = (1, 1)
    mesh2_2d = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=6,
                  n_cells_max=10_000
                  )
    coordinates_min = (-1, -1)
    coordinates_max = (1, 1)
    refinement_patches = (
    (type = "box", coordinates_min = (-1.0, -1.0), coordinates_max = (-0.5, -0.5)),
    (type = "box", coordinates_min = (-0.5, -0.5), coordinates_max = (0.0, 0.0)),
    (type = "box", coordinates_min = (0.0, 0.0), coordinates_max = (0.5, 0.5)),
    (type = "box", coordinates_min = (0.5, 0.5), coordinates_max = (1.0, 1.0)),
  )
    mesh3_2d = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=4,
                  refinement_patches=refinement_patches,
                  n_cells_max=10_000
                  )

    coordinates_min = (-1, -1)
    coordinates_max = (1, 1)
    refinement_patches = (
    (type = "box", coordinates_min = (-1.0, -1.0), coordinates_max = (-0.5, -0.5)),
    (type = "box", coordinates_min = (-0.5, -0.5), coordinates_max = (0.0, 0.0)),
    (type = "box", coordinates_min = (0.0, 0.0), coordinates_max = (0.5, 0.5)),
    (type = "box", coordinates_min = (0.5, 0.5), coordinates_max = (1.0, 1.0)),
    (type = "box", coordinates_min = (-1.0, -1.0), coordinates_max = (-0.75, -0.75)),
    (type = "box", coordinates_min = (-0.75, -0.75), coordinates_max = (-0.5, -0.5)),
    (type = "box", coordinates_min = (-0.5, -0.5), coordinates_max = (-0.25, -0.25)),
    (type = "box", coordinates_min = (-0.25, -0.25), coordinates_max = (0.0, 0.0)),
    (type = "box", coordinates_min = (0.0, 0.0), coordinates_max = (0.25, 0.25)),
    (type = "box", coordinates_min = (0.25, 0.25), coordinates_max = (0.5, 0.5)),
    (type = "box", coordinates_min = (0.5, 0.5), coordinates_max = (0.75, 0.75)),
    (type = "box", coordinates_min = (0.75, 0.75), coordinates_max = (1.0, 1.0)),
  )
    mesh4_2d = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=5,
                  refinement_patches=refinement_patches,
                  n_cells_max=10_000
                  )

    coordinates_min = (-2, -2)
    coordinates_max = (2, 2)
    refinement_patches = (
    (type = "box", coordinates_min = (-1.3, -1.3), coordinates_max = (-0.8, -0.8)),
    (type = "box", coordinates_min = (-1.1, -1.1), coordinates_max = (-0.9, -0.9)),
    (type = "box", coordinates_min = (0.8, 0.8), coordinates_max = (1.2, 1.2)),
    (type = "box", coordinates_min = (1.2, 1.2), coordinates_max = (1.8, 1.8)),
  )
    mesh5_2d = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=5,
                  refinement_patches=refinement_patches,
                  n_cells_max=10_000
                  )

  if i == 1
    return [mesh1_2d]
  elseif i == 2
    return[mesh2_2d]
  elseif i == 3
    return [mesh3_2d]
  elseif i == 4
    return [mesh4_2d]
  elseif i == 5
    return [mesh5_2d]
  elseif i == 6
    return [mesh1_2d,mesh2_2d,mesh3_2d,mesh4_2d,mesh5_2d]
  else
    error("Input not defined. Should be an Integer between 1-6")
  end

end

ext: jl


text:
#############################################################################
# relative_entropy_epicone.jl
# Constrains τ to
#   τ ⪰ e' * X^{1/2} * logm(X^{1/2}*Y^{-1}*X^{1/2}) * X^{1/2} * e
#
# This function implements the semidefinite programming approximation given in
# the reference below.  Parameters m and k control the accuracy of this
# approximation: m is the number of quadrature nodes to use and k the number
# of square-roots to take. See reference for more details.
#
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#
#REFERENCE
#   Ported from CVXQUAD which is based on the paper: "Semidefinite
#   approximations of matrix logarithm" by Hamza Fawzi, James Saunderson and
#   Pablo A. Parrilo (arXiv:1705.00812)
#############################################################################

struct RelativeEntropyEpiCone
    X::AbstractExpr
    Y::AbstractExpr
    m::Integer
    k::Integer
    e::AbstractMatrix
    size::Tuple{Int,Int}

    function RelativeEntropyEpiCone(
        X::AbstractExpr,
        Y::AbstractExpr,
        m::Integer = 3,
        k::Integer = 3,
        e::AbstractArray = Matrix(1.0 * I, size(X)),
    )
        if size(X) != size(Y)
            throw(DimensionMismatch("X and Y must be the same size"))
        end
        n = size(X)[1]
        if size(X) != (n, n)
            throw(DimensionMismatch("X and Y must be square"))
        end
        if length(size(e)) == 1
            e = reshape(e, (size(e)[1], 1))
        end
        erows, ecols = size(e)
        if ndims(e) != 2 || erows != n
            throw(DimensionMismatch("e matrix must have n rows"))
        end
        return new(X, Y, m, k, e, (ecols, ecols))
    end

    function RelativeEntropyEpiCone(
        X::Value,
        Y::AbstractExpr,
        m::Integer = 3,
        k::Integer = 3,
        e::AbstractArray = Matrix(1.0 * I, size(X)),
    )
        return RelativeEntropyEpiCone(Constant(X), Y, m, k, e)
    end
    function RelativeEntropyEpiCone(
        X::AbstractExpr,
        Y::Value,
        m::Integer = 3,
        k::Integer = 3,
        e::AbstractArray = Matrix(1.0 * I, size(X)),
    )
        return RelativeEntropyEpiCone(X, Constant(Y), m, k, e)
    end
    function RelativeEntropyEpiCone(
        X::Value,
        Y::Value,
        m::Integer = 3,
        k::Integer = 3,
        e::AbstractArray = Matrix(1.0 * I, size(X)),
    )
        return RelativeEntropyEpiCone(Constant(X), Constant(Y), m, k, e)
    end
end

struct RelativeEntropyEpiConeConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    τ::AbstractExpr
    cone::RelativeEntropyEpiCone

    function RelativeEntropyEpiConeConstraint(
        τ::AbstractExpr,
        cone::RelativeEntropyEpiCone,
    )
        if size(τ) != cone.size
            throw(DimensionMismatch("τ must be size $(cone.size)"))
        end
        id_hash = hash((
            cone.X,
            cone.Y,
            cone.m,
            cone.k,
            cone.e,
            :RelativeEntropyEpiCone,
        ))
        return new(:RelativeEntropyEpiCone, id_hash, τ, cone)
    end

    function RelativeEntropyEpiConeConstraint(
        τ::Value,
        cone::RelativeEntropyEpiCone,
    )
        return RelativeEntropyEpiConeConstraint(Constant(τ), cone)
    end
end

in(τ, cone::RelativeEntropyEpiCone) = RelativeEntropyEpiConeConstraint(τ, cone)

function AbstractTrees.children(constraint::RelativeEntropyEpiConeConstraint)
    return (constraint.τ, constraint.cone.X, constraint.cone.Y)
end

# This negative relative entropy function is matrix convex (arxiv:1705.00812).
# So if X and Y are convex sets, then τ ⪰ -D_op(X || Y) will be a convex set.
function vexity(constraint::RelativeEntropyEpiConeConstraint)
    X = vexity(constraint.cone.X)
    Y = vexity(constraint.cone.Y)
    τ = vexity(constraint.τ)

    # NOTE: can't say X == NotDcp() because the NotDcp constructor prints a warning message.
    if typeof(X) == ConcaveVexity || typeof(X) == NotDcp
        return NotDcp()
    end
    if typeof(Y) == ConcaveVexity || typeof(Y) == NotDcp
        return NotDcp()
    end
    # Copied from vexity(c::GtConstraint)
    vex = ConvexVexity() + (-τ)
    if vex == ConcaveVexity()
        vex = NotDcp()
    end
    return vex
end

function glquad(m)
    # Compute Gauss-Legendre quadrature nodes and weights on [0, 1]
    # Code below is from [Trefethen, "Is Gauss quadrature better than
    # Clenshaw-Curtis?", SIAM Review 2008] and computes the weights and
    # nodes on [-1, 1].
    beta = [0.5 ./ sqrt(1 - (2 * i)^-2) for i in 1:(m-1)] # 3-term recurrence coeffs
    T = diagm(1 => beta, -1 => beta) # Jacobi matrix
    s, V = eigen(T)
    w = 2 * V[1, :] .^ 2 # weights
    # Translate and scale to [0, 1]
    s = (s .+ 1) / 2
    w = w' / 2
    return s, w
end

function conic_form!(
    constraint::RelativeEntropyEpiConeConstraint,
    unique_conic_forms::UniqueConicForms,
)
    if !has_conic_form(unique_conic_forms, constraint)
        X = constraint.cone.X
        Y = constraint.cone.Y
        m = constraint.cone.m
        k = constraint.cone.k
        e = constraint.cone.e
        τ = constraint.τ
        n = size(X)[1]
        r = size(e)[2]

        s, w = glquad(m)

        is_complex =
            sign(X) == ComplexSign() ||
            sign(Y) == ComplexSign() ||
            sign(Constant(e)) == ComplexSign()
        if is_complex
            Z = ComplexVariable(n, n)
            T = [ComplexVariable(r, r) for i in 1:m]
        else
            Z = Variable(n, n)
            T = [Variable(r, r) for i in 1:m]
        end

        conic_form!(
            Z in GeomMeanHypoCone(X, Y, 1 // (2^k), false),
            unique_conic_forms,
        )

        for ii in 1:m
            # Note that we are dividing by w here because it is easier
            # to do this than to do sum w_i T(:,...,:,ii) later (cf. line that
            # involves τ)

            conic_form!(
                [
                    e'*X*e-s[ii]*T[ii]/w[ii] e'*X
                    X*e (1-s[ii])*X+s[ii]*Z
                ] ⪰ 0,
                unique_conic_forms,
            )
        end

        conic_form!((2^k) * sum(T) + τ ⪰ 0, unique_conic_forms)

        cache_conic_form!(
            unique_conic_forms,
            constraint,
            Array{Convex.ConicConstr,1}(),
        )
    end
    return get_conic_form(unique_conic_forms, constraint)
end


text:
############################################################################################
# types

############################################################################################
@testset "Types - Float" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = Float16(1.0)
            ReConstructor(val)
            reconstruct = ReConstructor(flatdefault, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == typeof(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test x_unflatAD isa eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test x_unflatAD2 isa eltype(x_flatAD)
################################################################################
# Constraint ~ Make Distribution as standard constraints are separately tested
            constraint_float = Distributions.Gamma(2,2)
            ReConstructor(constraint_float, val)
            reconstruct = ReConstructor(flatdefault, constraint_float, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == typeof(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test x_unflatAD isa eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test x_unflatAD2 isa eltype(x_flatAD)
################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint_float, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint_float, val)

            val_uncon = ModelWrappers.unconstrain(transformer, val)
            val_con = ModelWrappers.constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(x_inversetransformed)
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
        end
    end
end

@testset "Types - Float - Automatic Differentiation" begin
    val = Float32(1.0)
    constraint = Distributions.Gamma(2,2)
    reconstruct = ReConstructor(constraint, val)
    val_flat = flatten(reconstruct, val)

    check_AD = check_AD_closure(constraint, val)
    check_AD(val_flat)
    grad_mod_fd = ForwardDiff.gradient(check_AD, val_flat)
    grad_mod_rd = ReverseDiff.gradient(check_AD, val_flat)
    grad_mod_zy = Zygote.gradient(check_AD, val_flat)[1]
    @test sum(abs.(grad_mod_fd - grad_mod_rd)) ≈ 0 atol = _TOL
    @test sum(abs.(grad_mod_fd - grad_mod_zy)) ≈ 0 atol = _TOL

    grad_mod_fd = ForwardDiff.gradient(check_AD, flatten(reconstruct, val))
    grad_mod_rd = ReverseDiff.gradient(check_AD, flatten(reconstruct, val))
    grad_mod_zy = Zygote.gradient(check_AD, flatten(reconstruct, val))[1]
end

############################################################################################
@testset "Types - Vector Float" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = Float16.([1., 2.])
            ReConstructor(val)
            reconstruct = ReConstructor(flatdefault, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test eltype(x_unflatAD) == eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test eltype(x_unflatAD2) == eltype(x_flatAD)
################################################################################
# Constraint ~ Make Distribution as standard constraints are separately tested
            constraint = Distributions.MvNormal([1.,1.])
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test eltype(x_unflatAD) == eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test eltype(x_unflatAD2) == eltype(x_flatAD)

################################################################################
# Check if errors if flattened parameter not correct size
            @test_throws ArgumentError flatten(reconstruct, zeros(length(val)+1))
            @test_throws ArgumentError flatten(reconstruct, zeros(length(val)-1))
            @test_throws ArgumentError flattenAD(reconstruct, zeros(length(val)+1))
            @test_throws ArgumentError flattenAD(reconstruct, zeros(length(val)-1))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)-1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)-1))

################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(x_inversetransformed)
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            _check(_RNG, constraint, val)
        end
    end
end

@testset "Types - Vector Float - Automatic Differentiation" begin
    val = [1., 2.]
    constraint = Distributions.MvNormal([1., 1.])
    reconstruct = ReConstructor(constraint, val)
    val_flat = flatten(reconstruct, val)

    check_AD = check_AD_closure(constraint, val)
    check_AD(val_flat)
    grad_mod_fd = ForwardDiff.gradient(check_AD, val_flat)
    grad_mod_rd = ReverseDiff.gradient(check_AD, val_flat)
    grad_mod_zy = Zygote.gradient(check_AD, val_flat)[1]
    @test sum(abs.(grad_mod_fd - grad_mod_rd)) ≈ 0 atol = _TOL
#    @test sum(abs.(grad_mod_fd - grad_mod_zy)) ≈ 0 atol = _TOL

    grad_mod_fd = ForwardDiff.gradient(check_AD, flatten(reconstruct, val))
    grad_mod_rd = ReverseDiff.gradient(check_AD, flatten(reconstruct, val))
    grad_mod_zy = Zygote.gradient(check_AD, flatten(reconstruct, val))[1]
end

############################################################################################
@testset "Types - Array Float" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = Float16.([1. 0.3 ; .3 1.0])
            ReConstructor(val)
            reconstruct = ReConstructor(flatdefault, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test eltype(x_unflatAD) == eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test eltype(x_unflatAD2) == eltype(x_flatAD)

################################################################################
# Check if errors if flattened parameter not correct size
            @test_throws ArgumentError flatten(reconstruct, zeros(3,3))
            @test_throws ArgumentError flatten(reconstruct, zeros(1,1))
            @test_throws ArgumentError flattenAD(reconstruct, zeros(3,3))
            @test_throws ArgumentError flattenAD(reconstruct,  zeros(1,1))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)-1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)-1))

################################################################################
# Constraint ~ Make Distribution as standard constraints are separately tested
            constraint = Unconstrained()
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test eltype(x_unflatAD) == eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test eltype(x_unflatAD2) == eltype(x_flatAD)
################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(x_inversetransformed)
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            _check(_RNG, constraint, val)
        end
    end
end

@testset "Types - Array Float - Automatic Differentiation" begin
    val = [1. 0.3 ; .3 1.0]
    constraint = Distributions.Distributions.InverseWishart(10., [1. 0. ; 0. 1.])
    reconstruct = ReConstructor(constraint, val)
    val_flat = flatten(reconstruct, val)

    check_AD = check_AD_closure(constraint, val)
    check_AD(val_flat)
    grad_mod_fd = ForwardDiff.gradient(check_AD, val_flat)
    #grad_mod_rd = ReverseDiff.gradient(check_AD, val_flat)
    grad_mod_zy = Zygote.gradient(check_AD, val_flat)[1]
    #@test sum(abs.(grad_mod_fd - grad_mod_rd)) ≈ 0 atol = _TOL
    @test sum(abs.(grad_mod_fd - grad_mod_zy)) ≈ 0 atol = _TOL

    grad_mod_fd = ForwardDiff.gradient(check_AD, flatten(reconstruct, val))
    #grad_mod_rd = ReverseDiff.gradient(check_AD, flatten(reconstruct, val))
    grad_mod_zy = Zygote.gradient(check_AD, flatten(reconstruct, val))[1]
end

############################################################################################
@testset "Types - Integer" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = Int16(1.0)
            ReConstructor(val)
            reconstruct = ReConstructor(flatdefault, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            if flattentype isa FlattenAll
                @test eltype(x_flat) == output
            end
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test x_unflatAD isa eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test x_unflatAD2 isa eltype(x_flatAD)
        end
    end
end

############################################################################################
@testset "Types - Array Integer" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = Int16.([1 2 ; 3 4])
            ReConstructor(val)
            reconstruct = ReConstructor(flatdefault, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            if flattentype isa FlattenAll
                @test eltype(x_flat) == output
            end
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test eltype(x_unflatAD) == eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test eltype(x_unflatAD2) == eltype(x_flatAD)
        end
    end
end


text:
# FST passes/transforms

function flattenable(kind::Tokens.Kind)
    kind === Tokens.AND && return true
    kind === Tokens.OR && return true
    kind === Tokens.LAZY_AND && return true
    kind === Tokens.LAZY_OR && return true
    kind === Tokens.RPIPE && return true
    return false
end
flattenable(::Nothing) = false

"""
Flattens a binary operation call tree if the operation repeats 2 or more times.
"a && b && c" will be transformed while "a && b" will not.
"""
function flatten_binaryopcall(fst::FST; top = true)
    nodes = FST[]
    kind = op_kind(fst)

    lhs = fst[1]
    rhs = fst[end]
    lhs_kind = op_kind(lhs)
    rhs_kind = op_kind(rhs)
    lhs_same_op = lhs_kind === kind
    rhs_same_op = rhs_kind === kind
    idx = findlast(n -> n.typ === PLACEHOLDER, fst.nodes)

    if (top && !lhs_same_op && !rhs_same_op) || idx === nothing
        return nodes
    end

    if lhs_same_op
        push!(nodes, flatten_binaryopcall(lhs, top = false)...)
    else
        flatten_fst!(lhs)
        push!(nodes, lhs)
    end
    # everything except the indentation placeholder
    push!(nodes, fst[2:idx-1]...)

    if rhs_same_op
        push!(nodes, flatten_binaryopcall(rhs, top = false)...)
    else
        flatten_fst!(rhs)
        push!(nodes, rhs)
    end

    return nodes
end

function flatten_conditionalopcall(fst::FST)
    nodes = FST[]
    for n in fst.nodes
        if n.typ === Conditional
            push!(nodes, flatten_conditionalopcall(n)...)
        else
            push!(nodes, n)
        end
    end
    return nodes
end

function flatten_fst!(fst::FST)
    is_leaf(fst) && return
    for n in fst.nodes
        if is_leaf(n)
            continue
        elseif n.typ === Binary && flattenable(op_kind(n))
            # possibly convert Binary to Chain
            nnodes = flatten_binaryopcall(n)
            if length(nnodes) > 0
                n.typ = Chain
                n.nodes = nnodes
            else
                flatten_fst!(n)
            end
        else
            flatten_fst!(n)
        end
    end
end

"""
    pipe_to_function_call_pass!(fst::FST)

Rewrites `x |> f` to `f(x)`.
"""
function pipe_to_function_call_pass!(fst::FST)
    is_leaf(fst) && return

    if op_kind(fst) === Tokens.RPIPE
        fst.nodes = pipe_to_function_call(fst)
        fst.typ = Call
        return
    end

    for n in fst.nodes
        if is_leaf(n)
            continue
        elseif op_kind(n) === Tokens.RPIPE
            n.nodes = pipe_to_function_call(n)
            n.typ = Call
        else
            pipe_to_function_call_pass!(n)
        end
    end
end

function pipe_to_function_call(fst::FST)
    nodes = FST[]
    arg2 = fst[end]
    push!(nodes, arg2)
    paren = FST(PUNCTUATION, -1, arg2.endline, arg2.endline, "(")
    push!(nodes, paren)
    pipe_to_function_call_pass!(fst[1])
    arg1 = fst[1]
    push!(nodes, arg1)
    paren = FST(PUNCTUATION, -1, arg1.endline, arg1.endline, ")")
    push!(nodes, paren)
    return nodes
end

function import_to_usings(fst::FST, s::State)
    findfirst(n -> is_colon(n) || n.typ === As, fst.nodes) === nothing || return FST[]
    findfirst(n -> n.typ === PUNCTUATION && n.val == ".", fst[3].nodes) === nothing ||
        return FST[]

    usings = FST[]
    idxs = findall(n -> !is_leaf(n), fst.nodes)

    for i in idxs
        n = fst[i]
        sl = n.startline
        el = n.endline
        use = FST(Using, fst.indent)
        use.startline = n.startline
        use.endline = n.endline

        add_node!(use, FST(KEYWORD, -1, sl, el, "using"), s)
        add_node!(use, Whitespace(1), s)

        add_node!(use, n, s, join_lines = true)
        colon = FST(OPERATOR, -1, sl, el, ":")
        colon.opmeta = OpMeta(Tokens.COLON, false)
        add_node!(use, colon, s, join_lines = true)
        add_node!(use, Whitespace(1), s)
        add_node!(use, FST(IDENTIFIER, -1, sl, el, n[end].val), s, join_lines = true)

        push!(usings, use)
    end
    return usings
end

"""
    annotate_typefields_with_any!(fst::FST, s::State)

Annotates fields in a type definitions with `::Any` if
no type annotation is provided.
"""
function annotate_typefields_with_any!(fst::FST, s::State)
    is_leaf(fst) && return
    for (i, n) in enumerate(fst.nodes)
        if n.typ === IDENTIFIER
            nn = FST(Binary, n.indent)
            nn.startline = n.startline
            nn.endline = n.endline
            add_node!(nn, n, s)
            line_offset = n.line_offset + length(n)
            op = FST(OPERATOR, line_offset, n.startline, n.endline, "::")
            op.opmeta = OpMeta(Tokens.DECLARATION, false)
            add_node!(nn, op, s, join_lines = true)
            line_offset += 2
            add_node!(
                nn,
                FST(IDENTIFIER, line_offset, n.startline, n.endline, "Any"),
                s,
                join_lines = true,
            )
            fst[i] = nn
        else
            continue
        end
    end
end

"""
    short_to_long_function_def!(fst::FST, s::State)

Transforms a *short* function definition

```julia
f(arg1, arg2) = body
```

to a *long* function definition

```julia
function f(arg2, arg2)
    body
end
```
"""
function short_to_long_function_def!(fst::FST, s::State)
    (fst[1].typ !== Call && fst[1].typ !== Where) && return

    # 3 cases
    #
    # case 1
    #   func(a) = body
    #
    # case 2
    #   func(a::T) where T = body
    #
    # case 3
    #   func(a::T)::R where T = body
    funcdef = FST(FunctionN, fst.indent)
    # function
    kw = FST(KEYWORD, -1, fst[1].startline, fst[1].endline, "function")
    add_node!(funcdef, kw, s)
    add_node!(funcdef, Whitespace(1), s)

    # func(a)
    # OR
    # func(a) where T
    add_node!(funcdef, fst[1], s, join_lines = true)

    # body

    s.opts.always_use_return && prepend_return!(fst[end], s)
    if fst[end].typ === Block
        add_node!(funcdef, fst[end], s, max_padding = s.opts.indent)
    else
        # ```
        # function
        #     body
        # end
        # ```
        #
        # `body` is parsed wrapped block node. Wrapping it in
        # a `Block` node ensures the indent is correct.
        bl = FST(Block, fst[end].indent)
        add_node!(bl, fst[end], s)
        add_node!(funcdef, bl, s, max_padding = s.opts.indent)
    end
    add_indent!(funcdef[end], s, s.opts.indent)

    # end
    kw = FST(KEYWORD, -1, fst[end].startline, fst[end].endline, "end")
    add_node!(funcdef, kw, s)

    fst.typ = funcdef.typ
    fst.nodes = funcdef.nodes
    fst.len = funcdef.len
end

"""
    binaryop_to_whereop(fst::FST, s::State)

Handles the case of a function def defined
as:

```julia
foo(a::A)::R where A = body
```

In this case instead of it being parsed as (1):

```
Binary
 - Where
 - OP
 - RHS
```

It's parsed as (2):

```
Binary
 - Binary
  - LHS
  - OP
  - Where
   - R
   - ...
 - OP
 - RHS
```

(1) is preferrable since it's the same parsed result as:

```julia
foo(a::A) where A = body
```

This transformation converts (2) to (1).

ref https://github.com/julia-vscode/CSTParser.jl/issues/93
"""
function binaryop_to_whereop!(fst::FST, s::State)
    # transform fst[1] to a Where
    oldbinop = fst[1]
    oldwhereop = fst[1][end]
    binop = FST(Binary, fst[1].indent)

    # foo(a::A)
    add_node!(binop, oldbinop[1], s)
    # foo(a::A)::
    add_node!(binop, oldbinop[2], s, join_lines = true)
    # foo(a::A)::R
    add_node!(binop, oldwhereop[1], s, join_lines = true)

    whereop = FST(Where, fst[1].indent)
    add_node!(whereop, binop, s)

    # "foo(a::A)::R where A"
    for n in oldwhereop[2:end]
        add_node!(whereop, n, s, join_lines = true)
    end

    fst[1] = whereop
end

"""
    prepend_return!(fst::FST, s::State)

Prepends `return` to the last expression of a block if applicable.

```julia
function foo()
    a = 2 * 3
    a / 3
end
```

to

```julia
function foo()
    a = 2 * 3
    return a / 3
end
```
"""
function prepend_return!(fst::FST, s::State)
    fst.typ === Block || return
    ln = fst[end]
    is_block(ln) && return
    ln.typ === Return && return
    ln.typ === MacroCall && return
    ln.typ === MacroBlock && return
    ln.typ === MacroStr && return
    if length(fst.nodes) > 2 && (fst[end-2].typ === MacroStr || is_macrodoc(fst[end-2]))
        # The last node is has a docstring prior to it so a return should not be prepended
        # fst[end-1] is a newline
        return
    end

    ret = FST(Return, fst.indent)
    kw = FST(KEYWORD, -1, fst[end].startline, fst[end].endline, "return")
    add_node!(ret, kw, s)
    add_node!(ret, Whitespace(1), s)
    add_node!(ret, ln, s, join_lines = true)
    fst[end] = ret
    return
end

"""
    move_at_sign_to_the_end(fst::FST, s::State)

NOTE: Assumes `fst` is the caller name of a macrocall such as
`@macro` or `Module.@macro`.

Moves `@` to the last indentifier.

Example:

```julia
@Module.macro
```

to

```julia
Module.@macro
```
"""
function move_at_sign_to_the_end(fst::FST, s::State)
    t = FST[]
    f = (t) -> (n, s) -> begin
        if is_macrocall(n) || (n.typ === Quotenode && !is_leaf(n[1]))
            # 1. Do not move "@" in nested macro calls
            # 2. Do not move "@" if in the middle of a chain, i.e. "a.@b.c"
            # since it's semantically different to "@a.b.c" and "a.b.@c"
            push!(t, n)
            return false
        elseif is_leaf(n)
            push!(t, n)
        end
    end
    walk(f(t), fst, s)

    macroname = FST(Macroname, fst.indent)
    for (i, n) in enumerate(t)
        if n.val == "@"
            continue
        elseif n.typ === IDENTIFIER && i < length(t) && n.val[1] == '@'
            n.val = n.val[2:end]
            n.len -= 1
            add_node!(macroname, n, s, join_lines = true)
        elseif i < length(t) || n.typ == Quotenode
            add_node!(macroname, n, s, join_lines = true)
        else
            if n.typ === IDENTIFIER && n.val[1] != '@'
                n.val = "@" * n.val
                n.len += 1
                add_node!(macroname, n, s, join_lines = true)
            else
                add_node!(macroname, n, s, join_lines = true)
            end
        end
    end

    return macroname
end

function conditional_to_if_block!(fst::FST, s::State; top = true)
    t = FST(If, fst.indent)
    kw = FST(KEYWORD, -1, fst.startline, fst.startline, top ? "if" : "elseif")
    add_node!(t, kw, s, max_padding = 0)
    add_node!(t, Whitespace(1), s, join_lines = true)
    add_node!(t, fst[1], s, join_lines = true)

    idx1 = findfirst(n -> n.typ === OPERATOR && n.val == "?", fst.nodes)
    idx2 = findfirst(n -> n.typ === OPERATOR && n.val == ":", fst.nodes)

    block1 = FST(Block, fst.indent + s.opts.indent)
    for n in fst.nodes[idx1+1:idx2-1]
        if n.typ === PLACEHOLDER ||
           n.typ === WHITESPACE ||
           n.typ === NEWLINE ||
           is_comment(n)
            continue
        end
        add_node!(block1, n, s)
    end
    add_node!(t, block1, s, max_padding = s.opts.indent)

    block2 = FST(Block, fst.indent)
    padding = 0
    if fst[end].typ === Conditional
        conditional_to_if_block!(fst[end], s, top = false)
    else
        block2.indent += s.opts.indent
        padding = s.opts.indent
        kw = FST(KEYWORD, -1, -1, -1, "else")
        add_node!(t, kw, s, max_padding = 0)
    end
    add_node!(block2, fst[end], s)
    add_node!(t, block2, s, max_padding = 0)

    if top
        kw = FST(KEYWORD, -1, -1, -1, "end")
        add_node!(t, kw, s, max_padding = 0)
    end

    fst.typ = t.typ
    fst.nodes = t.nodes
    fst.len = t.len

    return nothing
end

"""
    separate_kwargs_with_semicolon!(fst::FST)

Ensures keyword arguments are separated with a ";".

### Examples

Replace "," with ";".

```julia
a = f(x, y = 3)

->

a = f(x; y = 3)
```

Move ";" to the prior to the first positional argument.

```julia
a = f(x = 1; y = 2)

->

a = f(; x = 1, y = 2)
```
"""
function separate_kwargs_with_semicolon!(fst::FST)
    kw_idx = findfirst(n -> n.typ === Kw, fst.nodes)
    kw_idx === nothing && return
    sc_idx = findfirst(n -> n.typ === SEMICOLON, fst.nodes)
    # first "," prior to a kwarg
    comma_idx = findlast(is_comma, fst.nodes[1:kw_idx-1])
    ph_idx = findlast(n -> n.typ === PLACEHOLDER, fst.nodes[1:kw_idx-1])
    # @info "" kw_idx sc_idx comma_idx ph_idx

    if sc_idx !== nothing && sc_idx > kw_idx
        # move ; prior to first kwarg
        fst[sc_idx].val = ","
        fst[sc_idx].typ = PUNCTUATION
        if comma_idx === nothing
            if ph_idx !== nothing
                fst[ph_idx] = Placeholder(1)
                insert!(fst, ph_idx, Semicolon())
            else
                insert!(fst, kw_idx, Placeholder(1))
                insert!(fst, kw_idx, Semicolon())
            end
        end
    elseif sc_idx === nothing && comma_idx === nothing
        if ph_idx !== nothing
            fst[ph_idx] = Placeholder(1)
            insert!(fst, ph_idx, Semicolon())
        else
            insert!(fst, kw_idx, Placeholder(1))
            insert!(fst, kw_idx, Semicolon())
        end
    elseif sc_idx === nothing
        fst[comma_idx].val = ";"
        fst[comma_idx].typ = SEMICOLON
    end
    return
end


text:
using Test

import Tar
import Pkg.GitTools: tree_hash

function gen_file(file::String, size::Int)
    open(file, write=true) do io
        for i = 1:size; write(io, i % UInt8); end
    end
end

function make_test_tarball()
    root = mktempdir()
    let i = 0, paths = String[]
        for a in [0, 10, 154, 155, 156, 255]
            dir = joinpath(root, "d"^a)
            push!(paths, normpath("..", dir))
            a > 0 && mkdir(dir)
            for b in [10, 99, 100, 101, 255]
                for s in [0, 511, 512, 513, 1000]
                    f = rpad("$s-", b, "f")
                    x = rpad("$s-", b, "x")
                    file = joinpath(dir, f)
                    exec = joinpath(dir, x)
                    gen_file(file, s)
                    gen_file(exec, s)
                    chmod(file, 0o644)
                    chmod(exec, 0o755)
                    push!(paths, file)
                    push!(paths, exec)
                end
                dir′ = joinpath(dir, "s"^b)
                mkpath(dir′)
                push!(paths, dir′)
                target = relpath(paths[i += 1], root)
                link = joinpath(dir, "l"^b)
                symlink(target, link)
                push!(paths, link)
                broken = joinpath(dir, "b"^b)
                if target != "."
                    symlink(chop(target), broken)
                    push!(paths, broken)
                end
            end
        end
    end
    hash = tree_hash(root)
    tarball = Tar.create(root)
    rm(root, force=true, recursive=true)
    return tarball, hash
end

@testset "empty tarball" begin
    dir = mktempdir()
    tarball = Tar.create(dir)
    rm(dir, recursive=true)
    @test Tar.list(tarball) == [Tar.Header(".", :directory, 0o755, 0, "")]
    dir = Tar.extract(tarball)
    @test isempty(readdir(dir))
    rm(dir, recursive=true)
    open(tarball, append=true) do io
        write(io, zeros(UInt8, 512))
    end
    dir = Tar.extract(tarball)
    @test isempty(readdir(dir))
    rm(dir, recursive=true)
end

function check_tree_hash(hash::Vector{UInt8}, root::AbstractString)
    @test tree_hash(root) == hash
    rm(root, recursive=true)
end

@testset "test tarball" begin
    tarball, hash = make_test_tarball()
    @testset "Tar.list & check properties" begin
        headers = Tar.list(tarball)
        @test issorted(headers, by = hdr -> hdr.path)
        for hdr in headers
            @test !isempty(hdr.path)
            @test hdr.type in (:file, :directory, :symlink)
            if hdr.type == :file
                @test hdr.path[end] != '/'
                @test hdr.mode in (0o644, 0o755)
                @test isempty(hdr.link)
            elseif hdr.type == :directory
                @test hdr.path[end] == '/'
                @test hdr.mode == 0o755
                @test hdr.size == 0
                @test isempty(hdr.link)
            elseif hdr.type == :symlink
                @test hdr.path[end] != '/'
                @test hdr.mode == 0o755
                @test hdr.size == 0
                @test !isempty(hdr.link)
            end
        end
        @testset "Tar.list from IO, process, pipeline" begin
            @test headers == open(Tar.list, tarball)
            @test headers == open(Tar.list, `cat $tarball`)
            @test headers == open(Tar.list, pipeline(`bzip2 -c -9 $tarball`, `bzcat`))
        end
    end
    @testset "extract with `tar` command" begin
        root = mktempdir()
        run(`tar -C $root -xf $tarball`)
        check_tree_hash(hash, root)
    end
    @testset "Tar.extract" begin
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
    @testset "Tar.extract from IO, process, pipeline" begin
        root = open(Tar.extract, tarball)
        check_tree_hash(hash, root)
        root = open(Tar.extract, `cat $tarball`)
        check_tree_hash(hash, root)
        root = open(Tar.extract, pipeline(`bzip2 -c -9 $tarball`, `bzcat`))
        check_tree_hash(hash, root)
    end
    open(tarball, append=true) do io
        write(io, zeros(UInt8, 512))
    end
    @testset "Tar.extract with trailing zeros" begin
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
    rm(tarball)
end

@testset "symlink attacks" begin
    # not dangerous but still not allowed
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("dir", :directory, 0o755, 0, ""))
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "dir"))
    Tar.write_header(io, Tar.Header("link/target", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
    # attempt to write through relative link out of root
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "../target"))
    Tar.write_header(io, Tar.Header("link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
    # attempt to write through absolute link
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "/tmp"))
    Tar.write_header(io, Tar.Header("link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
    # same attack with some obfuscation
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "/tmp"))
    Tar.write_header(io, Tar.Header("./link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
    # same attack with different obfuscation
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "/tmp"))
    Tar.write_header(io, Tar.Header("dir/../link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
end

@testset "symlink overwrite" begin
    # allowable and should work
    tarball₁, io = mktemp()
    Tar.write_header(io, Tar.Header("file", :symlink, 0o755, 0, "/tmp"))
    Tar.write_header(io, Tar.Header("file", :file, 0o644, 0, ""))
    close(io)
    tree₁ = Tar.extract(tarball₁)
    tarball₂, io = mktemp()
    Tar.write_header(io, Tar.Header("file", :file, 0o644, 0, ""))
    close(io)
    tree₂ = Tar.extract(tarball₂)
    @test tree_hash(tree₁) == tree_hash(tree₂)
    rm(tree₁, recursive=true)
    rm(tree₂, recursive=true)
    rm(tarball₁)
    rm(tarball₂)
end


text:
abstract type AbstractPointPose end
#NOTE names should match soft type names with, eg. Pose2 and Point2

# TODO look at replacing with traits

mutable struct Pose2{T<:AbstractFloat} <: AbstractPointPose
	x::T
	y::T
	θ::T
end

mutable struct Point2{T<:AbstractFloat} <: AbstractPointPose
	x::T
	y::T
end

mutable struct Point3{T<:AbstractFloat} <: AbstractPointPose
	x::T
	y::T
	z::T
end

"""
    $(SIGNATURES)
Set the draw transform for the meshcat root object.

Example:
```julia
# default is Z up
# make Z down
setGlobalDrawTransform!(vis, quat=UnitQuaternion(0.0,1.0,0.0,0.0))
```
"""
function setGlobalDrawTransform!(vis::Visualizer;trans=Translation(0.0,0.0,0.0), quat::Rotations.UnitQuaternion=UnitQuaternion(1.0,0.0,0.0,0.0))
	tform = trans ∘ LinearMap(quat)
	settransform!(vis, tform)
end

"""
    $(SIGNATURES)
Basic visualize function for a pose.
"""
function visPose!(vis::Visualizer,
				  tform::AbstractAffineMap,
				  updateonly::Bool=false;
				  scale::Float64=0.2,
				  sphereScale::Float64=0.05,
				  color::RGBA=RGBA(1., 1, 0, 0.5))
	if !updateonly
		setobject!(vis, Triad(scale))
		if sphereScale > 0.0
			sphere = GeometryBasics.HyperSphere(GeometryBasics.Point(0., 0, 0), sphereScale)
			matcolor = MeshPhongMaterial(color=color)
			setobject!(vis[:bal], sphere, matcolor)
		end
	end
	settransform!(vis, tform)
	nothing
end

"""
    $(SIGNATURES)
Basic visualize function for a point.
"""
function visPoint!(vis::Visualizer,
				   tform::AbstractAffineMap,
				   updateonly::Bool=false;
				   scale::Float64=0.1, #unused but needed for uniformity #TODO maybe replace params with dictionary?
				   sphereScale::Float64=0.1,
				   color::RGBA=RGBA(0., 1, 0, 0.5))

	if !updateonly
		sphere = GeometryBasics.HyperSphere(GeometryBasics.Point(0., 0, 0), sphereScale)
        matcolor = MeshPhongMaterial(color=color)
        setobject!(vis, sphere, matcolor)
	end
	settransform!(vis, tform)
	nothing
end

"""
    $(SIGNATURES)
Node Visualization functions.
"""
function visNode!(vis::Visualizer,
	              pose::Pose2,
				  updateonly::Bool=false;
				  zoffset::Float64=0.0,
				  kwargs...)::Nothing
				  # scale::Float64=0.2,
				  # sphereScale::Float64=0.05,
				  # color::RGBA=RGBA(1., 1, 0, 0.5))::Nothing

	tf = Translation(pose.x, pose.y, zoffset) ∘ LinearMap(Rotations.RotZ(pose.θ))
	visPose!(vis, tf, updateonly; kwargs...)
					  # scale=scale,
					  # sphereScale=sphereScale,
					  # color=color)
end

function visNode!(vis::Visualizer,
            	  point::Point2,
				  updateonly::Bool=false;
                  zoffset::Float64=0.0,
				  kwargs...)::Nothing
				  # scale::Float64=0.1, #TODO i tried with kwargs...), but missing something
				  # sphereScale::Float64=0.1,
				  # color::RGBA=RGBA(0., 1, 0, 0.5))::Nothing

	tf = Translation(point.x, point.y, zoffset)
	visPoint!(vis, tf, updateonly; kwargs...)
					   # scale=scale,
					   # sphereScale=sphereScale,
					   # color=color)
end

function visNode!(vis::Visualizer,
            	  point::Point3,
				  updateonly::Bool=false;
				  zoffset::Float64=0.0,
				  kwargs...)::Nothing
				  # scale::Float64=0.1, #unused but needed for uniformity #TODO maybe replace params with dictionary?
				  # sphereScale::Float64=0.1,
				  # color::RGBA=RGBA(0., 1, 0, 0.5))::Nothing

	tf = Translation(point...)
	visPoint!(vis, tf, updateonly; kwargs...)
					   # scale=scale,
					   # sphereScale=sphereScale,
					   # color=color)
end

max_stars_repo_name: Affie/Amphitheatre

max_issues_repo_name: Affie/Amphitheatre


text:
using DelimitedFiles
using Printf

function msmforward(data_list, T, emission, pi_i)
    ndata = length(data_list)
    nstate = length(T[1, :])
    logL = zeros(Float64, ndata)
    alpha_list = []
    factor_list = []
    for idata = 1:ndata
        data = data_list[idata]
        nframe = length(data)
        alpha  = zeros(Float64, (nframe, nstate))
        factor = zeros(Float64, nframe)
        alpha[1, :] = pi_i.*emission[:, data[1]]
        factor[1] = sum(alpha[1, :])
        alpha[1, :] = alpha[1, :]./factor[1]
        for iframe = 2:nframe
            alpha[iframe, :] = sum(alpha[iframe-1, :] .* T, dims=1)' .* emission[:, data[iframe]]
            factor[iframe] = sum(alpha[iframe, :])
            alpha[iframe, :] = alpha[iframe, :]./factor[iframe]
        end
        logL[idata] = sum(log.(factor))
        push!(alpha_list, alpha)
        push!(factor_list, factor)
    end
    logL, alpha_list, factor_list
end

function msmbackward(data_list, factor_list, T, emission, pi_i)
    ndata = length(data_list)
    nstate = length(T[1, :])
    logL = zeros(Float64, ndata)
    beta_list = []
    for idata = 1:ndata
        data   = data_list[idata]
        factor = factor_list[idata]
        nframe = length(data)
        beta   = zeros(Float64, (nframe, nstate))
        beta[nframe, :] .= 1.0
        for iframe = (nframe-1):-1:1
            beta[iframe, :] = sum((T .* (emission[:, data[iframe+1]] .* beta[iframe+1, :])'), dims=2) ./ factor[iframe+1]
        end
        logL[idata] = sum(log.(factor))
        push!(beta_list, beta)
    end
    logL, beta_list
end

function msmbaumwelch(data_list, T0, emission0, pi_i0)
    ## setup
    TOLERANCE = 10.0^(-4)
    check_convergence = Inf64
    count_iteration = 0
    logL_old = 1.0
    #if not isinstance(data_list, list):
    #    data_list = [data_list]
    ndata = length(data_list)
    nobs = length(emission0[1, :])
    nstate = length(T0[1, :])
    T = similar(T0)
    emission = similar(emission0)
    pi_i = similar(pi_i0)
    while check_convergence > TOLERANCE
        ## E-step
        logL, alpha_list, factor_list = msmforward(data_list, T0, emission0, pi_i0)
        print("1"); println(logL)
        logL2, beta_list = msmbackward(data_list, factor_list, T0, emission0, pi_i0)
        print("2"); println(logL2)
        log_alpha_list = []
        for a in alpha_list
            push!(log_alpha_list, log.(a))
        end
        log_beta_list = []
        for b in beta_list
            push!(log_beta_list, log.(b))
        end
        log_T0 = log.(T0)
        log_emission0 = log.(emission0)
        ## M-step
        # pi
        # pi = np.zeros(nstate, dtype=np.float64)
        # log_gamma_list = []
        # for idata in range(ndata):
        #     log_gamma_list.append(log_alpha_list[idata] + log_beta_list[idata])
        #     pi = pi + np.exp(log_gamma_list[idata][0, :])
        # pi = pi/np.sum(pi)
        pi_i = pi_i0
        # emission
        # emission = np.zeros((nstate, nobs), dtype=np.float64)
        # for idata in range(ndata):
        #     data = data_list[idata]
        #     for istate in range(nstate):
        #         for iobs in range(nobs):
        #             id = (data == iobs)
        #             if np.any(id):
        #                 emission[istate, iobs] = emission[istate, iobs] + np.sum(np.exp(log_gamma_list[idata][id, istate]))
        # emission[np.isnan(emission)] = 0.0
        # emission = emission / np.sum(emission, axis=1)[:, None]
        emission = emission0
        # T
        T = zeros(Float64, (nstate, nstate))
        for idata = 1:ndata
          data = data_list[idata]
          nframe = length(data)
          for iframe = 2:nframe
            #log_xi = bsxfun(@plus, log_alpha{idata}(iframe-1, :)', log_beta{idata}(iframe, :));
            log_xi = log_alpha_list[idata][iframe-1, :] .+ log_beta_list[idata][iframe, :]'
            #T = T .+ exp(bsxfun(@plus, log_xi, log_emission0(:, data(iframe))') + log_T0)./factor{idata}(iframe);
            T = T .+ exp.((log_xi .+ log_emission0[:, data[iframe]]') .+ log_T0) ./ factor_list[idata][iframe]
          end
        end
        #T[np.isnan(T)] = 0.0
        T = T ./ sum(T, dims=2)
        ## Check convergence
        count_iteration += 1
        logL = sum(logL)
        check_convergence = abs(logL_old - logL)
        Printf.@printf("%d iteration LogLikelihood = %e  delta = %e  tolerance = %e\n" , count_iteration, logL, check_convergence, TOLERANCE)
        logL_old = logL
        pi_i0 = pi_i
        emission0 = emission
        T0 = T
    end
    T, emission, pi_i
end

function run()
    # T = np.array([[0.1, 0.7, 0.2],
    #               [0.2, 0.1, 0.7],
    #               [0.7, 0.2, 0.1]], dtype=np.float64)
    # emission = np.array([[0.9, 0.1],
    #                      [0.6, 0.4],
    #                      [0.1, 0.9]], dtype=np.float64)
    # pi = np.array([8.0, 1.0, 1.0], dtype=np.float64)
    # nframe = 1000
    # state, data_list = msmgenerate(nframe, T, emission, pi)
    # T0 = np.array([[0.1, 0.1, 0.8],
    #               [0.8, 0.1, 0.1],
    #               [0.3, 0.2, 0.5]], dtype=np.float64)
    # emission0 = np.array([[0.9, 0.1],
    #                      [0.6, 0.4],
    #                      [0.1, 0.9]], dtype=np.float64)
    # pi0 = np.array([0.8, 0.1, 0.1], dtype=np.float64)
    #matlab = sio.loadmat('test_baumwelch.mat')
    #data = np.asarray(matlab['data'], dtype=np.int).flatten() - 1
    data = readdlm("data.dat", ',', Int32)
    #T = np.asarray(matlab['T'], dtype=np.float64)
    T = readdlm("T.dat", ',', Float64)
    #emission = np.asarray(matlab['emission'], dtype=np.float64)
    emission = readdlm("emission.dat", ',', Float64)
    #pi_i = np.asarray(matlab['pi_i'], dtype=np.float64).flatten()
    pi_i = readdlm("pi_i.dat", ',', Float64)
    pi_i = pi_i[:]
    #T0 = np.asarray(matlab['T0'], dtype=np.float64)
    T0 = readdlm("T0.dat", ',', Float64)
    data_list = [data]

    @time T1, emission1, pi1 = msmbaumwelch(data_list, T0, emission, pi_i);
    println(T)
    println(T0)
    println(T1)
end

run()


text:
using RootFinding
using Test: @test, @testset

function sqrt_r_f(x, r)
    return x*x - r
end

@testset "RootFinding tests" begin
    @testset "Bisection Tests" begin
        @test isapprox(bisect_root(x -> sqrt_r_f(x, 2.0), 1.0, 4.0, atol=1e-12, rtol=0), sqrt(2))

        r = rand()
        sr = sqrt(r)
        @test abs(bisect_root(x -> sqrt_r_f(x, r), 0.0, 1.0, atol=1e-6, rtol=0) - sr) < 1e-6
    end

    @testset "Secant Tests" begin
        @test isapprox(secant_root(x -> sqrt_r_f(x, 2.0), 1.0, 4.0, atol=1e-12, rtol=0), sqrt(2))

        r = rand()
        sr = sqrt(r)
        @test abs(secant_root(x -> sqrt_r_f(x, r), 0.0, 1.0, atol=1e-6, rtol=0) - sr) < 1e-6
    end
end


text:
# Autogenerated wrapper script for libclangex_jll for x86_64-w64-mingw32-cxx03-julia_version+1.7.0
export libclangex

JLLWrappers.@generate_wrapper_header("libclangex")
JLLWrappers.@declare_library_product(libclangex, "libclangex.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libclangex,
        "bin\\libclangex.dll",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
LetUsPretendThisIsAnIterator = 1:10

z = 0
for i in LetUsPretendThisIsAnIterator
	z += i * i
end

h = i + 1


text:
function get_QM_data(QM :: QuadraticModel)
    T = eltype(QM.meta.lvar)
    # constructs A and Q transposed so we can create K upper triangular.
    # As Q is symmetric (but lower triangular in QuadraticModels.jl) we leave its name unchanged.
    AT = sparse(QM.data.Acols, QM.data.Arows, QM.data.Avals, QM.meta.nvar, QM.meta.ncon)
    dropzeros!(AT)
    Q = sparse(QM.data.Hcols, QM.data.Hrows, QM.data.Hvals, QM.meta.nvar, QM.meta.nvar)
    dropzeros!(Q)
    id = QM_IntData([QM.meta.ilow; QM.meta.irng], [QM.meta.iupp; QM.meta.irng], QM.meta.irng, QM.meta.ifree,
                    QM.meta.ncon, QM.meta.nvar, 0, 0)
    id.nlow, id.nupp = length(id.ilow), length(id.iupp) # number of finite constraints
    @assert QM.meta.lcon == QM.meta.ucon # equality constraint (Ax=b)
    fd_T = QM_FloatData(Q, AT, QM.meta.lcon, QM.data.c, QM.data.c0, QM.meta.lvar, QM.meta.uvar)
    return fd_T, id, T
end

function initialize(fd :: QM_FloatData{T}, id :: QM_IntData, res :: Residuals{T}, iconf :: InputConfig{Tconf},
                    T0 :: DataType) where {T<:Real, Tconf<:Real}

    itd = IterData(zeros(T, id.nvar+id.ncon), # Δxy
                   zeros(T, id.nlow), # Δs_l
                   zeros(T, id.nupp), # Δs_u
                   zeros(T, id.nlow), # x_m_lvar
                   zeros(T, id.nupp), # uvar_m_x
                   zeros(T, id.nvar), # init Qx
                   zeros(T, id.nvar), # init ATy
                   zeros(T, id.ncon), # Ax
                   zero(T), #xTQx
                   zero(T), #cTx
                   zero(T), #pri_obj
                   zero(T), #dual_obj
                   zero(T), #μ
                   zero(T),#pdd
                   zeros(T, 6), #l_pdd
                   one(T), #mean_pdd
                   nnz(fd.Q) > 0
                   )

    dda_type = Symbol(:DescentDirectionAllocs, iconf.solve_method)
    dda = eval(dda_type)(id, T)

    pad = PreallocatedData(iconf.sp, fd, id, iconf)

    # init system
    # solve [-Q-D    A' ] [x] = [b]  to initialize (x, y, s_l, s_u)
    #       [  A     0  ] [y] = [0]
    itd.Δxy[id.nvar+1: end] = fd.b

    cnts = Counters(zero(Int), zero(Int), 0, 0, iconf.kc==-1 ? nb_corrector_steps(pad.K.colptr, id.ncon, id.nvar, T) : iconf.kc,
                    iconf.max_ref, zero(Int))

    pt0 = Point(zeros(T, id.nvar), zeros(T, id.ncon), zeros(T, id.nlow), zeros(T, id.nupp))
    out = solver!(pad, dda, pt0, itd, fd, id, res, cnts, T0, :init)
    pt0.x .= itd.Δxy[1:id.nvar]
    pt0.y .= itd.Δxy[id.nvar+1:end]

    return itd, dda, pad, pt0, cnts
end

function init_params(fd_T :: QM_FloatData{T}, id :: QM_IntData, ϵ :: Tolerances{T}, sc :: StopCrit{Tc},
                     iconf :: InputConfig{Tconf}, T0 :: DataType) where {T<:Real, Tc<:Real, Tconf<:Real}

    res = Residuals(zeros(T, id.ncon), zeros(T, id.nvar), zero(T), zero(T), zero(T))

    itd, dda, pad, pt, cnts = initialize(fd_T, id, res, iconf, T0)

    starting_points!(pt, fd_T, id, itd)

    # stopping criterion
    #     rcNorm, rbNorm = norm(rc), norm(rb)
    #     optimal = pdd < ϵ_pdd && rbNorm < ϵ_rb && rcNorm < ϵ_rc
    res.rb .= itd.Ax .- fd_T.b
    res.rc .= itd.ATy .-itd.Qx .- fd_T.c
    res.rc[id.ilow] .+= pt.s_l
    res.rc[id.iupp] .-= pt.s_u
    res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
    set_tol_residuals!(ϵ, res.rbNorm, res.rcNorm)

    sc.optimal = itd.pdd < ϵ.pdd && res.rbNorm < ϵ.tol_rb && res.rcNorm < ϵ.tol_rc
    sc.small_μ = itd.μ < ϵ.μ

    return itd, ϵ, dda, pad, pt, res, sc, cnts
end

function set_tol_residuals!(ϵ :: Tolerances{T}, rbNorm :: T, rcNorm :: T) where {T<:Real}
    if ϵ.normalize_rtol == true
        ϵ.tol_rb, ϵ.tol_rc = ϵ.rb*(one(T) + rbNorm), ϵ.rc*(one(T) + rcNorm)
    else
        ϵ.tol_rb, ϵ.tol_rc = ϵ.rb, ϵ.rc
    end
end

############ tools for sparse matrices ##############

function get_diag_sparseCSC(M_colptr, n; tri =:U)
    # get diagonal index of M.nzval
    # we assume all columns of M are non empty, and M triangular (:L or :U)
    @assert tri ==:U || tri == :L
    diagind = zeros(Int, n) # square matrix
    if tri == :U
        @inbounds @simd for i=1:n
            diagind[i] = M_colptr[i+1] - 1
        end
    else
        @inbounds @simd for i=1:n
            diagind[i] = M_colptr[i]
        end
    end
    return diagind
end

function get_diag_Q(Q_colptr, Q_rowval, Q_nzval :: Vector{T}, n) where {T<:Real}
    diagval = spzeros(T, n)
    @inbounds @simd for j=1:n
        for k=Q_colptr[j]:(Q_colptr[j+1]-1)
            if j == Q_rowval[k]
                diagval[j] = Q_nzval[k]
            end
        end
    end
    return diagval
end


text:
module MA
include("../src/common/ClassMacro.jl")
export A, benchmark, benchmark2
struct A
    x::Float32
    y::Integer
    z::Bool
    function A(x::Float32, y::Integer, z::Bool)
        this = new(x, y, z)
        return this
    end
end
@class A [:test, :test3, :benchmark, :benchmark2]

function test(this::A)
    println("------ test ", this)
end

function test3(this::A, b::Integer, c::Float32)
    println("------ test3 ", this, " ", b, " ", c)
end

function test3(this::A, b::Integer)
    println("------ test3 ", this, " ", b)
end

function benchmark(this::A)
    return 1
end

function benchmark2(this::A, b::Integer, c::Float32)
    return b + c
end

function benchmark2(this::A, b::Integer)
    return b - 1
end
end

module MB
include("../src/common/ClassMacro.jl")
export B
mutable struct B
    x::Float32
    y::Integer
    z::Bool
    function B(x::Float32, y::Integer, z::Bool)
        this = new(x, y, z)
        return this
    end
end
@class B [:funA, :funB]

function funA(this::B, b::Integer, c::Float32)
    println("------ funA ", this, " ", b, " ", c)
end

function funA(this::B, b::Integer)
    println("------ funA ", this, " ", b)
end

function funB(this::B)
    println("------ funB ", this)
end

end

using .MA, .MB
function testClassMethod()
    a = A(1.2f0, 3, false)
    println(Base.hasfield(A, :test3))
    println(Base.hasproperty(a, :test3))
    println(Base.hasproperty(a, :test33))
    a.test3(11, 1.f0)
    a.test3(22)
    a.test()
    println(a.x)

    b = B(4.5f0, 6, true)
    b.funA(33, 3.f0)
    b.funA(44)
    b.funB()

    return 0
end


text:
################################################################################
# This set of functions implements the Elston-Stewart algorithm for
# evaluating the likelihood of a pedigree.
################################################################################
#
# Required OpenMendel packages and modules.
#
# using DataStructures, ElstonStewartPreparation, ModelConstruction, ReadData
# using MendelSearch

export elston_stewart_loglikelihood

"""
Evaluate the loglikelihood of a collection of independent pedigrees
by the Elston-Stewart algorithm.
"""
function elston_stewart_loglikelihood(penetrance::Function, prior::Function,
  transmission::Function, pedigree::Pedigree, person::Person,
  locus::Locus, parameter::Parameter, instruction::Instruction,
  person_frame::DataFrame, keyword::Dict{AbstractString, Any})

  loglikelihood = 0.0
  fill!(pedigree.loglikelihood[:, 2], 0.0)
  par = parameter.par
  #
  # Given independence, loop over all pedigrees.
  #
  for ped = 1:pedigree.pedigrees
    #
    # Retrieve the population frequencies of the lumped alleles.
    #
    if length(locus.lumped_frequency) > 0
      for l = 1:locus.model_loci
        loc = locus.model_locus[l]
        locus.frequency[loc][:, end] = locus.lumped_frequency[ped, :, l]
      end
    end
    #
    # Compute the likelihood of the current pedigree.
    # Throw an error when the likelihood is 0.
    #
    (inconsistent, pedigree.loglikelihood[ped, 2]) =
      compute_likelihood(penetrance, prior, transmission, person,
      locus, instruction, par, person_frame, keyword, ped)
    if inconsistent
      name = pedigree.name[ped]
      throw(ArgumentError(
        "Error: in likelihood evaluation for pedigree $name.\n \n"))
    end
    loglikelihood = loglikelihood + pedigree.loglikelihood[ped, 2]
  end
  return loglikelihood
end # function elston_stewart_loglikelihood

"""
Compute the likelihood of a single pedigree
by the Elston-Stewart algorithm.
Various operations are performed on an array of arrays.
"""
function compute_likelihood(penetrance::Function, prior::Function,
  transmission::Function, person::Person, locus::Locus,
  instruction::Instruction, par::Vector{Float64},
  person_frame::DataFrame, keyword::Dict{AbstractString, Any}, ped::Int)
  #
  # Initialize the loglikelihood, the number of arrays, and the array of arrays.
  #
  inconsistent = false
  loglikelihood = 0.0
  work_arrays = instruction.instructfinish[ped] - instruction.instructstart[ped]
  work_array = Array{Array{Float64, 1}}(undef, work_arrays)
  #
  # Now proceed instruction by instruction.
  #
  extent = 0
  increment1 = 0
  increment2 = 0
  for n = instruction.instructstart[ped]:instruction.instructfinish[ped]
    #
    # Fetch the site, rank, and size of the new array created by the
    # current instruction. Then allocate the array. If a multiply
    # and add operation is indicated, then rank pertains to the
    # product array before the pivot index is summed over.
    #
    operation = instruction.operation[n]
    if operation != quit_processing_pedigree
      k = instruction.site[n]
      if operation != pure_add
        rankarray = instruction.rankarray[n]
        space = instruction.space[n]
        work_array[k] = zeros(space)
      end
    end
    #
    # Each instruction involves one of six operations. Two of these,
    # "multiply_and_add" and "pure_multiply" are treated simultaneously.
    #
    if operation == penetrance_and_prior_array
      #
      # Create a penetrance and prior array.
      #
      work_array[k] = construct_penetrance_prior!(penetrance, prior, person,
        locus, instruction, par, person_frame, keyword, work_array[k], n)
      (negative_entry, no_positive_entry) = check_array_entries(work_array[k])
      if negative_entry || no_positive_entry
        inconsistent = true
        println("Error: in penetrance or prior array.")
        break
      end
    #
    # Create a transmission array.
    #
    elseif operation == transmission_array
      work_array[k] = construct_transmission!(transmission, person, locus,
        instruction, par, person_frame, keyword, work_array[k], n)
      (negative_entry, no_positive_entry) = check_array_entries(work_array[k])
      if negative_entry || no_positive_entry
        inconsistent = true
        println("Error: in transmission array.")
        break
      end
    #
    # Perform a pure add.
    #
    elseif operation == pure_add
      (loglikelihood, zero_likelihood) =
        pure_add_operation(loglikelihood, work_array[k])
      if zero_likelihood
        inconsistent = true
        println("Error: in pure add operation.")
        break
      end
    #
    # Perform a multiply and add or a pure multiply.
    #
    elseif (operation == multiply_and_add || operation == pure_multiply)
      #
      # Fetch the positions of the multiplicand arrays. The number
      # pivot_range is 0 for a pure multiply and one less than the number
      # of entries to be added for a multiply and add. For a multiply
      # and add, record the summation strides for the muliplicand arrays.
      #
      i = instruction.extra[n][1]
      j = instruction.extra[n][2]
      pivot_range = instruction.extra[n][3]
      stride1 = instruction.extra[n][4]
      stride2 = instruction.extra[n][5]
      #
      # Insert into more easily accessed arrays, the extent of each dimension
      # of the product array and the increments to the positions within each
      # the multiplicand arrays.
      #
      if rankarray > 0
        extent = instruction.extra[n][6:5 + rankarray]
        increment1 = instruction.extra[n][6 + rankarray:5 + 2 * rankarray]
        increment2 = instruction.extra[n][6 + 2 * rankarray:5 + 3 * rankarray]
      end
      #
      # Perform the operation.
      #
      (loglikelihood, zero_likelihood, work_array[k]) =
        multiply_add_operation!(work_array[i], work_array[j], work_array[k],
        extent, increment1, increment2, rankarray, pivot_range,
        stride1, stride2, loglikelihood)
      if zero_likelihood
        inconsistent = true
        println("Error: in a multiply and add multiply operation.")
        break
      end
    elseif operation == quit_processing_pedigree
    #
    # Terminate processing of the current pedigree.
    #
      break
    end
  end
  return (inconsistent, loglikelihood)
end # function compute_likelihood

"""
Construct a penetrance array or the product of a penetrance and prior array.
"""
function construct_penetrance_prior!(penetrance::Function, prior::Function,
  person::Person, locus::Locus, instruction::Instruction, par::Vector{Float64},
  person_frame::DataFrame, keyword::Dict{AbstractString, Any},
  work_array::Vector{Float64}, n::Int)
  #
  # Fetch the starting & finishing loci and individual i.
  #
  startlocus = instruction.extra[n][1]
  finishlocus = instruction.extra[n][2]
  i = instruction.extra[n][3]
  #
  # Construct i's multiple locus genotypes.
  #
  genotypes = genotype_count(person, locus, i, startlocus, finishlocus)
  multi_genotype = construct_multigenotypes(person, locus, startlocus,
                                            finishlocus, genotypes, i)
  #
  # Compute the penetrance of each genotype of i.
  #
  for j = 1:genotypes
    a = penetrance(person, locus, multi_genotype[:, :, j], par,
                   person_frame, keyword, startlocus, finishlocus, i)
    if a < 0.0
      throw(ArgumentError(
        "Negative penetrance computed for $person.name[i].\n \n"))
    end
    #
    # Include the penetrance of each identical twin of i.
    #
    m = person.next_twin[i]
    while m > 0
      b = penetrance(person, locus, multi_genotype[:, :, j], par,
                     person_frame, keyword, startlocus, finishlocus, m)
      if b < 0.0
        throw(ArgumentError(
          "Negative penetrance computed for $person.name[m].\n \n"))
      end
      a = a * b
      m = next_twin[m]
    end
    #
    # Include the prior when i is a founder.
    #
    if person.mother[i] == 0
      b = prior(person, locus, multi_genotype[:, :, j], par,
                person_frame, keyword, startlocus, finishlocus, i)
      if b < 0.0
        throw(ArgumentError(
          "Negative prior computed for $person.name[i].\n \n"))
      end
      a = a * b
    end
    work_array[j] = a
  end
  return work_array
end # function construct_penetrance_prior!

"""
Construct a transmission array from parent i to child j.
"""
function construct_transmission!(transmission::Function, person::Person,
  locus::Locus, instruction::Instruction, par::Vector{Float64},
  person_frame::DataFrame, keyword::Dict{AbstractString, Any},
  work_array::Vector{Float64}, n::Int)
  #
  # Fetch the starting & finishing loci and i & j.
  #
  startlocus = instruction.extra[n][1]
  finishlocus = instruction.extra[n][2]
  i = instruction.extra[n][3]
  j = instruction.extra[n][4]
  #
  # Construct the parent's multiple locus genotypes.
  #
  i_genotypes = genotype_count(person, locus, i, startlocus, finishlocus)
  multi_genotype = construct_multigenotypes(person, locus, startlocus,
                                            finishlocus, i_genotypes, i)
  #
  # Construct the child's multiple locus gametes.
  #
  maternal = !person.male[i]
  j_genotypes = genotype_count(person, locus, j, startlocus, finishlocus)
  gamete = construct_gametes(person, locus, startlocus, finishlocus,
                             j_genotypes, j, maternal)
  #
  # Insert the transmission probabilities into the transmission array in the
  # correct reverse dictionary order.
  #
  m = 0
  for l = 1:j_genotypes
    for k = 1:i_genotypes
      m = m + 1
      work_array[m] = transmission(person, locus, gamete[:, l],
          multi_genotype[:, :, k], par, person_frame, keyword,
          startlocus, finishlocus, i, j)
      if work_array[m] < 0.0
        throw(ArgumentError("Negative transmission computed "*
          "from parent $person.name[i] to child $person.name[j].\n \n"))
      end
    end
  end
  return work_array
end # function construct_transmission!

"""
Count i's possible multilocus genotypes between loci start and finish.
"""
function genotype_count(person::Person, locus::Locus, i::Int,
  startlocus::Int, finishlocus::Int)

  count = 1
  for l = startlocus:finishlocus
    loc = locus.model_locus[l]
    count = count * length(person.genotype[i, loc])
  end
  return count
end # function genotype_count

"""
Check whether the array A has any negative entries and
any positive entries.
"""
function check_array_entries(array_a::Vector{Float64})

  negative_entry = minimum(array_a) < 0.0
  no_positive_entry = maximum(array_a) <= 0.0
  return (negative_entry, no_positive_entry)
end # function check_array_entries

"""
Carry out a pure add operation.
"""
function pure_add_operation(loglikelihood::Float64, work_array::Vector{Float64})

  array_sum = sum(work_array)
  if array_sum <= 0.0
    zero_likelihood = true
  else
    zero_likelihood = false
    loglikelihood = loglikelihood + log(array_sum)
  end
  return (loglikelihood, zero_likelihood)
end # function pure_add_operation

"""
Multiply work_array1 times work_array2 and sum on the pivot index if required.
Insert the result in work_array3 and rescale to avoid underflows and overflows.
?? Can this be recoded to call a LAPACK routine ??
"""
function multiply_add_operation!(work_array1::Vector{Float64},
  work_array2::Vector{Float64}, work_array3::Vector{Float64},
  extent::Vector{Int}, increment1::Vector{Int}, increment2::Vector{Int},
  rankarray::Int, pivot_range::Int, stride1::Int, stride2::Int,
  loglikelihood::Float64)

  address = ones(Int, rankarray)
  #
  # Initialize the positions in the multiplicand arrays
  # and the reverse dictionary address in the product array.
  #
  position1 = 1
  position2 = 1
  big = 0.0
  #
  # Compute the product array entry by entry.
  #
  for i = 1:length(work_array3)
    a = work_array1[position1] * work_array2[position2]
    #
    # Sum on the pivot index whenever pivot_range > 0.
    #
    for j = 1:pivot_range
      position1 = position1 + stride1
      position2 = position2 + stride2
      a = a + work_array1[position1] * work_array2[position2]
    end
    #
    # Save the result and keep a running maximum.
    #
    work_array3[i] = a
    big = max(big, a)
    #
    # Update the address in the product array
    # and the positions in the multiplicand arrays.
    #
    for j = 1:rankarray
      if address[j] < extent[j]
        address[j] = address[j] + 1
        position1 = position1 + increment1[j]
        position2 = position2 + increment2[j]
        break
      else
        address[j] = 1
      end
    end
  end
  #
  # Check for a zero-likelihood error. If none exists, rescale and
  # adjust the loglikelihood.
  #
  if big <= 0.0
    zero_likelihood = true
  else
    zero_likelihood = false
    loglikelihood = loglikelihood + log(big)
    work_array3 = work_array3 / big
  end
  return (loglikelihood, zero_likelihood, work_array3)
end # function multiply_add_operation!

"""
Create the multi-locus genotypes for person i from the single locus genotypes.
Only model loci between positions start and finish are considered.
"""
function construct_multigenotypes(person::Person, locus::Locus,
  startlocus::Int, finishlocus::Int, genotypes::Int, i::Int)

  multi_genotype = zeros(Int, 2, locus.model_loci, genotypes)
  #
  # Create the multiple locus genotypes by taking the appropriate
  # Cartesian product. Nested for loops here by computing the reverse
  # dictionary order location of each element of the product.
  #
  copies = 1
  increment = 1
  for l = startlocus:finishlocus
    loc = locus.model_locus[l]
    increment = increment * length(person.genotype[i, loc])
    initial_position = 1
    for (gm, gp) in person.genotype[i, loc]
      current = initial_position
      while current <= genotypes
        j = current
        for k = 1:copies
          multi_genotype[1, l, j] = gm
          multi_genotype[2, l, j] = gp
          j = j + 1
        end
        current = current + increment
      end
      initial_position = initial_position + copies
    end
    copies = increment
  end
  return multi_genotype
end # function construct_multigenotypes

"""
Create the maternal or paternal gametes for person i
from the single locus genotypes.
Only model loci between positions start and finish are considered.
"""
function construct_gametes(person::Person, locus::Locus,
  startlocus::Int, finishlocus::Int, genotypes::Int, i::Int, maternal::Bool)

  gamete = zeros(Int, locus.model_loci, genotypes)
  #
  # Create the multiple locus gametes by taking the appropriate Cartesian
  # product. Nested for loops are avoided here by computing the reverse
  # dictionary order location of each element of the product.
  #
  copies = 1
  increment = 1
  for l = startlocus:finishlocus
    loc = locus.model_locus[l]
    increment = increment * length(person.genotype[i, loc])
    initial_position = 1
    for (gm, gp) in person.genotype[i, loc]
      current = initial_position
      while current <= genotypes
        j = current
        for k = 1:copies
          if maternal
            gamete[l, j] = gm
          else
            gamete[l, j] = gp
          end
          j = j + 1
        end
        current = current + increment
      end
      initial_position = initial_position + copies
    end
    copies = increment
  end
  return gamete
end # function construct_gametes

max_stars_repo_name: OpenMendel/MendelBase

max_issues_repo_name: OpenMendel/MendelBase

max_forks_repo_name: OpenMendel/MendelBase


text:
using Documenter, SnpArrays

makedocs(
    format = Documenter.HTML(),
    sitename = "SnpArrays.jl",
    authors = "Hua Zhou",
    clean = true,
    debug = true,
    pages = [
        "index.md",
        "linalg.md"
    ]
)

deploydocs(
    repo   = "github.com/OpenMendel/SnpArrays.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
module CuCountMap

export countmap, cucountmap, cucountmap!

include("countmap.jl")
# include("gpuradixsort.jl")

end # module


text:
module Libtask

using CodeInfoTools
using FunctionWrappers: FunctionWrapper
using LRUCache

export TapedTask, consume, produce
export TArray, tzeros, tfill, TRef

include("tapedfunction.jl")
include("tapedtask.jl")

include("tarray.jl")
include("tref.jl")

end

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
#
# the API
#

"""
    initSOM(train, xdim, ydim = xdim;  norm = :zscore, topol = :hexagonal,
            toroidal = false)

Initialises a SOM.

# Arguments:
- `train`: training data
- `xdim, ydim`: geometry of the SOM
           If DataFrame, the column names will be used as attribute names.
           Codebook vectors will be sampled from the training data.
           for spherical SOMs ydim can be omitted.
- `norm`: optional normalisation; one of :`minmax, :zscore or :none`
- `topol`: topology of the SOM; one of `:rectangular, :hexagonal or :spherical`.
- `toroidal`: optional flag; if true, the SOM is toroidal.
"""
function initSOM( train, xdim, ydim = xdim;
             norm::Symbol = :none, topol = :hexagonal, toroidal = false)

    if typeof(train) == DataFrame
        colNames = [String(x) for x in names(train)]
    else
        colNames = ["x$i" for i in 1:ncol(train)]
    end
    train = convertTrainingData(train)

    if topol == :spherical
        toroidal = false
        nCodes = xdim
    else
        nCodes = xdim * ydim
    end

    som = initAll(train, colNames, norm,
                 xdim, ydim, nCodes,
                 topol, toroidal)
    return som
end

"""
    trainSOM(som::Som, train::Any, len;
             η = 0.2, kernelFun = gaussianKernel,
             r = 0.0, rDecay = true, ηDecay = true)

Train an initialised or pre-trained SOM.

# Arguments:
- `som`: object of type Som with a trained som
- `train`: training data
- `len`: number of single training steps (*not* epochs)
- `η`: learning rate
- `kernel::function`: optional distance kernel; one of (`bubbleKernel, gaussianKernel`)
            default is `gaussianKernel`
- `r`: optional training radius.
       If r is not specified, it defaults to √(xdim^2 + ydim^2) / 2
- `rDecay`: optional flag; if true, r decays to 0.0 during the training.
- `ηDecay`: optional flag; if true, learning rate η decays to 0.0
            during the training.

Training data must be convertable to Array{Float64,2} with
`convert()`. Training samples are row-wise; one sample per row.

An alternative kernel function can be provided to modify the distance-dependent
training. The function must fit to the signature fun(x, r) where
x is an arbitrary distance and
r is a parameter controlling the function and the return value is
between 0.0 and 1.0.
"""
function trainSOM(som::Som, train::Any, len;
                     η = 0.2, kernelFun::Function = gaussianKernel, r = 0.0,
                     rDecay = true, ηDecay = true)

    train = convertTrainingData(train)
    som = trainAll(som, train,
                 len, η, kernelFun, r,
                 rDecay, ηDecay)
    return som
end

#
#
# predict functions:
#
#
"""
    mapToSOM(som::Som, data)

Return a DataFrame with X-, Y-indices and index of winner neuron for
every row in data.

# Arguments
- `som`: a trained SOM
- `data`: Array or DataFrame with training data.

Data must have the same number of dimensions as the training dataset
and will be normalised with the same parameters.
"""
function mapToSOM(som::Som, data)

    data = convertTrainingData(data)

    if ncol(data) != ncol(som.codes)
        println("    data: $(ncol(data)), codes: $(ncol(som.codes))")
        error(SOM_ERRORS[:ERR_COL_NUM])
    end

    vis = visual(som.codes, data)
    x = [som.indices[i,:X] for i in vis]
    y = [som.indices[i,:Y] for i in vis]

    return DataFrame(X = x, Y = y, index = vis)
end

"""
    classFrequencies(som::Som, data, classes)

Return a DataFrame with class frequencies for all neurons.

# Arguments:
- `som`: a trained SOM
- `data`: data with row-wise samples and class information in each row
- `classes`: Name of column with class information.

Data must have the same number of dimensions as the training dataset.
The column with class labels is given as `classes` (name or index).
Returned DataFrame has the columns:
* X-, Y-indices and index: of winner neuron for every row in data
* population: number of samples mapped to the neuron
* frequencies: one column for each class label.
"""
function classFrequencies(som::Som, data, classes)

    if ncol(data) != ncol(som.codes) + 1
        println("    data: $(ncol(data)-1), codes: $(ncol(som.codes))")
        error(SOM_ERRORS[:ERR_COL_NUM])
    end

    x = deepcopy(data)
    deletecols!(x, classes)
    classes = data[classes]
    vis = visual(som.codes, x)

    df = makeClassFreqs(som, vis, classes)
    return df
end

#
#
# Plotting functions:
#
#
"""
    plotDensity(som::Som; predict = nothing,
                title = "Density of Self-Organising Map",
                paper = :a4r,
                colormap = "autumn_r",
                detail = 45,
                device = :display, fileName = "somplot")

Plot the population of neurons as colours.

# Arguments:
- `som`: the som of type `Som`; som is the only mandatory argument
- `predict`: DataFrame of mappings as outputed by `mapToSOM()`
- `title`: main title of plot
- `paper`: plot size; currentlx supported: `:a4, :a4r, :letter, :letterr`
- `colormap`: MatPlotLib colourmap (Python-style strings; e.g. `"Greys"`).
- `detail`: only relevant for 3D-plotting of spherical SOMs: higher
            values result in smoother display of the 3D-sphere
- `device`: one of `:display, :png, :svg, :pdf` or any file-type supported
            by MatPlotLib; default is `:display`
- `fileName`: name of image file. File extention overrides the setting of
              `device`.
"""
function plotDensity(som::Som; predict = nothing,
                     title = "Density of Self-Organising Map",
                     paper = :a4r,
                     colormap = "autumn_r",
                     detail = 45,
                     device = :display, fileName = "somplot")

    # do nothing, if matplotlib is not installed correctly:
    #
    if !MPL_INSTALLED
        error(SOM_ERRORS[:ERR_MPL])
    end

    # use population form the som itself, if no prediction is given as arg.
    #
    if predict == nothing
        population = som.population
    else
        population = makePopulation(som.nCodes, predict[:index])
    end

    if typeof(colormap) == Symbol
        colormap = string(colormap)
    end

    if som.topol == :spherical
        drawSpherePopulation(som, population, detail, title,
                             paper, colormap, device, fileName)
    else
        drawPopulation(som, population, title, paper, colormap, device, fileName)
    end
end

"""
    plotClasses(som::Som, frequencies;
                title = "Class Frequencies of Self-Organising Map",
                paper = :a4r,
                colors = "brg",
                detail = 45,
                device = :display, fileName = "somplot")

Plot the population of neurons as colours.

# Arguments:
- `som`: the som of type `Som`
- `frequencies`: DataFrame of frequencies as outputed by classFrequencies()
- `title`: main title of plot
- `paper`: plot size; currentlx supported: `:a4, :a4r, :letter, :letterr`
- `colors`: MatPlotLib colourmap (Python-style as string `"gray"` or
              Julia-style as Symbol `:gray`) *or* dictionary with
              classes as keys and colours as vals;
              keys can be provides as Strings or Symbols; colours must be
              valid coulour definitions (such as RGB, names, etc).
              Default: `brg`
- `detail`: only relevant for 3D-plotting of spherical SOMs: higher
            values result in smoother display of the 3D-sphere
- `device`: one of `:display, :png, :svg, :pdf` or any file-type supported
            by MatPlotLib; default is `:display`
- `fileName`: name of image file. File extention overrides the setting of
              `device`.
"""
function plotClasses(som::Som, frequencies;
                     title = "Class Frequencies of Self-Organising Map",
                     paper = :a4r,
                     colors = "brg",
                     detail = 45,
                     device = :display, fileName = "somplot")

    # do nothing, if matplotlib is not installed correctly:
    #
    if !MPL_INSTALLED
        error(SOM_ERRORS[:ERR_MPL])
    end

    # create dictionary of colours and classes if necessary:
    #
    if isa(colors, Symbol)
        colors = string(colors)
    end
    if isa(colors, String)
        numClasses = size(frequencies)[2] - 4
        classes = sort(names(frequencies)[5:end])
        cmap = get_cmap(colors)
        coloursRGB = cmap.(range(0.0, stop = 1.0, length = numClasses))

        colourDict = Dict((classes[i], coloursRGB[i]) for i in 1:numClasses)

    # else create Dict of colours with Symbols as keys:
    #
    elseif isa(colors, Dict)
        colourDict = Dict(Symbol(i) => colors[i] for i in keys(colors))
    else
        println(SOM_ERRORS[:ERR_COLOUR_DEF])
        return :ERR_COLOUR_DEF
    end

    if som.topol == :spherical
        drawSphereFreqs(som, frequencies, detail, title, paper, colourDict, device, fileName)
    else
        drawFrequencies(som, frequencies, title, paper, colourDict, device, fileName)
    end
end


text:
using Printf, Test

include("../../hybmc/termstructures/YieldCurve.jl")
include("../../hybmc/models/HullWhiteModel.jl")
include("../../hybmc/simulations/McSimulation.jl")
include("../../hybmc/simulations/Payoffs.jl")

function setup()
    curve = YieldCurve(0.03)
    mean = 0.03
    times = [ 1.0,    2.0,    5.0,    10.0     ]
    vols  = [ 0.0060, 0.0070, 0.0080,  0.0100  ]
    model = HullWhiteModel(curve,mean,times,vols)
    return model
end

@testset "test_zeroBondPrice" begin
    model = setup()
    obsTime = 10.0
    matTime = 20.0
    states  = [ -0.10, -0.05, 0.00, 0.03, 0.07, 0.12 ]
    refResult = [  # reference results checked against QuantLib
        1.7178994273756056,
        1.1153102854629025,
        0.7240918839816156,
        0.5587692049384843,
        0.39550396436404667,
        0.25677267968500767]
    for (x,res) in zip(states,refResult)
        pv = zeroBondPrice(model,obsTime,matTime,x)
        @test isapprox(pv,res,atol=1.0e-16, rtol=0.0)
        #println(pv)
    end
end

@testset "test_couponBondOption" begin
    model = setup()
    excTime = 10.0
    payTimes  = [ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 20.0 ]
    cashFlows = cat(0.025*ones(10), [1.0], dims=1) # bond cash flows
    callOrPut = 1.0 # call
    strikes   = [ 0.8, 0.9, 1.0, 1.1, 1.2 ]
    refResult = [  # reference results checked against QuantLib
    0.12572604970665557,
    0.07466919528750032,
    0.039969719617809936,
    0.019503405494683164,
    0.008804423283331253]
    for (strike,res) in zip(strikes,refResult)
        pv = couponBondOption(model,excTime,payTimes,cashFlows,strike,callOrPut)
        @test isapprox(pv,res,atol=1.0e-8, rtol=0.0)
        #println(pv)
    end
    return nothing
end

@testset "test_monteCarloSimulation" begin
    model = setup()
    times = Array(range(0.0, stop=10.0, length=11))
    nPaths = 2^11
    seed = 1234
    # risk-neutral simulation
    mcSim = McSimulation(model,times,nPaths,seed)
    discZeroBonds = [
        (1.0 / numeraire(model,times[i],mcSim.X[:,i,j]))
        for i = 1:size(times)[1], j = 1:nPaths ]
    mcZeroBondsRiskNeutral = mean(discZeroBonds, dims=2)
    zeroBonds = [ discount(model.yieldCurve,t) for t in times ]
    @printf("  T     ZeroRate    RiskNeutral\n")
    for k = 2:size(times)[1]
        t           = times[k]
        zeroRate    = -log(zeroBonds[k])/t
        riskNeutral = -log(mcZeroBondsRiskNeutral[k])/t
        @printf(" %4.1f   %8.5lf    %8.5lf\n", t, zeroRate, riskNeutral)
        @test isapprox(zeroRate,riskNeutral,atol=1.0e-8, rtol=2.0e-2)
    end
    return nothing
end

@testset "test_liborPayoff" begin
    model = setup()
    # Libor Rate payoff
    liborTimes = Array(range(0.0, stop=10.0, length=11))
    dT0 = 2.0 / 365
    dT1 = 0.5
    cfs = [ Pay(LiborRate(t,t+dT0,t+dT0+dT1),t+dT0+dT1) for t in liborTimes]
    times = observationTimes(cfs)
    #
    nPaths = 2^11
    seed = 4321
    # risk-neutral simulation
    mcSim = McSimulation(model,times,nPaths,seed)
    cfPVs = discountedAt(cfs, mcSim)
    cfPVs = mean(cfPVs, dims=2)
    discounts0 = [ discount(mcSim.model.yieldCurve,t+dT0) for t in liborTimes ]
    discounts1 = [ discount(mcSim.model.yieldCurve,t+dT0+dT1) for t in liborTimes ]
    liborsCv = (discounts0./discounts1 .- 1.0)./dT1
    liborsMC = cfPVs ./ discounts1
    println("  T     LiborRate  MnteCarlo")
    for (k,t) in enumerate(liborTimes)
        @printf(" %4.1f   %7.5f     %7.5f\n", t, liborsCv[k], liborsMC[k])
    end
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
include("arrviews.jl")
include("contrank.jl")
include("subviews.jl")
include("linalg.jl")
include("convenience.jl")


text:
struct RandCovariate <: AbstractCovariate{Float32}
    name::String
    length::Integer
end

function RandCovariate(length::Integer)
    RandCovariate("", length)
end

Base.length(covariate::RandCovariate) = covariate.length

function slice(covariate::RandCovariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    if fromobs > toobs
        EmptySeq{SubArray{Float32,1,Array{Float32,1},Tuple{UnitRange{Int64}},true}}()
    else
        fromobs = max(1, fromobs)
        toobs = min(toobs, length(covariate))
        slicelength = verifyslicelength(fromobs, toobs, slicelength)
        buffer = Vector{Float32}(undef, slicelength)
        map(Seq(Tuple{Int64, Int64}, (fromobs, toobs, slicelength), nextslice), SubArray{Float32,1,Array{Float32,1},Tuple{UnitRange{Int64}},true}) do rng
            rand!(buffer)
            view(buffer, 1:(rng[2] - rng[1] + 1))
        end
    end
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using MatchCore.Err: PatternCompilationError

PatternImpl = NamedTuple{
    (:and, :or, :literal, :wildcard, :decons, :guard, :effect),
}

PatternImpls{N} = NTuple{N, PatternImpl}

@nospecialize
and(args :: Any...) = and(collect(args))
and(ps::Vector) = function apply(impls::PatternImpls{N}) where N
    xs = [p(impls) for p in ps]
    me = Vector{Any}(undef, N)
    for i in 1:N
        @inbounds me[i] = impls[i].and(xs)
    end
    me
end

or(args :: Any...) = or(collect(args))
or(ps::Vector) = function apply(impls::PatternImpls{N}) where N
    xs = [p(impls) for p in ps]
    me = Vector{Any}(undef, N)
    for i in 1:N
        me[i] = impls[i].or(xs)
    end
    me
end

literal(val) = function apply(impls::PatternImpls{N}) where N
    me = Vector{Any}(undef, length(impls))
    for i in 1:N
        me[i] = impls[i].literal(val)
    end
    me
end

function wildcard(impls::PatternImpls{N}) where N
    me = Vector{Any}(undef, length(impls))
    for i in 1:N
        me[i] = impls[i].wildcard
    end
    me
end

guard(pred) = function apply(impls::PatternImpls{N}) where N
    me = Vector{Any}(undef, N)
    for i in 1:N
        me[i] = impls[i].guard(pred)
    end
    me
end

"""
abstract pure process
"""
abstract type APP end

struct NoncachablePre <: APP
    callable :: Any
end
(f::NoncachablePre)(target::Any) = f.callable(target)
struct NoPre <: APP end

"""composite pattern
"""
struct PComp
    repr :: AbstractString
    tcons :: Function
    guard1 :: APP
    view :: APP
    guard2 :: APP
end

invalid_extract(_, _) = error("impossible")

function PComp(
    repr :: AbstractString,
    tcons::Function;
    guard1::APP=NoPre(),
    view::APP=NoPre(),
    guard2::APP=NoPre()
)
    PComp(repr, tcons, guard1, view, guard2)
end

const any_type(_...) = Any
const as_is_comp = PComp("identity", any_type)

decons(extract, ps) = decons(as_is_comp, extract, ps)
decons(comp::PComp, ps; extract=invalid_extract) = decons(comp, extract, ps)
decons(comp::PComp, extract::Function, ps) = function apply(impls::PatternImpls{N}) where N
    xs = [p(impls) for p in ps]
    me = Vector{Any}(undef, N)
    for i in 1:N
        me[i] = impls[i].decons(comp, extract, xs)
    end
    me
end

effect(ctx_perf) = function apply(impls::PatternImpls{N}) where N
    me = Vector{Any}(undef, N)
    for i in 1:N
        me[i] = impls[i].effect(ctx_perf)
    end
    me
end

@specialize

const self = (
    and = and,
    or = or,
    literal = literal,
    wildcard = wildcard,
    decons = decons,
    guard = guard,
    effect = effect
)


text:
#!/usr/bin/env julia

using Dtree
if VERSION > v"0.5.0-dev"
    using Base.Threads
else
    # Pre-Julia 0.5 there are no threads
    nthreads() = 1
    threadid() = 1
    macro threads(x)
        x
    end
    SpinLock() = 1
    lock(l) = ()
    unlock(l) = ()
end

const cpu_hz = 2.3e9

const libdtree = joinpath(Pkg.dir("Dtree"), "deps", "Dtree",
        "libdtree.$(Libdl.dlext)")

@inline rdtsc() = ccall((:rdtsc, libdtree), Culonglong, ())
@inline cpu_pause() = ccall((:cpu_pause, libdtree), Void, ())
@inline secs2cpuhz(s) = s * cpu_hz::Float64

@inline function ntputs(nid,tid,s)
    ccall(:puts, Cint, (Cstring,), string("[", nid, "]<", tid, "> ", s))
    return
end

function threadfun(dt, ni, ci, li, ilock, rundt, dura)
    tid = threadid()
    if rundt && tid == 1
        ntputs(dt_nodeid, tid, "running tree")
        while runtree(dt)
            cpu_pause()
        end
    else
        ntputs(dt_nodeid, tid, string("begin, ", ni[], " items, ", length(dura), " available delays"))
        while ni[] > 0
            lock(ilock)
            if li[] == 0
                ntputs(dt_nodeid, tid, string("out of work"))
                unlock(ilock)
                break
            end
            if ci[] == li[]
                ntputs(dt_nodeid, tid, string("work consumed (last was ", li[], "); requesting more"))
                ni[], (ci[], li[]) = getwork(dt)
                ntputs(dt_nodeid, tid, string("got ", ni[], " work items (", ci[], " to ", li[], ")"))
                unlock(ilock)
                continue
            end
            item = ci[]
            ci[] = ci[] + 1
            unlock(ilock)

            # wait dura[item] seconds
            ticks = secs2cpuhz(dura[item])
            #ntputs(dt_nodeid, tid, string("item ", item, ", ", ticks, " ticks"))
            startts = rdtsc()
            while rdtsc() - startts < ticks
                cpu_pause()
            end
        end
    end
end

function bench(nwi, meani, stddevi, first_distrib, rest_distrib, min_distrib, fan_out = 1024)
    # create the tree
    dt, is_parent = DtreeScheduler(fan_out, nwi, true, 1.0, first_distrib, rest_distrib, min_distrib)

    # ---
    if dt_nodeid == 1
        println("dtreebench -- ", dt_nnodes, " nodes")
        println("  system clock speed is ", cpu_hz, " GHz")
    end

    # roughly how many work items will each node will handle?
    each, r = divrem(nwi, dt_nnodes)
    if r > 0
        each = each + 1
    end

    # ---
    if dt_nodeid == 1
        println("  ", nwi, " work items, ~", each, " per node")
    end

    # generate random numbers for work item durations
    dura = Float64[]
    mn = repmat([meani-0.5*stddevi, meani+0.5*stddevi], ceil(Int, dt_nnodes/2))
    mt = MersenneTwister(7777777)
    for i = 1:dt_nnodes
        r = randn(mt)*stddevi*0.25+mn[i]
        append!(dura, max(randn(mt, each)*stddevi+r, zero(Float64)))
    end

    # ---
    if dt_nodeid == 1
        println("  initializing...")
    end

    # get the initial allocation
    ilock = SpinLock()
    ni, (ci, li) = initwork(dt)

    # ---
    if dt_nodeid == 1
        println("  ...done.")
    end

    # start threads and run, or run single-threaded
    if VERSION > v"0.5.0-dev"
        tfargs = Core.svec(threadfun, dt, Ref(ni), Ref(ci), Ref(li), ilock, runtree(dt), dura)
        ccall(:jl_threading_run, Void, (Any,), tfargs)
    else
        threadfun(dt, Ref(ni), Ref(ci), Ref(li), ilock, runtree(dt), dura)
    end

    # ---
    if dt_nodeid == 1
        println("complete")
    end
    tic()
    finalize(dt)
    wait_done = toq()
    ntputs(dt_nodeid, 1, "wait for done: $wait_done secs")
end

#bench(80, 0.5, 0.125, 0.2, 0.5, nthreads())
bench(100, 0.5, 0.125, 0.5, 0.5, nthreads())