text:
module C_constraints
    include(dirname(@__FILE__)*"\\Utils_cons.jl")
    include(dirname(dirname(@__FILE__))*"\\Moments\\Moments.jl")
    using .Utils_cons
    using .Moments
    using LinearAlgebra
    const la = LinearAlgebra
    const uc = Utils_cons
    const mom = Moments

    export make_mon_expo_keys,
           make_PSD_con,
           make_ord4_con,
           make_loc_cons_S_inf,
           make_loc_cons_S₂,
           make_loc_cons_S₂₁,
           make_Gᴿ_con

    make_mon_expo_keys(d,t::Int) = Moments.make_mon_expo(d,t*2)

    """L([xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜᵀ) ⪰ 0"""
    function make_PSD_con(d,t,Lx;noBlock=false)
        MMexᴿ,MMCoefᴿ = Moments.get_ℂ_block_diag(d,t;noBlock=noBlock)
        println(keys(MMexᴿ))
        return  Utils_cons.idx2var_dic(Lx,MMexᴿ,MMCoefᴿ)
    end

    """L(xx*⊗yy*) ⪰ 0 ⟺ """
    make_ord4_con(d,Lx) = uc.idx2varxx̄ᵀtyȳᵀ(Lx,mom.make_xx̄ᵀ_tens_yȳᵀ(d))

    """ Lᴿ(gᴿ⋅[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁ᵀ) ⪰ 0
        gᴿ = √(ρₘₐₓ) - ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ²) , i ∈ [d_1],
        or   √(ρₘₐₓ) - ((yᵣₑ)ⱼ² + (yᵢₘ)ⱼ²) , j ∈ [d_2] """
    function make_loc_con(Lx,sqρ,eᵤ,eᵥ,B,C)
        Lxᵢxⱼ = uc.idx2var_arr(Lx,B,C,[2*eᵤ]) + uc.idx2var_arr(Lx,B,C,[2*eᵥ]) #  L((xᵣₑ)ᵢ²+ (xᵢₘ)ᵢ²⋅ηₜ₋₁)
        return  sqρ*uc.idx2var_arr(Lx,B,C) - Lxᵢxⱼ  # L((√ρₘₐₓ - ((xᵣₑ)ᵢ²+(xᵢₘ)ᵢ²))⋅ηₜ₋₁)
    end
    function make_loc_cons_S_inf(ρ,d,t,Lx;noBlock=false)
        d₁,d₂ = d ; n = sum(2 .*d) ; sqρ   = sqrt(maximum(norm.(ρ)))
        MMexᴿ,MMCoefᴿ = mom.get_ℂ_block_diag(d,t .- 1;noBlock=noBlock)
        loc_con = Dict()
        for b in keys(MMCoefᴿ)
            for k in 1:d₁ # Constraint:  Lᴿ( (√ρₘₐₓ-((xᵣₑ)ᵢ²-(xᵢₘ)ᵢ²))⋅ηₜ₋₁)) ⪰ 0 for k ∈ [d₁]
                loc_con[(b,"x²_$k")] = make_loc_con(Lx,sqρ,uc.eᵢ(n,k),uc.eᵢ(n,k+d₁),MMexᴿ[b],MMCoefᴿ[b])
            end
            for k in 1:d₂ # Constraint:  Lᴿ( (√ρₘₐₓ-((yᵣₑ)ᵢ²+(yᵢₘ)ᵢ²))⋅ηₜ₋₁) ⪰ 0 for k ∈ [d₂]
                loc_con[(b,"y²_$k")] = make_loc_con(Lx,sqρ,uc.eᵢ(n,k+2*d₁),uc.eᵢ(n,k+2*d₁+d₂),MMexᴿ[b],MMCoefᴿ[b])
            end
        end
        return loc_con
    end

    """Lᴿ(gᴿ⋅[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁ᵀ) ⪰ 0
        gᴿ   = √Tr(ρ) - ∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ²),
                √Tr(ρ) - ∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ²)"""
    tmp(Lx,B,C,n,s,l) = sum([uc.idx2var_arr(Lx,B,C,[2* uc.eᵢ(n,k)] ) for k in s:(s+l)]) #### SOmething is wrong
    function make_loc_cons_S₂(ρ,d,t,Lx;noBlock=false)
        d₁,d₂ = d ; n = sum(2 .*d) ; sqrt_tr_ρ = sqrt(real(tr(ρ)))
        MMexᴿ,MMCoefᴿ = mom.get_ℂ_block_diag(d, t.- 1,noBlock=noBlock)

        loc_con = Dict(); g₂ = Dict()
        for b in keys(MMCoefᴿ)
            xRterm = tmp(Lx,MMexᴿ[b],MMCoefᴿ[b],n,1,d₁)    + tmp(Lx,MMexᴿ[b],MMCoefᴿ[b],n,d₁,d₁)       #  Lᴿ(∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            yRterm = tmp(Lx,MMexᴿ[b],MMCoefᴿ[b],n,2*d₁,d₂) + tmp(Lx,MMexᴿ[b],MMCoefᴿ[b],n,2*d₁+d₂,d₂)  #  Lᴿ(∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            loc_con[b,"x"] = sqrt_tr_ρ*uc.idx2var_arr(Lx,MMexᴿ[b],MMCoefᴿ[b]) - xRterm #  √Tr(ρ) ⋅ Lᴿ(ηₜ₋₁) - Lᴿ(∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            # loc_con[b,"y"] = sqrt_tr_ρ*uc.idx2var_arr(Lx,MMexᴿ[b],MMCoefᴿ[b]) - yRterm #  √Tr(ρ) ⋅ Lᴿ(ηₜ₋₁) - Lᴿ(∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            g₂[b] = xRterm - yRterm
        end
        return loc_con, g₂
    end

    """Lᴿ(gᴿ⋅[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁ᵀ) ⪰ 0
        g^ℝ = Tr(ρ) - ∑ᵈᵢ(uₓᵢ² + vₓᵢ²),
                ∑ᵈᵢ(u_yⱼ² + v_yⱼ²)  - 1"""
    function make_loc_cons_S₂₁(ρ,d,t,Lx;noBlock=false)
        d₁,d₂ = d
        n     = sum(2 .*d)
        # MB     = mom.get_ℂ_block_diag(d,t .- 1)
        MMexᴿ,MMCoefᴿ = mom.get_ℂ_block_diag(d, t.- 1,noBlock=noBlock)
        tr_ρ  = real(tr(ρ))

        loc_con    = Dict()
        loc_con_eq = Dict()
        for b in keys(MMCoefᴿ)
            xRterm = tmp(Lx,MMexᴿ[b],MMCoefᴿ[b],n,1,d₁)    + tmp(Lx,MMexᴿ[b],MMCoefᴿ[b],n,d₁,d₁)        #  L(∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            yRterm = tmp(Lx,MMexᴿ[b],MMCoefᴿ[b],n,2*d₁,d₂) + tmp(Lx,MMexᴿ[b],MMCoefᴿ[b],n,2*d₁+d₂,d₂)   #  L(∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            loc_con[b]    = tr_ρ*uc.idx2var_arr(Lx,MMexᴿ[b],MMCoefᴿ[b]) - xRterm       #  Tr(ρ)⋅L(ηₜ₋₁) - L(∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ ) ⪰ 0
            loc_con_eq[b] = uc.idx2var_arr(Lx, MMexᴿ[b],MMCoefᴿ[b]) - yRterm            #  1⋅L(ηₜ₋₁) -  L(∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ ) = 0
        end
        return loc_con, loc_con_eq
    end

    """ρ⊗L(η) - L( (xx*⊗yy*) ⊗ η )⪰ 0 ∀ η ∈ blocks of ([x,̄x,y,̄y]ₜ₋₂[x,̄x,y,̄y]*ₜ₋₂)
     L( ρℝ ⊗ η)  - L( Gℝ[k] ⊗ η)⪰ 0 ∀ η ∈ blocks of [uₓ,vₓ,u_y,v_y]ₜ₋₂[uₓ,vₓ,u_y,v_y]^Tₜ₋₂  """

    function make_Gᴿ_con(ρ,d,t,Lx;noBlock=false)
        n = sum(2 .*d) ; D = prod(d) ; ρᴿ = [real(ρ) -imag(ρ); imag(ρ) real(ρ)]
        Gᴿ,sm        = uc.get_Gᴿ(mom.make_xx̄ᵀ_tens_yȳᵀ(d)) #
        MMexᴿ,MMCoefᴿ = mom.get_ℂ_block_diag(d, t.- 2,noBlock=noBlock)

        LGℝη         = Dict()
        tmp2(B,C,k) = la.kron(sm[k],ones(D .*size(B))) .*
            uc.idx2var_arr(
                Lx,
                mom.var_kron_C(Gᴿ[k],B),
                mom.var_kron_C(fill(0.0,size(Gᴿ[k])...),C)) # L(Gᴿ[k]⋅ηₜ₋₂)
        for b in keys(MMCoefᴿ)
            TEMP1 = sum([tmp2(MMexᴿ[b],MMCoefᴿ[b],k) for k in 1:8]) # ∑ₖL(Gᴿ[k]⋅ηₜ₋₂)
            TEMP2 = la.kron(ρᴿ,uc.idx2var_arr(Lx,MMexᴿ[b],MMCoefᴿ[b]))            #
            LGℝη[b]   = TEMP2 - TEMP1 # L(ρᴿ ⊗ ηₜ₋₂)  - ∑ₖL(Gᴿ[k] ⊗ ηₜ₋₂)
        end
        return  LGℝη
    end
end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
# ======================================================================================================
# operations for balancing trees

# -----------------------------------------------------------------------------------------------------------------------
# nodes operations

opposite(dir::Symbol) = dir == :left ? (:right) : (:left)
opposite(::Nothing) = nothing ###

isnull(Node::AbstractNode) = false
isnull(Node::NullNode) = true

height(Node::AbstractNode) = 0
height(Node::HBN) = Node.height
height(Node::NullNode) = -1

"""
    link!(<parent>, <child>, <direction>)

Link parent and child in specific direction.
"""
function link!(parent::AbstractNode, child::AbstractNode, dir::Symbol)
    setfield!(parent, dir, child)
    setfield!(child, :parent, parent)
end

link!(parent::AbstractNode, child::NullNode, dir::Symbol) = setfield!(parent, dir, child)
link!(parent::NullNode, child::AbstractNode, ::Symbol) = setfield!(child, :parent, parent)
link!(parent::NullNode, child::AbstractNode, ::Nothing) = setfield!(child, :parent, parent)
link!(::NullNode, ::NullNode, ::Symbol) = nothing
link!(::NullNode, ::NullNode, ::Nothing) = nothing

"""
    cut!(<parent>, <child>, <direction>)
    cut!(<parent>, <child>)

Cut parent and child.
"""
function cut!(parent::AbstractNode, child::AbstractNode, dir::Symbol)
    setfield!(parent, dir, NullNode())
    setfield!(child, :parent, NullNode())
end

cut!(parent::AbstractNode, child::AbstractNode) = cut!(parent, child, getdir(parent, child))
cut!(::AbstractNode, ::NullNode, ::Symbol) = nothing
cut!(::NullNode, ::AbstractNode, ::Symbol) = nothing

"""
    getdir(<parent>, <child>)

Get direction of parent to child.
"""
getdir(parent::AbstractBinaryNode, child::AbstractBinaryNode) = parent.left == child ? (:left) : (:right)
getdir(::NullNode, ::AbstractNode) = nothing
getdir(::AbstractNode, ::NullNode) = nothing
getdir(::NullNode, ::NullNode) = nothing

# -------------------------------------------------------------------------------------------------
# rotations

function SingleRotation!(grandparent::HBN,
                        parent::HBN,
                        child::HBN,
                        dir1::Symbol)
    dir2 = opposite(dir1)
    sister = getproperty(parent, dir2)
    ggparent = grandparent.parent
    link!(grandparent, sister, dir1)
    link!(parent, grandparent, dir2)
    dir3 = getdir(ggparent, grandparent)
    link!(ggparent, parent, dir3)
    grandparent.height -= 1
    ggparent, parent, child, dir3
end

function DoubleRotation!(grandparent::HBN,
                        parent::HBN,
                        child::HBN,
                        dir1::Symbol,
                        dir2::Symbol)
    ggparent = grandparent.parent
    outerchild = getproperty(child, dir1)
    innerchild = getproperty(child, dir2)
    link!(parent, outerchild, dir2)
    link!(grandparent, innerchild, dir1)
    link!(child, parent, dir1)
    link!(child, grandparent, dir2)
    dir3 = getdir(ggparent, grandparent)
    link!(ggparent, child, dir3)
    child.height += 1
    parent.height -= 1
    grandparent.height -= 1
    ggparent, child, parent, dir3
end

function Zig!(grandparent::SBN,
            parent::SBN,
            dir1::Symbol)
    dir2 = opposite(dir1)
    child = getproperty(parent, dir2)
    link!(grandparent, child, dir1)
    link!(parent, grandparent, dir2)
    parent.parent = NullNode()
    parent
end

function ZigZig!(grandparent::SBN,
                parent::SBN,
                child::SBN,
                dir1::Symbol,
                dir2::Symbol)
    sister = getproperty(parent, dir2)
    grandchild = getproperty(child, dir2)
    ggparent = grandparent.parent
    dir3 = getdir(ggparent, grandparent)
    link!(grandparent, sister, dir1)
    link!(parent, grandparent, dir2)
    link!(parent, grandchild, dir1)
    link!(child, parent, dir2)
    link!(ggparent, child, dir3)
    ggparent, dir3
end

function ZigZag!(grandparent::SBN,
                parent::SBN,
                child::SBN,
                dir1::Symbol,
                dir2::Symbol)
    grandchild1 = getproperty(child, dir1)
    grandchild2 = getproperty(child, dir2)
    ggparent = grandparent.parent
    dir3 = getdir(ggparent, grandparent)
    link!(parent, grandchild1, dir2)
    link!(child, parent, dir1)
    link!(grandparent, grandchild2, dir1)
    link!(child, grandparent, dir2)
    link!(ggparent, child, dir3)
    ggparent, dir3
end

function TopDownZig!(parent::SBN,
                    child::SBN,
                    dir::Symbol)
    cut!(parent, child, dir)
    dir == :left ? (child, NullNode(), parent) : (child, parent, NullNode())
end

function TopDownZigZig!(grandparent::SBN,
                    parent::SBN,
                    child::SBN,
                    dir1::Symbol,
                    dir2::Symbol)
    sister = getproperty(parent, dir2)
    cut!(grandparent, parent, dir1)
    cut!(parent, child, dir1)
    link!(parent, grandparent, dir2)
    link!(grandparent, sister, dir1)
    dir1 == :left ? (child, NullNode(), parent) : (child, parent, NullNode())
end

function TopDownZigZag!(grandparent::SBN,
                    parent::SBN,
                    child::SBN,
                    dir1::Symbol,
                    dir2::Symbol)
    cut!(grandparent, parent, dir1)
    cut!(parent, child, dir2)
    dir1 == :left ? (child, parent, grandparent) : (child, grandparent, parent)
end

# -----------------------------------------------------------------------------------------------
# BinarySearchTree traverse

function traverse!(tree::BST,
                node::HBN,
                nodes...)
    value = max(-1, height(node.left), height(node.right)) + 1
    if node.height != value
        node.height = value
        return traverse!(tree, node.parent, node)
    end
    tree
end

traverse!(tree::BST,
        ::SBN,
        nodes...) = tree

function traverse!(tree::BST,
                ::NullNode,
                nodes...)
    tree.height = first(nodes).height
    tree
end

# --------------------------------------------------------------------------------------------------
# AVL tree traverse

## Root case
function traverse!(tree::AVLTree,
                ::NullNode,
                parent::HBN,
                nodes...)
    tree.height = parent.height
    tree
end

# traverse back to root
function AVLtraverse!(tree::AVLTree,
                    ::NullNode,
                    parent::HBN,
                    nodes...)
    tree.root = parent
    tree.height = parent.height
    tree
end

## insert! case
function traverse!(tree::AVLTree,
                parent::HBN,
                child::HBN)
    if parent.height == 0
        parent.height = 1
        return AVLtraverse!(tree, parent.parent, parent, child, getdir(parent.parent, parent), getdir(parent, child))
    end
    tree
end

function AVLtraverse!(tree::AVLTree,
                    grandparent::HBN,
                    parent::HBN,
                    child::HBN,
                    dir1::Symbol,
                    dir2::Symbol)
    delta = [child.height for child in grandparent]
    append!(delta, [-1, -1])
    if abs(delta[1] - delta[2]) > 1
        # AVL property violation
        if dir1 == dir2
            grandparent, parent, child, dir3 = SingleRotation!(grandparent, parent, child, dir1)
        else
            grandparent, parent, child, dir3 = DoubleRotation!(grandparent, parent, child, dir1, dir2)
        end
        AVLtraverse!(tree, grandparent, parent, child, dir3, dir1)
    else
        value = maximum(delta) + 1
        if grandparent.height != value
            # upper node is affected
            grandparent.height = value
            ggparent = grandparent.parent
            return AVLtraverse!(tree, ggparent, grandparent, parent, getdir(ggparent, grandparent), dir1)
        end
        tree
    end
end

## delete! case
function traverse!(tree::AVLTree, grandparent::HBN)
    parent = grandparent.left
    aunt = grandparent.right
    if height(parent) > height(aunt)
        AVLtraverse!(tree, grandparent, aunt, :left, :right)
    else
        AVLtraverse!(tree, grandparent, parent, :right, :left)
    end
end

function AVLtraverse!(tree::AVLTree,
                    grandparent::HBN,
                    aunt::Union{NullNode, HBN},
                    dir1::Symbol,
                    dir2::Symbol)
    parent, child1, child2 = findchild(grandparent, dir1)
    if abs(height(parent) - height(aunt)) > 1
        # AVL property violation
        grandparent.height -= 1
        if height(child1) > height(child2)
            grandparent, parent, child, dir3 = SingleRotation!(grandparent, parent, child1, dir1)
        elseif height(child1) < height(child2)
            grandparent, parent, child, dir3 = DoubleRotation!(grandparent, parent, child2, dir1, dir2)
        else
            grandparent, parent, child, dir3 = DoubleRotation!(grandparent, parent, child2, dir1, dir2)
            parent.height += 1
            child.height += 1
            if isnull(grandparent)
                tree.root = parent
                tree.height = parent.height
            end
            return tree
        end
        AVLtraverse!(tree, grandparent, parent, opposite(dir3), dir3)
    else
        value = max(height(parent), height(aunt)) + 1
        if grandparent.height != value
            # upper node is affected
            grandparent.height = value
            ggparent = grandparent.parent
            dir3 = getdir(ggparent, grandparent)
            return AVLtraverse!(tree, ggparent, grandparent, opposite(dir3), dir3)
        end
        tree
    end
end

# ------------------------------------------------------------------------------------
# Splay tree traverse
function Splaytraverse!(tree::SplayTree,
                    parent::SBN,
                    child::SBN,
                    dir2::Symbol)
    grandparent = parent.parent
    if isnull(grandparent)
        # to the root
        Zig!(parent, child, dir2)
        tree.root = child
        tree
    else
        dir1 = getdir(grandparent, parent)
        if dir1 == dir2
            parent, dir3 = ZigZig!(grandparent, parent, child, dir1, opposite(dir1))
        else
            parent, dir3 = ZigZag!(grandparent, parent, child, dir1, dir2)
        end
        Splaytraverse!(tree, parent, child, dir3)
    end
end

# Root case
function Splaytraverse!(tree::SplayTree,
                    ::NullNode,
                    child::SBN,
                    ::Nothing)
    tree.root = child
    tree
end

ext: jl


text:
using DifferentialEquations, thesis

function f′(du, u, p, t)
    du .= -p.α .* u
    du[1] += p.input(t)
    du[2:end] .+= p.α*u[1:end-1]
end

p = (
    input=t-> 0.5*(sin(t)) + 0.3*(cos(0.7+1.5*t)),
    α = 2.0,
    δ = zeros(5),
    spiketimes = [Float64[] for i ∈ 1:5]
)
tspan = (0.0,10.0)
u₀ = copy(p.δ)
cb = CallbackSet((ContinuousCallback((u,t,integrator) -> u[i]-integrator.p.δ[i], integrator->push!(integrator.p.spiketimes[i], integrator.t)) for i ∈ eachindex(p.δ))...)

sol = solve(ODEProblem(f′, u₀, tspan.-(10.0,0), p); callback=cb)

p_in = plot(sol, tspan=tspan, xlims=tspan.+(0,0.2), xlabel=false, ylabel="signal embedding")
plot!(p.input, tspan..., color=:black, label="s(t)")
hline!([p.δ], color=thesis_colors[end], linestyle=:dash, label="threshold")

p_spike = plot(xlabel="time", ylabel="population code")
for (i,spiketrain) ∈ enumerate(p.spiketimes)
    spiketrain = filter(s->tspan[1]<=s<=tspan[2], spiketrain)
    scatter!(p_in, spiketrain, fill(p.δ[i],length(spiketrain)), markersize=3, markercolor=:black, markerstrokecolor=i, markerstrokewidth=2, label="")
    plot!(p_spike, [spiketrain'; spiketrain'], [fill(i,length(spiketrain))';fill(i+0.75,length(spiketrain))'], linewidth=2, color=i, label="")
    scatter!(p_spike, spiketrain, fill(i+0.75,length(spiketrain)), markersize=3,  markercolor=:black, markerstrokecolor=i, markerstrokewidth=2, label="")
end
plot!(p_spike, xlims=tspan.+(0,0.2), ylims=[1.0,6.25], yticks=1:5, label=false)

plot(p_in, p_spike, layout=grid(2,1, heights=[0.7,0.3]))

savefig("figures/embedding.pdf")
savefig("figures/embedding.svg")

max_stars_repo_name: jleugeri/phd

max_issues_repo_name: jleugeri/phd

max_forks_repo_name: jleugeri/phd


text:
using LeastSquaresOptim, LinearAlgebra

tests = ["nonlinearsolvers.jl", "nonlinearleastsquares.jl", "nonlinearfitting.jl", "bounds.jl"]

println("Running tests:")

for test in tests
	try
		include(test)
		println("\t\033[1m\033[32mPASSED\033[0m: $(test)")
	 catch e
	 	println("\t\033[1m\033[31mFAILED\033[0m: $(test)")
	 	showerror(stdout, e, backtrace())
	 	rethrow(e)
	 end
end

# test README
function rosenbrock(x)
	[1 - x[1], 100 * (x[2]-x[1]^2)]
end
x0 = zeros(2)
optimize(rosenbrock, x0, Dogleg())
optimize(rosenbrock, x0, LevenbergMarquardt())

function rosenbrock_f!(out, x)
 out[1] = 1 - x[1]
 out[2] = 100 * (x[2]-x[1]^2)
end
optimize!(LeastSquaresProblem(x = zeros(2), f! = rosenbrock_f!, output_length = 2, autodiff = :central), Dogleg())

# if you want to use gradient
function rosenbrock_g!(J, x)
    J[1, 1] = -1
    J[1, 2] = 0
    J[2, 1] = -200 * x[1]
    J[2, 2] = 100
end
optimize!(LeastSquaresProblem(x = zeros(2), f! = rosenbrock_f!, g! = rosenbrock_g!, output_length = 2), Dogleg())

max_stars_repo_name: matthieugomez/LeastSquares

max_issues_repo_name: matthieugomez/LeastSquares

max_forks_repo_name: matthieugomez/LeastSquares


text:
const MOD_REPORT_CHANNEL_NAME = "mod-report"
const MOD_REPORT_CHANNEL = Ref{DiscordChannel}()

const MOD_BAD_WORDS = Set{BadWord}()
const MOD_BAD_WORD_REGEXES = Dict{BadWord,Regex}()

function handler(
    c::Client,
    e::Union{MessageCreate, MessageUpdate},
    ::Val{:mod}
)
    if ismissing(e.message.content)
        return # Possible e.g. only embed is updated by Discord (flipping discourse)
    end

    if e.message.author.id == c.state.user.id
        return # No need to check if the message came from the bot itself
    end

    mod_report_channel = mod_get_report_channel(c, e.message.guild_id)
    if mod_report_channel === nothing
        return # config error; nothing we can do about...
    end
    if e.message.channel_id == mod_report_channel.id
        return # Do not moderate the mod channel itself
    end

    mod_init() # lazy initialization
    result = mod_check_message(e.message.content)
    if !isempty(result)
        user_id = e.message.author.id
        report_message_id = e.message.id
        # Censor message
        new_content = mod_censor_message(e.message.content, result)
        if new_content != e.message.content
            delete_message(c, e.message.channel_id, e.message.id)
            censored_message = @discord create_message(c, e.message.channel_id;
                content = "<@!$(user_id)> said: $new_content")
            report_message_id = censored_message.id
        end
        # Update the mod-report channel
        report = mod_report(user_id, e.message.content, result,
            report_message_id, e.message.channel_id, e.message.guild_id)
        if mod_report_channel !== nothing
            create_message(c, mod_report_channel.id; content = report)
        end
    end
    return nothing
end

"Get the mod-report channel."
function mod_get_report_channel(c::Client, guild_id::Integer)
    if !isassigned(MOD_REPORT_CHANNEL)
        channels = @discord get_guild_channels(c, guild_id)
        channel = filter(c -> c.name == MOD_REPORT_CHANNEL_NAME, channels)
        if length(channel) == 1
            MOD_REPORT_CHANNEL[] = channel[1]
            return MOD_REPORT_CHANNEL[]
        else
            @error "Unable to find mod report channel: $MOD_REPORT_CHANNEL_NAME"
            return nothing
        end
    end
    return MOD_REPORT_CHANNEL[]
end

"Initialize mod function e.g. read bad words list."
function mod_init(; force = false)
    if force
        empty!(MOD_BAD_WORDS)
    end
    if isempty(MOD_BAD_WORDS)
        bad_words_file = joinpath(@__DIR__, "..", "..", "config", "mod", "bad-words.txt")
        bad_words = readlines(bad_words_file)
        for word in bad_words
            if startswith(word, '?')
                mod_add_bad_word(BadWord(Questionable, word[2:end]))
            elseif startswith(word, '-')
                mod_add_bad_word(BadWord(Overridden, word[2:end]))
            else
                mod_add_bad_word(BadWord(Bad, word))
            end
        end
    end
end

"Create a regex that matches a word."
function mod_make_regex(word::AbstractString)
    return Regex(string(
        "\\b",           # any word boundary
        "(",
        "\\Q",           # parse symbols as-is with \Q and \E marker
        word,
        "\\E",
        ")",
        "\\b",           # any word boundary
    ), "i")
end

"Add a new bad word to the global list."
function mod_add_bad_word(w::BadWord)
    try
        MOD_BAD_WORD_REGEXES[w] = mod_make_regex(w.word)
        push!(MOD_BAD_WORDS, w)
    catch ex
        @warn "Cannot compile regex for word" w.word regex_str
    end
end

"Return true if `content` contains `word` considering word boundaries."
function mod_contains(content::AbstractString, w::BadWord)
    regex = MOD_BAD_WORD_REGEXES[w]
    return match(regex, content) !== nothing
end

"""
    mod_check_message(content::AbstractString)

Check the message content and return a set
of identified bad words. An empty set is returned if
nothing is found.
"""
function mod_check_message(content::AbstractString)
    # split by spoiler marker
    # after this, every even-number token is within the spoiler markers
    tokens = split(content, "||")

    result = Set{BadWord}()
    for (i, t) in enumerate(tokens)
        if isodd(i)
            bad_word_set = mod_check(t)
            union!(result, bad_word_set)
        end
    end
    return result
end

"""
    mod_check(content::AbstractString)

Check the content against a profanity set of words.
"""
function mod_check(content::AbstractString)
    return Set{BadWord}(w for w in MOD_BAD_WORDS if mod_contains(content, w))
end

"""
    mod_report

Return a string about the issue that will be sent to the mod-report channel.
"""
function mod_report(
    user_id::Integer,
    content::AbstractString,
    bad_words::AbstractSet{BadWord},
    message_id::Integer,
    channel_id::Integer,
    guild_id::Integer,
)
    str = join([w.word for w in bad_words], ", ")
    badness = mod_badness(bad_words)
    return "$badness message from <@!$user_id> with `$str`: $(content)\n" *
        "Ref: https://discord.com/channels/$guild_id/$channel_id/$message_id"
end

"""
    mod_badness(bad_words::AbstractSet{BadWord})

Return a string that describe how bad the message is based
upon the set of bad words already found from the message.
"""
function mod_badness(bad_words::AbstractSet{BadWord})
    return if any(x -> x.class === Bad, bad_words)
        "Bad"
    elseif any(x -> x.class === Questionable, bad_words)
        "Questionable"
    elseif any(x -> x.class === Overridden, bad_words)
        "Overridden"
    else
        "Unknown"
    end
end

"""
    mod_censor_message(content, bad_words)

Returned a censored string for the message. For example,
bad words may be hidden using spoiler tags. Text that are within
Discord spoiler markers are ignored. Only `bad_words` are
considered. See `mod_check_message`[@ref] about how to gather
the list.
"""
function mod_censor_message(
    content::AbstractString,
    bad_words::AbstractSet{BadWord}
)
    original_content = content

    tokens = split(content, "||")
    buf = IOBuffer()
    for (i, token) in enumerate(tokens)
        if isodd(i)
            token = mod_censor(token, bad_words)
        else
            token = "||" * token * "||"
        end
        print(buf, token)
    end
    content = String(take!(buf))

    @info "Censored: original=$original_content new=$content"
    return content
end

"""
    mod_censor(content, bad_words)

If the `content` contains any bad word from the list, then it is
hidden by Discord spoiler markers. See `mod_check_message`[@ref]
about how to gather the list.
"""
function mod_censor(content::AbstractString, bad_words::AbstractSet{BadWord})
    for w in bad_words
        regex = MOD_BAD_WORD_REGEXES[w]
        if w.class in (Bad, Questionable)
            content = replace(content, regex => s"||\1||")
        end
    end
    return content
end


text:
export iauCp
"""
Copy a p-vector.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   p        double[3]     p-vector to be copied

Returned:
   c        double[3]     copy

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauCp(double p[3], double c[3])

function iauCp(p::AbstractVector{<:Real})

   # Allocate return value
   c = zeros(Float64, 3)

   ccall((:iauCp, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}),
         convert(Array{Float64, 1}, p), c)

   return SVector{3}(c)
end


text:
# This file is a part of BayesianBlocks.jl, licensed under the MIT License (MIT).

import Test

Test.@testset "Package BayesianBlocks" begin
    include("test.jl")
end # testset


text:
# # using SymPy

# x,y,t = @vars x y t
# u = exp(-t)*x*(1-x)*y*(1-y)
# v = u
# ux = diff(u, x)
# uy = diff(u, y)
# vx = diff(v, x)
# vy = diff(v, y)
# e = -[ux
#     vy
#     uy+vx
#     ]
# σ = Cs * e
# f1 = -(diff(σ[1], x) + diff(σ[3], y)) + u
# print(replace(replace(sympy.julia_code(simplify(f1)), ".^"=>"^"), ".*"=>"*"))

# f2 = -(diff(σ[3], x) + diff(σ[2], y)) + u
# print(replace(replace(sympy.julia_code(simplify(f2)), ".^"=>"^"), ".*"=>"*"))

@testset "viscous" begin
    domain = example_domain(10, 10, 0.1)

    EBC = zeros(Int64, domain.nnodes, 2)
    FBC = zeros(Int64, domain.nnodes, 2)
    g = zeros(domain.nnodes, 2)
    f = zeros(domain.nnodes, 2)
    for i = 1:domain.nnodes
        if domain.nodes[i,1]<1e-5 || domain.nodes[i,2]<1e-5 || domain.nodes[i,1]>1-1e-5 || domain.nodes[i,2]>1-1e-5
            EBC[i,:] .= -1
        end
    end

    Cs = [3. 2. 0.0
        2. 5. 0.0
        0.0 0.0 1.0]
    Hs = zeros(3, 3)

    domain = Domain(domain.nodes, domain.elements, 2, EBC, g, FBC, f)
    globdat = example_global_data(domain)
    globdat.Body_func = (x, y, t)->begin
        f1 = @. 1.0*(1.0*x^2*y^2 - 1.0*x^2*y + 2.0*x^2 - 1.0*x*y^2 + 13.0*x*y - 8.0*x + 6.0*y^2 - 12.0*y + 3.0)*exp(-t)
        f2 = @. 1.0*(1.0*x^2*y^2 - 1.0*x^2*y + 10.0*x^2 - 1.0*x*y^2 + 13.0*x*y - 16.0*x + 2.0*y^2 - 8.0*y + 3.0)*exp(-t)
        out = [f1 f2]
    end
    Δt = 0.01
    NT = 100

    ts = GeneralizedAlphaSolverTime(Δt, NT)
    ubd, abd = compute_boundary_info(domain, globdat, ts)
    Fext = compute_external_force(domain, globdat, ts)

    # xy = reshape(domain.nodes[domain.dof_to_eq], sum(domain.dof_to_eq)÷2, 2)
    x, y = domain.nodes[:,1], domain.nodes[:,2]
    d0 = @. x*(1-x)*y*(1-y)
    v0 = - @. x*(1-x)*y*(1-y)
    a0 = @. x*(1-x)*y*(1-y)
    d0 = [d0;d0]
    v0 = [v0;v0]
    a0 = [a0;a0]

    d, v, a = GeneralizedAlphaSolver(globdat, domain, d0, v0, a0, Δt, NT, Cs, Hs, Fext, ubd, abd)

    sess = Session(); init(sess)
    d_ = run(sess, d)

    uexact = zeros(NT+1, domain.nnodes*2)
    for i = 1:NT+1
        t = (i-1)*Δt
        uexact[i,:] = [(@. exp(-t) * x*(1-x)*y*(1-y));(@. exp(-t) * x*(1-x)*y*(1-y))]
    end

    uexact - d_

end


text:
using Test
using UltraDark: Grids, PencilGrids, output_grids, OutputConfig

for grid_type in [Grids, ]
    grids = grid_type(1.0, 16)

    for (name, output_config) in [
              ("default", OutputConfig(mktempdir(), [])),
              ("nothing", OutputConfig(mktempdir(), []; box=false, slice=false, psi=true, rho=true)),
              ("box", OutputConfig(mktempdir(), []; box=true, slice=false, psi=true, rho=true)),
              ("slice", OutputConfig(mktempdir(), []; box=true, slice=true, psi=true, rho=true)),
              ("box & slice", OutputConfig(mktempdir(), []; box=true, slice=true, psi=true, rho=true)),
                         ]
        mkpath(output_config.directory)
        @test output_grids(grids, output_config, 0) == nothing
    end
end

for grid_type in [PencilGrids, ]
    grids = grid_type(1.0, 16)

    for (name, output_config) in [
              ("default", OutputConfig(mktempdir(), [])),
              ("nothing", OutputConfig(mktempdir(), []; box=false, slice=false, psi=true, rho=true)),
              ("box", OutputConfig(mktempdir(), []; box=true, slice=false, psi=true, rho=true)),
                         ]
        mkpath(output_config.directory)
        @test output_grids(grids, output_config, 0) == nothing
    end

    for (name, output_config) in [
              ("slice", OutputConfig(mktempdir(), []; box=true, slice=true, psi=true, rho=true)),
              ("box & slice", OutputConfig(mktempdir(), []; box=true, slice=true, psi=true, rho=true)),
                         ]
        mkpath(output_config.directory)
        @test_broken output_grids(grids, output_config, 0) == nothing
    end
end


text:
using StatsBase
using Roots

# include("functions_align.jl")

# reduced_is_same restricts comparisons to sites that are different
# from the reference in either sequence being compared.
# NB this could be changed to the consensus sequence for the alignment to avoid
# reliance on an another file
function reduced_is_same(nuc_bit_array, ref_array)
    same_array_size = size(nuc_bit_array, 2)

    same_array = trues(same_array_size, same_array_size)

    difference_list = get_list_of_differences(nuc_bit_array, ref_array)

    for a in 1:size(nuc_bit_array,2) - 1

        for b in (a + 1):size(nuc_bit_array,2)

            start_again = true

            for r in difference_list[a]
                @inbounds x = nuc_bit_array[r,a]
                @inbounds y = nuc_bit_array[r,b]
                different = (x & y) < 16

                # store the different True/False result in an array:
                if different
                    same_array[a,b] = !different
                    same_array[b,a] = !different
                    start_again = false
                    break
                end
            end

            # need to check the second seq's sites too, if we get this far:
            if start_again
                for r in difference_list[b]
                    if r in difference_list[a]
                        continue
                    end

                    @inbounds x = nuc_bit_array[r,a]
                    @inbounds y = nuc_bit_array[r,b]
                    different = (x & y) < 16

                    # store the different True/False result in an array:
                    if different
                        same_array[a,b] = !different
                        same_array[b,a] = !different
                        break
                    end
                end
            end
        end
    end

    return same_array
end

# # compare sites (rows) until a difference is found, for all pairs of sequences (columns)
# # - break as soon as you know they're different - don't need any more information from this pair
# function is_same(nuc_bit_array)
#     same_array_size = size(nuc_bit_array, 2)
#
#     same_array = trues(same_array_size, same_array_size)
#
#     println("number of threads for sequence comparisons: ", Threads.nthreads())
#
#     for a in 1:size(nuc_bit_array,2) - 1
#         for b in (a + 1):size(nuc_bit_array,2)
#             for r in Iterators.reverse(1:size(nuc_bit_array,1))
#                 @inbounds x = nuc_bit_array[r,a]
#                 @inbounds y = nuc_bit_array[r,b]
#                 different = (x & y) < 16
#
#                 # store the different True/False result in an array:
#                 if different
#                     same_array[a,b] = !different
#                     same_array[b,a] = !different
#                     break
#                 end
#             end
#         end
#     end
#
#     return same_array
# end

# Functions for parallelising
# ---------------------------

# This to provide the chunks to chop up the n*n matrix into approximately equal
# amounts of work per thread when we parallelise it.

# Because we are effectively only working on one triangle of a square matrix,
# the first rows have many more columns to fill, so that giving threads
# equal numbers of rows to populate means very unequal amounts of work, per thread.

# Instead, we chop up the triangle into equal-area trapezoids.

# For a given sample size (n) and quantile (q), generate an expression that
# we can find the root for (solve for x) - this will provide the row number
# that splits a right-angled triangle of area n^2 / 2 into quantiles q & 1-q by
# area:
function make_equation(n, q)
    f = x-> x * n  - ((x^2 + x) / 2) - ((n^2 - n) / 2) * q
    return f
end

# get many splits:
function get_chunk_splits(n, q_Array)
    chunk_splits = Array{Int, 1}(undef, length(q_Array))
    for (i, q) in enumerate(q_Array)
        f = make_equation(n, q)
        split = find_zero(f, (0, n), Bisection())
        chunk_splits[i] = floor(split)
    end

    return chunk_splits
end

# turn the splits into chunks that are ranges to iterate over:
function get_ranges(n, N_threads)
    q_Array = Array{Float64, 1}(undef, N_threads - 1)
    for i in Iterators.reverse(1:N_threads - 1)
        q_Array[i] = i/N_threads
    end

    chunk_splits = get_chunk_splits(n, q_Array)

    ranges = Array{Tuple{Int, Int}, 1}(undef, N_threads)
    previous_finish = 0
    for (i,split) in enumerate(chunk_splits)
        start = previous_finish + 1
        finish = split
        ranges[i] = (start, finish)
        previous_finish = finish
    end
    ranges[lastindex(ranges)] = (previous_finish + 1, n - 1)

    return ranges
end

function parallel_reduced_is_same(nuc_bit_array, ref_array)

    same_array_size = size(nuc_bit_array, 2)

    same_array = Array{Bool,2}(undef, same_array_size, same_array_size)

    for i in 1:size(same_array, 2)
        same_array[i,i] = true
    end

    # make as many chunks as we have threads
    n_Chunks = Threads.nthreads()

    # ...unless we have too many threads for this many sequences, in which case, down-size.
    # A chunk size less than or equal to half the total number of sequences makes valid ranges
    if n_Chunks > size(nuc_bit_array, 2) / 2
        n_Chunks::Int = floor(size(nuc_bit_array, 2) / 2)
    end

    ranges = get_ranges(size(nuc_bit_array, 2), n_Chunks)

    difference_list = get_list_of_differences(nuc_bit_array, ref_array)

    Threads.@threads for range in ranges
        for a in range[1]:range[2]
            for b in (a + 1):size(nuc_bit_array,2)

                start_again = true
                all_same = true

                for r in difference_list[a]
                    @inbounds x = nuc_bit_array[r,a]
                    @inbounds y = nuc_bit_array[r,b]
                    different = (x & y) < 16

                    # store the different True/False result in an array:
                    if different
                        all_same = false
                        start_again = false
                        break
                    end
                end

                # need to check the second seq's sites too, if we get this far:
                if start_again
                    for r in difference_list[b]
                        if r in difference_list[a]
                            continue
                        end

                        @inbounds x = nuc_bit_array[r,a]
                        @inbounds y = nuc_bit_array[r,b]
                        different = (x & y) < 16

                        # store the different True/False result in an array:
                        if different
                            all_same = false
                            break
                        end
                    end
                end

                same_array[a,b] = all_same
                same_array[b,a] = all_same
            end
        end
    end

    return same_array
end

# # A parallel version of is_same()
# function parallel_is_same(nuc_bit_array)
#     same_array_size = size(nuc_bit_array, 2)
#
#     # initialises as falses so have to change this:
#     # same_array = SharedArray{Bool,2}((same_array_size,same_array_size))
#     same_array = Array{Bool,2}(undef, same_array_size, same_array_size)
#
#     for i in 1:size(same_array, 2)
#         same_array[i,i] = true
#     end
#
#     # make as many chunks as we have threads
#     n_Chunks = Threads.nthreads()
#
#     # ...unless we have too many threads for this many sequences, in which case, down-size.
#     # A chunk size less than or equal to half the total number of sequences makes valid ranges
#     if n_Chunks > size(nuc_bit_array, 2) / 2
#         n_Chunks::Int = floor(size(nuc_bit_array, 2) / 2)
#     end
#
#     ranges = get_ranges(size(nuc_bit_array, 2), n_Chunks)
#
#     println("number of threads for sequence comparisons: ", Threads.nthreads())
#     # println("number of chunks: ", length(ranges))
#
#     Threads.@threads for range in ranges
#         for a in range[1]:range[2]
#             for b in (a + 1):size(nuc_bit_array,2)
#                 # comb_test+=1
#                 all_same = true
#                 for r in Iterators.reverse(1:size(nuc_bit_array,1))
#                     @inbounds x = nuc_bit_array[r,a]
#                     @inbounds y = nuc_bit_array[r,b]
#                     different = (x & y) < 16
#
#                     # store the different True/False result in an array:
#                     if different
#                         all_same = false
#                         break
#                     end
#                 end
#                 same_array[a,b] = all_same
#                 same_array[b,a] = all_same
#             end
#         end
#     end
#
#     return same_array
# end

# end of functions for parallelising
# ---------------------------

function get_one_set_from_view(bool_view)
    # one subset of rows of the view
    A = Array{Int64,1}()
    original_row_numbers = bool_view.indices[1]

    for (index, value) in enumerate(bool_view)
        if value
            push!(A, original_row_numbers[index])
        end
    end
    return Set(A)
end

function get_sets_from_view(bool_view)
    S = Set{Set}()
    for i in 1:size(bool_view, 2)
        push!(S, get_one_set_from_view(view(bool_view, :, i)))
    end
    return collect(S)
end

function check_the_view_is_good(a_bool_array)
    # because these logical arrays have symmetry only need
    # to check one triangle, but haven't here
    for i in a_bool_array
        if !i
            return false
        end
    end
    return true
end

# function is_subset(s1, s2)
#     test = true
#     for x in s1
#         if ! (x in s2)
#             test = false
#             break
#         end
#     end
#     return test
# end

# # test all sets for being supersets, and throw them out if they are
# # NB - this is 1/10 of the run time and allocates a lot of memory in total
# function get_subsets(array_of_sets)
#
#     A = Array{BitSet,1}()
#     sizehint!(A, length(array_of_sets))
#
#     for i in 1:size(array_of_sets, 1)
#         valid = true
#         for j in 1:size(array_of_sets, 1)
#             if i == j
#                 continue
#             end
#             # if array_of_sets[i] is a superset of comparison set, then break, else push to A
#             # (note the reverse logic here, because there is no superset function)
#             if issubset(array_of_sets[j], array_of_sets[i]) && array_of_sets[i] != array_of_sets[j]
#                 valid = false
#                 break
#             end
#         end
#
#         if valid
#             push!(A, array_of_sets[i])
#         end
#     end
#
#     return A
# end

# a parallel version of get_subsets()
function parallel_get_subsets(array_of_sets)

    A = Array{BitSet,1}()
    sizehint!(A, length(array_of_sets))

    tests = Array{Bool, 1}(undef, length(array_of_sets))

    for i in eachindex(tests)
        tests[i] = true
    end

    # NB - dont need to compare all the sets with all the other sets because the only sets that can
    # possibly be subsets of the tested set are those that share an element with it - so just iterate
    # over the focal set's rows. But can't do this as it stands because array_of_sets is a view of the
    # sets so its indices don't match the numbers in each set
    Threads.@threads for i in 1:size(array_of_sets, 1)
        for j in 1:size(array_of_sets, 1)
            if i == j
                continue
            end
            # if array_of_sets[i] is a superset of comparison set, then break, else push to A
            # (note the reverse logic here, because there is no superset function)
            if issubset(array_of_sets[j], array_of_sets[i]) && array_of_sets[i] != array_of_sets[j]
                tests[i] = false
                break
            end
        end
    end

    for (i, test) in enumerate(tests)
        if test
            push!(A, array_of_sets[i])
        end
    end

    return A
end

function get_one_truefalse_set_from_vector(bool_vector)
    # get the indices of Trues in bool_vector as a BitSet
    # bool vector must be 1-dimensional
    T = Array{Int64,1}()
    F = Array{Int64,1}()

    for i in 1:length(bool_vector)
        if bool_vector[i]
            push!(T, i)
        else
            push!(F, i)
        end
    end

    BT = BitSet(T)
    BF = BitSet(F)

    return BT, BF
end

function get_truefalse_row_indices(bool_array)
    # will return two arrays of BitSets which are:
    # 1) the rows in bool_array that are true, for every column in bool_array - matches between pairs of sequences
    # 2) the rows in bool_array that are false, for every column in bool_array - mismatches between pairs of sequences

    AT = Array{BitSet,1}(undef, size(bool_array, 2))
    AF = Array{BitSet,1}(undef, size(bool_array, 2))

    for i in 1:size(bool_array, 2)
        AT[i], AF[i] = get_one_truefalse_set_from_vector(view(bool_array, :, i))
    end

    return AT, AF
end

function get_sets_in_one_go(the_whole_bool_array)
    # returns an array of bitsets which are the sets

    # println("number of threads for set comparisons: ", Threads.nthreads())

    good_views = []
    bad_views = []

    the_whole_thing_is_fine = check_the_view_is_good(the_whole_bool_array)

    if the_whole_thing_is_fine
        return [BitSet(collect(1:size(the_whole_bool_array, 2)))]
    else
        push!(bad_views, the_whole_bool_array)
    end

    while length(bad_views) > 0
        new_bad_views = []
        for bv in bad_views
            new_bitsets = parallel_get_subsets(get_sets_from_view(bv))
            for bitset in new_bitsets
                set = collect(bitset)
                new_view = @view the_whole_bool_array[set, set]
                this_view_is_good = check_the_view_is_good(new_view)
                if this_view_is_good
                    push!(good_views, new_view)
                else
                    push!(new_bad_views, new_view)
                end
            end
        end
        bad_views = new_bad_views
    end

    A = Array{BitSet,1}()

    for gv in good_views
        inds = gv.indices[2]
        push!(A, BitSet(inds))
    end

    return Set(A)
end

# Functions for choosing haplotypes to represent sets
# ---------------------------------------------------

# returns the INDICES of a vector in sorted order. Like sortperm()
# but with ties broken by a second vector
function double_sort_perm(main_vector, secondary_vector;
                          main_increasing = false, secondary_increasing = false)

    @assert length(main_vector) == length(secondary_vector)

    final_order = Array{Int, 1}(undef, length(main_vector))

    sp_1 = sortperm(main_vector, rev = !main_increasing)

    last_item = nothing
    chunk_size = 0
    start = 0
    for (i,item) in enumerate(main_vector[sp_1])
        if i == 1
            last_item = item
            continue
        end

        if item == last_item
            if chunk_size == 0
                start = i - 1
            end

            last_item = item
            chunk_size += 1

        elseif chunk_size > 0
            new_chunk_order = sortperm(secondary_vector[sp_1][start:start + chunk_size], rev = !secondary_increasing)
            sp_1[start:start + chunk_size] = sp_1[start:start + chunk_size][new_chunk_order]
            chunk_size = 0
            last_item = item

        elseif i < length(main_vector[sp_1])
            last_item = item

        end

        if chunk_size > 0
            new_chunk_order = sortperm(secondary_vector[sp_1][start:start + chunk_size], rev = !secondary_increasing)
            sp_1[start:start + chunk_size] = sp_1[start:start + chunk_size][new_chunk_order]
        end

    end

    final_order = sp_1

    return final_order
end

# score sequences based on the number of other sequences
# they're different from
function score_differences(pairwise_identity)

    different_scores = Array{Int, 1}(undef, size(pairwise_identity, 2))

    for i in 1:size(pairwise_identity, 2)
        different_scores[i] = size(pairwise_identity, 1) - sum(pairwise_identity[:,i])
    end

    return different_scores
end

function get_output_sets(pairwise_identity, final_sets, retained, whole_nuc_bit_array, whole_ID_array)

    completeness_scores = score_alignment2(whole_nuc_bit_array)
    difference_scores = score_differences(pairwise_identity)

    # convert the set to an Array
    final_sets_as_array = collect(final_sets)

    # split up singletons and others to relieve the burden on sorting?
    singletons = Array{BitSet,1}()
    non_single = Array{BitSet,1}()

    for set in final_sets_as_array
        if length(set) == 1
            push!(singletons, set)
        else
            push!(non_single, set)
        end
    end

    # sort non-singletons by length in increasing order
    sort!(non_single, by = length)

    # and here are singletons + sorted longer sets, concatenated together
    sorted_final_sets = vcat(singletons, non_single)

    # initiate a ditionary that we will use as a check for whether a particular sequence
    # has been used (to represent a set) yet
    output_set_dict = Dict{String, Array{String, 1}}()

    for (i, set) in enumerate(sorted_final_sets)

        # cs is the columns in the alignment that this set represents
        cs = collect(set)

        # If it's a singleton it can't have been so we can just get on with things.
        if length(cs) == 1

            output_set_dict[whole_ID_array[cs[1]]] = []

        # otherwise get the most-different-from-all-other-sequences sequence
        # (breaking ties with genome completeness)
        else

            set_completeness_scores = completeness_scores[cs]
            set_difference_scores = difference_scores[cs]
            set_nuc_array = whole_nuc_bit_array[:,cs]
            set_IDs = whole_ID_array[cs]

            max_indices = double_sort_perm(set_difference_scores, set_completeness_scores,
                                           main_increasing = false, secondary_increasing = false)

            indx = 1
            best_column = max_indices[indx]
            best_ID = set_IDs[best_column]

            # A while loop to check that the highest scoring sequence hasn't
            # already been used to represent a set - if it has, go through
            # each next best sequence in turn.
            while in(best_ID, keys(output_set_dict)) && indx < length(max_indices)
                indx += 1
                best_column = max_indices[indx]
                best_ID = set_IDs[best_column]
            end

            # if, after the iteration above, all of the IDs in this set are already
            # representing sets (are in used_names), then we can just skip this set (because
            # all its members are assigned to a set anyway).
            if in(best_ID, keys(output_set_dict))
                continue
            end

            best_seq = set_nuc_array[:,best_column]
            other_IDs = setdiff(set_IDs, [best_ID])

            output_set_dict[best_ID] = collect(other_IDs)

        end
    end

    # add retained sequences to the alignment if they aren't included already,
    # and add them to used names so that don't get duplicated in the mapping
    for seqname in retained
        if !in(seqname, keys(output_set_dict))

            indx = findfirst(isequal(seqname), whole_ID_array)

            if indx == nothing
                println(stderr, "warning: " * seqname * " can't be retained because it isn't in the input alignment")
                continue
            end

            if seqname != whole_ID_array[indx]
                throw("bad indexing when forced to retain " * seqname)
            end

            output_set_dict[seqname] = []
        end
    end

    return output_set_dict
end

function get_id_to_column_dict(ID_array)
    ID_to_cl_dict = Dict{String, Int}()
    for (i, ID) in enumerate(ID_array)
        ID_to_cl_dict[ID] = i
    end

    return ID_to_cl_dict
end

function write_fasta(output_set_dict, whole_nuc_bit_array, whole_ID_array, alignment_out, append_IDs)

    ID_to_cl_dict = get_id_to_column_dict(whole_ID_array)

    open(alignment_out, "w") do io

        for tip in keys(output_set_dict)

            clm = ID_to_cl_dict[tip]

            if !append_IDs
                id = ">" * string(tip)
            else
                id = ">" * join(vcat(tip, output_set_dict[tip]), "|")
            end

            seq = get_seq_from_1D_byte_array(whole_nuc_bit_array[:,clm])

            # write fasta header and sequence:
            println(io, id)
            println(io, seq)

        end

        close(io)
    end
end

function get_redundant_seq_to_tip_relationships(representative_to_set)
    #= representative_to_set is a dict with representative
       haplotypes to be written to file as keys, and the other
       members of the set as values.

       Want to get a map from each possible member => representatives
    =#

    D = Dict{String, Array{String, 1}}()

    for (key, value) in representative_to_set
        # key is the representative sequence for this set
        # value is an array of other sequences in the set
        for v in value
            # if this guy is already a representative of a set,
            # then skip it
            if in(v, keys(representative_to_set))
                continue
            end

            if in(v, keys(D))
                D[v] = vcat(D[v], [key])
            else
                D[v] = [key]
            end
        end
    end

    return D
end

function write_redundant_seq_to_tip_relationships(D, filepath)
    # write a table of sequences that are not included in the
    # collapsed alignment, with information about their place(s)
    # on the tree i.e., for each excluded seq, what tip it can
    # be placed at
    open(filepath, "w") do io
        println(io, "sequence,tips")
        for (key, value) in D
            println(io, key * "," * join(value, "|"))
        end
    close(io)
    end
end

function write_tip_to_redundant_seq_relationships(D, filepath, retained)
    # write the converse file to write_redundant_seq_to_tip_relationships()
    # i.e., tips of the tree and what other seqs they represent (if any)
    open(filepath, "w") do io
        print(io, "tip,sequences")
        for (key, value) in D
            # subtract retained sequences from the set that is represented

            print(io, "\n")
            print(io, key * ",")

            first = true
            for v in value
                # if this seq is in the vector of retained tips, skip
                if v in retained
                    continue
                end
                # if this seq is a tip itself, skip
                if v in keys(D)
                    continue
                end
                if first
                    print(io, v)
                    first = false
                    continue
                end
                print(io, "|" * v)
            end
        end
        print(io, "\n")
    close(io)
    end
end

function TEST_get_levels(collection_of_bitsets)
    #=
    return the set of all rows
    that make it into the output - if things
    have worked, then the length of this should
    be the same as the number of sequences - all
    sequences are included in the dataset (at least once)
    =#
    levels = []
    for x in collection_of_bitsets
        for y in x
            push!(levels, y)
        end
    end
    return Set(levels)
end

function TEST_compare_within_sets(array_of_bitsets, nuc_bit_array)
    #=
    compare sequences within sets using the original data -
    =#

    function TEST_all_same(my_small_array)
        for a in 1:size(my_small_array,2) - 1
            for b in (a + 1):size(my_small_array,2)
                for r in 1:size(my_small_array,1)
                    @inbounds x = my_small_array[r,a]
                    @inbounds y = my_small_array[r,b]
                    different = (x & y) < 16

                    # store the different True/False result in an array:
                    if different
                        return false
                    end
                end
            end
        end

        return true
    end

    samples = sample(collect(1:length(array_of_bitsets)), length(array_of_bitsets), replace = false)
    fails = []
    for smp in samples
        bitset = array_of_bitsets[smp]
        columns = collect(bitset)
        if length(columns) > 1
            sub_nuc_bit_array = nuc_bit_array[:,columns]
            test = TEST_all_same(sub_nuc_bit_array)

            if !test
                push!(fails, columns)
            end
        end
    end

    if length(fails) > 0
        return false
    end

    return true
end

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
using Libz

# TODO: make the type below concrete (abstract InputHandle and
# concrete handles for each file extension?)

type InputHandle{N}
    handles::NTuple{N,FASTQ.Reader}
    name::String
    maxreads::Float64
    curreads::Int
end

function InputHandle{N}(reads::NTuple{N,String}; maxreads=Inf)
    name = basename(reads[1])
    # this removes the whole extension .fastq.gz instead of just .gz
    name = split(name,".")[1]
    handles = map(tryopen,reads)
    InputHandle{N}((handles...), name, maxreads, 0)
end

function Base.read!{N}(ih::InputHandle{N},records::NTuple{N,FASTQ.Record})
    for i in 1:N
        read!(ih.handles[i],records[i])
    end
    ih.curreads+=1
end

function Base.eof{N}(ih::InputHandle{N})
    return ih.curreads >= ih.maxreads || any(map(eof,ih.handles))
end

function Base.close(ih::InputHandle)
    map(close,ih.handles)
end

function tryopen(f)
    name, ext = splitext(f)
    if ext == ".fastq"
        FASTQ.Reader(open(f,"r"))
    elseif ext == ".gz"
        # TODO: some files are too long to use with Libz
        FASTQ.Reader(ZlibInflateInputStream(open(f,"r")))
        # FASTQ.Reader(open(pipeline(f,`zcat`))[1])
        # zcatfifo(f)
    else
        error("Unrecognized extension: $ext of the file $f")
    end
end

"""

TODO: This is a terrible hack, but it seems much faster then any of
the Libz or pipeline solutions.

"""
function zcatfifo(f)
    fifo, handle = mktemp()
    close(handle)
    zcat=joinpath(Pkg.dir("FASTQDemultiplexer"),"src","zcat.sh")
    run(`$zcat $f $fifo`)
    FASTQ.Reader(open(fifo))
end

"""

Lists FASTQ files and collects them into tuples of different types of
reads.  The reads can then be opened with InputHandle

"""

function listfastq{N}(inputdir, protocol::Protocol{N})
    # TODO: improve the pattern matching, this approach is probably not good enough
    patterns = map(Regex,protocol.readnames)
    filenames = map(f->joinpath(inputdir,f),readdir(inputdir))

    # TODO: here we assume that the filelist is not malformed, in such
    # case the different reads will not be assigned properly.  I
    # should fix this.
    readnames = map(patterns) do pat
        filtered = filter(filenames) do f
            ismatch(pat,f)
        end
        sort!(filtered)
    end

    # check if all the lengths are the same
    lengths = map(length,readnames)
    if length(unique(lengths)) != 1
        error("Uneven number of files: $lengths")
    end

    return collect(zip(readnames...))
end

max_stars_repo_name: pwl/FASTQDemultiplexer

max_issues_repo_name: pwl/FASTQDemultiplexer

max_forks_repo_name: pwl/FASTQDemultiplexer


text:
module Parallel

function map(count, body::Body, args...) where {Body}
    ccall(:jl_enter_threaded_region, Cvoid, ())
    n_threads = Base.Threads.nthreads()
    n_gen = min(n_threads, count)
    tasks = Vector{Task}(undef, n_gen)
    cnt = (count + n_gen - 1) ÷ n_gen
    for i = 0:(n_gen-1)
        let start = i * cnt, endv = min(count, (i+1) * cnt)-1
        t = Task() do
           for j in start:endv
              body(j+1, args...)
           end
           nothing
        end
        t.sticky = true
        ccall(:jl_set_task_tid, Cint, (Any, Cint), t, i)
        @inbounds tasks[i+1] = t
        schedule(t)
        end
    end
    try
        for t in tasks
            wait(t)
        end
    finally
        ccall(:jl_exit_threaded_region, Cvoid, ())
    end
end

macro parallel(args...)
  captured = args[1:end-1]
  ex = args[end]
  if !(isa(ex, Expr) && ex.head === :for)
    throw(ArgumentError("@threads requires a `for` loop expression"))
  end
  if !(ex.args[1] isa Expr && ex.args[1].head === :(=))
        throw(ArgumentError("nested outer loops are not currently supported by @threads"))
   end
   iter = ex.args[1]
   lidx = iter.args[1]         # index
   range = iter.args[2]
   body = ex.args[2]
   @show body, body.head, body.args
   quote
     let range = $(esc(range))
       function bodyf(idx, iter, $(Base.map(esc, captured)...))
         local $(esc(lidx)) = iter[idx]
         $(esc(body))
       end
       lenr = length(range)
       map(lenr, bodyf, range, $(captured...))
     end
   end
end

end


text:
# Accessors
getLabelDict(dfg::MetaGraphsDFG) = copy(dfg.g.metaindex[:label])
getDescription(dfg::MetaGraphsDFG) = dfg.description
setDescription(dfg::MetaGraphsDFG, description::String) = dfg.description = description
getAddHistory(dfg::MetaGraphsDFG) = dfg.addHistory
getSolverParams(dfg::MetaGraphsDFG) = dfg.solverParams

# setSolverParams(dfg::MetaGraphsDFG, solverParams) = dfg.solverParams = solverParams
function setSolverParams(dfg::MetaGraphsDFG, solverParams::P) where P <: AbstractParams
  dfg.solverParams = solverParams
end

"""
    $(SIGNATURES)
True if the variable or factor exists in the graph.
"""
function exists(dfg::MetaGraphsDFG, node::N) where N <: DFGNode
    return haskey(dfg.g.metaindex[:label], node.label)
end
exists(dfg::MetaGraphsDFG, nId::Symbol) = haskey(dfg.g.metaindex[:label], nId)

"""
    $(SIGNATURES)
Add a DFGVariable to a DFG.
"""
function addVariable!(dfg::MetaGraphsDFG, variable::DFGVariable)::Bool
    if haskey(dfg.g.metaindex[:label], variable.label)
        error("Variable '$(variable.label)' already exists in the factor graph")
    end

    #NOTE Internal ID always set to zero as it is not needed?
    # variable._internalId = 0

    #If other properties are needed in the graph itself, maybe :tags
    # props = Dict{:Symbol, Any}()
    # props[:tags] = variable.tags
    # props[:variable] = variable
    props = Dict(:variable=>variable)

    MetaGraphs.add_vertex!(dfg.g, :label, variable.label) || return false
    MetaGraphs.set_props!(dfg.g, nv(dfg.g), props) || return false

    # Track insertion
    push!(dfg.addHistory, variable.label)

    return true
end

"""
    $(SIGNATURES)
Add a DFGFactor to a DFG.
"""
function addFactor!(dfg::MetaGraphsDFG, variables::Vector{DFGVariable}, factor::DFGFactor)::Bool
    if haskey(dfg.g.metaindex[:label], factor.label)
        error("Factor '$(factor.label)' already exists in the factor graph")
    end
    for v in variables
        if !(v.label in keys(dfg.g.metaindex[:label]))
            error("Variable '$(v.label)' not found in graph when creating Factor '$(factor.label)'")
        end
    end

    factor._variableOrderSymbols = map(v->v.label, variables)

    #NOTE something like this or the next props definition
    # props = Dict{:Symbol, Any}()
    # props[:tags] = factor.tags
    # props[:factor] = factor

    props = Dict(:factor=>factor)

    MetaGraphs.add_vertex!(dfg.g, :label, factor.label) || return false
    set_props!(dfg.g, nv(dfg.g), props) || return false

    # Add index
    # push!(dfg.labels, factor.label)
    # Add the edges...
    for variable in variables
        MetaGraphs.add_edge!(dfg.g, dfg.g[variable.label,:label], dfg.g[factor.label,:label]) || return false
    end

    return true
end

"""
    $(SIGNATURES)
Add a DFGFactor to a DFG.
"""
function addFactor!(dfg::MetaGraphsDFG, variableIds::Vector{Symbol}, factor::DFGFactor)::Bool
    variables = map(vId -> getVariable(dfg, vId), variableIds)
    return addFactor!(dfg, variables, factor)
end

"""
    $(SIGNATURES)
Get a DFGVariable from a DFG using its underlying integer ID.
"""
function getVariable(dfg::MetaGraphsDFG, variableId::Int64)::DFGVariable
    return get_prop(dfg.g, variableId, :variable)
end

function getVariable(g::MetaGraph, variableId::Int64)::DFGVariable
    return get_prop(g, variableId, :variable)
end

"""
    $(SIGNATURES)
Get a DFGVariable from a DFG using its label.
"""
function getVariable(dfg::MetaGraphsDFG, label::Union{Symbol, String})::DFGVariable
    if typeof(label) == String
        label = Symbol(label)
    end
    if !haskey(dfg.g.metaindex[:label], label)
        error("Variable label '$(label)' does not exist in the factor graph")
    end
    return get_prop(dfg.g, dfg.g[label,:label], :variable)
end

"""
    $(SIGNATURES)
Get a DFGFactor from a DFG using its underlying integer ID.
"""
function getFactor(dfg::MetaGraphsDFG, factorId::Int64)::DFGFactor
    # if !(factorId in values(dfg.g.metaindex[:label]))
    #     error("Factor ID '$(factorId)' does not exist in the factor graph")
    # end
    return get_prop(dfg.g, factorId, :factor)
end

function getFactor(g::MetaGraph, factorId::Int64)::DFGFactor
    # if !(factorId in values(dfg.g.metaindex[:label]))
    #     error("Factor ID '$(factorId)' does not exist in the factor graph")
    # end
    return get_prop(g, factorId, :factor)
end

"""
    $(SIGNATURES)
Get a DFGFactor from a DFG using its label.
"""
function getFactor(dfg::MetaGraphsDFG, label::Union{Symbol, String})::DFGFactor
    if typeof(label) == String
        label = Symbol(label)
    end
    if !haskey(dfg.g.metaindex[:label], label)
        error("Factor label '$(label)' does not exist in the factor graph")
    end
    return get_prop(dfg.g, dfg.g[label,:label], :factor)
end

"""
    $(SIGNATURES)
Update a complete DFGVariable in the DFG.
"""
function updateVariable!(dfg::MetaGraphsDFG, variable::DFGVariable)::DFGVariable
    if !haskey(dfg.g.metaindex[:label], variable.label)
        error("Variable label '$(variable.label)' does not exist in the factor graph")
    end
    set_prop!(dfg.g, dfg.g[variable.label,:label], :variable, variable)
    return variable
end

"""
    $(SIGNATURES)
Update a complete DFGFactor in the DFG.
"""
function updateFactor!(dfg::MetaGraphsDFG, factor::DFGFactor)::DFGFactor
    if !haskey(dfg.g.metaindex[:label], factor.label)
        error("Factor label '$(factor.label)' does not exist in the factor graph")
    end
    set_prop!(dfg.g, dfg.g[factor.label,:label], :factor, factor)
    return factor
end

"""
    $(SIGNATURES)
Delete a DFGVariable from the DFG using its label.
"""
function deleteVariable!(dfg::MetaGraphsDFG, label::Symbol)::DFGVariable
    if !haskey(dfg.g.metaindex[:label], label)
        error("Variable label '$(label)' does not exist in the factor graph")
    end
    variable = get_prop(dfg.g, dfg.g[label,:label], :variable)
    rem_vertex!(dfg.g, dfg.g[label,:label])

    return variable
end

#Alias
"""
    $(SIGNATURES)
Delete a referenced DFGVariable from the DFG.
"""
deleteVariable!(dfg::MetaGraphsDFG, variable::DFGVariable)::DFGVariable = deleteVariable!(dfg, variable.label)

"""
    $(SIGNATURES)
Delete a DFGFactor from the DFG using its label.
"""
function deleteFactor!(dfg::MetaGraphsDFG, label::Symbol)::DFGFactor
    if !haskey(dfg.g.metaindex[:label], label)
        error("Factor label '$(label)' does not exist in the factor graph")
    end
    factor = get_prop(dfg.g, dfg.g[label,:label], :factor)
    MetaGraphs.rem_vertex!(dfg.g, dfg.g[label,:label])
    return factor
end

# Alias
"""
    $(SIGNATURES)
Delete the referened DFGFactor from the DFG.
"""
deleteFactor!(dfg::MetaGraphsDFG, factor::DFGFactor)::DFGFactor = deleteFactor!(dfg, factor.label)

"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
"""
function getVariables(dfg::MetaGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{DFGVariable}

    # variables = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGVariable, vertices(dfg.g)))
    variableIds = collect(filter_vertices(dfg.g, :variable))
    variables = map(vId->getVariable(dfg, vId), variableIds)
    if regexFilter != nothing
        variables = filter(v -> occursin(regexFilter, String(v.label)), variables)
    end
    if length(tags) > 0
        mask = map(v -> length(intersect(v.tags, tags)) > 0, variables )
        return variables[mask]
    end
    return variables
end

"""
    $(SIGNATURES)
Get a list of IDs of the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.

Example
```julia
getVariableIds(dfg, r"l", tags=[:APRILTAG;])
```

Related

ls
"""
function getVariableIds(dfg::MetaGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{Symbol}
  vars = getVariables(dfg, regexFilter, tags=tags)
  # mask = map(v -> length(intersect(v.tags, tags)) > 0, vars )
  map(v -> v.label, vars)
end

# Alias
"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
"""
ls(dfg::MetaGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{Symbol} = getVariableIds(dfg, regexFilter, tags=tags)

"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function getFactors(dfg::MetaGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{DFGFactor}
    # factors = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGFactor, vertices(dfg.g)))
    factorIds = collect(filter_vertices(dfg.g, :factor))
    factors = map(vId->getFactor(dfg, vId), factorIds)
    if regexFilter != nothing
        factors = filter(f -> occursin(regexFilter, String(f.label)), factors)
    end
    return factors
end

"""
    $(SIGNATURES)
Get a list of the IDs of the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
getFactorIds(dfg::MetaGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol} = map(f -> f.label, getFactors(dfg, regexFilter))

"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
# Alias
lsf(dfg::MetaGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol} = getFactorIds(dfg, regexFilter)

"""
    $(SIGNATURES)
Alias for getNeighbors - returns neighbors around a given node label.
"""
function lsf(dfg::MetaGraphsDFG, label::Symbol)::Vector{Symbol}
  return getNeighbors(dfg, label)
end

"""
    $(SIGNATURES)
Checks if the graph is fully connected, returns true if so.
"""
function isFullyConnected(dfg::MetaGraphsDFG)::Bool
    return length(LightGraphs.connected_components(dfg.g)) == 1
end

#Alias
"""
    $(SIGNATURES)
Checks if the graph is not fully connected, returns true if it is not contiguous.
"""
hasOrphans(dfg::MetaGraphsDFG)::Bool = !isFullyConnected(dfg)

function _isready(dfg::MetaGraphsDFG, idx::Int, ready::Int)::Bool
    p = props(dfg.g, idx)
    haskey(p, :variable) && (return p[:variable].solvable >= ready)
    haskey(p, :factor) && (return p[:factor].solvable >= ready)

    #TODO should this be an error?
    @warn "Node not a factor or variable"
    return false
end

function _isbackendset(dfg::MetaGraphsDFG, idx::Int, backendset::Int)::Bool
    p = props(dfg.g, idx)
    haskey(p, :variable) && (return p[:variable].solveInProgress == backendset)
    haskey(p, :factor) && (return p[:factor].solveInProgress == backendset)

    #TODO should this be an error?
    @warn "Node not a factor or variable"
    return false
end
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor.
"""
function getNeighbors(dfg::MetaGraphsDFG, node::T; ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol}  where T <: DFGNode
    if !haskey(dfg.g.metaindex[:label], node.label)
        error("Variable/factor with label '$(node.label)' does not exist in the factor graph")
    end

    neighbors = map(idx->get_prop(dfg.g, idx, :label),  LightGraphs.neighbors(dfg.g, dfg.g[node.label,:label]))
    # Additional filtering
    neighbors = ready != nothing ? filter(lbl -> _isready(dfg, dfg.g[lbl,:label], ready), neighbors) : neighbors
    neighbors = backendset != nothing ? filter(lbl -> _isbackendset(dfg, dfg.g[lbl,:label], backendset), neighbors) : neighbors

    # Variable sorting (order is important)
    if node isa DFGFactor
        order = intersect(node._variableOrderSymbols, neighbors)#map(v->v.dfgNode.label, neighbors))
        return order
    end

    return neighbors
end
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor specified by its label.
"""
function getNeighbors(dfg::MetaGraphsDFG, label::Symbol; ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol}  where T <: DFGNode
    if !haskey(dfg.g.metaindex[:label], label)
        error("Variable/factor with label '$(label)' does not exist in the factor graph")
    end

    # neighbors = in_neighbors(vert, dfg.g) #Don't use out_neighbors! It enforces directiveness even if we don't want it
    neighbors = map(idx->get_prop(dfg.g, idx, :label),  LightGraphs.neighbors(dfg.g, dfg.g[label,:label]))
    # Additional filtering
    neighbors = ready != nothing ? filter(lbl -> _isready(dfg, dfg.g[lbl,:label], ready), neighbors) : neighbors
    neighbors = backendset != nothing ? filter(lbl -> _isbackendset(dfg, dfg.g[lbl,:label], backendset), neighbors) : neighbors

    # Variable sorting when using a factor (function order is important)
    if has_prop(dfg.g, dfg.g[label,:label], :factor)
        node = get_prop(dfg.g, dfg.g[label,:label], :factor)
        order = intersect(node._variableOrderSymbols, neighbors)
        return order
    end

    return neighbors

end

# Aliases
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor.
"""
function ls(dfg::MetaGraphsDFG, node::T)::Vector{Symbol} where T <: DFGNode
    return getNeighbors(dfg, node)
end
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor specified by its label.
"""
function ls(dfg::MetaGraphsDFG, label::Symbol)::Vector{Symbol} where T <: DFGNode
    return getNeighbors(dfg, label)
end

function _copyIntoGraph!(sourceDFG::MetaGraphsDFG, destDFG::MetaGraphsDFG, ns::Vector{Int}, includeOrphanFactors::Bool=false)::Nothing
    # Split into variables and factors
    subgraph = sourceDFG.g[ns]
    sourceVariableIds = collect(filter_vertices(subgraph, :variable))
    sourceFactorIds = collect(filter_vertices(subgraph, :factor))
    # or copy out of sourceDFG.
    # sourceFactorIds = intersect(ns, collect(filter_vertices(sourceDFG.g, :factor)))

    #get factor and variable nodes
    variables = map(vId->getVariable(subgraph, vId), sourceVariableIds)
    factors = map(fId->getFactor(subgraph, fId), sourceFactorIds)

    # Now we have to add all variables first,
    for v in variables
        if !haskey(destDFG.g.metaindex[:label], v.label)
            addVariable!(destDFG, deepcopy(v))
        end
    end

    # And then all factors to the destDFG.
    for f in factors
        if !haskey(destDFG.g.metaindex[:label], f.label)
            # Get the original factor variables (we need them to create it)
            # variables = in_neighbors(factor, sourceDFG.g)
            # variables = getNeighbors(sourceDFG, f)
            varIds = LightGraphs.neighbors(sourceDFG.g, sourceDFG.g[f.label,:label])
            variables = map(vId->getVariable(sourceDFG.g, vId), varIds)

            # Find the labels and associated variables in our new subgraph
            factVariables = DFGVariable[]
            for v in variables
                if haskey(destDFG.g.metaindex[:label], v.label)
                    push!(factVariables, getVariable(destDFG, v.label))
                    #otherwise ignore
                end
            end

            # Only if we have all of them should we add it (otherwise strange things may happen on evaluation)
            if includeOrphanFactors || length(factVariables) == length(variables)
                addFactor!(destDFG, factVariables, deepcopy(f))
            end
        end
    end
    return nothing
end

"""
    $(SIGNATURES)
Retrieve a deep subgraph copy around a given variable or factor.
Optionally provide a distance to specify the number of edges should be followed.
Optionally provide an existing subgraph addToDFG, the extracted nodes will be copied into this graph. By default a new subgraph will be created.
Note: By default orphaned factors (where the subgraph does not contain all the related variables) are not returned. Set includeOrphanFactors to return the orphans irrespective of whether the subgraph contains all the variables.
"""
function getSubgraphAroundNode(dfg::MetaGraphsDFG{P}, node::T, distance::Int64=1, includeOrphanFactors::Bool=false, addToDFG::MetaGraphsDFG=MetaGraphsDFG{P}())::MetaGraphsDFG where {P <: AbstractParams, T <: DFGNode}
    if !haskey(dfg.g.metaindex[:label], node.label)
        error("Variable/factor with label '$(node.label)' does not exist in the factor graph")
    end

    # Get a list of all unique neighbors inside 'distance'
    ns = neighborhood(dfg.g, dfg.g[node.label,:label], distance)

    # Copy the section of graph we want
    # _copyIntoGraph!(dfg, addToDFG, map(n->get_prop(dfg.g, n, :label), ns), includeOrphanFactors)
    _copyIntoGraph!(dfg, addToDFG, ns, includeOrphanFactors)
    return addToDFG
end

function getSubgraphAroundNode(dfg::MetaGraphsDFG{<:AbstractParams}, node::DFGNode, distance::Int64, includeOrphanFactors::Bool, addToDFG::AbstractDFG)::AbstractDFG
    if !haskey(dfg.g.metaindex[:label], node.label)
        error("Variable/factor with label '$(node.label)' does not exist in the factor graph")
    end

    # Get a list of all unique neighbors inside 'distance'
    ns = neighborhood(dfg.g, dfg.g[node.label,:label], distance)

    # Copy the section of graph we want
    _copyIntoGraph!(dfg, addToDFG, map(n->get_prop(dfg.g, n, :label), ns), includeOrphanFactors)
    # _copyIntoGraph!(dfg, addToDFG, ns, includeOrphanFactors)
    return addToDFG
end

"""
    $(SIGNATURES)
Get a deep subgraph copy from the DFG given a list of variables and factors.
Optionally provide an existing subgraph addToDFG, the extracted nodes will be copied into this graph. By default a new subgraph will be created.
Note: By default orphaned factors (where the subgraph does not contain all the related variables) are not returned. Set includeOrphanFactors to return the orphans irrespective of whether the subgraph contains all the variables.
"""
function getSubgraph(dfg::MetaGraphsDFG, variableFactorLabels::Vector{Symbol}, includeOrphanFactors::Bool=false, addToDFG::MetaGraphsDFG=MetaGraphsDFG{AbstractParams}())::MetaGraphsDFG
    for label in variableFactorLabels
        if !haskey(dfg.g.metaindex[:label], label)
            error("Variable/factor with label '$(label)' does not exist in the factor graph")
        end
    end

    _copyIntoGraph!(dfg, addToDFG, variableFactorLabels, includeOrphanFactors)
    return addToDFG
end

"""
    $(SIGNATURES)
Get an adjacency matrix for the DFG, returned as a Matrix{Union{Nothing, Symbol}}.
Rows are all factors, columns are all variables, and each cell contains either nothing or the symbol of the relating factor.
The first row and first column are factor and variable headings respectively.
"""
function getIncidenceMatrix(dfg::MetaGraphsDFG)::Matrix{Union{Nothing, Symbol}}
    varLabels = map(v->v.label, getVariables(dfg))
    factLabels = map(f->f.label, getFactors(dfg))
    vDict = Dict(varLabels .=> [1:length(varLabels)...].+1)

    adjMat = Matrix{Union{Nothing, Symbol}}(nothing, length(factLabels)+1, length(varLabels)+1)
    # Set row/col headings
    adjMat[2:end, 1] = factLabels
    adjMat[1, 2:end] = varLabels
    for (fIndex, factLabel) in enumerate(factLabels)
        factVars = getNeighbors(dfg, getFactor(dfg, factLabel))
        map(vLabel -> adjMat[fIndex+1,vDict[vLabel]] = factLabel, factVars)
    end
    return adjMat
end

function getIncidenceMatrixSparse(dfg::MetaGraphsDFG)::Tuple{LightGraphs.SparseMatrixCSC, Vector{Symbol}, Vector{Symbol}}
    adj = LightGraphs.adjacency_matrix(dfg.g)
    v_labels = getVariableIds(dfg)
    f_labels = getFactorIds(dfg)
    v_index = [dfg.g[s,:label] for s in v_labels]
    f_index = [dfg.g[s,:label] for s in f_labels]
    adjvf = adj[f_index, v_index]
    return adjvf, v_labels, f_labels
end

#FIXME remove _ when issue #86 is resolved, for now force dispach to generic toDot
"""
    $(SIGNATURES)
Produces a dot-format of the graph for visualization.
"""
function _toDot(dfg::MetaGraphsDFG)::String
    @error "toDot(dfg::MetaGraphsDFG) is not sopported yet, see https://github.com/JuliaGraphs/MetaGraphs.jl/issues/86"
    m = PipeBuffer()
    MetaGraphs.savedot(m, dfg.g)
    data = take!(m)
    close(m)
    return String(data)
end

"""
    $(SIGNATURES)
Produces a dot file of the graph for visualization.
Download XDot to see the data

Note
- Default location "/tmp/dfg.dot" -- MIGHT BE REMOVED
- Can be viewed with the `xdot` system application.
- Based on graphviz.org
"""
function _toDotFile(dfg::MetaGraphsDFG, fileName::String="/tmp/dfg.dot")::Nothing
    @error "toDotFile(dfg::MetaGraphsDFG,filename) is not supported yet, see https://github.com/JuliaGraphs/MetaGraphs.jl/issues/86"
    open(fileName, "w") do fid
        MetaGraphs.savedot(fid, dfg.g)
    end
    return nothing
end

"""
    $(SIGNATURES)
Gets an empty and unique CloudGraphsDFG derived from an existing DFG.
"""
function _getDuplicatedEmptyDFG(dfg::MetaGraphsDFG{P})::MetaGraphsDFG where {P <: AbstractParams}
    newDfg = MetaGraphsDFG{P}(;
        userId=dfg.userId, robotId=dfg.robotId, sessionId=dfg.sessionId,
        params=deepcopy(dfg.solverParams))
    newDfg.description ="(Copy of) $(dfg.description)"
    return newDfg
end


text:
export SExponents

"""
    SExponents{N}

Store the exponents of the terms of a polynomial (the support) with `N` terms as an type. This results
in an **unique** type for each possible support.
"""
struct SExponents{N}
    exponents::NTuple{N, UInt8}
    size::Tuple{Int,Int} # nvars, nterms
end

function SExponents(exponents::Matrix{<:Integer})
    E = ntuple(i -> convert(UInt8, exponents[i]), length(exponents))

    return SExponents(E, size(exponents))
end

function SExponents{N}(;exponents::NTuple{N, UInt8}=nothing, size::Tuple{Int, Int}=nothing) where {N}
    SExponents{N}(exponents, size)
end

Base.isbits(::Type{<:SExponents}) = true
Base.length(::SExponents{N}) where N = N
function Base.:(==)(f::SExponents{N}, g::SExponents{N}) where {N}
    f.exponents == g.exponents && f.size == g.size
end
Base.hash(f::SExponents, h) = hash(f.exponents, hash(f.size, h))
Base.size(SE::SExponents) = SE.size
Base.size(SE::SExponents, i) = size(SE)[1]

"""
    exponents(::SExponents)

Converts exponents stored in a `SExponents` to a matrix.
"""
function exponents(SE::SExponents)
    nvars, nterms = SE.size
    exps = fill(0, nvars, nterms)
    for k=1:nvars*nterms
        exps[k] = SE.exponents[k]
    end
    exps
end
exponents(::Type{Nothing}) = nothing

exponents(::Type{Nothing}, SE::SExponents) = exponents(SE)
exponents(A::SExponents, B::SExponents) = [exponents(A); exponents(B)]

@static if VERSION ≥ v"0.7-"
    function Base.show(io::IO, SE::SExponents{N}) where {N}
        n = hash(SE.exponents)
        exps_hash = string(n, base=16, pad=sizeof(n) * 2)
        print(io, "SExponents{$N}($(exps_hash))")
    end
else
    function Base.show(io::IO, SE::SExponents{N}) where {N}
        print(io, "SExponents{$N}($(num2hex(hash(SE.exponents))))")
    end
end


text:
debugser(x, level) = nothing;#println(" "^level, typeof(x))
debugser(x::Observable, level) = nothing# println(" "^level, x.id)

function debugser(node::Node, level)
    # return println(" "^level, "DOM")
end

function debugser(x::Vector{T}, level) where {T<:Number}
    # return println(" "^level, "TypedVector: ", T)
end

function debugser(x::Union{AbstractArray, Tuple}, level)
    return foreach(x-> debugser(x, level+1), x)
end

function debugser(dict::AbstractDict, level)
    if haskey(dict, "__javascript_type__")
        print(" "^level, dict["__javascript_type__"])
        if dict["__javascript_type__"] == "Reference"
            println(dict["payload"])
        else
            println()
        end
    else
        for (k, v) in dict
            # println(" "^level, k, ":")
            debugser(v, level + 1)
        end
    end
end

debugser(jss::JSString, level) = nothing#println(" "^level, "JSCode")

function debugser(jsc::JSCode, level)
    # println(" "^level, "JSCode")
end

function debugser(asset::Asset, level)
    # println(" "^level, "Asset")
end

function printsize(name, x)
    sz = Base.summarysize(x)
    if sz/10^6 > 0.01
        @info(name, " > ", Base.format_bytes(sz))
    end
end

debug_size(x, level) = printsize(typeof(x), x);#println(" "^level, typeof(x))
debug_size(x::Observable, level) = nothing# println(" "^level, x.id)

function debug_size(node::Node, level)
    # return println(" "^level, "DOM")
end

function debug_size(x::Vector{T}, level) where {T<:Number}
    printsize(typeof(x), x)
end

function debug_size(x::Union{AbstractArray, Tuple}, level)
    return foreach(x-> debug_size(x, level+1), x)
end

function debug_size(dict::AbstractDict, level)
    if haskey(dict, "__javascript_type__")
        printsize(dict["__javascript_type__"], dict["payload"])
    else
        for (k, v) in dict
            # println(" "^level, k, ":")
            debug_size(v, level + 1)
        end
    end
end

debug_size(jss::JSString, level) = nothing#println(" "^level, "JSCode")

function debug_size(jsc::JSCode, level)
    # println(" "^level, "JSCode")
end

function debug_size(asset::Asset, level)
    # println(" "^level, "Asset")
end


text:
module MultiTest

# Override @testset to use task-heritable storage
import Test

@eval Test using TaskHeritableStorage
@eval Test begin
    #-----------------------------------------------------------------------
    # Various helper methods for test sets

    # Update Test helper methods to store the stack of active testsets in Task Heritable
    # Storage instead of in task_local_storage, which makes Test support concurrent testsets.

    # This is not enough to make them *thread-safe*, since the `results` array in
    # DefaultTestSet is written to in parallel. The next block of code, below, addresses
    # that by adding a lock around setters. A nicer fix would be to use a lock-free accumulator.

    # Key for Task Local Storage
    struct TestSetStackKey end

    # Fork the stack of testsets when new tasks are spawned, so that child tasks don't
    # clobber eachother's stacks.
    TaskHeritableStorage.fork_task_local_storage(k::TestSetStackKey, v) = (k, copy(v))

    # Get the stack of currently executing testsets, inherited from parent tasks
    function get_current_testsets()
        get!(task_local_storage(), TestSetStackKey(), AbstractTestSet[])
    end

    function get_testset()
        testsets = get_current_testsets()
        return isempty(testsets) ? fallback_testset : testsets[end]
    end

    function push_testset(ts::AbstractTestSet)
        testsets = get_current_testsets()
        push!(testsets, ts)
    end

    function pop_testset()
        testsets = get_current_testsets()
        ret = isempty(testsets) ? fallback_testset : pop!(testsets)
        return ret
    end

    function get_testset_depth()
        testsets = get_current_testsets()
        return length(testsets)
    end

    # ------------------------------------------
    # Make Test _thread-safe_ by locking when writing to the results of a parent testset:

    if VERSION >= v"1.3-"

        test_lock_ = ReentrantLock()

        record(ts::DefaultTestSet, t::Broken) = begin
            lock(test_lock_); push!(ts.results, t); unlock(test_lock_);
            t
        end
        # For a passed result, do not store the result since it uses a lot of memory
        record(ts::DefaultTestSet, t::Pass) = begin
            lock(test_lock_); ts.n_passed += 1; unlock(test_lock_);
            t
        end

        # For the other result types, immediately print the error message
        # but do not terminate. Print a backtrace.
        function record(ts::DefaultTestSet, t::Union{Fail, Error})
            if myid() == 1
                printstyled(ts.description, ": ", color=:white)
                # don't print for interrupted tests
                if !(t isa Error) || t.test_type != :test_interrupted
                    print(t)
                    # don't print the backtrace for Errors because it gets printed in the show
                    # method
                    if !isa(t, Error)
                        Base.show_backtrace(stdout, scrub_backtrace(backtrace()))
                    end
                    println()
                end
            end
            lock(test_lock_); push!(ts.results, t); unlock(test_lock_);
            t, isa(t, Error) || backtrace()
        end

        # When a DefaultTestSet finishes, it records itself to its parent
        # testset, if there is one. This allows for recursive printing of
        # the results at the end of the tests
        record(ts::DefaultTestSet, t::AbstractTestSet) = begin
            lock(test_lock_); push!(ts.results, t); unlock(test_lock_);
        end

        function record(ts::DefaultTestSet, t::LogTestFailure)
            if myid() == 1
                printstyled(ts.description, ": ", color=:white)
                print(t)
                Base.show_backtrace(stdout, scrub_backtrace(backtrace()))
                println()
            end
            # Hack: convert to `Fail` so that test summarization works correctly
            lock(test_lock_);
            push!(ts.results, Fail(:test, t.orig_expr, t.logs, nothing, t.source))
            unlock(test_lock_);
            t
        end

    end  # VERSION >= v1.3
end

end


text:
function KnuthBendix.Alphabet(S::AbstractVector{<:GSymbol})
    S = unique!([S; inv.(S)])
    inversions = [findfirst(==(inv(s)), S) for s in S]
    return Alphabet(S, inversions)
end

struct AutomorphismGroup{G<:Group,T,R,S} <: AbstractFPGroup
    group::G
    gens::Vector{T}
    rws::R
    domain::S
end

object(G::AutomorphismGroup) = G.group
rewriting(G::AutomorphismGroup) = G.rws

function equality_data(f::AbstractFPGroupElement{<:AutomorphismGroup})
    imf = evaluate(f)
    # return normalform!.(imf)

    tmp = one(first(imf))
    for g in imf
        normalform!(tmp, g)
        copyto!(g, tmp)
    end
    return imf
end

function Base.:(==)(g::A, h::A) where {A<:AbstractFPGroupElement{<:AutomorphismGroup}}
    @assert parent(g) === parent(h)

    if _isvalidhash(g) && _isvalidhash(h)
        hash(g) != hash(h) && return false
    end

    length(word(g)) > 8 && normalform!(g)
    length(word(h)) > 8 && normalform!(h)

    word(g) == word(h) && return true

    img_computed, imh_computed = false, false

    if !_isvalidhash(g)
        img = equality_data(g)
        _update_savedhash!(g, img)
        img_computed = true
    end
    if !_isvalidhash(h)
        imh = equality_data(h)
        _update_savedhash!(h, imh)
        imh_computed = true
    end

    @assert _isvalidhash(g)
    @assert _isvalidhash(h)

    hash(g) != hash(h) && return false

    # words are different, but hashes agree
    if !img_computed
        img = equality_data(g)
    end
    if !imh_computed
        imh = equality_data(h)
    end

    equal = img == imh
    equal || @warn "hash collision in == :" g h

    return equal
end

function Base.isone(g::AbstractFPGroupElement{<:AutomorphismGroup})
    if length(word(g)) > 8
        normalform!(g)
    end
    return evaluate(g) == parent(g).domain
end

# eye-candy

Base.show(io::IO, ::Type{<:FPGroupElement{<:AutomorphismGroup{T}}}) where {T} =
    print(io, "Automorphism{$T, …}")

Base.show(io::IO, A::AutomorphismGroup) = print(io, "automorphism group of ", object(A))

function Base.show(io::IO, ::MIME"text/plain", a::AbstractFPGroupElement{<:AutomorphismGroup})
    println(io, " ┌ $(a):")
    d = domain(a)
    im = evaluate(a)
    for (x, imx) in zip(d, im[1:end-1])
        println(io, " │   $x ↦ $imx")
    end
    println(io, " └   $(last(d)) ↦ $(last(im))")
end

## Automorphism Evaluation

domain(f::AbstractFPGroupElement{<:AutomorphismGroup}) = deepcopy(parent(f).domain)
# tuple(gens(object(parent(f)))...)

evaluate(f::AbstractFPGroupElement{<:AutomorphismGroup}) = evaluate!(domain(f), f)

function evaluate!(
    t::NTuple{N,T},
    f::AbstractFPGroupElement{<:AutomorphismGroup{<:Group}},
    tmp = one(first(t)),
) where {N, T<:FPGroupElement}
    A = alphabet(f)
    for idx in word(f)
        t = @inbounds evaluate!(t, A[idx], tmp)::NTuple{N,T}
    end
    return t
end

evaluate!(t::NTuple{N, T}, s::GSymbol, tmp=nothing) where {N, T} = throw("you need to implement `evaluate!(::$(typeof(t)), ::$(typeof(s)), ::Alphabet, tmp=one(first(t)))`")

# forward evaluate by substitution

struct LettersMap{T, A}
    indices_map::Dict{Int, T}
    A::A
end

function LettersMap(a::FPGroupElement{<:AutomorphismGroup})
    dom = domain(a)
    @assert all(isone ∘ length ∘ word, dom)
    A = alphabet(first(dom))
    first_letters = first.(word.(dom))
    img = evaluate!(dom, a)

    # (dom[i] → img[i] is a map from domain to images)
    # we need a map from alphabet indices → (gens, gens⁻¹) → images
    # here we do it for elements of the domain
    # (trusting it's a set of generators that define a)
    @assert length(dom) == length(img)

    indices_map = Dict(A[A[fl]] => word(im) for (fl, im) in zip(first_letters, img))
    # inverses of generators are dealt lazily in getindex

    return LettersMap(indices_map, A)
end

function Base.getindex(lm::LettersMap, i::Integer)
    # here i is an index of an alphabet
    @boundscheck 1 ≤ i ≤ length(KnuthBendix.letters(lm.A))

    if !haskey(lm.indices_map, i)
        img = if haskey(lm.indices_map, inv(lm.A, i))
            inv(lm.A, lm.indices_map[inv(lm.A, i)])
        else
            @warn "LetterMap: neither $i nor its inverse has assigned value"
            one(valtype(lm.indices_map))
        end
        lm.indices_map[i] = img
    end
    return lm.indices_map[i]
end

function (a::FPGroupElement{<:AutomorphismGroup})(g::FPGroupElement)
    @assert object(parent(a)) === parent(g)
    img_w = evaluate(word(g), LettersMap(a))
    return parent(g)(img_w)
end

evaluate(w::AbstractWord, lm::LettersMap) = evaluate!(one(w), w, lm)

function evaluate!(res::AbstractWord, w::AbstractWord, lm::LettersMap)
    for i in w
        append!(res, lm[i])
    end
    return res
end

# compile automorphisms

compiled(a) = eval(generated_evaluate(a))

function generated_evaluate(a::FPGroupElement{<:AutomorphismGroup})
    lm = Groups.LettersMap(a)
    d = Groups.domain(a)
    @assert all(length.(word.(d)) .== 1)
    A = alphabet(first(d))
    first_ltrs = first.(word.(d))

    args = [Expr(:call, :*) for _ in first_ltrs]

    for (idx, letter) in enumerate(first_ltrs)
        for l in lm[letter]
            k = findfirst(==(l), first_ltrs)
            if k !== nothing
                push!(args[idx].args, :(d[$k]))
                continue
            end
            k = findfirst(==(inv(A, l)), first_ltrs)
            if k !== nothing
                push!(args[idx].args, :(inv(d[$k])))
                continue
            end
            throw("Letter $l doesn't seem to be mapped anywhere!")
        end
    end
    locals = Dict{Expr, Symbol}()
    locals_counter = 0
    for (i,v) in enumerate(args)
        @assert length(v.args) >= 2
        if length(v.args) > 2
            for (j, a) in pairs(v.args)
                if a isa Expr &&  a.head == :call "$a"
                    @assert a.args[1] == :inv
                    if !(a in keys(locals))
                        locals[a] = Symbol("var_#$locals_counter")
                        locals_counter += 1
                    end
                    v.args[j] = locals[a]
                end
            end
        else
            args[i] = v.args[2]
        end
    end

    q = quote
        $([:(local $v = $k) for (k,v) in locals]...)
    end

    # return args, locals

    return :(d -> begin
        @boundscheck @assert length(d) == $(length(d))
        $q
        @inbounds tuple($(args...))
    end)
end


text:
struct McCraryBinwidth <: BinwidthSelector end

function binwidth(::McCraryBinwidth, Zs::RunningVariable)
    2 * std(Zs) / sqrt(length(Zs))
end

function StatsDiscretizations.samplehull_discretizer(Zs::RunningVariable, b)
    b = binwidth(b, Zs)
    rmin, rmax = extrema(Zs)
    cutoff = Zs.cutoff

    leftend = floor(rmin/b - eps(rmin/b))*b
    leftgrid_length =  round(Int,(cutoff-leftend)/b) + 1

    rightend = (1+ceil(rmax/b + eps(rmax/b)))*b
    rightgrid_length =  round(Int,(rightend-cutoff)/b)

    if Zs.treated ∈ (:≥, :<)
        L, R = :open, :closed
    elseif Zs.treated ∈ (:≤, :>)
        L, R = :closed, :open
    end

    discr_left =  BoundedIntervalDiscretizer{L,R}(
        range(leftend; stop=cutoff, length=leftgrid_length)
    )
    discr_right = BoundedIntervalDiscretizer{L,R}(
        range(cutoff; stop=rightend, length=rightgrid_length + 1)
    )

    if Zs.treated ∈ (:≥, :>)
        discr = (untreated = discr_left, treated = discr_right)
    elseif Zs.treated ∈ (:<, :≤)
        discr = (untreated = discr_rightt, treated = discr_left)
    end
    discr
end