text:
function getweight(gradient::LossGradient{T}, λ::T) where {T<:AbstractFloat}
    -gradient.∂𝑙 / (gradient.∂²𝑙 + λ)
end

function getloss(∂𝑙::T, ∂²𝑙::T, λ::T) where {T<:AbstractFloat}
    a::T = -0.5
    a * ∂𝑙 * ∂𝑙 / (∂²𝑙 + λ)
end

function getloss(node::LeafNode{T}, λ::T) where {T<:AbstractFloat}
    ∂𝑙 = node.gradient.∂𝑙
    ∂²𝑙 = node.gradient.∂²𝑙
    getloss(∂𝑙, ∂²𝑙, λ)
end

function getloss(node::SplitNode{T}, λ::T) where {T<:AbstractFloat}
    node.loss
end

function getgain(node::LeafNode{T}) where {T<:AbstractFloat}
    zero(T)
end

function getgain(node::SplitNode{T}) where {T<:AbstractFloat}
    node.gain
end

function sumgradientslice!(∂𝑙sum0, ∂²𝑙sum0, nodeids::Vector{<:Integer}, nodecansplit::Vector{Bool}, factor::AbstractFactor,
                           ∂𝑙covariate::AbstractCovariate{T}, ∂²𝑙covariate::AbstractCovariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat}

    nodeslices = slice(nodeids, fromobs, toobs, slicelength)
    factorslices = slice(factor, fromobs, toobs, slicelength)
    ∂𝑙slices = slice(∂𝑙covariate, fromobs, toobs, slicelength)
    ∂²𝑙slices = slice(∂²𝑙covariate, fromobs, toobs, slicelength)
    zipslices = zip(nodeslices, factorslices, ∂𝑙slices, ∂²𝑙slices)
    fold((∂𝑙sum0, ∂²𝑙sum0), zipslices) do gradsum, zipslice
        nodeslice, factorslice, ∂𝑙slice, ∂²𝑙slice = zipslice
        ∂𝑙sum, ∂²𝑙sum = gradsum
        @inbounds for i in 1:length(nodeslice)
            nodeid = nodeslice[i]
            if nodecansplit[nodeid]
                levelindex = factorslice[i] + 1
                ∂𝑙sum[nodeid][levelindex] += ∂𝑙slice[i]
                ∂²𝑙sum[nodeid][levelindex] += ∂²𝑙slice[i]
            end
        end
        (∂𝑙sum, ∂²𝑙sum)
    end
end

function sumgradient(nodeids::Vector{<:Integer}, nodecansplit::Vector{Bool}, factor::AbstractFactor, partitions::Vector{LevelPartition},
                     ∂𝑙covariate::AbstractCovariate{T}, ∂²𝑙covariate::AbstractCovariate{T}, slicelength::Integer, singlethread::Bool) where {T<:AbstractFloat}

    nodecount = length(nodecansplit)
    levelcounts = [length(p.mask) + 1 for p in partitions]
    fromobs = 1
    toobs = length(nodeids)

    nthreads = singlethread ? 1 : Threads.nthreads()
    threadspace = map((x -> Int64(floor(x))), range(fromobs, toobs, length = nthreads + 1))
    ∂𝑙sum = [[(nodecansplit[node] ? [zero(Float64) for i in 1:(levelcounts[node])] : Vector{Float64}()) for node in 1:nodecount] for i in 1:nthreads]
    ∂²𝑙sum = [[(nodecansplit[node] ? [zero(Float64) for i in 1:(levelcounts[node])] : Vector{Float64}()) for node in 1:nodecount] for i in 1:nthreads]

    if nthreads > 1
        Threads.@threads for i in 1:nthreads
            sumgradientslice!(∂𝑙sum[i], ∂²𝑙sum[i], nodeids, nodecansplit, factor,
                            ∂𝑙covariate, ∂²𝑙covariate, (i == 1 ? threadspace[i] : threadspace[i] + 1),
                            threadspace[i + 1], slicelength)
        end
        ∂𝑙sum = reduce(+, ∂𝑙sum)
        ∂²𝑙sum = reduce(+, ∂²𝑙sum)
        [(nodecansplit[node] ? [LossGradient{T}(∂𝑙sum[node][i], ∂²𝑙sum[node][i]) for i in 1:(levelcounts[node])] : Vector{LossGradient{T}}()) for node in 1:nodecount]
    else
        sumgradientslice!(∂𝑙sum[1], ∂²𝑙sum[1], nodeids, nodecansplit, factor,
                          ∂𝑙covariate, ∂²𝑙covariate, fromobs, toobs, slicelength)
        [(nodecansplit[node] ? [LossGradient{T}(∂𝑙sum[1][node][i], ∂²𝑙sum[1][node][i]) for i in 1:(levelcounts[node])] : Vector{LossGradient{T}}()) for node in 1:nodecount]
    end
end

function splitnodeidsslice!(nodeids::Vector{<:Integer}, factors, issplitnode::Vector{Bool}, nodemap::Vector{Int64},
                            leftpartitions::Vector{Vector{Bool}}, factorindex::Vector{Int64},
                            fromobs::Integer, toobs::Integer, slicelength::Integer)
    if length(factors) > 0
        factorslices = zip(Tuple([slice(factor, fromobs, toobs, slicelength) for factor in factors]))
        nodeslices = slice(nodeids, fromobs, toobs, slicelength)
        foreach(zip(nodeslices, factorslices)) do x
            nodeslice, fslices = x
            @inbounds for i in 1:length(nodeslice)
                nodeid = nodeslice[i]
                if issplitnode[nodeid]
                    levelindex = fslices[factorindex[nodeid]][i]
                    nodeslice[i] = (leftpartitions[nodeid][levelindex + 1]) ? nodemap[nodeslice[i]] : nodemap[nodeslice[i]] + 1
                else
                    nodeslice[i] = nodemap[nodeslice[i]]
                end
            end
        end
    end
end

function splitnodeids!(nodeids::Vector{<:Integer}, layer::TreeLayer{T}, slicelength::Integer, singlethread::Bool) where {T<:AbstractFloat}
    nodes = layer.nodes
    nodecount = length(nodes)
    len = length(nodeids)
    fromobs = 1
    toobs = len
    issplitnode = [isa(n, SplitNode) && n.isactive for n in nodes]
    nodemap = Vector{Int64}()
    splitnodecount = 0
    for (i, x) in enumerate(issplitnode)
        push!(nodemap, i + splitnodecount)
        if x
            splitnodecount += 1
        end
    end
    factors = Vector{AbstractFactor}()
    factorindex = zeros(Int64, nodecount)
    for i in 1:nodecount
         if issplitnode[i]
             factor = nodes[i].factor
             index = findfirst((x -> x == factor), factors)
             if index === nothing
                 push!(factors, factor)
             end
             factorindex[i] = findfirst((x -> x == factor), factors)
         end
    end
    leftpartitions = [isa(n, SplitNode) && n.isactive ? [n.leftnode.partitions[n.factor].inclmissing; n.leftnode.partitions[n.factor].mask] : Vector{Bool}() for n in nodes]

    nthreads = singlethread ? 1 : Threads.nthreads()
    if nthreads > 1
        threadspace = map((x -> Int64(floor(x))), range(fromobs, toobs, length = nthreads + 1))
        Threads.@threads for j in 1:nthreads
             splitnodeidsslice!(nodeids, factors, issplitnode, nodemap, leftpartitions, factorindex,
                                j == 1 ? threadspace[j] : threadspace[j] + 1,
                                threadspace[j + 1], slicelength)
        end
    else
        splitnodeidsslice!(nodeids, factors, issplitnode, nodemap, leftpartitions, factorindex,
                           fromobs, toobs, slicelength)
    end
    nodeids
end

function getsplitnode(factor::AbstractFactor, leafnode::LeafNode{T}, gradient::Vector{LossGradient{T}},
                      λ::T, min∂²𝑙::T, ordstumps::Bool, optsplit::Bool) where {T<:AbstractFloat}

    partition = leafnode.partitions[factor]
    isord = isordinal(factor)
    gradstart = findfirst(partition.mask) + 1
    ∂𝑙sum0 = sum((grad -> grad.∂𝑙), gradient[gradstart:end])
    ∂²𝑙sum0 = sum((grad -> grad.∂²𝑙), gradient[gradstart:end])
    k = length(gradient) - gradstart + 1
    f = isord ? collect(1:k) :
        (optsplit ? sortperm([gradient[i].∂𝑙 / gradient[i].∂²𝑙 for i in gradstart:length(gradient)]) : collect(1:k))
    miss∂𝑙 = gradient[1].∂𝑙
    miss∂²𝑙 = gradient[1].∂²𝑙
    currloss = getloss(∂𝑙sum0 + miss∂𝑙, ∂²𝑙sum0 + miss∂²𝑙, λ)
    bestloss = typemax(T)
    levelcount = length(partition.mask)
    leftnode = LeafNode{T}(LossGradient{T}(∂𝑙sum0 + miss∂𝑙, ∂²𝑙sum0 + miss∂²𝑙), true, copy(leafnode.partitions))
    rightnode = LeafNode{T}(LossGradient{T}(zero(T), zero(T)), true, copy(leafnode.partitions))
    split = SplitNode{T}(factor, leftnode, rightnode, bestloss, false, zero(T))
    leftpartition = deepcopy(leftnode.partitions[factor])
    rightpartition = deepcopy(rightnode.partitions[factor])

    left∂𝑙sum = isord ? gradient[gradstart].∂𝑙 : gradient[gradstart + f[1] - 1].∂𝑙
    left∂²𝑙sum = isord ? gradient[gradstart].∂²𝑙 : gradient[gradstart + f[1] - 1].∂²𝑙

    firstlevelwithmiss = getloss(left∂𝑙sum + miss∂𝑙, left∂²𝑙sum + miss∂²𝑙, λ) + getloss(∂𝑙sum0 - left∂𝑙sum, ∂²𝑙sum0 - left∂²𝑙sum, λ)
    firstlevelwitouthmiss = getloss(left∂𝑙sum, left∂²𝑙sum, λ) + getloss(∂𝑙sum0 - left∂𝑙sum + miss∂𝑙, ∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙, λ)

    if firstlevelwithmiss < bestloss && (left∂²𝑙sum + miss∂²𝑙 >= min∂²𝑙) && (∂²𝑙sum0 - left∂²𝑙sum >= min∂²𝑙)
        if firstlevelwitouthmiss < firstlevelwithmiss && (left∂²𝑙sum >= min∂²𝑙) && (∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙 >= min∂²𝑙)
            split.leftnode.gradient.∂𝑙 = left∂𝑙sum
            split.leftnode.gradient.∂²𝑙 = left∂²𝑙sum
            split.rightnode.gradient.∂𝑙 = ∂𝑙sum0 - left∂𝑙sum + miss∂𝑙
            split.rightnode.gradient.∂²𝑙 = ∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙
            for j in (gradstart - 1):levelcount
                leftpartition.mask[j] = j == (gradstart + f[1] - 2)
                rightpartition.mask[j] = j == (gradstart + f[1] - 2) ? false : partition.mask[j]
            end
            leftpartition.inclmissing = false
            rightpartition.inclmissing = partition.inclmissing
            split.loss = firstlevelwitouthmiss
        else
            split.leftnode.gradient.∂𝑙 = left∂𝑙sum + miss∂𝑙
            split.leftnode.gradient.∂²𝑙 = left∂²𝑙sum + miss∂²𝑙
            split.rightnode.gradient.∂𝑙 = ∂𝑙sum0 - left∂𝑙sum
            split.rightnode.gradient.∂²𝑙 = ∂²𝑙sum0 - left∂²𝑙sum
            for j in (gradstart - 1):levelcount
                leftpartition.mask[j] = j == (gradstart + f[1] - 2)
                rightpartition.mask[j] = j == (gradstart + f[1] - 2) ? false : partition.mask[j]
            end
            leftpartition.inclmissing = partition.inclmissing
            rightpartition.inclmissing = false
            split.loss = firstlevelwithmiss
        end
    end

    @inbounds for i in (gradstart + 1):(levelcount + 1)
        fi = isord ? i : (gradstart - 1) + f[(i - gradstart) + 1]
        if !partition.mask[fi - 1]
            continue
        end
        ∂𝑙 = gradient[fi].∂𝑙
        ∂²𝑙 = gradient[fi].∂²𝑙

        singlelevelwithmisstotal = getloss(∂𝑙 + miss∂𝑙, ∂²𝑙 + miss∂²𝑙, λ) + getloss(∂𝑙sum0 - ∂𝑙, ∂²𝑙sum0 - ∂²𝑙, λ)
        singlelevelwitouthmisstotal = getloss(∂𝑙, ∂²𝑙, λ) + getloss(∂𝑙sum0 - ∂𝑙 + miss∂𝑙, ∂²𝑙sum0 - ∂²𝑙 + miss∂²𝑙, λ)

        left∂𝑙sum += ∂𝑙
        left∂²𝑙sum += ∂²𝑙

        leftwithmisstotal = getloss(left∂𝑙sum + miss∂𝑙, left∂²𝑙sum + miss∂²𝑙, λ) + getloss(∂𝑙sum0 - left∂𝑙sum, ∂²𝑙sum0 - left∂²𝑙sum, λ)
        leftwithoutmisstotal = getloss(left∂𝑙sum, left∂²𝑙sum, λ) + getloss(∂𝑙sum0 - left∂𝑙sum + miss∂𝑙, ∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙, λ)

        if !isord
            if singlelevelwithmisstotal < split.loss && (∂²𝑙 + miss∂²𝑙 >= min∂²𝑙) && (∂²𝑙sum0 - ∂²𝑙 >= min∂²𝑙)
                if singlelevelwitouthmisstotal < singlelevelwithmisstotal && (∂²𝑙 >= min∂²𝑙) && (∂²𝑙sum0 - ∂²𝑙 + miss∂²𝑙 >= min∂²𝑙)
                    split.leftnode.gradient.∂𝑙 = ∂𝑙
                    split.leftnode.gradient.∂²𝑙 =  ∂²𝑙
                    split.rightnode.gradient.∂𝑙 = ∂𝑙sum0 - ∂𝑙 + miss∂𝑙
                    split.rightnode.gradient.∂²𝑙 = ∂²𝑙sum0 - ∂²𝑙 + miss∂²𝑙
                    leftpartition.inclmissing = false
                    rightpartition.inclmissing = partition.inclmissing
                    split.loss = singlelevelwitouthmisstotal
                else
                    split.leftnode.gradient.∂𝑙 = ∂𝑙 + miss∂𝑙
                    split.leftnode.gradient.∂²𝑙 =  ∂²𝑙 + miss∂²𝑙
                    split.rightnode.gradient.∂𝑙 = ∂𝑙sum0 - ∂𝑙
                    split.rightnode.gradient.∂²𝑙 = ∂²𝑙sum0 - ∂²𝑙
                    leftpartition.inclmissing = partition.inclmissing
                    rightpartition.inclmissing = false
                    split.loss = singlelevelwithmisstotal
                end
                fi = isord ? (i - 1) : (fi - 1)
                for j in (gradstart - 1):levelcount
                    leftpartition.mask[j] = j == fi
                    rightpartition.mask[j] = j == fi ? false : partition.mask[j]
                end
            end
        end

        if isord
            if leftwithmisstotal < split.loss && (left∂²𝑙sum + miss∂²𝑙 >= min∂²𝑙) && (∂²𝑙sum0 - left∂²𝑙sum >= min∂²𝑙)
                if leftwithoutmisstotal < leftwithmisstotal && (left∂²𝑙sum >= min∂²𝑙) && (∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙 >= min∂²𝑙)
                    split.leftnode.gradient.∂𝑙 = left∂𝑙sum
                    split.leftnode.gradient.∂²𝑙 = left∂²𝑙sum
                    split.rightnode.gradient.∂𝑙 = ∂𝑙sum0 - left∂𝑙sum + miss∂𝑙
                    split.rightnode.gradient.∂²𝑙 = ∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙
                    leftpartition.inclmissing = false
                    rightpartition.inclmissing = partition.inclmissing
                    split.loss = leftwithoutmisstotal
                else
                    split.leftnode.gradient.∂𝑙 = left∂𝑙sum + miss∂𝑙
                    split.leftnode.gradient.∂²𝑙 = left∂²𝑙sum + miss∂²𝑙
                    split.rightnode.gradient.∂𝑙 = ∂𝑙sum0 - left∂𝑙sum
                    split.rightnode.gradient.∂²𝑙 = ∂²𝑙sum0 - left∂²𝑙sum
                    leftpartition.inclmissing = partition.inclmissing
                    rightpartition.inclmissing = false
                    split.loss = leftwithmisstotal
                end
                for j in (gradstart - 1):levelcount
                    fj = isord ? j : (gradstart - 2) + f[(j - gradstart) + 2]
                    if j <= i - 1
                        leftpartition.mask[fj] = partition.mask[fj]
                        rightpartition.mask[fj] = false
                    else
                        leftpartition.mask[fj] = false
                        rightpartition.mask[fj] = partition.mask[fj]
                    end
                end
            end
        end
    end

    if split.loss < typemax(T)
        split.gain = currloss - split.loss
        leftnode.partitions[factor] = leftpartition
        leftnode.cansplit = true
        rightnode.partitions[factor] = rightpartition
        rightnode.cansplit = true
        split
    else
        nothing
    end
end

function getnewsplit(gradient::Vector{Vector{LossGradient{T}}}, nodes::Vector{TreeNode{T}}, factor::AbstractFactor,
                     λ::T, γ::T, min∂²𝑙::T, ordstumps::Bool, optsplit::Bool, singlethread::Bool) where {T<:AbstractFloat}
    newsplit = Vector{Union{SplitNode{T}, Nothing}}(undef, length(gradient))
    if !singlethread && length(gradient) > 2 * Threads.nthreads()
        Threads.@threads for i in 1:length(gradient)
            grad = gradient[i]
            if isa(nodes[i], LeafNode) && nodes[i].cansplit
                partition = nodes[i].partitions[factor]
                if count(partition.mask) > 1
                    newsplit[i] = getsplitnode(factor, nodes[i],  grad, λ, min∂²𝑙, ordstumps, optsplit)
                else
                    newsplit[i] = nothing
                end
            else
                newsplit[i] = nothing
            end
        end
    else
        for i in 1:length(gradient)
            grad = gradient[i]
            if isa(nodes[i], LeafNode) && nodes[i].cansplit
                partition = nodes[i].partitions[factor]
                if count(partition.mask) > 1
                    newsplit[i] = getsplitnode(factor, nodes[i],  grad, λ, min∂²𝑙, ordstumps, optsplit)
                else
                    newsplit[i] = nothing
                end
            else
                newsplit[i] = nothing
            end
        end
    end
    newsplit
end

function findbestsplit(state::TreeGrowState{T}) where {T<:AbstractFloat}

    nodecansplit = [isa(n, LeafNode) && n.cansplit for n in state.nodes]
    mingain = T(0.5) * state.γ
    currloss = [getloss(n, state.λ) for n in state.nodes]
    res = foldl(enumerate(state.factors); init = state.nodes) do currsplit, nfactor
        n, factor = nfactor
        partitions = [isa(node, LeafNode) && node.cansplit ? node.partitions[factor] : LevelPartition(Vector{Bool}(), false) for node in state.nodes]

        gradient = sumgradient(state.nodeids, nodecansplit, factor, partitions, state.∂𝑙covariate, state.∂²𝑙covariate, state.slicelength, state.singlethread)

        newsplit = getnewsplit(gradient, state.nodes, factor, state.λ, state.γ, state.min∂²𝑙, state.ordstumps, state.optsplit, state.singlethread)

        res = Vector{TreeNode{T}}(undef, length(newsplit))
        @inbounds for i in 1:length(newsplit)
            if newsplit[i] !== nothing
               newloss = newsplit[i].loss
               newgain = newsplit[i].gain
               if newgain > mingain && newloss < getloss(currsplit[i], state.λ)
                   res[i] = newsplit[i]
               else
                   res[i] = currsplit[i]
               end
            else
               res[i] = currsplit[i]
            end
        end
        res

        # if state.leafwise
        #     gain = map(getgain, res)
        #     (maxgain, imax) = findmax(gain)
        #     if maxgain > mingain

        #     end

        #     gain = [isnull(newsplit[i]) ? T(NaN32) : currloss[i] - get(newsplit[i]).loss for i in 1:length(newsplit)]
        #     (maxgain, imax) = findmax(gain)
        #     if maxgain > bestgain
        #         bestgain = maxgain
        #         @inbounds for i in 1:length(newsplit)
        #             if i == imax
        #                 res[i] = get(newsplit[i])
        #             else
        #                 res[i] = state.nodes[i]
        #             end
        #         end
        #     else
        #         @inbounds for i in 1:length(newsplit)
        #                 res[i] = currsplit[i]
        #         end
        #     end
        # else
        #     mingain = T(0.5) * state.γ
        #     @inbounds for i in 1:length(newsplit)
        #         if !isnull(newsplit[i])
        #            newloss = get(newsplit[i]).loss
        #            if currloss[i] - newloss > mingain && newloss < getloss(currsplit[i], state.λ)
        #                res[i] = get(newsplit[i])
        #            else
        #                res[i] = currsplit[i]
        #            end
        #         else
        #            res[i] = currsplit[i]
        #         end
        #    end
        #end
    end
    if state.leafwise
        gain = map(getgain, res)
        (maxgain, imax) = findmax(gain)
        if maxgain > mingain
            @inbounds for i in 1:length(res)
                if i == imax && isa(res[i], SplitNode)
                    res[i].isactive = true
                elseif isa(res[i], SplitNode)
                    res[i].isactive = false
                end
            end
        end
    else
        @inbounds for i in 1:length(res)
            if isa(res[i], SplitNode)
                res[i].isactive = true
            end
        end
    end
    res
end

function updatestate(state::TreeGrowState{T}, layer::TreeLayer{T}) where {T<:AbstractFloat}
    splitnodeids!(state.nodeids, layer, state.slicelength, state.singlethread)
    factors = state.factors
    newnodes = Vector{TreeNode{T}}()
    @inbounds for (i, n) in enumerate(layer.nodes)
        if isa(n, SplitNode) && n.isactive
            push!(newnodes, n.leftnode)
            push!(newnodes, n.rightnode)
            # if isa(state.nodes[i], LeafNode)
            #     leftpartitions = map(state.nodes[i].partitions) do x
            #         f, p = x
            #         if f == n.factor
            #             f => n.leftpartition
            #         else
            #             x
            #         end
            #     end
            #     rightpartitions = map(state.nodes[i].partitions) do x
            #         f, p = x
            #         if f == n.factor
            #             f => n.rightpartition
            #         else
            #             x
            #         end
            #     end
            #     push!(newnodes, LeafNode{T}(n.leftgradient,
            #                                 n.leftgradient.∂²𝑙 >= state.min∂²𝑙,
            #                                 leftpartitions))
            #     push!(newnodes, LeafNode{T}(n.rightgradient,
            #                                 n.rightgradient.∂²𝑙 >= state.min∂²𝑙,
            #                                 rightpartitions))
        else
            push!(newnodes, n)
        end
    end
    # activefactors = filter(factors) do f
    #     any(map((n -> count(n.partitions[f].mask) > 1), newnodes))
    # end
    # state.factors = activefactors
    # for n in newnodes
    #     n.partitions = filter(n.partitions) do f, p
    #         f in activefactors
    #     end
    # end
    state.nodes = newnodes
    state
end

function nextlayer(state::TreeGrowState{T}) where {T<:AbstractFloat}
    layernodes = findbestsplit(state)
    layer = TreeLayer{T}(layernodes)
    updatestate(state, layer)
    layer, state
end

function predict(treelayer::TreeLayer{T}, nodeids::Vector{<:Integer}, λ::T) where {T<:AbstractFloat}
    weights = Vector{T}()
    @inbounds for (i, node) in enumerate(treelayer.nodes)
        if isa(node, SplitNode) && node.isactive
            push!(weights, getweight(node.leftnode.gradient, λ))
            push!(weights, getweight(node.rightnode.gradient, λ))
        elseif isa(node, SplitNode)
            push!(weights, getweight(node.leftnode.gradient + node.rightnode.gradient, λ))
        else
            push!(weights, getweight(node.gradient, λ))
        end
    end
    (nodeid -> nodeid > 0 ? weights[nodeid] : T(NaN32)).(nodeids)
end

function getlevelmap(fromfactor::AbstractFactor, tofactor::AbstractFactor)
    fromlevels = getlevels(fromfactor)
    tolevels = getlevels(tofactor)
    levelmap = Dict{Int64, Int64}()
    for (i, level) in enumerate(fromlevels)
        j = findfirst((x -> x == level), tolevels)
        if j !== nothing
            levelmap[i] = j
        end
    end
    levelmap
end

function getnewindices(fromfactor::AbstractFactor, tofactor::AbstractFactor)
    fromlevels = getlevels(fromfactor)
    tolevels = getlevels(tofactor)
    newind = Set{Int64}()
    for (i, level) in enumerate(tolevels)
        j = findfirst((x -> x == level), fromlevels)
        if j === nothing
            push!(newind, i)
        end
    end
    newind
end

function Base.map(node::SplitNode{T}, dataframe::AbstractDataFrame,
                  factormap::Dict{AbstractFactor, Tuple{AbstractFactor, Dict{Int64, Int64}, Set{Int64}, Int64}}) where {T<:AbstractFloat}

    leftdict = Dict{AbstractFactor, LevelPartition}()
    rightdict = Dict{AbstractFactor, LevelPartition}()

    for f in keys(node.leftnode.partitions)
        if f in keys(factormap)
            factor, levelmap, newind, levelcount = factormap[f]
            leftmask = Vector{Bool}(undef, levelcount)
            for (i, j) in levelmap
                leftmask[j] = node.leftnode.partitions[f].mask[i]
            end
            for i in newind
                leftmask[i] = false
            end
            leftdict[factor] = LevelPartition(leftmask, node.leftnode.partitions[f].inclmissing)
        end
    end

    for f in keys(node.rightnode.partitions)
        if f in keys(factormap)
            factor, levelmap, newind, levelcount = factormap[f]
            rightmask = Vector{Bool}(undef, levelcount)
            for (i, j) in levelmap
                rightmask[j] = node.rightnode.partitions[f].mask[i]
            end
            for i in newind
                rightmask[i] = true
            end
            rightdict[factor] = LevelPartition(rightmask, node.rightnode.partitions[f].inclmissing)
        end
    end

    leftnode = LeafNode(node.leftnode.gradient, node.leftnode.cansplit, leftdict)
    rightnode = LeafNode(node.rightnode.gradient, node.rightnode.cansplit, rightdict)

    factor, levelmap, newind, levelcount = factormap[node.factor]
    SplitNode{T}(factor, leftnode, rightnode, node.loss, node.isactive, node.gain)
end

function Base.map(node::LeafNode{T}, dataframe::AbstractDataFrame,
                  factormap::Dict{AbstractFactor, Tuple{AbstractFactor, Dict{Int64, Int64}, Set{Int64}, Int64}}) where {T<:AbstractFloat}
    node
end

function predict(tree::XGTree{T}, dataframe::AbstractDataFrame) where {T<:AbstractFloat}
    len = length(dataframe)
    maxnodecount = tree.leafwise ? tree.maxleaves : 2 ^ tree.maxdepth
    nodeids = maxnodecount <= typemax(UInt8) ? ones(UInt8, len) : (maxnodecount <= typemax(UInt16) ? ones(UInt16, len) : ones(UInt32, len))
    for layer in tree.layers
        splitnodeids!(nodeids, layer, tree.slicelength, tree.singlethread)
    end
    predict(tree.layers[end], nodeids, tree.λ)
end

function growtree(factors::Vector{<:AbstractFactor}, ∂𝑙covariate::AbstractCovariate{T},
                  ∂²𝑙covariate::AbstractCovariate{T}, maxdepth::Integer, λ::T, γ::T, leafwise::Bool, maxleaves::Integer,
                  min∂²𝑙::T, ordstumps::Bool, optsplit::Bool, pruning::Bool, slicelength::Integer, singlethread::Bool) where {T<:AbstractFloat}

    len = length(∂𝑙covariate)
    maxnodecount = leafwise ? maxleaves : 2 ^ maxdepth
    maxsteps = leafwise ? (maxleaves - 1) : maxdepth
    nodeids = maxnodecount <= typemax(UInt8) ? ones(UInt8, len) : (maxnodecount <= typemax(UInt16) ? ones(UInt16, len) : ones(UInt32, len))
    intercept = ConstFactor(len)
    grad0 = sumgradient(nodeids, [true], intercept, [LevelPartition([true], false)], ∂𝑙covariate, ∂²𝑙covariate, slicelength, singlethread)[1][2]

    nodes0 = Vector{TreeNode{T}}()
    push!(nodes0, LeafNode{T}(grad0, true, Dict([f => LevelPartition(ones(Bool, length(getlevels(f))), true) for f in factors])))
    state0 = TreeGrowState{T}(nodeids, nodes0, factors, ∂𝑙covariate, ∂²𝑙covariate, λ, γ, min∂²𝑙, ordstumps, optsplit, pruning, leafwise, slicelength, singlethread)
    layers = collect(Iterators.take(Seq(TreeLayer{T}, state0, nextlayer), maxsteps))
    xgtree = XGTree{T}(layers, λ, γ, min∂²𝑙, maxdepth, leafwise, maxleaves, slicelength, singlethread)
    if pruning
        tree = convert(Tree{TreeNode{T}}, xgtree)
        pruned = prune(tree, λ, γ)
        prunedlayers = map((nodes -> TreeLayer{T}(nodes)) , convert(Vector{Vector{TreeNode{T}}}, convert(List{List{TreeNode{T}}}, rebalance(pruned, maxdepth))))
        xgtree = XGTree{T}(prunedlayers, λ, γ, min∂²𝑙, maxdepth, leafwise, maxleaves, slicelength, singlethread)
        pred = predict(xgtree.layers[end], nodeids, λ)
        xgtree, pred
    else
        pred = predict(xgtree.layers[end], nodeids, λ)
        xgtree, pred
    end
end

function Base.convert(::Type{Tree{TreeNode{T}}}, xgtree::XGTree{T}) where {T<:AbstractFloat}
    layers = xgtree.layers
    maxdepth = xgtree.maxdepth
    gettree  = (depth::Integer, nodeid::Integer) ->
        begin
            node = layers[depth].nodes[nodeid]
            if depth == maxdepth
                ConsTree{TreeNode{T}}(node)
            else
                lefttree = gettree(depth + 1, 2 * nodeid - 1)
                righttree = gettree(depth + 1, 2 * nodeid)
                ConsTree{TreeNode{T}}(node, lefttree, righttree)
            end
        end
    gettree(1, 1)
end

function Base.convert(::Type{List{List{TreeNode{T}}}}, tree::EmptyTree{TreeNode{T}}) where {T<:AbstractFloat}
    EmptyList{EmptyList{TreeNode{T}}}()
end

function Base.convert(::Type{List{List{TreeNode{T}}}}, tree::ConsTree{TreeNode{T}}) where {T<:AbstractFloat}
    node = tree.value
    left = convert(List{List{TreeNode{T}}}, tree.lefttree)
    right = convert(List{List{TreeNode{T}}}, tree.righttree)
    ConsList{List{TreeNode{T}}}(ConsList{TreeNode{T}}(node), map((x -> x[1] + x[2]), zip(left, right), List{TreeNode{T}}))
end

function prune(node::SplitNode{T}, λ::T, γ::T) where {T<:AbstractFloat}
    sumgrad = LossGradient(node.leftgradient.∂𝑙 + node.rightgradient.∂𝑙, node.leftgradient.∂²𝑙 + node.rightgradient.∂²𝑙)
    totloss = getloss(sumgrad.∂𝑙, sumgrad.∂²𝑙, λ)
    if totloss - node.loss > γ * T(0.5)
        node::TreeNode{T}
    else
        LeafNode(sumgrad, false, Dict{AbstractFactor, LevelPartition}())::TreeNode{T}
    end
end

function prune(node::LeafNode{T}, λ::T, γ::T) where {T<:AbstractFloat}
    node
end

function prune(tree::ConsTree{<:TreeNode{T}}, λ::T, γ::T) where {T<:AbstractFloat}
    node = tree.value
    if isempty(tree.lefttree) && isempty(tree.righttree)
        ConsTree{TreeNode{T}}(prune(node, λ, γ))
    else
        prunednode = prune(node, λ, γ)
        if isa(prunednode, LeafNode{T})
            ConsTree{TreeNode{T}}(prunednode)
        else
            left = prune(tree.lefttree, λ, γ)
            right = prune(tree.righttree, λ, γ)
            ConsTree{TreeNode{T}}(prunednode, left, right)
        end
    end
end

function rebalance(prunedtree::ConsTree{TreeNode{T}}, maxdepth::Integer) where {T<:AbstractFloat}
    if maxdepth == 1
        prunedtree
    else
        if isempty(prunedtree.lefttree) && isempty(prunedtree.righttree)
            left = rebalance(ConsTree{TreeNode{T}}(prunedtree.value), maxdepth - 1)
            right = rebalance(ConsTree{TreeNode{T}}(prunedtree.value), maxdepth - 1)
            ConsTree{TreeNode{T}}(prunedtree.value, left, right)
        else
            left = rebalance(prunedtree.lefttree, maxdepth - 1)
            right = rebalance(prunedtree.righttree, maxdepth - 1)
            ConsTree{TreeNode{T}}(prunedtree.value, left, right)
        end
    end
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
####################
# Simplex bijector #
####################
struct SimplexBijector{T} <: Bijector{1} where {T} end
SimplexBijector() = SimplexBijector{true}()

function _clamp(x::T, b::Union{SimplexBijector, Inverse{<:SimplexBijector}}) where {T}
    bounds = (zero(T), one(T))
    clamped_x = clamp(x, bounds...)
    DEBUG && _debug("x = $x, bounds = $bounds, clamped_x = $clamped_x")
    return clamped_x
end

(b::SimplexBijector)(x::AbstractVector) = _simplex_bijector(x, b)
(b::SimplexBijector)(y::AbstractVector, x::AbstractVector) = _simplex_bijector!(y, x, b)
function _simplex_bijector(x::AbstractVector, b::SimplexBijector)
    return _simplex_bijector!(similar(x), x, b)
end
function _simplex_bijector!(y, x::AbstractVector, ::SimplexBijector{proj}) where {proj}
    K = length(x)
    @assert K > 1 "x needs to be of length greater than 1"
    T = eltype(x)
    ϵ = _eps(T)
    sum_tmp = zero(T)
    @inbounds z = x[1] * (one(T) - 2ϵ) + ϵ # z ∈ [ϵ, 1-ϵ]
    @inbounds y[1] = StatsFuns.logit(z) + log(T(K - 1))
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k - 1]
        # z ∈ [ϵ, 1-ϵ]
        # x[k] = 0 && sum_tmp = 1 -> z ≈ 1
        z = (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        y[k] = StatsFuns.logit(z) + log(T(K - k))
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if proj
        y[K] = zero(T)
    else
        y[K] = one(T) - sum_tmp - x[K]
    end

    return y
end

# Vectorised implementation of the above.
function (b::SimplexBijector)(X::AbstractMatrix)
    _simplex_bijector(X, b)
end
function (b::SimplexBijector)(
    Y::AbstractMatrix,
    X::AbstractMatrix,
)
    _simplex_bijector!(Y, X, b)
end
function _simplex_bijector(X::AbstractMatrix, b::SimplexBijector)
    _simplex_bijector!(similar(X), X, b)
end
function _simplex_bijector!(Y, X::AbstractMatrix, ::SimplexBijector{proj}) where {proj}
    K, N = size(X, 1), size(X, 2)
    @assert K > 1 "x needs to be of length greater than 1"
    T = eltype(X)
    ϵ = _eps(T)
    @inbounds @simd for n in 1:size(X, 2)
        sum_tmp = zero(T)
        z = X[1, n] * (one(T) - 2ϵ) + ϵ
        Y[1, n] = StatsFuns.logit(z) + log(T(K - 1))
        for k in 2:(K - 1)
            sum_tmp += X[k - 1, n]
            z = (X[k, n] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
            Y[k, n] = StatsFuns.logit(z) + log(T(K - k))
        end
        sum_tmp += X[K-1, n]
        if proj
            Y[K, n] = zero(T)
        else
            Y[K, n] = one(T) - sum_tmp - X[K, n]
        end
    end

    return Y
end

function (ib::Inverse{<:SimplexBijector{proj}})(y::AbstractVector{T}) where {T, proj}
    _simplex_inv_bijector(y, ib.orig)
end
function (ib::Inverse{<:SimplexBijector})(
    x::AbstractVector{T},
    y::AbstractVector{T},
) where {T}
    _simplex_inv_bijector!(x, y, ib.orig)
end
function _simplex_inv_bijector(y::AbstractVector, b::SimplexBijector)
    return _simplex_inv_bijector!(similar(y), y, b)
end
function _simplex_inv_bijector!(x, y::AbstractVector, b::SimplexBijector{proj}) where {proj}
    K = length(y)
    @assert K > 1 "x needs to be of length greater than 1"
    T = eltype(y)
    ϵ = _eps(T)
    @inbounds z = StatsFuns.logistic(y[1] - log(T(K - 1)))
    @inbounds x[1] = _clamp((z - ϵ) / (one(T) - 2ϵ), b)
    sum_tmp = zero(T)
    @inbounds @simd for k = 2:(K - 1)
        z = StatsFuns.logistic(y[k] - log(T(K - k)))
        sum_tmp += x[k-1]
        x[k] = _clamp(((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ, b)
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if proj
        x[K] = _clamp(one(T) - sum_tmp, b)
    else
        x[K] = _clamp(one(T) - sum_tmp - y[K], b)
    end

    return x
end

# Vectorised implementation of the above.
function (ib::Inverse{<:SimplexBijector})(Y::AbstractMatrix)
    _simplex_inv_bijector(Y, ib.orig)
end
function (ib::Inverse{<:SimplexBijector})(
    X::AbstractMatrix{T},
    Y::AbstractMatrix{T},
) where {T <: Real}
    _simplex_inv_bijector!(X, Y, ib.orig)
end
function _simplex_inv_bijector(Y::AbstractMatrix, b::SimplexBijector)
    _simplex_inv_bijector!(similar(Y), Y, b)
end
function _simplex_inv_bijector!(X, Y::AbstractMatrix, b::SimplexBijector{proj}) where {proj}
    K, N = size(Y, 1), size(Y, 2)
    @assert K > 1 "x needs to be of length greater than 1"
    T = eltype(Y)
    ϵ = _eps(T)
    @inbounds @simd for n in 1:size(X, 2)
        sum_tmp, z = zero(T), StatsFuns.logistic(Y[1, n] - log(T(K - 1)))
        X[1, n] = _clamp((z - ϵ) / (one(T) - 2ϵ), b)
        for k in 2:(K - 1)
            z = StatsFuns.logistic(Y[k, n] - log(T(K - k)))
            sum_tmp += X[k - 1, n]
            X[k, n] = _clamp(((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ, b)
        end
        sum_tmp += X[K - 1, n]
        if proj
            X[K, n] = _clamp(one(T) - sum_tmp, b)
        else
            X[K, n] = _clamp(one(T) - sum_tmp - Y[K, n], b)
        end
    end

    return X
end

function logabsdetjac(b::SimplexBijector, x::AbstractVector{T}) where {T}
    ϵ = _eps(T)
    lp = zero(T)

    K = length(x)

    sum_tmp = zero(eltype(x))
    @inbounds z = x[1]
    lp += log(z + ϵ) + log((one(T) + ϵ) - z)
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k-1]
        z = x[k] / ((one(T) + ϵ) - sum_tmp)
        lp += log(z + ϵ) + log((one(T) + ϵ) - z) + log((one(T) + ϵ) - sum_tmp)
    end

    return -lp
end
function simplex_logabsdetjac_gradient(x::AbstractVector)
    T = eltype(x)
    ϵ = _eps(T)
    K = length(x)
    g = similar(x)
    g .= 0
    sum_tmp = zero(eltype(x))
    @inbounds z = x[1]
    #lp += log(z + ϵ) + log((one(T) + ϵ) - z)
    g[1] = -1/(z + ϵ) + 1/((one(T) + ϵ) - z)
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k-1]
        temp = ((one(T) + ϵ) - sum_tmp)
        z = x[k] / temp
        #lp += log(z + ϵ) + log((one(T) + ϵ) - z) + log(temp)
        dzdx = 1 / temp
        dldz = (1/(z + ϵ) - 1/((one(T) + ϵ) - z))
        dldx = dldz * dzdx
	    g[k] -= (1/(z + ϵ) - 1/((one(T) + ϵ) - z)) * 1 / ((one(T) + ϵ) - sum_tmp)
        for i in 1:k-1
	        dzdxp = x[k] * dzdx^2
	        dldxp = dldz * dzdxp - 1 / temp
	        g[i] -= dldxp
	    end
    end
    return g
end
function logabsdetjac(b::SimplexBijector, x::AbstractMatrix{<:Real})
    mapvcat(eachcol(x)) do c
        logabsdetjac(b, c)
    end
end

#=
function simplex_link_val_adjoint(
    x::AbstractVector{T},
    Δ::AbstractVector,
    proj::Bool = true,
) where {T <: Real}
    K = length(x)
    @assert K > 1 "x needs to be of length greater than 1"
    y = similar(x)
    ϵ = _eps(T)
    sum_tmp = zero(T)
    nΔ = zeros(T, length(Δ))
    @inbounds z = x[1] * (one(T) - 2ϵ) + ϵ # z ∈ [ϵ, 1-ϵ]
    @inbounds y[1] = StatsFuns.logit(z) + log(T(K - 1))
    @inbounds nΔ[1] = Δ[1] * (1/z + 1/(1-z)) * (one(T) - 2ϵ)
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k - 1]
        # z ∈ [ϵ, 1-ϵ]
        # x[k] = 0 && sum_tmp = 1 -> z ≈ 1
        z = (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        y[k] = StatsFuns.logit(z) + log(T(K - k))
        nΔ[k] += Δ[k] * (1/z + 1/(1-z)) * (one(T) - 2ϵ) / ((one(T) + ϵ) - sum_tmp)
        for i in 1:k-1
            nΔ[i] += Δ[k] * (1/z + 1/(1-z)) * (x[k] + ϵ) * (one(T) - 2ϵ) / ((one(T) + ϵ) - sum_tmp)^2
        end
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if proj
        y[K] = zero(T)
    else
        y[K] = one(T) - sum_tmp - x[K]
        @simd for i in 1:K
            nΔ[i] += -Δ[K]
        end
    end

    return y, nΔ
end

function simplex_link_val_jacobian(
    x::AbstractVector{T},
    proj::Bool = true,
) where {T <: Real}
    K = length(x)
    @assert K > 1 "x needs to be of length greater than 1"
    y = similar(x)
    dydxt = similar(x, length(x), length(x))
    @inbounds dydxt .= 0
    ϵ = _eps(T)
    sum_tmp = zero(T)

    @inbounds z = x[1] * (one(T) - 2ϵ) + ϵ # z ∈ [ϵ, 1-ϵ]
    @inbounds y[1] = StatsFuns.logit(z) + log(T(K - 1))
    @inbounds dydxt[1,1] = (1/z + 1/(1-z)) * (one(T) - 2ϵ)
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k - 1]
        # z ∈ [ϵ, 1-ϵ]
        # x[k] = 0 && sum_tmp = 1 -> z ≈ 1
        z = (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        y[k] = StatsFuns.logit(z) + log(T(K - k))
        dydxt[k,k] = (1/z + 1/(1-z)) * (one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        for i in 1:k-1
            dydxt[i,k] = (1/z + 1/(1-z)) * (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)^2
        end
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if proj
        y[K] = zero(T)
    else
        y[K] = one(T) - sum_tmp - x[K]
        @simd for i in 1:K
            dydxt[i,K] = -1
        end
    end

    return y, UpperTriangular(dydxt)'
end
=#
function simplex_link_jacobian(
    x::AbstractVector{T},
    proj::Bool = true,
) where {T <: Real}
    K = length(x)
    @assert K > 1 "x needs to be of length greater than 1"
    dydxt = similar(x, length(x), length(x))
    @inbounds dydxt .= 0
    ϵ = _eps(T)
    sum_tmp = zero(T)

    @inbounds z = x[1] * (one(T) - 2ϵ) + ϵ # z ∈ [ϵ, 1-ϵ]
    @inbounds dydxt[1,1] = (1/z + 1/(1-z)) * (one(T) - 2ϵ)
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k - 1]
        # z ∈ [ϵ, 1-ϵ]
        # x[k] = 0 && sum_tmp = 1 -> z ≈ 1
        z = (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        dydxt[k,k] = (1/z + 1/(1-z)) * (one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        for i in 1:k-1
            dydxt[i,k] = (1/z + 1/(1-z)) * (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)^2
        end
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if !proj
        @simd for i in 1:K
            dydxt[i,K] = -1
        end
    end

    return UpperTriangular(dydxt)'
end

#=
function simplex_invlink_val_adjoint(
    y::AbstractVector{T},
    Δ::AbstractVector,
    proj::Bool = true,
) where {T <: Real}
    K = length(y)
    @assert K > 1 "x needs to be of length greater than 1"
    x = similar(y)
    nΔ = similar(Δ)
    @inbounds nΔ .= 0
    dxdy = similar(y, length(y), length(y))
    @inbounds dxdy .= 0

    ϵ = _eps(T)
    @inbounds z = StatsFuns.logistic(y[1] - log(T(K - 1)))
    unclamped_x = (z - ϵ) / (one(T) - 2ϵ)
    clamped_x = _clamp(unclamped_x, SimplexBijector())
    @inbounds x[1] = clamped_x
    @inbounds if unclamped_x == clamped_x
        dxdy[1,1] = z * (1 - z) / (one(T) - 2ϵ)
        nΔ[1] = Δ[1] * dxdy[1,1]
    end
    sum_tmp = zero(T)
    @inbounds for k = 2:(K - 1)
        z = StatsFuns.logistic(y[k] - log(T(K - k)))
        sum_tmp += clamped_x
        unclamped_x = ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ
        clamped_x = _clamp(unclamped_x, SimplexBijector())
        x[k] = clamped_x
        if unclamped_x == clamped_x
            dxdy[k,k] = z * (1 - z) * ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ)
            nΔ[k] = Δ[k] * dxdy[k,k]
            for i in 1:k-1
                for j in i:k-1
                    temp = -dxdy[j,i] * z / (one(T) - 2ϵ)
                    dxdy[k,i] += temp
                    nΔ[i] += Δ[k] * temp
                end
            end
        end
    end
    @inbounds sum_tmp += clamped_x
    @inbounds if proj
    	unclamped_x = one(T) - sum_tmp
        clamped_x = _clamp(unclamped_x, SimplexBijector())
    else
    	unclamped_x = one(T) - sum_tmp - y[K]
        clamped_x = _clamp(unclamped_x, SimplexBijector())
        if unclamped_x == clamped_x
            nΔ[K] = -Δ[K]
        end
    end
    x[K] = clamped_x
    @inbounds if unclamped_x == clamped_x
        for i in 1:K-1
            @simd for j in i:K-1
                nΔ[i] += -Δ[K] * dxdy[j,i]
            end
        end
    end
    return x, nΔ
end

function simplex_invlink_val_jacobian(
    y::AbstractVector{T},
    proj::Bool = true,
) where {T <: Real}
    K = length(y)
    @assert K > 1 "x needs to be of length greater than 1"
    x = similar(y)
    dxdy = similar(y, length(y), length(y))
    @inbounds dxdy .= 0

    ϵ = _eps(T)
    @inbounds z = StatsFuns.logistic(y[1] - log(T(K - 1)))
    unclamped_x = (z - ϵ) / (one(T) - 2ϵ)
    clamped_x = _clamp(unclamped_x, SimplexBijector())
    @inbounds x[1] = clamped_x
    @inbounds if unclamped_x == clamped_x
        dxdy[1,1] = z * (1 - z) / (one(T) - 2ϵ)
    end
    sum_tmp = zero(T)
    @inbounds for k = 2:(K - 1)
        z = StatsFuns.logistic(y[k] - log(T(K - k)))
        sum_tmp += clamped_x
        unclamped_x = ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ
        clamped_x = _clamp(unclamped_x, SimplexBijector())
        x[k] = clamped_x
        if unclamped_x == clamped_x
            dxdy[k,k] = z * (1 - z) * ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ)
            for i in 1:k-1
                for j in i:k-1
                    dxdy[k,i] += -dxdy[j,i] * z / (one(T) - 2ϵ)
                end
            end
        end
    end
    @inbounds sum_tmp += clamped_x
    @inbounds if proj
    	unclamped_x = one(T) - sum_tmp
        clamped_x = _clamp(unclamped_x, SimplexBijector())
    else
    	unclamped_x = one(T) - sum_tmp - y[K]
        clamped_x = _clamp(unclamped_x, SimplexBijector())
        if unclamped_x == clamped_x
            dxdy[K,K] = -1
        end
    end
    x[K] = clamped_x
    @inbounds if unclamped_x == clamped_x
        for i in 1:K-1
            @simd for j in i:K-1
                dxdy[K,i] += -dxdy[j,i]
            end
        end
    end
    return x, LowerTriangular(dxdy)
end
=#
function simplex_invlink_jacobian(
    y::AbstractVector{T},
    proj::Bool = true,
) where {T <: Real}
    K = length(y)
    @assert K > 1 "x needs to be of length greater than 1"
    dxdy = similar(y, length(y), length(y))
    @inbounds dxdy .= 0

    ϵ = _eps(T)
    @inbounds z = StatsFuns.logistic(y[1] - log(T(K - 1)))
    unclamped_x = (z - ϵ) / (one(T) - 2ϵ)
    clamped_x = _clamp(unclamped_x, SimplexBijector())
    @inbounds if unclamped_x == clamped_x
        dxdy[1,1] = z * (1 - z) / (one(T) - 2ϵ)
    end
    sum_tmp = zero(T)
    @inbounds for k = 2:(K - 1)
        z = StatsFuns.logistic(y[k] - log(T(K - k)))
        sum_tmp += clamped_x
        unclamped_x = ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ
        clamped_x = _clamp(unclamped_x, SimplexBijector())
        if unclamped_x == clamped_x
            dxdy[k,k] = z * (1 - z) * ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ)
            for i in 1:k-1
                for j in i:k-1
                    dxdy[k,i] += -dxdy[j,i] * z / (one(T) - 2ϵ)
                end
            end
        end
    end
    @inbounds sum_tmp += clamped_x
    @inbounds if proj
    	unclamped_x = one(T) - sum_tmp
        clamped_x = _clamp(unclamped_x, SimplexBijector())
    else
    	unclamped_x = one(T) - sum_tmp - y[K]
        clamped_x = _clamp(unclamped_x, SimplexBijector())
        if unclamped_x == clamped_x
            dxdy[K,K] = -1
        end
    end
    @inbounds if unclamped_x == clamped_x
        for i in 1:K-1
            @simd for j in i:K-1
                dxdy[K,i] += -dxdy[j,i]
            end
        end
    end
    return LowerTriangular(dxdy)
end

# jacobian
function jacobian(b::SimplexBijector{proj}, x::AbstractVector{T}) where {proj, T}
    K = length(x)
    dydxt = similar(x, length(x), length(x))
    @inbounds dydxt .= 0
    ϵ = _eps(T)
    sum_tmp = zero(T)

    @inbounds z = x[1] * (one(T) - 2ϵ) + ϵ # z ∈ [ϵ, 1-ϵ]
    @inbounds dydxt[1,1] = (1/z + 1/(1-z)) * (one(T) - 2ϵ)
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k - 1]
        # z ∈ [ϵ, 1-ϵ]
        # x[k] = 0 && sum_tmp = 1 -> z ≈ 1
        z = (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        dydxt[k,k] = (1/z + 1/(1-z)) * (one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        for i in 1:k-1
            dydxt[i,k] = (1/z + 1/(1-z)) * (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)^2
        end
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if !proj
        @simd for i in 1:K
            dydxt[i,K] = -1
        end
    end

    return UpperTriangular(dydxt)'
end

function jacobian(ib::Inverse{<:SimplexBijector{proj}}, y::AbstractVector{T}) where {proj, T}
    b = ib.orig

    K = length(y)
    dxdy = similar(y, length(y), length(y))
    @inbounds dxdy .= 0

    ϵ = _eps(T)
    @inbounds z = StatsFuns.logistic(y[1] - log(T(K - 1)))
    unclamped_x = (z - ϵ) / (one(T) - 2ϵ)
    clamped_x = _clamp(unclamped_x, b)
    @inbounds if unclamped_x == clamped_x
        dxdy[1,1] = z * (1 - z) / (one(T) - 2ϵ)
    end
    sum_tmp = zero(T)
    @inbounds for k = 2:(K - 1)
        z = StatsFuns.logistic(y[k] - log(T(K - k)))
        sum_tmp += clamped_x
        unclamped_x = ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ
        clamped_x = _clamp(unclamped_x, b)
        if unclamped_x == clamped_x
            dxdy[k,k] = z * (1 - z) * ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ)
            for i in 1:k-1
                for j in i:k-1
                    dxdy[k,i] += -dxdy[j,i] * z / (one(T) - 2ϵ)
                end
            end
        end
    end
    @inbounds sum_tmp += clamped_x
    @inbounds if proj
    	unclamped_x = one(T) - sum_tmp
        clamped_x = _clamp(unclamped_x, b)
    else
    	unclamped_x = one(T) - sum_tmp - y[K]
        clamped_x = _clamp(unclamped_x, b)
        if unclamped_x == clamped_x
            dxdy[K,K] = -1
        end
    end
    @inbounds if unclamped_x == clamped_x
        for i in 1:K-1
            @simd for j in i:K-1
                dxdy[K,i] += -dxdy[j,i]
            end
        end
    end
    return LowerTriangular(dxdy)
end

ext: jl


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("MLJTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# ## Initial data processing## In this example, we consider the [UCI "wine" dataset](https://archive.ics.uci.edu/ml/datasets/wine)## > These data are the results of a chemical analysis of wines grown in the same region in Italy but derived from three different cultivars. The analysis determined the quantities of 13 constituents found in each of the three types of wines.## ### Getting the data# Let's download the data thanks to the [HTTP.jl](HTTP.get("https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data")) package and load it into a DataFrame via the [CSV.jl](https://github.com/JuliaData/CSV.jl) package:
using HTTP, CSV, MLJ, StatsBase, PyPlot

req = HTTP.get("https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data")
data = CSV.read(req.body,
                header=["Class", "Alcool", "Malic acid",
                        "Ash", "Alcalinity of ash", "Magnesium",
                        "Total phenols", "Flavanoids",
                        "Nonflavanoid phenols", "Proanthcyanins",
                        "Color intensity", "Hue",
                        "OD280/OD315 of diluted wines", "Proline"])
# the target is the Class column, everything else is a feature
y, X = unpack(data, ==(:Class), colname->true);

# ### Setting the scientific type## Let's explore the scientific type attributed by default to the target and the features
scitype(y)

# this should be changed as it should be considered as an ordered factor
yc = coerce(y, OrderedFactor);

# Let's now consider the features
scitype(X)

# So there are `Continuous` values (encoded as floating point) and `Count` values (integer).# Note also that there are no missing value (otherwise one of the scientific type would have been a `Union{Missing,*}`).# Let's check which column is what:
sch = schema(X)
println(rpad(" Name", 28), "| Scitype")
println("-"^45)
for (name, scitype) in zip(sch.names, sch.scitypes)
    println(rpad("$name", 30), scitype)
end

# The two variable that are encoded as `Count` can  probably be re-interpreted; let's have a look at the `Proline` one to see what it looks like
X[1:5, :Proline]

# It can likely be interpreted as a Continuous as well (it would be better to know precisely what it is but for now let's just go with the hunch).# We'll do the same with `:Magnesium`:
Xc = coerce(X, :Proline=>Continuous, :Magnesium=>Continuous);

# Finally, let's have a quick look at the mean and standard deviation of each feature to get a feel for their amplitude:
describe(Xc, :mean, :std)

# Right so it varies a fair bit which would invite to standardise the data.## **Note**: to complete such a first step, one could explore histograms of the various features for instance, check that there is enough variation among the continuous features and that there does not seem to be problems in the encoding, we cut this out to shorten the tutorial. We could also have checked that the data was balanced.## ## Getting a baseline## It's a multiclass classification problem with continuous inputs so a sensible start is  to test two very simple classifiers to get a baseline.# We'll train two simple pipelines:# - a Standardizer + KNN classifier and# - a Standardizer + Multinomial classifier (logistic regression).
@load KNNClassifier pkg="NearestNeighbors"
@load MultinomialClassifier pkg="MLJLinearModels";

@pipeline KnnPipe(std=Standardizer(), clf=KNNClassifier()) is_probabilistic=true
@pipeline MnPipe(std=Standardizer(), clf=MultinomialClassifier()) is_probabilistic=true

# We can now fit this on a train split of the data setting aside 20% of the data for eventual testing.
train, test = partition(eachindex(yc), 0.8, shuffle=true, rng=111)
Xtrain = selectrows(Xc, train)
Xtest = selectrows(Xc, test)
ytrain = selectrows(yc, train)
ytest = selectrows(yc, test);

# Let's now wrap an instance of these models with data (all hyperparameters are set to default here):
knn = machine(KnnPipe(), Xtrain, ytrain)
multi = machine(MnPipe(), Xtrain, ytrain)

# Let's train a KNNClassifier with default hyperparameters and get a baseline misclassification rate using 90% of the training data to train the model and the remaining 10% to evaluate it:
opts = (resampling=Holdout(fraction_train=0.9), measure=cross_entropy)
res = evaluate!(knn; opts...)
round(res.measurement[1], sigdigits=3)

# Now we do the same with a MultinomialClassifier
res = evaluate!(multi; opts...)
round(res.measurement[1], sigdigits=3)

# Both methods seem to offer comparable levels of performance.# Let's check the misclassification over the full training set:
mcr_k = misclassification_rate(predict_mode(knn, Xtrain), ytrain)
mcr_m = misclassification_rate(predict_mode(multi, Xtrain), ytrain)
println(rpad("KNN mcr:", 10), round(mcr_k, sigdigits=3))
println(rpad("MNC mcr:", 10), round(mcr_m, sigdigits=3))

# So here we have done no hyperparameter training and already have a misclassification rate below 5%.# Clearly the problem is not very difficult.## ## Visualising the classes## One way to get intuition for why the dataset is so easy to classify is to project it onto a 2D space using the PCA and display the two classes to see if they are well separated; we use the arrow-syntax here (if you're on Julia <= 1.2, use the commented-out lines as you won't be able to use the arrow-syntax)
# @pipeline PCAPipe(std=Standardizer(), t=PCA(maxoutdim=2))
# pca = machine(PCAPipe(), Xtrain)
# fit!(pca, Xtrain)
# W = transform(pca, Xtrain)

@load PCA

pca = Xc |> Standardizer() |> PCA(maxoutdim=2)
fit!(pca)
W = pca(rows=train);

# Let's now display this using different colours for the different classes:
x1 = W.x1
x2 = W.x2

mask_1 = ytrain .== 1
mask_2 = ytrain .== 2
mask_3 = ytrain .== 3

figure(figsize=(8, 6))
plot(x1[mask_1], x2[mask_1], linestyle="none", marker="o", color="red")
plot(x1[mask_2], x2[mask_2], linestyle="none", marker="o", color="blue")
plot(x1[mask_3], x2[mask_3], linestyle="none", marker="o", color="magenta")

xlabel("PCA dimension 1", fontsize=14)
ylabel("PCA dimension 2", fontsize=14)
legend(["Class 1", "Class 2", "Class 3"], fontsize=12)
xticks(fontsize=12)
yticks(fontsize=12)

# \figalt{PCA}{EX-wine-pca.svg}## On that figure it now becomes quite clear why we managed to achieve such high scores with very simple classifiers.# At this point it's a bit pointless to dig much deaper into parameter tuning etc.## As a last step, we can report performances of the models on the test set which we set aside earlier:
perf_k = misclassification_rate(predict_mode(knn, Xtest), ytest)
perf_m = misclassification_rate(predict_mode(multi, Xtest), ytest)
println(rpad("KNN mcr:", 10), round(perf_k, sigdigits=3))
println(rpad("MNC mcr:", 10), round(perf_m, sigdigits=3))

# Pretty good for so little work!

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
macro export_all_files(cur_item, package_name)
  all_files = get_all_files(eval(cur_item))

  cur_module = getfield(Main, Symbol(package_name))

  cur_expression = quote end
  for file in all_files
    cur_export = split(file, "/")[end]

    !endswith(cur_export, ".jl") && continue

    cur_export = replace(cur_export, ".jl", "")

    in(cur_export, julia_reserved_words) && continue

    for added_string in ["", "!"]
      cur_export *= added_string
      exported_object = Symbol(cur_export)

      !isdefined(cur_module, exported_object) && continue

      cur_block = quote
        export $(esc(exported_object))
      end

      append!(cur_expression.args, cur_block.args)
    end
  end

  cur_expression
end


text:
# when new brainflow binaries are released, then this script needs to be re-executed.

using SHA
using Pkg
using Pkg.Artifacts
using Tar

function sha256sum(tarball_path)
    return open(tarball_path, "r") do io
        return bytes2hex(sha256(io))
    end
end

function add_brainflow_artifact!(
    url,
    artifact_toml_path = "Artifacts.toml",
    artifact_name = "brainflow",
    )

    download_path = "$(tempname())-download.tar"
    download(url, download_path)
    tar_hash_sha256 = sha256sum(download_path)

    brainflow_hash = create_artifact() do artifact_dir
        # Pkg.PlatformEngines.unpack() gives errors for some users on Windows
        Tar.extract(download_path, artifact_dir)
    end

    rm(download_path)

    Pkg.Artifacts.bind_artifact!(
        artifact_toml_path,
        artifact_name,
        brainflow_hash;
        download_info=[(url, tar_hash_sha256)],
        force=true,
        lazy=true,
    )

    return brainflow_hash
end

cd(@__DIR__)
include("src/brainflow_url.jl")
add_brainflow_artifact!(brainflow_url())

max_stars_repo_name: p-riebs/brainflow

max_issues_repo_name: p-riebs/brainflow

max_forks_repo_name: p-riebs/brainflow


text:
function getInputString()
    open("day17input.txt") do f
        return strip(read(f, String))
    end
end

function getMinMaxXY(lines)
    xs = []
    ys = []
    for line in lines
        fixedCoord, coordRange = split(line)

        fixedValue = parse(Int, fixedCoord[3:end-1])
        if fixedCoord[1] == 'x'
            push!(xs, fixedValue)
        elseif fixedCoord[1] == 'y'
            push!(ys, fixedValue)
        end

        rangeStart, rangeEnd = map(s -> parse(Int, s), split(coordRange[3:end], ".."))
        if coordRange[1] == 'x'
            push!(xs, rangeStart:rangeEnd...)
        elseif coordRange[1] == 'y'
            push!(ys, rangeStart:rangeEnd...)
        end
    end
    return (min(xs...), min(ys...), max(xs...), max(ys...))
end

function parseInput(s)
    lines = filter(!isempty, split(s, '\n'))
    xMin, yMin, xMax, yMax = getMinMaxXY(lines)
    grid = Array{Char}(undef, xMax+1, yMax)
    fill!(grid, '.')
    for line in lines
        fixedCoord, coordRange = split(line)

        fixedValue = parse(Int, fixedCoord[3:end-1])

        rangeStart, rangeEnd = map(s -> parse(Int, s), split(coordRange[3:end], ".."))

        if fixedCoord[1] == 'x'
            grid[fixedValue, rangeStart:rangeEnd] .= '#'
        elseif fixedCoord[1] == 'y'
            grid[rangeStart:rangeEnd, fixedValue] .= '#'
        end
     end
     return grid, yMin, yMax
end

function printGrid(grid)
    println("Grid:")
    minX = Inf
    for y in 1:size(grid, 2)
        for x in 1:size(grid, 1)
            if x < minX && grid[x,y] != '.'
                minX = x
            end
        end
    end

    for y in 1:size(grid, 2)
        for x in minX:size(grid, 1)
            print(grid[x,y])
        end
        print("\n")
    end
end

function isPermeable(c::Char)
    c != '#' && c != '~'
end

function floodRow!(loc, grid)
    grid[loc...] = '~'

    left = loc .+ (-1, 0)
    while grid[left...] == '.' || grid[left...] == '|'
        grid[left...] = '~'
        left = left .+ (-1, 0)
    end

    right = loc .+ (1, 0)
    while grid[right...] == '.' || grid[right...] == '|'
        grid[right...] = '~'
        right = right .+ (1, 0)
    end

    return loc .+ (0, -1)
end

function findPourOvers(loc, grid)
    pourOvers = []
    grid[loc...] = '|'

    left = loc .+ (-1, 0)
    while grid[left...] == '.' || grid[left...] == '|'
        grid[left...] = '|'
        if isPermeable(grid[(left .+ (0, 1))...])
            push!(pourOvers, left)
            break
        end
        left = left .+ (-1, 0)
    end

    right = loc .+ (1, 0)
    while grid[right...] == '.' || grid[right...] == '|'
        grid[right...] = '|'
        if isPermeable(grid[(right .+ (0, 1))...])
            push!(pourOvers, right)
            break
        end
        right = right .+ (1, 0)
    end
    return pourOvers
end

function flowToPourOver(loc, grid)
    if grid[loc...] == '~'
        return []
    end

    # Mark this square as hit with water
    grid[loc...] = '|'

    # Flow down as far as possible
    down = loc
    while isPermeable(grid[down...])
        grid[down...] = '|'
        loc = down
        down = down .+ (0, 1)
        if down[2] > size(grid, 2)
            return []
        end
    end

    # we hit something, so find downspouts where we pour over
    pourOvers = findPourOvers(loc, grid)
    if length(pourOvers) > 0
        return pourOvers
    end

    # We didn't find any down spouts, so we're in a well. Flood it until we find downspouts
    while length(pourOvers) == 0
        loc = floodRow!(loc, grid)
        pourOvers = findPourOvers(loc, grid)
    end

    return pourOvers
end

function fillGrid!(grid)
    # Keep track of the places where water is traveling
    q = deque(Tuple{Int,Int})
    processed = Set()

    push!(q, (500, 1))

    while length(q) > 0
        square = popfirst!(q)
        if square in processed
            continue
        end

        toFlow = flowToPourOver(square, grid)
        if length(toFlow) > 0
            push!(q, toFlow...)
        end
        push!(processed, square)
        n = length(q)
        println("$n pourovers to process")
    end
end

function countWetSquares(grid, minY, maxY)
    length(filter(c -> c == '|' || c == '~', grid[:, minY:maxY]))
end

function countFilledSquares(grid, minY, maxY)
    length(filter(c -> c == '~', grid[:, minY:maxY]))
end

function solvePartA()
    grid, yMin, yMax = parseInput(getInputString())
    fillGrid!(grid)
    countWetSquares(grid, yMin, yMax)
end

function solvePartB()
    grid, yMin, yMax = parseInput(getInputString())
    fillGrid!(grid)
    countFilledSquares(grid, yMin, yMax)
end


text:
importall Bukdu

type WelcomeController <: ApplicationController
end

index(::WelcomeController) = "Hello Bukdu"

Router() do
    get("/", WelcomeController, index)
end

Bukdu.start(8080)

Endpoint() do
    plug(Plug.Logger)
    plug(Router)
end

# wait()

# Bukdu.stop()


text:
export relatedness, merge_relatedness

"""
    bootstrap_summary(::Vector{Union{Missing, Float64}}, width::Tuple{Float64, Float64})

Return the mean, median, standard error, and quantiles (given by `witdth`) of relatedness resampling.
"""
@inline function bootstrap_summary(boot_out::Vector{Union{Missing, Float64}}, width::Tuple{Float64, Float64})
    all(ismissing.(boot_out)) == true && return missing, missing, missing, missing
    boot_skipmissing = collect(skipmissing(boot_out))
    n_nonmiss = length(boot_skipmissing)
    Mean = mean(boot_skipmissing)
    Median = median(boot_skipmissing)
    SE = sqrt(sum((boot_skipmissing - (boot_skipmissing / n_nonmiss)).^2) / (n_nonmiss - 1))
    quants = quantile(boot_skipmissing, width)

    return Mean, Median, SE, quants
end

"""
    bootstrap_genos_all(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}, n_per_loc::Vector{Int}, alleles::NamedTuple; method::Function, iterations::Int)

Perform `iterations` number of bootstrap resampling iterations of all genotypes between pair (`ind1` `ind2`). Returns a vector of length `interatotions`
of the relatedness estimate given by method `method`. This is an internal function with `locus_names`, `n_per_loc`, and `alleles` supplied by `relatedness_boot_all`.
"""
@inline function bootstrap_genos_all(ind1::T, ind2::T, locus_names::Vector{Symbol}, n_per_loc::Vector{Int}, alleles::U; method::F, iterations::Int, inbreeding::Bool, n::Int) where T <: GenoArray where U <: NamedTuple where F
    relate_vec_boot = Vector{Union{Missing,Float64}}(undef, iterations)
    n_loc = length(locus_names)
    @sync for iter in 1:iterations
        Base.Threads.@spawn begin
            # bootstrap the indices
            boot_idx = rand(Xoroshiro128Star(), 1:n_loc, n_loc)
            # sample the source vectors with the resampled/bootstrapped indices
            ind1_boot, ind2_boot, loc_boot, n_per_loci = map(i -> getindex(i, boot_idx), [ind1, ind2, locus_names, n_per_loc])
            # get index for genotype appearing missing in at least one individual in the pair
            keep_idx = nonmissings(ind1_boot, ind2_boot)
            relate_vec_boot[iter] = method(ind1_boot[keep_idx], ind2_boot[keep_idx], loc_boot[keep_idx], alleles, loc_n = n_per_loci[keep_idx], n_samples = n_loc, inbreeding = inbreeding)
        end
    end
    return relate_vec_boot
end

"""
    bootstrap_genos_nonmissing(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}, n_per_loc::Vector{Int}, alleles::NamedTuple; method::Function, iterations::Int)

Perform `iterations` number of bootstrap resampling iterations of only shared (nonmissing) genotypes between pair (`ind1` `ind2`). Returns a vector of length `interatotions`
of the relatedness estimate given by method `method`. This is an internal function with `locus_names`, `n_per_loc`, and `alleles` supplied by `relatedness_boot_nonmissing`.
"""
@inline function bootstrap_genos_nonmissing(ind1::T, ind2::T, locus_names::Vector{Symbol}, n_per_loc::Vector{Int}, alleles::U; method::F, iterations::Int, inbreeding::Bool) where T <: GenoArray where U <: NamedTuple where F
    relate_vec_boot = Vector{Union{Missing,Float64}}(undef, iterations)
    n_loc = length(locus_names)
    @sync for iter in 1:iterations
        Base.Threads.@spawn begin
            # bootstrap the indices
            boot_idx = rand(Xoroshiro128Star(), 1:n_loc, n_loc)
            # sample the source vectors with the resampled/bootstrapped indices
            ind1_boot, ind2_boot, loc_boot, n_per_loci = map(i -> getindex(i, boot_idx), [ind1, ind2, locus_names, n_per_loc])
            # faster/cheaper n counting
            relate_vec_boot[iter] = method(ind1_boot, ind2_boot, loc_boot, alleles, loc_n = n_per_loci, n_samples = n_loc, inbreeding = inbreeding)
        end
    end
    return relate_vec_boot
end

#FEATURE namedtuple output
"""
    relatedness_boot_all(::PopData, sample_names::Vector{String}; method::Function, iterations::Int, interval::Tuple{Float64, Float64})

Calculate pairwise relatedness between all combinations of the provided `sample_names` for each `method` provided. Bootstrapping resamples using
the `all` method, where resampling occurs over all loci. This is an internal function with all arguments provided by `relatedness`.
"""
function relatedness_boot_all(data::PopData, sample_names::Vector{String}; method::F, iterations::Int = 100, interval::Tuple{Float64, Float64} = (0.025, 0.975), inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    uniq_pops = unique(data.meta.population)
    if first(uniq_pops) ∈ ["fullsib", "halfsib", "unrelated", "parent_offspring"]
        sample_pairs = sim_pairs(sample_names)
    else
        sample_pairs = pairwise_pairs(sample_names)
    end
    n_pairs = length(sample_pairs)
    n_samples = length(samples(data))
    allele_frequencies = allele_freq(data)
    n_per_loci = DataFrames.combine(groupby(data.loci, :locus), :genotype => nonmissing => :n)[:, :n]
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, n_pairs), 1:length(method))
    boot_means, boot_medians, boot_ses = map(i -> deepcopy(relate_vecs), 1:3)
    boot_CI = map(i -> Vector{Union{Missing,Tuple{Float64,Float64}}}(undef, n_pairs), 1:length(method))
    shared_loci = Vector{Int}(undef, n_pairs)
    p = Progress(n_pairs*length(method), dt = 1, color = :blue, barglyphs = BarGlyphs("|=> |"), barlen = 30)
    popdata_idx = groupby(data.loci, :name)
    @inbounds for i in 1:n_pairs
        @inbounds geno1 = popdata_idx[(sample_pairs[i][1],)].genotype
        @inbounds geno2 = popdata_idx[(sample_pairs[i][2],)].genotype
        # get index for genotype appearing missing in at least one individual in the pair
        keep_idx = nonmissings(geno1, geno2)
        # generate nonmissing genotype data
        gen1, gen2, loc, n_per_loc = (i[keep_idx] for i in (geno1, geno2, loci_names, n_per_loci))
        @inbounds shared_loci[i] = length(keep_idx)
        @sync @inbounds for (j, mthd) in enumerate(method)
            Base.Threads.@spawn begin
                @inbounds relate_vecs[j][i] = mthd(gen1, gen2, loc, allele_frequencies, loc_n = n_per_loci, n_samples = n_samples, inbreeding = inbreeding)
                boot_out = bootstrap_genos_all(geno1, geno2, loci_names, n_per_loci, allele_frequencies, method = mthd, iterations = iterations, inbreeding = inbreeding, n = j+1)
                @inbounds boot_means[j][i], boot_medians[j][i], boot_ses[j][i], boot_CI[j][i] = bootstrap_summary(boot_out, interval)
                pair_text = sample_pairs[i][1] * " × " * sample_pairs[i][2] * "  ($i" * " of " * "$(n_pairs)" * ")"
                ProgressMeter.next!(p; showvalues = [(:Pair, pair_text), (:Method, mthd)])
            end
        end
    end
    method_colnames = [Symbol("$i") for i in method]
    boot_mean_colnames = [Symbol("$i"*"_mean") for i in method]
    boot_median_colnames = [Symbol("$i"*"_median") for i in method]
    boot_se_colnames = [Symbol("$i"*"_SE") for i in method]
    CI_percent = convert(Int64, round(interval[2] - interval[1], digits = 2) * 100)
    boot_CI_colnames = [Symbol("$i"*"_CI_"*"$CI_percent") for i in method]

    out_dfs = map(method_colnames) do mthod
        DataFrame(:sample_1 => map(i -> i[1], sample_pairs), :sample_2 => map(i -> i[2], sample_pairs), :n_loci => shared_loci)
    end
    @inbounds for (i, mth) in enumerate(method_colnames)
        out_dfs[i][:, mth] = relate_vecs[i]
        out_dfs[i][:, boot_mean_colnames[i]] = boot_means[i]
        out_dfs[i][:, boot_median_colnames[i]] = boot_medians[i]
        out_dfs[i][:, boot_se_colnames[i]] = boot_ses[i]
        out_dfs[i][:, boot_CI_colnames[i]] = boot_CI[i]
    end
    return (; (method_colnames .=> out_dfs)...)
end

"""
    relatedness_boot_nonmissing(::PopData, sample_names::Vector{String}; method::F, iterations::Int, interval::Tuple{Float64, Float64}) where F

Calculate pairwise relatedness between all combinations of the provided `sample_names` for each `method` provided. Bootstrapping resamples using
the `nonmissing` method, where resampling occurs over only shared non-missing loci. This is an internal function with all arguments provided by `relatedness`.
"""
function relatedness_boot_nonmissing(data::PopData, sample_names::Vector{String}; method::F, iterations::Int, interval::Tuple{Float64, Float64} = (0.025, 0.975), inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    uniq_pops = unique(data.meta.population)
    if first(uniq_pops) ∈ ["fullsib", "halfsib", "unrelated", "parent_offspring"]
        sample_pairs = sim_pairs(sample_names)
    else
        sample_pairs = pairwise_pairs(sample_names)
    end
    n_pairs = length(sample_pairs)
    n_samples = length(samples(data))
    n_per_loci = DataFrames.combine(groupby(data.loci, :locus), :genotype => nonmissing => :n)[:, :n]
    allele_frequencies = allele_freq(data)
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, n_pairs), 1:length(method))
    boot_means, boot_medians, boot_ses = map(i -> deepcopy(relate_vecs), 1:3)
    boot_CI = map(i -> Vector{Union{Missing,Tuple{Float64,Float64}}}(undef, n_pairs), 1:length(method))
    shared_loci = Vector{Int}(undef, n_pairs)
    p = Progress(n_pairs * length(method), dt = 1, color = :blue, barglyphs = BarGlyphs("|=> |"), barlen = 30)
    popdata_idx = groupby(data.loci, :name)
    @inbounds for i in 1:n_pairs
        @inbounds geno1 = popdata_idx[(sample_pairs[i][1],)].genotype
        @inbounds geno2 = popdata_idx[(sample_pairs[i][2],)].genotype
        # get index for genotype appearing missing in at least one individual in the pair
        keep_idx = nonmissings(geno1, geno2)
        # generate nonmissing genotype data
        gen1, gen2, loc, n_per_loc = (i[keep_idx] for i in (geno1, geno2, loci_names, n_per_loci))
        @inbounds shared_loci[i] = length(keep_idx)
        @inbounds @sync for (j, mthd) in enumerate(method)
            Base.Threads.@spawn begin
                @inbounds relate_vecs[j][i] = mthd(gen1, gen2, loc, allele_frequencies, loc_n = n_per_loc, n_samples = n_samples, inbreeding = inbreeding)
                boot_out = bootstrap_genos_nonmissing(gen1, gen2, loc, n_per_loc, allele_frequencies, method = mthd, iterations = iterations, inbreeding = inbreeding)
                @inbounds boot_means[j][i], boot_medians[j][i], boot_ses[j][i], boot_CI[j][i] = bootstrap_summary(boot_out, interval)
                pair_text = sample_pairs[i][1] * " × " * sample_pairs[i][2] * "  ($i" * " of " * "$n_pairs" * ")"
                ProgressMeter.next!(p; showvalues = [(:Pair, pair_text), (:Method, mthd)])
            end
        end
    end
    method_colnames = [Symbol("$i") for i in method]
    boot_mean_colnames = [Symbol("$i"*"_mean") for i in method]
    boot_median_colnames = [Symbol("$i"*"_median") for i in method]
    boot_se_colnames = [Symbol("$i"*"_SE") for i in method]
    CI_percent = convert(Int64, round(interval[2] - interval[1], digits = 2) * 100)
    boot_CI_colnames = [Symbol("$i"*"_CI_"*"$CI_percent") for i in method]

    out_dfs = map(method_colnames) do mthod
        DataFrame(:sample_1 => map(i -> i[1], sample_pairs), :sample_2 => map(i -> i[2], sample_pairs), :n_loci => shared_loci)
    end
    @inbounds for (i, mth) in enumerate(method_colnames)
        out_df[:, mth] = relate_vecs[i]
        out_df[:, boot_mean_colnames[i]] = boot_means[i]
        out_df[:, boot_median_colnames[i]] = boot_medians[i]
        out_df[:, boot_se_colnames[i]] = boot_ses[i]
        out_df[:, boot_CI_colnames[i]] = boot_CI[i]
    end
    return (; (method_colnames .=> out_dfs)...)
end

"""
    relatedness_no_boot(::PopData, sample_names::Vector{String}; method::F) where F

Calculate pairwise relatedness between all combinations of the provided `sample_names` for each `method` provided.
This is an internal function with arguments provided by `relatedness`.
"""
function relatedness_no_boot(data::PopData, sample_names::Vector{String}; method::F, inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    n_samples = length(samples(data))
    uniq_pops = unique(data.meta.population)
    if first(uniq_pops) ∈ ["fullsib", "halfsib", "unrelated", "parent_offspring"]
        sample_pairs = sim_pairs(sample_names)
    else
        sample_pairs = pairwise_pairs(sample_names)
    end
    n_pairs = length(sample_pairs)
    n_per_loci = DataFrames.combine(groupby(data.loci, :locus), :genotype => nonmissing => :n)[:, :n]
    allele_frequencies = allele_freq(data)
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, n_pairs), 1:length(method))
    shared_loci = Vector{Int}(undef, n_pairs)
    p = Progress(n_pairs* length(method), dt = 1, color = :blue, barglyphs = BarGlyphs("|=> |"), barlen = 30)
    popdata_idx = groupby(data.loci, :name)
    @inbounds @sync for i in 1:n_pairs
        Base.Threads.@spawn begin
            @inbounds geno1 = popdata_idx[(sample_pairs[i][1],)].genotype
            @inbounds geno2 = popdata_idx[(sample_pairs[i][2],)].genotype
            keep_idx = nonmissings(geno1, geno2)
            @inbounds shared_loci[i] = length(keep_idx)
            @inbounds for (j, mthd) in enumerate(method)
                @inbounds relate_vecs[j][i] = mthd(geno1[keep_idx], geno2[keep_idx], loci_names[keep_idx], allele_frequencies, loc_n = n_per_loci[keep_idx], n_samples = n_samples, inbreeding = inbreeding)
                pair_text = sample_pairs[i][1] * " × " * sample_pairs[i][2] * "  ($i" * " of " * "$(n_pairs)" * ")"
                ProgressMeter.next!(p; showvalues = [(:Pair, pair_text), (:Method, mthd)])
            end
        end
    end
    method_colnames = [Symbol("$i") for i in method]
    out_df = DataFrame(:sample_1 => getindex.(sample_pairs, 1), :sample_2 => getindex.(sample_pairs, 2), :n_loci => shared_loci)
    [out_df[:, mth] = relate_vecs[i] for (i, mth) in enumerate(method_colnames)]
    return out_df
end

"""
    # compare all samples
    relatedness(::PopData; method::Function, iterations::Int64, interval::Tuple{Float64, Float64}, resample::String, inbreeding::Bool = false)

```
# to compare specific samples
relatedness(::PopData, samples::Vector{String}; method::F, iterations::Int64, interval::Tuple{Float64, Float64}, resample::String, inbreeding::Bool = false)
```
Return a dataframe of pairwise relatedness estimates for all or select pairs of `samples` in a `PopData` object using
method(s) `F` where `F` is one or several of the methods listed below. If no bootstrapping is required, then the only
necessary keyword to provide is `method = ` and `inbreeding = ` for the `dyadicLikelihood` method (see examples below).
**Note:** samples must be diploid.

### Estimator methods
The available estimators are listed below and are functions themselves. `relatedness` takes the
function names as arguments (**case sensitive**), therefore do not use quotes or colons
in specifying the methods. Multiple methods can be supplied as a vector. All of these methods will tab-autocomplete.
For more information on a specific method, please see the respective docstring (e.g. `?Loiselle`).

- `Blouin`
- `LiHorvitz`
- `Loiselle`
- `Lynch`
- `LynchLi`
- `LynchRitland`
- `Moran`
- `QuellerGoodnight`
- `Ritland`

### Simulated siblingship comparison
If validating the estimators using `PopGenSims.jl` to simulate sibship relationships, `relatedness`
will recognize `PopData` generated in that manner (the `population` column) and only compare siblingship
pairs.

### Inbreeding
Use the `inbreeding` keyword to specify whether to allow inbreeding (`true`) or not (`false`, default).
This is only relevant for the `dyadicLikelihood` method (not yet released)

### Bootstrapping
To calculate means, medians, standard errors, and confidence intervals using bootstrapping,
set `iterations = n` where `n` is an integer greater than `0` (the default) corresponding to the number
of bootstrap iterations you wish to perform for each pair. The default confidence interval is `(0.05, 0.95)` (i.e. 90%),
however that can be changed by supplying the keyword `interval = (low, high)` where `low` and `high` are the intervals you want
(as `AbstractFloat`). Performing bootstrapping will return a NamedTuple of DataFrames, with each field in the NamedTuple
corresponding to the estimator `Method` it describes, which can be merged into one large dataframe using `merge_relatedness()`.

#### Resampling methods
There are two available resampling methods, `"all"` (default  & recommended) and `"nonmissing"`.
- `"all"` : resamples all loci for a pair of individuals and then drops missing loci between them
    - speed: slower
    - pro: better resampling variation
    - con: by chance some iterations may have a lot of missing loci that have to be dropped
- `"nonmissing"` : resamples only the shared non-missing loci between the pair
    - speed: faster
    - pro: every iteration guarantees the same number of loci compared between the pair
    - con: too-tight confidence intervals due to less possible variation

**Examples**
```
julia> cats = @nancycats;

julia> relatedness(cats, method = Ritland)
27966×4 DataFrame
   Row │ sample_1  sample_2  n_loci  Ritland
       │ String    String    Int64   Float64?
───────┼─────────────────────────────────────────
     1 │ N215      N216           8   0.258824
     2 │ N215      N217           8   0.193238
     3 │ N215      N218           8   0.127497
     4 │ N215      N219           8   0.0453471
   ⋮   │    ⋮         ⋮        ⋮          ⋮
 27963 │ N297      N290           7   0.189647
 27964 │ N281      N289           8   0.0892068
 27965 │ N281      N290           7   0.104614
 27966 │ N289      N290           7   0.0511663
                               27958 rows omitted

julia> relatedness(cats, ["N7", "N111", "N115"], method = [Ritland, Loiselle])
3×5 DataFrame
 Row │ sample_1  sample_2  n_loci  Ritland     Loiselle
     │ String    String    Int64   Float64?    Float64?
─────┼────────────────────────────────────────────────────
   1 │ N7        N111           9  -0.129432   -0.101618
   2 │ N7        N115           9  -0.0183925  -0.0428898
   3 │ N111      N115           9   0.0240152   0.13681

julia> rel_out = relatedness(cats, ["N7", "N111", "N115"], method = [Loiselle, Moran], iterations = 100, interval = (0.025, 0.975));

julia> rel_out.Loiselle
3×8 DataFrame
 Row │ sample_1  sample_2  n_loci  Loiselle    Loiselle_mean  Loiselle_median  Loiselle_SE  Loiselle_CI_95
     │ String    String    Int64   Float64?    Float64?       Float64?         Float64?     Tuple…?
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ N7        N111           9  -0.101618       0.0164593        0.0168983    0.0511942  (-0.0670707, 0.104008)
   2 │ N7        N115           9  -0.0428898      0.0771664        0.0880589    0.0982462  (-0.0428068, 0.183625)
   3 │ N111      N115           9   0.13681        0.255124         0.24991      0.267319   (0.103864, 0.421979)

# merge results into one big dataframe

julia> merge_relatedness(rel_out)
3×13 DataFrame
 Row │ sample_1  sample_2  n_loci  Loiselle    Loiselle_mean  Loiselle_ ⋯
     │ String    String    Int64   Float64?    Float64?       Float64?  ⋯
─────┼───────────────────────────────────────────────────────────────────
   1 │ N7        N111           9  -0.101618       0.0164593        0.0 ⋯
   2 │ N7        N115           9  -0.0428898      0.0771664        0.0
   3 │ N111      N115           9   0.13681        0.255124         0.2
                                                        8 columns omitted
```
"""
function relatedness(data::PopData, sample_names::Vector{String}; method::F, iterations::Int64 = 0, interval::Tuple{Float64, Float64} = (0.025, 0.975), resample::String = "all", inbreeding::Bool = false) where F
    all(data.meta[data.meta.name .∈ Ref(sample_names), :ploidy] .== 2) == false && error("Relatedness analyses currently only support diploid samples")
    errs = ""
    all_samples = samples(data)
    if sample_names != all_samples
        [errs *= "$i," for i in sample_names if i ∉ all_samples]
        errs != "" && error("Samples not found in the PopData: " * errs)
    end
    if eltype(method) != Function
        method = [method]
    end
    relate_mthds = [:QuellerGoodnight, :Ritland, :Lynch, :LynchLi, :LynchRitland, :Wang, :Loiselle, :Blouin, :Moran, :LiHorvitz, :dyadicLikelihood]
    [errs *= "$i is not a valid method\n" for i in Symbol.(method) if i ∉ relate_mthds]
    errs != "" && throw(ArgumentError(errs * "Methods are case-sensitive. Please see the docstring (?relatedness) for additional help."))
    if iterations > 0
        if resample == "all"
            relatedness_boot_all(data, sample_names, method = method, iterations = iterations, interval = interval, inbreeding = inbreeding)
        elseif resample == "nonmissing"
            relatedness_boot_nonmissing(data, sample_names, method = method, iterations = iterations, interval = interval, inbreeding = inbreeding)
        else
            throw(ArgumentError("Invalid resample method. Please choose from resample methods \"all\" or \"nonmissing\""))
        end
    else
        relatedness_no_boot(data, sample_names, method = method, inbreeding = inbreeding)
    end
end

function relatedness(data::PopData; method::F, iterations::Int64 = 0, interval::Tuple{Float64, Float64} = (0.025, 0.975), resample::String = "all", inbreeding::Bool = false) where F
    sample_names = samples(data) |> collect
    relatedness(data, sample_names, method = method, iterations = iterations, interval = interval, resample = resample, inbreeding = inbreeding)
end

"""
    merge_relatedness(data::NamedTuple)
A convenience function that takes the `NamedTuple` output from `relatedness` performed with bootstrapping
and returns one large DataFrame.
"""
function merge_relatedness(data::NamedTuple)
    k = keys(data)
    k1 = k[1]
    outdf = deepcopy(data[Symbol(k1)])
    for key in k[2:end]
        outdf = innerjoin(
            outdf,
            select(data[Symbol(key)], Not(:n_loci)),
            on = [:sample_1, :sample_2]
        )
    end
    return outdf
end


text:
# Example Configuration for nearly neutral simuation
# Command line to run:  julia n_neutral.jl examples/nn_example
# Creates output file examples/nn_example.csv

const nn_simtype = 1    # nearly neutral infinite alleles
mu_list_flag =true
type_str="I1"
@everywhere const N_list     = [100,200]   # popsize N
const N_mu_list= [1.0,2.0]   # Mutation rate as a multiple of 1.0/N
const mu_list = [0.1, 0.02]
const L = 1           # number of loci # not used in infinite alleles
const ngens = 1000      # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of N
const use_poplist=false   # save a list of populations as a check of correctness
dfe=dfe_neutral
dfe_str = "dfe_neutral"

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using Trixi
# Testumgebung, um die Korrektheit von q=∇u zu testen
# Test mit Polynomen funktioniert -> Achte aber auf stetige RB
# test mit Trigonometrischen Funktionen Funktioniert
# Krumme Gitter sehen jetzt auch gut aus.
# Fehler für Konstante Ableitung = 0 wird immer größer mit steigendem c, N... alles andere wird besser

######## notwendige funktionen die geladen werden müssen ##############################################
# mapping O-mesh
function mapping(xi_, eta_)
    ξ = xi_
    η = eta_
    x = 5 * (2 + ξ ) * cos(π * (η + 1))
    y = 5 * (2 + ξ ) * sin(π * (η + 1))
    return SVector(x, y)
end
function mapping1zu1(xi_, eta_)
  x = xi_
  y = eta_
  return SVector(x, y)
end
function mappingTri(xi_, eta_)
    ξ = xi_
    η = eta_
    x = cos(pi* ξ)
    y = η
    return SVector(x, y)
  end
  function mappingLin(xi_, eta_)
    ξ = xi_
    η = eta_
    x = 1.2* ξ
    y = η
    return SVector(x, y)
  end
function WR2_initial_condition_constant(x, t, equations::Trixi.AbstractEquations)
    rho = 1.0
    rho_v1 = 0.1
    rho_v2 = -0.2
    rho_e = 10.0
    return SVector(rho, rho_v1, rho_v2, rho_e)
end
function WR2_initial_condition_polynomial(x, t, equations::Trixi.AbstractEquations)
    return SVector(2, x[1] + x[2], x[1]^2, x[2]^2+x[1]^3)
end
function ABLx_initial_condition_polynomial(x, t, equations::Trixi.AbstractEquations)
    return SVector(0, 1, x[1]*2, 3*x[1]^2) end
function ABLy_initial_condition_polynomial(x, t, equations::Trixi.AbstractEquations)
    return SVector(0, 1, 0, 2*x[2]) end

function WR2_initial_condition_trigonometric(x, t, equations::Trixi.AbstractEquations)
    return SVector(sin(x[1]*π)/π, cos(x[1]*π)/π, sin(x[2]*π)/π + sin(x[1]*π)/π, sin(x[2]*π)/π + cos(x[2]*π)/π + sin(x[1]*π)/π+cos(x[1]*π)/π)
end
    # braucht hohe Auflösung !
function ABLx_initial_condition_trigonometric(x, t, equations::Trixi.AbstractEquations)
    return SVector(cos(x[1]*π), -sin(x[1]*π), cos(x[1]*π), cos(x[1]*π)-sin(x[1]*π)) end
function ABLy_initial_condition_trigonometric(x, t, equations::Trixi.AbstractEquations)
    return SVector(0, 0, cos(x[2]*π), cos(x[2]*π)-sin(x[2]*π)) end

function WR2_initial_condition_convergence_test(x, t, equations::Trixi.AbstractEquations)
    rho = sin(pi * x[1])
    rho_v1 = sin(pi * x[1])
    rho_v2 = sin(pi * x[1])
    rho_e = sin(pi * x[1])
    return SVector(rho, rho_v1, rho_v2, rho_e)
end
function ABLx_initial_condition_convergence_test(x, t, equations::Trixi.AbstractEquations)
    rho = pi * cos(pi * x[1])
    rho_v1 = pi * cos(pi * x[1])
    rho_v2 = pi * cos(pi * x[1])
    rho_e = pi * cos(pi * x[1])
    return SVector(rho, rho_v1, rho_v2, rho_e)
end
function ABLy_initial_condition_convergence_test(x, t, equations::Trixi.AbstractEquations)
    rho = 0
    rho_v1 = 0
    rho_v2 = 0
    rho_e = 0
    return SVector(rho, rho_v1, rho_v2, rho_e)
end
function boundary_condition_constant( u_inner, q1_inner, q2_inner, orientation, direction, x, t,
    surface_flux_function,
    equations::Trixi.AbstractEquations)
    # Far Field Conditions
    u_boundary = initial_condition(x , t, equations)

    # Calculate boundary flux
    if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
    else # direction == 4 # u_boundary is "left" of boundary, u_inner is "right" of boundary
        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
    end
    return flux
  end
function wrap_array(u_ode::AbstractVector, mesh::Union{TreeMesh{2},CurvedMesh{2},UnstructuredQuadMesh}, equations, dg::DG, cache)
    @boundscheck begin
        @assert length(u_ode) == nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)
    end
    unsafe_wrap(Array{eltype(u_ode),ndims(mesh) + 2}, pointer(u_ode),
                (nvariables(equations), nnodes(dg), nnodes(dg), nelements(dg, cache)))
end

###########################################################################################

# initial_condition = WR2_initial_condition_constant
# initial_condition = WR2_initial_condition_polynomial
initial_condition = WR2_initial_condition_trigonometric
# initial_condition =   WR2_initial_condition_convergence_test
# Zum Vergleich mit exakter Ableitung !
initial_condition2 = ABLx_initial_condition_trigonometric # ABLx_initial_condition_polynomial # ABLx_initial_condition_convergence_test #
initial_condition3 = ABLy_initial_condition_trigonometric # ABLy_initial_condition_polynomial # ABLy_initial_condition_convergence_test #
# boundary_conditions = boundary_condition_constant
# boundary_conditions = boundary_condition_periodic
boundary_conditions = (x_neg=boundary_condition_constant,
                       x_pos=boundary_condition_constant,
                       y_neg=boundary_condition_periodic,
                       y_pos=boundary_condition_periodic)

N = 4
Nq = 64
# coordinates_min = (0.0, 0.0)
# coordinates_max = (2.0,  2.0)
Zellen = [8, 16, 32, 64]

eq = Trixi.AuxiliaryEquation()
equations = CompressibleEulerEquations2D(1.4, viscous = true)

# for Nq in Zellen
cells_per_dimension = (Nq, Nq)
# mesh = CurvedMesh(cells_per_dimension, mapping1zu1, periodicity = false)
mesh = CurvedMesh(cells_per_dimension, mapping, periodicity = (false, true))
# mesh = CurvedMesh(cells_per_dimension, mappingLin, periodicity = true) #mappingTri

volume_integral = VolumeIntegralWeakForm()
surface_flux = flux_lax_friedrichs
basis = LobattoLegendreBasis(N)
dg = DGSEM(basis, surface_flux, volume_integral)

solver = DGSEM(basis, surface_flux, volume_integral)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)
ode = semidiscretize(semi, (0.0, 0.0))

u = wrap_array(ode.u0, mesh, equations, solver, semi.cache)

@unpack q1 = semi.cache
@unpack q2 = semi.cache
t = 0
Trixi.calc_nabla_u!(
    q1, q2, u, t, mesh, eq,
    boundary_conditions,
    dg, semi.cache)

# Zum vergleich mit exakter Ableitung !
solver = DGSEM(basis, surface_flux, volume_integral)
semi_x = SemidiscretizationHyperbolic(mesh, equations, initial_condition2, solver, boundary_conditions=boundary_conditions)
ode_x = semidiscretize(semi_x, (0.0, 0.0))
semi_y = SemidiscretizationHyperbolic(mesh, equations, initial_condition3, solver, boundary_conditions=boundary_conditions)
ode_y = semidiscretize(semi_y, (0.0, 0.0))

u_ABLx = wrap_array(ode_x.u0, mesh, equations, solver, semi.cache)
u_ABLy = wrap_array(ode_y.u0, mesh, equations, solver, semi.cache)

println("N = $N und Nq = $Nq")
for vgl = 1:4
    println("q1 cons.var. $vgl max error: ", maximum(abs.(u_ABLx[vgl, :, :, :]-q1[vgl, :, :, :])))
end
for vgl = 1:4
    println("q2 cons.var. $vgl max error: ", maximum(abs.(u_ABLy[vgl, :, :, :]-q2[vgl, :, :, :])))
end
# end

max_stars_repo_path: src/WR2/MS5 Auxiliary Convergence Test.jl


text:
export skipindex, @symmetrical

using SpecialFunctions
using Rocket

import Base: show, similar
import Base: IteratorSize, HasLength
import Base: IteratorEltype, HasEltype
import Base: eltype, length, size, sum
import Base: IndexStyle, IndexLinear, getindex

import Rocket: similar_typeof

struct SkipIndexIterator{T, I} <: AbstractVector{T}
    iterator :: I
    skip     :: Int
end

skip(iter::SkipIndexIterator) = iter.skip

function skipindex(iterator::I, skip::Int) where I
    @assert skip >= 1
    @assert length(iterator) >= 1
    return SkipIndexIterator{eltype(I), I}(iterator, skip)
end

Base.IteratorSize(::Type{<:SkipIndexIterator})   = HasLength()
Base.IteratorEltype(::Type{<:SkipIndexIterator}) = HasEltype()
Base.IndexStyle(::Type{<:SkipIndexIterator})     = IndexLinear()

Base.eltype(::Type{<:SkipIndexIterator{T}}) where T = T
Base.length(iter::SkipIndexIterator)                = length(iter.iterator) - 1
Base.size(iter::SkipIndexIterator)                  = (length(iter), )

Base.getindex(iter::SkipIndexIterator, i::Int)               = i < skip(iter) ? @inbounds(iter.iterator[i]) : @inbounds(iter.iterator[i + 1])
Base.getindex(iter::SkipIndexIterator, i::CartesianIndex{1}) = Base.getindex(iter, first(i.I))

Rocket.similar_typeof(::SkipIndexIterator, ::Type{L}) where L = Vector{L}

reduce_with_sum(array) = reduce(+, array)

##

import Base: +, -, *, /, convert, float, isfinite, isinf, zero, eltype

struct InfCountingReal{ T <: Real }
    value :: T
    infs  :: Int
end

InfCountingReal(value::T)            where { T <: Real } = InfCountingReal{T}(value, 0)
InfCountingReal(::Type{T}, inf::Int) where { T <: Real } = InfCountingReal{T}(zero(T), inf)

Infinity(::Type{T}) where { T <: Real } = InfCountingReal(T, 1)

const ∞ = Infinity(Float64)

value(a::InfCountingReal) = a.value
infs(a::InfCountingReal)  = a.infs

Base.isfinite(a::InfCountingReal) = infs(a) === 0
Base.isinf(a::InfCountingReal)    = !(isfinite(a))

Base.eltype(::Type{ InfCountingReal{T} }) where T = T
Base.eltype(::Type{ InfCountingReal })            = Real

Base.eltype(::T) where { T <: InfCountingReal } = eltype(T)

Base.:+(a::InfCountingReal) = InfCountingReal(+value(a), +infs(a))
Base.:-(a::InfCountingReal) = InfCountingReal(-value(a), -infs(a))

Base.:+(a::InfCountingReal, b::Real) = InfCountingReal(value(a) + b, infs(a))
Base.:-(a::InfCountingReal, b::Real) = InfCountingReal(value(a) - b, infs(a))
Base.:+(b::Real, a::InfCountingReal) = InfCountingReal(b + value(a), +infs(a))
Base.:-(b::Real, a::InfCountingReal) = InfCountingReal(b - value(a), -infs(a))

Base.:*(::InfCountingReal, ::Real) = error("InfCountingReal multiplication is dissalowed")
Base.:/(::InfCountingReal, ::Real) = error("InfCountingReal division is dissalowed")
Base.:*(::Real, ::InfCountingReal) = error("InfCountingReal multiplication is dissalowed")
Base.:/(::Real, ::InfCountingReal) = error("InfCountingReal division is dissalowed")

Base.:+(a::InfCountingReal, b::InfCountingReal) = InfCountingReal(value(a) + value(b), infs(a) + infs(b))
Base.:-(a::InfCountingReal, b::InfCountingReal) = InfCountingReal(value(a) - value(b), infs(a) - infs(b))

Base.convert(::Type{ InfCountingReal },    v::T)                  where { T <: Real }            = InfCountingReal(v)
Base.convert(::Type{ InfCountingReal{T} }, v::T)                  where { T <: Real }            = InfCountingReal(v)
Base.convert(::Type{ InfCountingReal{T} }, v::R)                  where { T <: Real, R <: Real } = InfCountingReal(convert(T, v))
Base.convert(::Type{ InfCountingReal{T} }, v::InfCountingReal{R}) where { T <: Real, R <: Real } = InfCountingReal{T}(convert(T, value(v)), infs(v))

Base.float(a::InfCountingReal) = isfinite(a) ? value(a) : Inf

Base.zero(::Type{InfCountingReal{T}}) where { T <: Real } = InfCountingReal(zero(T))

Base.show(io::IO, a::InfCountingReal{T}) where T = print(io, "InfCountingReal($(value(a)), $(infs(a))∞)")

Base.promote_rule(::Type{ InfCountingReal{T1} }, ::Type{ T2 }) where { T1 <: Real, T2 <: Real } = InfCountingReal{ promote_type(T1, T2) }
Base.promote_rule(::Type{ InfCountingReal },     ::Type{ T })  where { T <: Real }              = InfCountingReal{ T }

Base.:(==)(left::InfCountingReal{T}, right::InfCountingReal{T}) where { T } = (left.value == right.value) && (left.infs == right.infs)

# Union helpers

union_types(x::Union) = (x.a, union_types(x.b)...)
union_types(x::Type)  = (x,)

# Symbol helpers

__extract_val_type(::Type{ Type{ Val{ S } } }) where S = S
__extract_val_type(::Type{ Val{ S } })         where S = S

@generated function split_underscored_symbol(symbol_val)
    S = __extract_val_type(symbol_val)
    R = tuple(map(Symbol, split(string(S), "_"))...)
    return :(Val{ $R })
end

# NamedTuple helpers

fields(::NamedTuple{ F }) where F  = F
hasfield(field::Symbol, ntuple::NamedTuple) = field ∈ fields(ntuple)

function swapped(tuple::Tuple, i, j)
    @assert j > i
    return (tuple[1:i-1]..., tuple[j], tuple[i+1:j-1]..., tuple[i], tuple[j+1:end]...)
end

function swapped(array::AbstractArray, i, j)
    array = copy(array)
    array[i], array[j] = array[j], array[i]
    return array
end

##

__check_all(fn::Function, iterator)     = all(fn, iterator)
__check_all(fn::Function, tuple::Tuple) = TupleTools.prod(map(fn, tuple))
__check_all(fn::Function, ::Nothing)    = true

##

is_clamped_or_initial(something) = is_clamped(something) || is_initial(something)

## Meta utils

"""
    @symmetrical `function_definition`
Duplicate a method definition with the order of the first two arguments swapped.
This macro is used to duplicate methods that are symmetrical in their first two input arguments,
but require explicit definitions for the different argument orders.
Example:
    @symmetrical function prod!(x, y, z)
        ...
    end
"""
macro symmetrical(fn::Expr)
    # Check if macro is applied to a function definition
    # Valid function definitions include:
    # 1. foo([ args... ]) [ where ... [ where ... [ ... ] ] ] = :block
    # 2. function foo([ args... ]) [ where ... [ where ... [ ... ] ] ]
    #        :block
    #    end
    if (fn.head === :(=) || fn.head === :function) &&
        (fn.args[1] isa Expr && fn.args[2] isa Expr) &&
        (fn.args[2].head === :block)
        return esc(quote
            $fn
            $(swap_arguments(fn))
        end)
    else
        error("@symmetrical macro can be applied only to function definitions")
    end
end

function swap_arguments(fn::Expr)
    swapped = copy(fn)

    if swapped.args[1].head === :where
        swapped.args[1] = swap_arguments(swapped.args[1])
    elseif swapped.args[1].head === :call && length(fn.args[1].args) >= 3 # Note: >= 3, because the first argument is a function name
        swapped.args[1].args[2] = fn.args[1].args[3]
        swapped.args[1].args[3] = fn.args[1].args[2]
    else
        error("Function method passed for @symmetrical macro must have more than 2 arguments")
    end

    return swapped
end

## Other helpers

# We override this function for some specific types
function is_typeof_equal(left, right)
    _isequal = typeof(left) === typeof(right)
    if !_isequal
        @warn "typeof($left) !== typeof($right)"
    end
    return _isequal
end

function custom_isapprox(left, right; kwargs...)
    _isapprox = isapprox(left, right; kwargs...)
    if !_isapprox
        @warn "$left !≈ $right"
    end
    return _isapprox
end

custom_isapprox(left::NamedTuple, right::NamedTuple; kwargs...) = false

function custom_isapprox(left::NamedTuple{K}, right::NamedTuple{K}; kwargs...) where { K }
    _isapprox = true
    for key in keys(left)
        _isapprox = _isapprox && custom_isapprox(left[key], right[key]; kwargs...)
    end
    if !_isapprox
        @warn "$left !≈ $right"
    end
    return _isapprox
end

##

deep_eltype(::Type{ T })  where { T <: Number } = T
deep_eltype(::Type{ T })  where T               = deep_eltype(eltype(T))
deep_eltype(::T)          where T               = deep_eltype(T)

##

# See: https://github.com/JuliaLang/julia/issues/42795
function fill_bitarray!(V::SubArray{Bool, <:Any, <:BitArray, <:Tuple{AbstractUnitRange{Int}}}, x)
    B = V.parent
    I0 = V.indices[1]
    l0 = length(I0)
    l0 == 0 && return V
    Base.fill_chunks!(B.chunks, Bool(x), first(I0), l0)
    return V
end

##

default_if_nothing(::Nothing, default) = default
default_if_nothing(any::Any, default)  = any

##

forward_range(range::OrdinalRange) = step(range) > 0 ? range : last(range):first(range)


text:
y1 = zeros(5)
y2 = zeros(5,2)
print("Testing byrow: ")
y1,ny1,nt1 = byrow(y1)
y2,ny2,nt2 = byrow(y2)
if (ny1 == 1) && (ny2 == 2)
	println("OK")
end