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

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

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

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

function getloss(tree::ConsTree{TreeNode{T}}, λ::T, γ::T) where {T<:AbstractFloat}
    node = tree.value
    if isempty(tree.lefttree) && isempty(tree.righttree)
        getloss(node, λ, γ)
    else
        getloss(tree.lefttree, λ, γ) + getloss(tree.righttree, λ, γ)
    end
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 = zip4(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))), LinSpace(fromobs, toobs, nthreads + 1))
    ∂𝑙sum = [[(nodecansplit[node] ? [zero(T) for i in 1:(levelcounts[node])] : Vector{T}()) for node in 1:nodecount] for i in 1:nthreads]
    ∂²𝑙sum = [[(nodecansplit[node] ? [zero(T) for i in 1:(levelcounts[node])] : Vector{T}()) 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},
                            leftpartitions::Vector{Vector{Bool}}, factorindex::Vector{Int64},
                            fromobs::Integer, toobs::Integer, slicelength::Integer)
    if length(factors) == 0
        for i in fromobs:toobs
            nodeids[i] = 2 * nodeids[i] - 1
        end
    else
        factorslices = zipn([slice(factor, fromobs, toobs, slicelength) for factor in factors])
        nodeslices = slice(nodeids, fromobs, toobs, slicelength)
        foreach(zip2(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]) ? (2 * nodeslice[i] - 1) : (2 * nodeslice[i])
                else
                    nodeslice[i] = 2 * nodeslice[i] - 1
                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) for n in nodes]
    factors = Vector{AbstractFactor}()
    factorindex = Vector{Int64}(nodecount)
    for i in 1:nodecount
         if issplitnode[i]
             factor = nodes[i].factor
             index = findfirst(factors, factor)
             if index == 0
                 push!(factors, factor)
             end
             factorindex[i] = findfirst(factors, factor)
         end
    end
    factors = widenfactors(factors)
    leftpartitions = [isa(n, SplitNode) ? [n.leftpartition.inclmissing; n.leftpartition.mask] : Vector{Bool}() for n in nodes]

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

function getsplitnode(factor::AbstractFactor, partition::LevelPartition, gradient::Vector{LossGradient{T}},
                      λ::T, γ::T, min∂²𝑙::T) where {T<:AbstractFloat}

    isord = isordinal(factor)
    gradstart = 2
    ∂𝑙sum0 = sum((grad -> grad.∂𝑙), gradient[gradstart:end])
    ∂²𝑙sum0 = sum((grad -> grad.∂²𝑙), gradient[gradstart:end])
    miss∂𝑙 = gradient[1].∂𝑙
    miss∂²𝑙 = gradient[1].∂²𝑙
    bestloss = typemax(T) #getloss(∂𝑙sum0 + miss∂𝑙, ∂²𝑙sum0 + miss∂²𝑙, λ, γ)
    levelcount = length(partition.mask)
    split = SplitNode{T}(factor, LevelPartition(copy(partition.mask), partition.inclmissing), LevelPartition(zeros(Bool, levelcount), !partition.inclmissing),
                         LossGradient{T}(∂𝑙sum0 + miss∂𝑙, ∂²𝑙sum0 + miss∂²𝑙), LossGradient{T}(zero(T), zero(T)),
                         bestloss)

    left∂𝑙sum = gradient[gradstart].∂𝑙
    left∂²𝑙sum = gradient[gradstart].∂²𝑙

    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.leftgradient.∂𝑙 = left∂𝑙sum
            split.leftgradient.∂²𝑙 = left∂²𝑙sum
            split.rightgradient.∂𝑙 = ∂𝑙sum0 - left∂𝑙sum + miss∂𝑙
            split.rightgradient.∂²𝑙 = ∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙
            for j in 1:levelcount
                split.leftpartition.mask[j] = j == 1
                split.rightpartition.mask[j] = j == 1 ? false : partition.mask[j]
            end
            split.leftpartition.inclmissing = false
            split.rightpartition.inclmissing = partition.inclmissing
            split.loss = firstlevelwitouthmiss
        else
            split.leftgradient.∂𝑙 = left∂𝑙sum + miss∂𝑙
            split.leftgradient.∂²𝑙 = left∂²𝑙sum + miss∂²𝑙
            split.rightgradient.∂𝑙 = ∂𝑙sum0 - left∂𝑙sum
            split.rightgradient.∂²𝑙 = ∂²𝑙sum0 - left∂²𝑙sum
            for j in 1:levelcount
                split.leftpartition.mask[j] = j == 1
                split.rightpartition.mask[j] = j == 1 ? false : partition.mask[j]
            end
            split.leftpartition.inclmissing = partition.inclmissing
            split.rightpartition.inclmissing = false
            split.loss = firstlevelwithmiss
        end
    end

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

        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 leftwithmisstotal < split.loss && (left∂²𝑙sum + miss∂²𝑙 >= min∂²𝑙) && (∂²𝑙sum0 - left∂²𝑙sum >= min∂²𝑙)
                if leftwithoutmisstotal < leftwithmisstotal && (left∂²𝑙sum >= min∂²𝑙) && (∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙 >= min∂²𝑙)
                    split.leftgradient.∂𝑙 = left∂𝑙sum
                    split.leftgradient.∂²𝑙 = left∂²𝑙sum
                    split.rightgradient.∂𝑙 = ∂𝑙sum0 - left∂𝑙sum + miss∂𝑙
                    split.rightgradient.∂²𝑙 = ∂²𝑙sum0 - left∂²𝑙sum + miss∂²𝑙
                    for j in 1:levelcount
                        split.leftpartition.mask[j] = j <= i ? partition.mask[j] : false
                        split.rightpartition.mask[j] = j <= i ? false : partition.mask[j]
                    end
                    split.leftpartition.inclmissing = false
                    split.rightpartition.inclmissing = partition.inclmissing
                    split.loss = leftwithoutmisstotal
                else
                    split.leftgradient.∂𝑙 = left∂𝑙sum + miss∂𝑙
                    split.leftgradient.∂²𝑙 = left∂²𝑙sum + miss∂²𝑙
                    split.rightgradient.∂𝑙 = ∂𝑙sum0 - left∂𝑙sum
                    split.rightgradient.∂²𝑙 = ∂²𝑙sum0 - left∂²𝑙sum
                    for j in 1:levelcount
                        split.leftpartition.mask[j] = j <= i ? partition.mask[j] : false
                        split.rightpartition.mask[j] = j <= i ? false : partition.mask[j]
                    end
                    split.leftpartition.inclmissing = partition.inclmissing
                    split.rightpartition.inclmissing = false
                    split.loss = leftwithmisstotal
                end
            end
        else
            if singlelevelwithmisstotal < split.loss && (∂²𝑙 + miss∂²𝑙 >= min∂²𝑙) && (∂²𝑙sum0 - ∂²𝑙 >= min∂²𝑙)
                if singlelevelwitouthmisstotal < singlelevelwithmisstotal && (∂²𝑙 >= min∂²𝑙) && (∂²𝑙sum0 - ∂²𝑙 + miss∂²𝑙 >= min∂²𝑙)
                    split.leftgradient.∂𝑙 = ∂𝑙
                    split.leftgradient.∂²𝑙 =  ∂²𝑙
                    split.rightgradient.∂𝑙 = ∂𝑙sum0 - ∂𝑙 + miss∂𝑙
                    split.rightgradient.∂²𝑙 = ∂²𝑙sum0 - ∂²𝑙 + miss∂²𝑙
                    for j in 1:levelcount
                        split.leftpartition.mask[j] = j == i
                        split.rightpartition.mask[j] = j == i ? false : partition.mask[j]
                    end
                    split.leftpartition.inclmissing = false
                    split.rightpartition.inclmissing = partition.inclmissing
                    split.loss = singlelevelwitouthmisstotal
                else
                    split.leftgradient.∂𝑙 = ∂𝑙 + miss∂𝑙
                    split.leftgradient.∂²𝑙 =  ∂²𝑙 + miss∂²𝑙
                    split.rightgradient.∂𝑙 = ∂𝑙sum0 - ∂𝑙
                    split.rightgradient.∂²𝑙 = ∂²𝑙sum0 - ∂²𝑙
                    for j in 1:levelcount
                        split.leftpartition.mask[j] = j == i
                        split.rightpartition.mask[j] = j == i ? false : partition.mask[j]
                    end
                    split.leftpartition.inclmissing = partition.inclmissing
                    split.rightpartition.inclmissing = false
                    split.loss = singlelevelwithmisstotal
                end
            end
        end
    end
    if count(split.rightpartition.mask) > 0 && split.loss < typemax(T)
        Nullable{SplitNode{T}}(split)
    else
        Nullable{SplitNode{T}}()
    end
end

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

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

    nodecansplit = [n.cansplit for n in state.nodes]
    # factorcount = length(state.factors)
    # nthreads = state.singlethread ? 1 : Threads.nthreads()
    # newsplits = Vector{Vector{Nullable{SplitNode{T}}}}(factorcount)

    # if nthreads == 1
    #     for i in 1:factorcount
    #         factor = state.factors[i]
    #         partitions = [node.partitions[factor] for node in state.nodes]
    #         gradient = sumgradient(state.nodeids, nodecansplit, factor, partitions, state.∂𝑙covariate, state.∂²𝑙covariate, state.slicelength, true)
    #         newsplits[i] = getnewsplit(gradient, state.nodes, factor, state.λ, state.γ, state.min∂²𝑙, true)
    #     end
    # else
    #     Threads.@threads for i in 1:factorcount
    #         factor = state.factors[i]
    #         partitions = [node.partitions[factor] for node in state.nodes]
    #         gradient = sumgradient(state.nodeids, nodecansplit, factor, partitions, state.∂𝑙covariate, state.∂²𝑙covariate, state.slicelength, true)
    #         newsplits[i] = getnewsplit(gradient, state.nodes, factor, state.λ, state.γ, state.min∂²𝑙, true)
    #     end
    # end

    # foldl(state.nodes, newsplits) do currsplit, newsplit
    #     res = Vector{TreeNode{T}}(length(newsplit))
    #     for i in 1:length(newsplit)
    #          if !isnull(newsplit[i]) &&
    #             (isa(currsplit[i], LeafNode{T}) || get(newsplit[i]).loss < getloss(currsplit[i], state.λ, state.γ))
    #             res[i] = get(newsplit[i])
    #          else
    #             res[i] = currsplit[i]
    #          end
    #     end
    #     res
    # end

    foldl(state.nodes, enumerate(state.factors)) do currsplit, nfactor
        n, factor = nfactor
        partitions = [node.partitions[factor] 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.singlethread)

        res = Vector{TreeNode{T}}(length(newsplit))
        @inbounds for i in 1:length(newsplit)
             if !isnull(newsplit[i]) &&
                (isa(currsplit[i], LeafNode{T}) || get(newsplit[i]).loss < getloss(currsplit[i], state.λ, state.γ))
                res[i] = get(newsplit[i])
             else
                res[i] = currsplit[i]
             end
        end
        res
    end
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{LeafNode{T}}(2 * length(state.nodes))
    @inbounds for (i, n) in enumerate(layer.nodes)
        if isa(n, SplitNode)
            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
            newnodes[2 * i - 1] = LeafNode{T}(n.leftgradient,
                                              n.leftgradient.∂²𝑙 >= state.min∂²𝑙,
                                              leftpartitions)
            newnodes[2 * i] = LeafNode{T}(n.rightgradient,
                                          n.rightgradient.∂²𝑙 >= state.min∂²𝑙,
                                          rightpartitions)
        else
            newnodes[2 * i - 1] = LeafNode{T}(n.gradient, false, n.partitions)
            newnodes[2 * i] = LeafNode{T}(n.gradient, false, n.partitions)
        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)
    Nullable{TreeLayer{T}}(layer), state
end

function predict(treelayer::TreeLayer{T}, nodeids::Vector{<:Integer}, λ::T) where {T<:AbstractFloat}
    weights = Vector{T}(2 * length(treelayer.nodes))
    @inbounds for (i, node) in enumerate(treelayer.nodes)
        if isa(node, SplitNode)
            weights[2 * i - 1] = getweight(node.leftgradient, λ)
            weights[2 * i] = getweight(node.rightgradient, λ)
        else
            weights[2 * i - 1] = getweight(node.gradient, λ)
            weights[2 * i] = 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(tolevels, level)
        if j > 0
            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(fromlevels, level)
        if j == 0
            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}

    factor, levelmap, newind, levelcount = factormap[node.factor]
    leftmask = Vector{Bool}(levelcount)
    rightmask = Vector{Bool}(levelcount)
    for (i, j) in levelmap
        leftmask[j] = node.leftpartition.mask[i]
        rightmask[j] = node.rightpartition.mask[i]
    end
    for i in newind
        leftmask[i] = false
        rightmask[i] = true
    end
    SplitNode{T}(factor, LevelPartition(leftmask, node.leftpartition.inclmissing),
                 LevelPartition(rightmask, node.rightpartition.inclmissing),
                 node.leftgradient, node.rightgradient, node.loss)
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 = 2 ^ tree.maxdepth
    nodeids = maxnodecount <= typemax(UInt8) ? ones(UInt8, len) : (maxnodecount <= typemax(UInt16) ? ones(UInt16, len) : ones(UInt32, len))
    nodes = Vector{TreeNode{T}}()
    factormap = Dict{AbstractFactor, Tuple{AbstractFactor, Dict{Int64, Int64}, Set{Int64}, Int64}}()
    for layer in tree.layers
        for node in layer.nodes
            if isa(node, SplitNode) && !(node.factor in keys(factormap))
                mappedfactor = map(node.factor, dataframe) |> cache
                levelmap = getlevelmap(node.factor, mappedfactor)
                newind = getnewindices(node.factor, mappedfactor)
                levelcount = length(getlevels(mappedfactor))
                factormap[node.factor] = (mappedfactor, levelmap, newind, levelcount)
            end
        end
    end

    for layer in tree.layers
        nodes = [map(n, dataframe, factormap) for n in layer.nodes]
        #nodes = [isa(n, SplitNode) ? SplitNode{T}(factormap[n.factor], n.leftpartition, n.rightpartition, n.leftgradient, n.rightgradient, n.loss) : n for n in layer.nodes]
        splitnodeids!(nodeids, TreeLayer{T}(nodes), tree.slicelength, tree.singlethread)
    end
    predict(TreeLayer{T}(nodes), nodeids, tree.λ)
end

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

    len = length(∂𝑙covariate)
    maxnodecount = 2 ^ 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][1]
    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∂²𝑙, slicelength, singlethread)
    layers = collect(Iterators.take(Seq(TreeLayer{T}, state0, nextlayer), maxdepth))
    xgtree = XGTree{T}(layers, λ, γ, min∂²𝑙, maxdepth, slicelength, singlethread)
    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, slicelength, singlethread)
    pred = predict(xgtree.layers[end], nodeids, λ)
    xgtree, pred
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]), zip2(left, right), List{TreeNode{T}}))
end

function prune(tree::ConsTree{<:TreeNode{T}}, λ::T, γ::T) where {T<:AbstractFloat}
    node = tree.value
    if isempty(tree.lefttree) && isempty(tree.righttree)
        tree
    else
        left = prune(tree.lefttree, λ, γ)
        right = prune(tree.righttree, λ, γ)
        if getloss(left, λ, γ) + getloss(right, λ, γ) < getloss(node, λ, γ)
            ConsTree{TreeNode{T}}(node, left, right)
        else
            ConsTree{TreeNode{T}}(node)
        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: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
"""
    MaternKernel(; ν = 1.0)

The matern kernel is a Mercer kernel given by the formula:
```
    κ(x,y) = 2^{1-ν}/Γ(ν)*(√(2ν)‖x-y‖)^ν K_ν(√(2ν)‖x-y‖)
```
For `ν=n+1/2, n=0,1,2,...` it can be simplified and you should instead use
[`ExponentialKernel`](@ref) for `n=0`, [`Matern32Kernel`](@ref), for `n=1`,
[`Matern52Kernel`](@ref) for `n=2` and [`SqExponentialKernel`](@ref) for `n=∞`.
"""
struct MaternKernel{Tν<:Real} <: SimpleKernel
    ν::Vector{Tν}
    function MaternKernel(;nu::T=1.5, ν::T=nu) where {T<:Real}
        @check_args(MaternKernel, ν, ν > zero(T), "ν > 0")
        return new{T}([ν])
    end
end

@functor MaternKernel

@inline function kappa(κ::MaternKernel, d::Real)
    result = _matern(first(κ.ν), d)
    return ifelse(iszero(d), one(result), result)
end

function _matern(ν::Real, d::Real)
    y = sqrt(2ν) * d
    return exp((one(d) - ν) * logtwo - loggamma(ν) + ν * log(y) + log(besselk(ν, y)))
end

metric(::MaternKernel) = Euclidean()

Base.show(io::IO, κ::MaternKernel) = print(io, "Matern Kernel (ν = ", first(κ.ν), ")")

"""
    Matern32Kernel()

The matern 3/2 kernel is a Mercer kernel given by the formula:
```
    κ(x,y) = (1+√(3)‖x-y‖)exp(-√(3)‖x-y‖)
```
"""
struct Matern32Kernel <: SimpleKernel end

kappa(κ::Matern32Kernel, d::Real) = (1 + sqrt(3) * d) * exp(-sqrt(3) * d)

metric(::Matern32Kernel) = Euclidean()

Base.show(io::IO, ::Matern32Kernel) = print(io, "Matern 3/2 Kernel")

"""
    Matern52Kernel()

The matern 5/2 kernel is a Mercer kernel given by the formula:
```
    κ(x,y) = (1+√(5)‖x-y‖ + 5/3‖x-y‖^2)exp(-√(5)‖x-y‖)
```
"""
struct Matern52Kernel <: SimpleKernel end

kappa(κ::Matern52Kernel, d::Real) = (1 + sqrt(5) * d + 5 * d^2 / 3) * exp(-sqrt(5) * d)

metric(::Matern52Kernel) = Euclidean()

Base.show(io::IO, ::Matern52Kernel) = print(io, "Matern 5/2 Kernel")

ext: jl


text:
"""
my_wong(;kwargs)
Description of ```my_wong```
------------------------------
sets custom theme for Plots.

    If no keyword arguments are provided, the default has the following attributes:

+ :wong2,
+ legend =:false,
+ left_margin = 8Plots.mm,
+ right_margin = 8Plots.mm,
+ top_margin  = 8Plots.mm,
+ bottom_margin  = 8Plots.mm,
+ framestyle = :box,
+ grid = :off,
+ size = (1200,600),
+ dpi = 200,
+ color = "blue"
+ titlefontsize = 24
+ guidefontsize = 16
+ tickfontsize  = 16
+ legendfontsize= 16

-------------------------------
    To change any of the default arguments, simply use the standard <keyword = value>
(i.e.) my_wong(color = red)
"""
function  my_wong(;kwargs...)
    ## default values for the plotting style
    plot_args = Dict(
        :legend       => :false,
        :left_margin   => 8Plots.mm,
        :right_margin  => 8Plots.mm,
        :top_margin    => 8Plots.mm,
        :bottom_margin => 8Plots.mm,
        :framestyle    => :box,
        :grid          => :off,
        :size          => (1200,600),
        :dpi           => 200,
        :color         => "blue",
        :titlefontsize => 24,
        :guidefontsize => 16,
        :tickfontsize  => 16,
        :legendfontsize=> 16,
        :minorgrid     => :true,
        :minorticks    => 10,
        )

    ## checks if kwargs are provided
    if length(kwargs) == 0                ## if no kwargs are provided, use my default theme
        Plots.theme(:wong2; plot_args...)
    else                                  ## if there are kwargs, they are unpacked using keys()
        kwarg_names = keys(kwargs)        ## and iterated over to change the default dictionary <plot_args>
        for key in kwarg_names
            plot_args[key] = kwargs[key]
        end
        Plots.theme(:wong2; plot_args...)
    end
end


text:
using Documenter, ModelingToolkit, DiffEqBase

makedocs(
    sitename="ModelingToolkit.jl",
    modules=[ModelingToolkit],
    pages=[
        "Home" => "index.md",
        "api.md",
    ]
)


text:
using Ledger
using Test

@testset "Ledger.jl" begin

    # Extending the interface for my own types/valuation basis

    struct MyValuationBasis <: Ledger.ValuationBasis end

    struct MyBond <: Ledger.Asset
        par
    end

    Ledger.value(::MyValuationBasis,b::MyBond) = b.par

    gl = GeneralLedger()

    addAsset!(gl,MyBond(100))

    bal = Balance(MyValuationBasis(),gl)

    @test sum(bal.assets) == 100.0
    @test length(bal.liabilities) == 0
    @test sum(bal.equity) == 100.0

    struct myDebt <: Ledger.Liability
        balance
        int_rate
        years_due
    end

    Ledger.value(::MyValuationBasis,b::myDebt) = b.balance / (1+ b.int_rate) ^ b.years_due

    addLiability!(gl,myDebt(50,0.05,3))

    bal = Balance(MyValuationBasis(),gl)

    @test sum(bal.assets) == 100.0
    @test sum(bal.liabilities) == 50 / 1.05 ^ 3
    @test sum(bal.equity) == 100.0 - 50 / 1.05 ^ 3

end

max_stars_repo_name: alecloudenback/Ledger

max_issues_repo_name: alecloudenback/Ledger

max_forks_repo_name: alecloudenback/Ledger


text:
module GlassCatalog

include("parse-library.jl")
include("dispersion.jl")

end # module GlassCatalog


text:
## -------------------------------------------------------------------
# Collect
DAT = ChU.MagicDict()
let

    DATfile = joinpath(iJR.MODEL_PROCESSED_DATA_DIR, "2.1_DAT.jls")
    # CACHE
    if isfile(DATfile)
        global DAT = deserialize(DATfile)
        @info("DAT CACHE LOADED")
        return
    end

    objider = iJR.BIOMASS_IDER
    DAT[:EXCH_FLX_IDERS] = EXCH_FLX_IDERS
    DAT[:EXPS] = []

    # Find exps
    for exp in 1:4
        ok = false
        for method in ALL_METHODS
            ok = haskey(INDEX, method, :DFILE, exp) &&
                INDEX[method, :DFILE, exp] != :unfeasible
            !ok && break
        end
        !ok && continue
        push!(DAT[:EXPS], exp)
    end
    max_model = iJR.load_model("max_model"; uncompress = false)

    for exp in DAT[:EXPS], method in ALL_METHODS

        # !haskey(INDEX, method, :DFILE, exp) && continue
        datfile = INDEX[method, :DFILE, exp]
        dat = deserialize(datfile)

        model = dat[:model]
        objidx = ChU.rxnindex(model, objider)
        epouts = dat[:epouts]
        exp_beta = maximum(keys(epouts))
        epout = epouts[exp_beta]
        exp_xi = Fd.val(:xi, exp)
        fva_model = iJR.load_model("fva_models", exp; uncompress = false)

        println()
        @info("Doing", exp, method, length(dat[:epouts]), epout.iter);

        # Biomass
        ep_biom = ChU.av(model, epout, objidx)
        ep_std = sqrt(ChU.va(model, epout, objidx))
        Fd_biom = Fd.val("D", exp)
        Fd_biom_err = Fd.err("D", exp)
        max_lb, max_ub = ChU.bounds(max_model, objidx)
        fva_lb, fva_ub = ChU.bounds(fva_model, objidx)
        lb = max(max_lb, fva_lb)
        ub = min(max_ub, fva_ub)

        # store
        DAT[method, :ep   , :flx, objider, exp] = ep_biom
        DAT[method, :eperr, :flx, objider, exp] = ep_std
        DAT[method, :Fd   , :flx, objider, exp] = Fd_biom
        DAT[method, :Fderr, :flx, objider, exp] = Fd_biom_err
        DAT[:Fd   , :flx, objider, exp] = Fd_biom
        DAT[:Fderr, :flx, objider, exp] = Fd_biom_err
        DAT[method, :bounds, :flx, objider, exp] = (lb, ub)

        # fluxes
        for Fd_met in EXCH_FLX_IDERS

                model_met = Fd_mets_map[Fd_met]
                model_exch = Fd_rxns_map[Fd_met]
                model_exchi = ChU.rxnindex(model, model_exch)

                proj = ChLP.projection2D(model, objider, model_exchi; l = 50)
                ep_av = ChU.av(model, epout, model_exchi)
                ep_std = sqrt(ChU.va(model, epout, model_exchi))
                Fd_flx = Fd.val("u$Fd_met", exp)
                Fd_err = Fd.err("u$Fd_met", exp)

                max_lb, max_ub = ChU.bounds(max_model, Fd_rxns_map[Fd_met])
                fva_lb, fva_ub = ChU.bounds(fva_model, Fd_rxns_map[Fd_met])
                lb = max(max_lb, fva_lb)
                ub = min(max_ub, fva_ub)

                DAT[method, :ep, :proj, Fd_met, exp] = proj
                DAT[method, :Fd, :flx, Fd_met, exp] = Fd_flx
                DAT[method, :Fderr, :flx, Fd_met, exp] = Fd_err
                DAT[:Fd, :flx, Fd_met, exp] = Fd_flx
                DAT[:Fderr, :flx, Fd_met, exp] = Fd_err
                DAT[method, :ep, :flx, Fd_met, exp] = ep_av
                DAT[method, :eperr, :flx, Fd_met, exp] = ep_std

                DAT[method, :bounds, :flx, Fd_met, exp] = (lb, ub)

        end # for Fd_met
    end # for exp in EXPS, for method

    DAT[:EXPS] |> unique! |> sort!
    serialize(DATfile, DAT)
end;

# -------------------------------------------------------------------
# Inter project comunication
let
    CORR_DAT = isfile(iJR.CORR_DAT_FILE) ? ChU.load_data(iJR.CORR_DAT_FILE) : Dict()
    CORR_DAT[:MAXENT_EP] = DAT
    ChU.save_data(iJR.CORR_DAT_FILE, CORR_DAT)
end


text:
module TestJSON
import JSON

mutable struct S
	a
	b
end

mutable struct T
	a
	b
end

s = S(1,T(10, 25))

println("JSON.print(s):")
JSON.print(s)
println()

n = (a=1, b=(a=10, b=25))

println("JSON.print(n):")
JSON.print(n)
println()

println("JSON.print(s) == JSON.print(n):")
println(JSON.print(s) == JSON.print(n))

e = :(x = sin(time)*(a+1))

println()
JSON.print(e)
println()

println()
ej = JSON.json(e)
@show ej

println()
ee = JSON.parse(ej)
@show eval(ee)

end


text:
using Compat

if VERSION.minor == 7
  @eval using LinearAlgebra, SparseMatrices
end

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

A = Float64[
  1/2 1/3 1/4 1/5 1/6;
  1/3 1/4 1/5 1/6 1/7;
  1/4 1/5 1/6 1/7 1/8;
  1/5 1/6 1/7 1/8 1/9;
  1/6 1/7 1/8 1/9 1/10;
]
b = ones(5)

(Al, Au, Ap) = lu(A)

@show y = Al \ b[Ap]
println()

@show x = Au \ y
println()

@show As = sparse(A)
println()

@show As\b
println()

@show F = lufact(As)
println()

if VERSION.minor == 7
  @show F.L * F.U == F.Rs .* A[F.p, :]
  println()

  @show F.L*F.U
  println()

  @show sparse(F.Rs.* A[F.p, F.q])
  println()
else
  @show ys = F[:L] \ (F[:Rs] .* b)
  println()

  @show x = F[:U] \ ys
  println()

  @show F[:L]*F[:U] == (F[:Rs] .* A)[F[:p], F[:q]]
  println()

  @show F[:L]*F[:U]
  println()

  @show sparse((F[:Rs] .* A)[F[:p], F[:q]])
  println()
end

@show lufact(A)
println()

end


text:
module Nettle

include( "../deps/deps.jl")
include( "hash.jl" )
include( "hmac.jl" )

# similar to Python's hmac.HMAC.hexdigest
hexdigest!(state::Union(HMACState,HashState)) = bytes2hex(digest!(state))

end


text:
using Test, DataFrames, PairsMacros
df = DataFrame(x = [1, 2], y = [3, 4], z = [5, 6])

# one arg
@test (@cols(z = sum(x))) == (:x => sum => :z)
@test (@cols(sum(x))) == (:x => sum)
@test (@cols(z = x)) == (:x => :z)
@test (@cols(x)) == :x

@test select(df, @cols(x)).x == df.x
@test (@cols(z = sum(skipmissing(x)))) == (:x => sum ∘ skipmissing => :z)
@test transform(df, @cols(z = exp.(x))).z == exp.(df.x)
@test transform(df, @cols(z = 1)).z == fill(1, size(df, 1))
@test transform(df, @cols(z = sum([1, 2, 3]))).z == fill(sum([1, 2, 3]), size(df, 1))
@test size(filter(@cols(x > 1), df), 1) == 1
@test size(filter(@cols((x > 1) & (y < 3)), df), 1) == 0

# multiple args
@test transform(df, @cols sum(x) first(x)).x_sum == [3, 3]
@test transform(df, @cols xfirst = first(x) xlast = last(x)).xfirst == [1, 1]

@test transform(df, @rows(z = exp.(x))).z == exp.(df.x)
@test isequal(transform(df, @rows(z2 = x >= 2 ? missing : z)).z2, [5, missing])

# rows
@test (@rows(z = rand())) == (Any[] => (DataFrames.ByRow{typeof(rand)}(rand) => :z))
@test transform(df, @rows(z = x - y)).z  == df.x .- df.y
@test transform(df, @rows(z = y - x)).z  == df.y .- df.x
@test transform(df, @rows(z = x - x)).z  == df.x .- df.x
@test transform(df, @rows abs(x) x^2).x_abs == abs.(df.x)
@test transform(df, @rows(z = x^2)).z  == df.x.^2
@test transform(df, @rows abs(x) x^2).x_abs  == abs.(df.x)

df = DataFrame(x = [[-1, -2], [-3, -4]])
@test transform(df, @rows(z = abs.(x))).z == [abs.(x) for x in df.x]

# test $
u = :y
@test (@cols(z = sum($u))) == (:y => sum => :z)
@test (@cols(z = sum($u))) == (:y => sum => :z)
@test (@rows($u = x)) == (:x => :y)
@test (@cols(z = sum(skipmissing($u)))) == (:y => sum ∘ skipmissing => :z)
u = [:y]
@test (@cols(z = sum($(u[1])))) == (:y => sum => :z)

# test obj
u = [3, 4]
@test combine(df, @cols(z = sum(esc(u)))).z == [sum(u)]

# test hygiene
mean = x -> x
@test transform(df, @cols(z = mean(x))).z == df.x

# weird symbols
df = DataFrame(x = [1, 2])
@test transform!(df, @cols($("ok now") = x))."ok now" == df.x
@test transform(df, @cols(z = sum($("ok now")))).z == [3, 3]


text:
module SQS

import URIParser

using Requests
using LightXML
using AWS.AWSEnv
using AWS

type SQSError
    typ::AbstractString
    code::AbstractString
    msg::AbstractString
    detail::AbstractString
    request_id::Union{AbstractString, Void}
end
export SQSError

sqs_error_str(o::SQSError) = "type: $(o.typ), code: $(o.code), msg : $(o.msg), $(o.request_id)"
export sqs_error_str

type SQSResponse
    http_code::Int
    headers
    body::Union{AbstractString, Void}
    pd::Union{LightXML.XMLElement, Void}
    obj::Any

    SQSResponse() = new(0, Dict{Any, Any}(), "", nothing, nothing)
end
export SQSResponse

function sqs_execute(env_::AWSEnv, action::AbstractString, ep, params_in, use_post)
    # Adjust endpoint.
    env = env_
    if ep != nothing
        env = AWSEnv(env_; ep=ep)
    end

    # Prepare the standard params
    params = copy(params_in)

    push!(params, ("Action", action))
    push!(params, ("AWSAccessKeyId", env.aws_id))
    push!(params, ("Timestamp", get_utc_timestamp()))
    push!(params, ("Version", "2012-11-05"))

    amz_headers, data, signed_querystr = canonicalize_and_sign(env, "sqs", use_post, params)

    complete_url = "https://" * ep_host(env, "sqs") * env.ep_path
    if length(signed_querystr) > 0
        complete_url *= "?" * signed_querystr
    end
    if (env.dbg) || (env.dry_run)
        println("URL:\n$complete_url\n")
    end

	 headers = Dict(
        "User-Agent" => "Requests.jl/0.0.0",
        "Accept" => "text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8"
    )
    for (k,v) in amz_headers
        setindex!(headers, v, k)
    end

    #make the request
    sqsresp = SQSResponse()
    if !(env.dry_run)
        ## ro = HTTPC.RequestOptions(headers = amz_headers, request_timeout = env.timeout)
		## TODO timout value to be passed to Requests !
        if use_post
            ## resp = HTTPC.post(complete_url, data, ro)
			resp = Requests.post(URIParser.parse_url(complete_url), data; headers=headers)
        else
            ## resp = HTTPC.get(complete_url, ro)
			resp = Requests.get(complete_url; headers = headers)
        end

        sqsresp.http_code = resp.status
        sqsresp.headers = resp.headers
        sqsresp.body = bytestring(resp.data)

        if (env.dbg)
            print("HTTPCode: ", sqsresp.http_code, "\nHeaders: ", sqsresp.headers, "\nBody : ", sqsresp.body, "\n")
        end

        if (sqsresp.http_code >= 200) && (sqsresp.http_code <= 299)
             if (search(Base.get(resp.headers, "Content-Type", [""]), "/xml") != 0:-1)
#            if  haskey(resp.headers, "Content-Type") && (resp.headers["Content-Type"] == "application/xml")
                sqsresp.pd = LightXML.root(LightXML.parse_string(bytestring(sqsresp.body)))
            end
        elseif (sqsresp.http_code >= 400) && (sqsresp.http_code <= 599)
            if length(sqsresp.body) > 0
                xom = LightXML.root(LightXML.parse_string(bytestring(sqsresp.body)))
                epd = LightXML.find_element(xom, "Error")
                if epd == nothing
                    error("HTTP error : $(resp.status), $(sqsresp.body)")
                end
                sqsresp.obj = SQSError(LightXML.content(LightXML.find_element(epd, "Type")), LightXML.content(LightXML.find_element(epd, "Code")), LightXML.content(LightXML.find_element(epd, "Message")), LightXML.content(LightXML.find_element(epd, "Detail")), LightXML.content(LightXML.find_element(xom, "RequestId")))
            else
                error("HTTP error : $(resp.http_code)")
            end
        else
            error("HTTP error : $(resp.http_code), $(sqsresp.body)")
        end
    end

    sqsresp
end

include("sqs_typed.jl")

end


text:
module StatusPage

using InfluxDB, TOML, DataFrames, Dates, OrderedCollections, Artifacts

function run(config="config.toml")
    config = abspath(config)
    # Read config
    toml = TOML.parsefile(config)
    # TODO: Defaults
    server = InfluxServer(toml["influxdb"]["server"])
    measurement = toml["influxdb"]["measurement"]
    db = toml["influxdb"]["db"]
    output_directory = abspath(dirname(config), toml["output_directory"])
    mkpath(output_directory)
    interval = toml["interval"]
    ## Build the output structure
    group_service_map = OrderedDict()
    for service in toml["services"]
        vect = get!(Vector, group_service_map, get(service, "group", nothing))
        push!(vect, service)
    end
    ## Sort by group array
    if (groups = get(toml, "groups", nothing); groups !== nothing)
        # TODO: ...
    end
    sort!(group_service_map; lt = (x, y) -> x === nothing ? false : y === nothing ? true : x < y)
    ## All urls we are interested in
    all_urls = Set{String}(url for s in toml["services"] for url in s["checks"])

    while true
        t0 = time()
        # Make the InfluxDB query
        q = SELECT(measurements=[measurement], condition="time > now() - 91d")
        dfs = InfluxDB.query(server, db, q)
        if dfs === nothing
            df = DataFrame(time=String[], result=String[], server=String[])
        else
            df = dfs[1]
        end
        ## Keep only what we need
        select!(df, [:time, :result, :server])
        ## Convert to DateTime Date
        conv_dt(X) = map(x -> DateTime(x, dateformat"yyyy-mm-dd\THH:MM:SS\Z"), X)
        transform!(df, :time => conv_dt => :time)
        transform!(df, :time => (X -> Date.(X)) => :date)
        # Filter out from the config
        filter!(:server => (x -> x in all_urls), df)
        ## Add service and group from the config
        function assign_name_group(X)
            n = String[]
            g = Union{String,Nothing}[]
            for x in X
                for (group, services) in group_service_map, service in services
                    if x in service["checks"]
                        push!(n, service["name"])
                        push!(g, group)
                        break
                    end
                end
            end
            return (name=n, group=g) # ??
        end
        transform!(df, :server => assign_name_group => [:config_name, :config_group])
        # Convert result to binary
        transform!(df, :result => (X -> map(x -> x == "success", X)) => :result)

        # Generate html
        index = joinpath(output_directory, "index.html")
        index_tmp = index * ".tmp"
        try
            open(index_tmp, "w") do io
                generate_html(io, df, group_service_map, toml)
            end
            mv(index_tmp, index; force=true)
        finally
            rm(index_tmp; force=true)
        end
        ## Copy asset files
        copy_assets(output_directory)

        sleeptime = interval - (time() - t0)
        @info "Generated status page, regenerating in $(round(Int, sleeptime)) seconds."
        if sleeptime > 0
            sleep(sleeptime)
        end
    end
end

function copy_safe(src, dst)
    if filesize(src) == filesize(dst)
        return # Assume file has not changed...
    end
    tmp = dst * ".tmp"
    try
        open(tmp, "w") do io; open(src, "r") do io2
            write(io, io2)
        end end
        mv(tmp, dst; force=true)
    catch err
        @error "Something went wrong when copying asset..." src dst exception=(err, catch_backtrace())
    finally
        rm(tmp; force=true)
    end
end

function copy_assets(out)
    # CSS file
    src_files = joinpath.(@__DIR__, ["status-page.css"])
    # JuliaMono webfonts
    jlmono = joinpath.(artifact"JuliaMono", "juliamono-0.030", "webfonts", [
            "JuliaMono-RegularLatin.woff2",
            "JuliaMono-BoldLatin.woff2",
            "JuliaMono-Regular.woff2",
            "JuliaMono-Bold.woff2"
        ])
    append!(src_files, jlmono)
    dst = joinpath(out, "assets")
    mkpath(dst)
    for src_file in src_files
        dst_file = joinpath(dst, basename(src_file))
        copy_safe(src_file, dst_file)
    end
end

function generate_html(io, df, group_service_map, toml)
    title = get(toml, "title", "Status page")
    desc = get(toml, "desc", "")
    print(io, """
        <!doctype html>

        <meta charset="UTF-8">
        <meta name="viewport" content="width=device-width, initial-scale=1">
        <meta name="description" content="$(desc)">
        <title>$(title)</title>
        <link rel="stylesheet" href="/assets/status-page.css">

$(title)

$(desc)

        """)
    # Groups
    dfg = groupby(df, [:config_group, :config_name])
    for (group_name, services) in group_service_map
        # df_group = df_groups[(cg,)]
        # gname = group_name === nothing ? "" : group_name
        print(io, """

$(something(group_name, ""))

            """)
        # config_services = [x["name"] for x in toml["services"] if get(x, "group", nothing) == cg]
        # df_services = groupby(df_group, :config_name)
        for service in services
            df_service = get(dfg, (group_name, service["name"]), DataFrame(date=Date[], result=Bool[], server=String[]))
            print_service(io, df_service, service)
        end
        println(io, """
""")
    end
    last_update = Dates.format(unix2datetime(time()), dateformat"yyyy-mm-dd\THH:MM:SS\Z")
    print(io, """

Status page built with StatusPage.jl. Last update: $(last_update).

        """)
end

const OPERATIONAL_SYMBOL = """<svg version="1.1" viewBox="0 0 512 512"><path pid="0" _fill="currentColor" d="M504 256c0 136.967-111.033 248-248 248S8 392.967 8 256 119.033 8 256 8s248 111.033 248 248zM227.314 387.314l184-184c6.248-6.248 6.248-16.379 0-22.627l-22.627-22.627c-6.248-6.249-16.379-6.249-22.628 0L216 308.118l-70.059-70.059c-6.248-6.248-16.379-6.248-22.628 0l-22.627 22.627c-6.248 6.248-6.248 16.379 0 22.627l104 104c6.249 6.249 16.379 6.249 22.628.001z"></path></svg>"""
const PARTIAL_SYMBOL = """<svg version="1.1" viewBox="0 0 512 512"><path pid="0" _fill="currentColor" d="M256 8C119 8 8 119 8 256s111 248 248 248 248-111 248-248S393 8 256 8zM124 296c-6.6 0-12-5.4-12-12v-56c0-6.6 5.4-12 12-12h264c6.6 0 12 5.4 12 12v56c0 6.6-5.4 12-12 12H124z"></path></svg>"""
const MAJOR_SYMBOL = """<svg version="1.1" viewBox="0 0 512 512"><path pid="0" _fill="currentColor" d="M256 8C119 8 8 119 8 256s111 248 248 248 248-111 248-248S393 8 256 8zm121.6 313.1c4.7 4.7 4.7 12.3 0 17L338 377.6c-4.7 4.7-12.3 4.7-17 0L256 312l-65.1 65.6c-4.7 4.7-12.3 4.7-17 0L134.4 338c-4.7-4.7-4.7-12.3 0-17l65.6-65-65.6-65.1c-4.7-4.7-4.7-12.3 0-17l39.6-39.6c4.7-4.7 12.3-4.7 17 0l65 65.7 65.1-65.6c4.7-4.7 12.3-4.7 17 0l39.6 39.6c4.7 4.7 4.7 12.3 0 17L312 256l65.6 65.1z"></path></svg>"""

function print_service(io::IO, df_service, service)
    service_name = service["name"]
    ## Compute average for last 90 days
    day_avgs = combine(groupby(df_service, :date), :result => (X -> sum(X)/length(X)) => :daily_avg)
    dict = Dict(day_avgs.date .=> day_avgs.daily_avg)
    today = Date(unix2datetime(time()))
    avgs = [get(dict, day, nothing) for day in (today - Day(89)):Day(1):today]
    history = mapreduce(*, avgs) do uptime
        st = uptime === nothing ? "unknown" :
             uptime > 0.9965 ? "operational" :   # <  5 minutes/day
             uptime > 0.993 ? "partial-outage" : # < 10 minutes/day
                            "major-outage"
        return """

\n"""
    end
    ## Compute current
    gdf = groupby(df_service, :server)
    if length(gdf) < length(service["checks"]) # no measurement for some checks
        status_text = "Unknown"
        status_css = "unknown"
        status_sym = PARTIAL_SYMBOL
    else
        g_last = Bool[last(g).result for g in gdf]
        current = sum(g_last) / length(service["checks"])
        if current > 0.99
            status_text = "Operational"
            status_css = "operational"
            status_sym = OPERATIONAL_SYMBOL
        elseif current < 0.1
            status_text = "Major outage"
            status_css = "major-outage"
            status_sym = MAJOR_SYMBOL
        else
            status_text = "Partial outage"
            status_css = "partial-outage"
            status_sym = PARTIAL_SYMBOL
        end
    end
    ## Print it out
    print(io, """

        <input id="$(service_name)" type="checkbox" checked="false" style="display: none;">

        $(service_name)

        $(status_text) $(status_sym)

        $(history)

        90 days ago
        Today

        """)
end

end # module


text:
### Utility Functions needed for the construction of the Power System, mostly used for consistency checking ####

## Check that all the buses have a type defintion ##

function buscheck(buses)
    for b in buses
        if isnothing(b.bustype)
            @warn "Bus/Nodes data does not contain information to build an a network"
        end
    end
end

## Slack Bus Definition ##

function slack_bus_check(buses)
    slack = -9
    for b in buses
        if b.bustype == BusTypes.REF
            slack = b.number
            break
        end
    end
    if slack == -9
        @error "Model doesn't contain a slack bus"
    end
end

# TODO: Check for islanded Buses

# check for minimum timediff
function minimumtimestep(time_series::Array{T}) where {T <: TimeSeriesData}
    if length(time_series[1].data) > 1
        timeseries = time_series[1].data
        n = length(timeseries) - 1
        ts = []
        for i in 1:n
            push!(
                ts,
                TimeSeries.timestamp(timeseries)[n + 1] -
                TimeSeries.timestamp(timeseries)[n],
            )
        end
        return minimum(ts)
    else
        ts = Dates.Dates.Minute(1)
        return ts
    end
end

function critical_components_check(sys::System)
    critical_component_types = [Bus, Generator, ElectricLoad]
    for component_type in critical_component_types
        components = get_components(component_type, sys)
        if length(components) == 0
            @warn "There are no $(component_type) Components in the System"
        end
    end
end

"""
    adequacy_check(sys::System)

Checks the system for sum(generator ratings) >= sum(load ratings).

# Arguments
- `sys::System`: system
"""
function adequacy_check(sys::System)
    gen = total_capacity_rating(sys)
    load = total_load_rating(sys)
    load > gen && @warn "System peak load ($load) exceeds total capacity capability ($gen)."
end

"""
    total_load_rating(sys::System)

Checks the system for sum(generator ratings) >= sum(load ratings).

# Arguments
- `sys::System`: system
"""
function total_load_rating(sys::System)
    base_power = get_base_power(sys)
    controllable_loads = get_components(ControllableLoad, sys)
    cl =
        isempty(controllable_loads) ? 0.0 :
        sum(get_max_active_power.(controllable_loads)) * base_power
    @debug "System has $cl MW of ControllableLoad" _group = IS.LOG_GROUP_SYSTEM_CHECKS
    static_loads = get_components(StaticLoad, sys)
    sl = isempty(static_loads) ? 0.0 : sum(get_max_active_power.(static_loads)) * base_power
    @debug "System has $sl MW of StaticLoad" _group = IS.LOG_GROUP_SYSTEM_CHECKS
    # Total load calculation assumes  P = Real(V^2/Y) assuming V=1.0
    fa_loads = get_components(FixedAdmittance, sys)
    fa =
        isempty(fa_loads) ? 0.0 :
        sum(real.(get_base_voltage.(get_bus.(fa_loads)) .^ 2 ./ get_Y.(fa_loads)))
    @debug "System has $fa MW of FixedAdmittance assuming admittance values are in P.U." _group =
        IS.LOG_GROUP_SYSTEM_CHECKS
    total_load = cl + sl + fa
    @debug "Total System Load: $total_load" _group = IS.LOG_GROUP_SYSTEM_CHECKS
    return total_load
end

"""
    total_capacity_rating(sys::System)

Sum of system generator and storage ratings.

# Arguments
- `sys::System`: system
"""
function total_capacity_rating(sys::System)
    total = 0
    for component_type in (Generator, Storage)
        components = get_components(Generator, sys)
        if !isempty(components)
            component_total = sum(get_rating.(components)) * get_base_power(sys)
            @debug "total rating for $component_type = $component_total" _group =
                IS.LOG_GROUP_SYSTEM_CHECKS
            total += component_total
        end
    end

    @debug "Total System capacity: $total" _group = IS.LOG_GROUP_SYSTEM_CHECKS
    return total
end


text:
using Correlate
using Test
using Statistics

# We have run correlate and want to see if the result is within tolerance.

function check_correlation(zdata::Array{Float64,2}, corm::Array{Float64,2}, tol::Float64)
  (nrows,ncols) = size(zdata)

  cor_after = cor(zdata)

  for c in 1:ncols
    for r in (c+1):ncols
      if abs(corm[r,c] - cor_after[r,c]) > tol
        return false
      end
    end
  end

  return true
end

function test1()
  data = rand(1000,3)
  corm = [1 0.2 -0.3; 0.2 1 0.85; -0.3 .85 1]
  tol = 0.01
  rc = correlate!(data, corm, tol)
  if ! rc
    return false
  end

  if ! check_correlation(data, corm, tol)
    return false
  end

  (nrows,ncols) = size(data)

# now reverse it

  for c in 1:ncols
    for r in 1:ncols
      if r != c
        corm[r,c] = - corm[r,c]
      end
    end
  end

# write(STDERR, "CORM $(corm)\n")

  if ! correlate!(data, corm, tol)
    return false
  end

  if ! check_correlation(data, corm, tol)
    return false
  end

# now send back to uncorrelated

  for c in 1:ncols
    for r in 1:ncols
      if r != c
        corm[r,c] = 0.0
      end
    end
  end

  if ! correlate!(data, corm, tol)
    return false
  end

  return check_correlation(data, corm, tol)
end

# write your own tests here
@test test1()

max_stars_repo_name: IanAWatson/Correlate

max_issues_repo_name: IanAWatson/Correlate

max_forks_repo_name: IanAWatson/Correlate


text:
"""
Showing how to use Term.jl to create styled console representation for your types.
"""

using Term: Panel, vLine, rvstack, lvstack, RenderableText, @with_repr, termshow

struct myType
    name::String
    height::Int
    width::Int
    mass::Float64
end

""" Custom show method """
function Base.show(io::IO, ::MIME"text/plain", obj::myType)
    fields = (:height, :width, :mass)

    info = map(f -> RenderableText(string(f); style = "bold"), fields)
    vals = map(
        f -> RenderableText(" " * string(getfield(obj, f)); style = "bright_blue"),
        fields,
    )
    obj_details = rvstack(info...) * vLine(3; style = "dim") * lvstack(vals...)

    return print(
        io,
        Panel(
            obj_details;
            title = obj.name,
            style = "red dim",
            title_style = "default bright_red bold",
            fit = true,
            padding = (2, 2, 1, 1),
        ),
    )
end

obj = myType("Rocket", 10, 10, 99.9)

# or just use our macro!

@with_repr mutable struct Rocket
    width::Int
    height::Int
    mass::Float64
    engines_activity::Vector

    manufacturer::String
end

obj = Rocket(10, 50, 5000, [1.1, 123123, 1], "NASA")


text:
using Test
using DECAES
using DECAES.NormalHermiteSplines
using StaticArrays
using Random
using LinearAlgebra

function test_mock_surrogate_search_problem(
        opts::T2mapOptions = DECAES.mock_t2map_opts(;
            MatrixSize = (1, 1, 1),
            nRefAngles = 8,
        )
    )
    function A(α, β)
        theta = DECAES.EPGOptions((; α = α, TE = opts.TE, T2 = 0.0, T1 = opts.T1, β = β), Val(32), Float64)
        T2_times = DECAES.logrange(opts.T2Range..., opts.nT2)
        DECAES.epg_decay_basis(theta, T2_times)
    end

    function f!(work, prob, α, β)
        DECAES.solve!(work, A(α, β), prob.b)
        return log(max(DECAES.chi2(work), eps()))
    end

    function fg_approx!(work, prob, α, β; h)
        l   = f!(work, prob, α, β)
        lα⁺ = f!(work, prob, α + h, β)
        lα⁻ = f!(work, prob, α - h, β)
        lβ⁺ = f!(work, prob, α, β + h)
        lβ⁻ = f!(work, prob, α, β - h)
        ∂l_∂α = (lα⁺ - lα⁻) / 2h
        ∂l_∂β = (lβ⁺ - lβ⁻) / 2h
        ∇l  = SA[∂l_∂α, ∂l_∂β]
        return l, ∇l
    end

    function fg_surrogate!(prob, I)
        l  = DECAES.loss!(prob, I)
        ∇l = DECAES.∇loss!(prob, I)
        return l, ∇l
    end

    prob = DECAES.mock_surrogate_search_problem(Val(2), Val(opts.nTE), opts)
    work = DECAES.lsqnonneg_work(zeros(opts.nTE, opts.nT2), zeros(opts.nTE))

    for I in CartesianIndices(prob.αs)
        α, β = prob.αs[I]
        l′, ∇l′ = fg_approx!(work, prob, α, β; h = 1e-6)
        l , ∇l  = fg_surrogate!(prob, I)
        @test l ≈ l′
        @test ∇l ≈ ∇l′ rtol = 1e-6 atol = 1e-8
    end
end

function test_bounding_box()
    for bounds in [
            ((2,3),),
            ((2,3), (5,7)),
            ((2,3), (5,7), (-3,-1)),
        ]
        box = DECAES.BoundingBox(bounds)
        for I in box.corners
            Iopp = CartesianIndex(ntuple(d -> ifelse(I[d] == bounds[d][1], bounds[d][2], bounds[d][1]), length(bounds)))
            @test DECAES.opposite_corner(box, I) == Iopp
        end
    end
end

function test_discrete_searcher()
    grid  = DECAES.meshgrid(SVector{2,Float64}, 1:5, 1:10)
    surr  = DECAES.HermiteSplineSurrogate(I -> (1.0, zero(SVector{2,Float64})), grid)
    state = DECAES.DiscreteSurrogateSearcher(surr; mineval = 0, maxeval = 0)
    state.seen[[1,5], [1,10]] .= true

    @test DECAES.sorted_corners(state, DECAES.BoundingBox(((2,4), (3,7))), SA[3.5, 6.5]) == SMatrix{2,2}((CartesianIndex(4, 7), CartesianIndex(2, 7), CartesianIndex(4, 3), CartesianIndex(2, 3)))
    @test DECAES.sorted_corners(state, DECAES.BoundingBox(((2,4), (3,7))), SA[2.5, 6.5]) == SMatrix{2,2}((CartesianIndex(2, 7), CartesianIndex(4, 7), CartesianIndex(2, 3), CartesianIndex(4, 3)))
    @test DECAES.sorted_corners(state, DECAES.BoundingBox(((2,4), (3,7))), SA[2.5, 4.5]) == SMatrix{2,2}((CartesianIndex(2, 3), CartesianIndex(4, 3), CartesianIndex(2, 7), CartesianIndex(4, 7)))
    @test DECAES.sorted_corners(state, DECAES.BoundingBox(((2,4), (3,7))), SA[3.5, 4.5]) == SMatrix{2,2}((CartesianIndex(4, 3), CartesianIndex(2, 3), CartesianIndex(4, 7), CartesianIndex(2, 7)))

    #  1  -------  1  ----------  1
    #  |  0  0  0  |  0  0  0  0  |
    #  |  0  0  0  |  0  0  0  0  |
    #  |  0  0  0  |  0  0  0  0  |
    #  1  -------  1  ----------  1
    @test DECAES.minimal_bounding_box(state, SA[1.5, 2.0]) == DECAES.BoundingBox(((1,5), (1,5)))
    @test DECAES.minimal_bounding_box(state, SA[4.5, 6.5]) == DECAES.BoundingBox(((1,5), (5,10)))
    state.seen[[1,5], 5] .= true

    #  1  -------  1  -- 1 -----  1
    #  |  0  0  0  |  0  |  0  0  |
    #  |  -------  |  0  |  0  0  |
    #  |  0  0  0  |  0  |  0  0  |
    #  1  -------  1  -- 1 -----  1
    @test DECAES.minimal_bounding_box(state, SA[1.5, 2.0]) == DECAES.BoundingBox(((1,3), (1,5)))
    @test DECAES.minimal_bounding_box(state, SA[4.5, 2.0]) == DECAES.BoundingBox(((3,5), (1,5)))
    @test DECAES.minimal_bounding_box(state, SA[1.5, 6.5]) == DECAES.BoundingBox(((1,5), (5,7)))
    @test DECAES.minimal_bounding_box(state, SA[4.5, 6.5]) == DECAES.BoundingBox(((1,5), (5,7)))
    @test DECAES.minimal_bounding_box(state, SA[1.5, 8.5]) == DECAES.BoundingBox(((1,5), (7,10)))
    @test DECAES.minimal_bounding_box(state, SA[4.5, 8.5]) == DECAES.BoundingBox(((1,5), (7,10)))
end

@testset "Splines" begin
    @testset "mock surrogate search problem" begin
        test_mock_surrogate_search_problem()
    end
    @testset "bounding box" begin
        test_bounding_box()
        test_discrete_searcher()
    end
end


text:
using ClimateTools
using AxisArrays
# const axes = Base.axes
using PyPlot
using Shapefile
using NetCDF
using NCDatasets
using Dates
using Statistics
using Random
using ArgCheck
# using Unitful: K, °C, m, mm, s, kg, Pa
# using Unitful: @u_str, ustrip, uconvert
# import Unitful
using Test

include("interface_test.jl")
include("utils_test.jl")
include("extract_test.jl")
include("biascorrect_test.jl")
include("indices_test.jl")
include("mapping_test.jl")
include("functions_test.jl")
include("spatial_test.jl")


text:
using Test
using OpenQuantumSystems
using Random, SparseArrays, LinearAlgebra, StableRNGs
import QuantumOpticsBase

import DelayDiffEq

@testset "master" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    mode1 = Mode(0.2, 1.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 2, Energy)
    mol2 = Molecule([mode1], 2, Energy)
    agg = Aggregate([mol1, mol2])
    aggInds = getIndices(agg)
    aggIndsLen = length(aggInds)
    basis = GenericBasis([aggIndsLen])
    FCFact = getFranckCondonFactors(agg, aggInds)
    Ham = getAggHamiltonian(agg, aggInds, FCFact)

    Ham_I = getAggHamInteraction(agg, aggInds, FCFact)
    Ham_0 = getAggHamSystemBath(agg, aggInds, FCFact)
    Ham = getAggHamiltonian(agg, aggInds, FCFact)

    data = Matrix(Hermitian(rand(ComplexF64, aggIndsLen, aggIndsLen)))
    rho0 = DenseOperator(basis, basis, data)
    normalize!(rho0)
    # tests have to be quick enough
    tspan = [0.0:0.01:0.1;]

    # TODO: chenge back reltol, abstol after solving QuadGK
    T, rho_t = master_int(
        rho0,
        tspan,
        Ham_0,
        Ham_I;
        reltol = 1e-6,
        abstol = 1e-6,
        int_reltol = 1e-8,
        int_abstol = 0.0,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()),
    )#, alg=OrdinaryDiffEq.Vern7())
    rho_prev = deepcopy(rho0)
    for t_i = 2:length(rho_t)
        t = T[t_i]
        rho_I = rho_t[t_i]
        U_op_S = evolutionOperator(Ham_0, t)
        rho = U_op_S * rho_I * U_op_S'
        U_op = evolutionOperator(Ham, t)
        rho_ref = U_op * rho0 * U_op'
        @test 1e-4 > D(rho_ref, rho)
    end

    T, rho_t = master(
        rho0,
        tspan,
        Ham;
        reltol = 1e-6,
        abstol = 1e-6,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()),
    )#, alg=OrdinaryDiffEq.Vern7())
    rho_prev = deepcopy(rho0)
    for t_i = 2:length(rho_t)
        t = T[t_i]
        rho = rho_t[t_i]
        U_op = evolutionOperator(Ham, t)
        rho_ref = U_op * rho0 * U_op'
        @test 1e-4 > D(rho_ref, rho)
    end

end # testset


text:
include("../src/julia/AlgamesRoundabout.jl")
using .AlgamesRoundabout
using Plots

# game settings

uc = 1
xpoly = [[ 9.92483328e+02,  1.00992169e+00, -1.06373371e-02,
         5.44482755e-03, -1.56124376e-03,  2.96003128e-04,
        -3.97313295e-05,  3.90312145e-06, -2.85274315e-07,
         1.56024980e-08, -6.35690001e-10,  1.89076217e-11,
        -3.87873308e-13,  4.48568557e-15,  1.04254642e-17,
        -1.56423667e-18,  3.29567455e-20, -3.92492131e-22,
         2.88391940e-24, -1.22655411e-26,  2.32697067e-29],
       [ 1.06021128e+03, -3.94695552e-01, -2.65480724e-01,
         7.72331895e-02, -1.28764151e-02,  1.35288108e-03,
        -9.63032818e-05,  4.88162591e-06, -1.81700838e-07,
         5.04535619e-09, -1.04705323e-10,  1.59722674e-12,
        -1.69571241e-14,  1.03019971e-16,  9.91898387e-20,
        -9.43413842e-21,  1.05428864e-22, -6.63261341e-25,
         2.57206072e-27, -5.76876538e-30,  5.76518416e-33],
       [ 1.03956886e+03, -6.66404535e-01,  1.68722243e-01,
        -9.61339363e-02,  2.85736090e-02, -4.98160515e-03,
         5.74110190e-04, -4.64814345e-05,  2.72600951e-06,
        -1.17131710e-07,  3.66123146e-09, -8.02534848e-11,
         1.09485575e-12, -4.37048296e-15, -1.60842484e-16,
         3.75509773e-18, -3.66648589e-20,  1.13749036e-22,
         1.09091124e-24, -1.18463934e-26,  3.49468666e-29],
       [ 1.04368180e+03, -1.00002693e+00,  1.62315284e-02,
        -3.67536866e-03,  5.53573396e-04, -6.29877665e-05,
         5.78415323e-06, -4.25140558e-07,  2.40307950e-08,
        -1.01864918e-09,  3.20099297e-11, -7.33625975e-13,
         1.17000629e-14, -1.09904618e-16,  1.15040614e-20,
         1.78120874e-20, -3.04212971e-22,  2.83825404e-24,
        -1.61941073e-26,  5.32615567e-29, -7.79110941e-32]]

ypoly = [[ 9.85056258e+02, -7.33797447e-03, -1.60582416e-02,
         7.69597502e-03, -2.24558622e-03,  4.25337577e-04,
        -5.54541321e-05,  5.15766542e-06, -3.50257135e-07,
         1.75989111e-08, -6.54629343e-10,  1.76918782e-11,
        -3.26218029e-13,  3.19403241e-15,  1.82888213e-17,
        -1.28191809e-18,  2.40451512e-20, -2.64439070e-22,
         1.81562474e-24, -7.25616139e-27,  1.29732722e-29],
       [ 9.72977375e+02,  1.68013496e-01,  2.62859673e-01,
        -8.19509164e-02,  1.45446629e-02, -1.62040306e-03,
         1.21181800e-04, -6.37599460e-06,  2.43629872e-07,
        -6.88482976e-09,  1.44532818e-10, -2.22237103e-12,
         2.37577213e-14, -1.45865218e-16, -1.22469420e-19,
         1.32130012e-20, -1.49434739e-22,  9.49787654e-25,
        -3.72296646e-27,  8.44740565e-30, -8.54771388e-33],
       [ 9.71432643e+02,  1.09644346e+00, -3.16524750e-01,
         2.00055092e-01, -6.62052012e-02,  1.33578100e-02,
        -1.79359783e-03,  1.69149381e-04, -1.15739527e-05,
         5.84432938e-07, -2.18342229e-08,  5.93275610e-10,
        -1.10412944e-11,  1.11264523e-13,  5.14067848e-16,
        -4.13379623e-17,  7.86212065e-19, -8.68704866e-21,
         5.97414273e-23, -2.38779614e-25,  4.26590812e-28],
       [ 1.00765571e+03,  1.40122934e-01,  4.66769186e-03,
        -9.84994942e-04,  3.54607296e-04, -5.61687052e-05,
         4.94425778e-06, -2.64201889e-07,  7.95049087e-09,
        -6.11160141e-11, -5.32885165e-12,  2.53520405e-13,
        -5.65667160e-15,  6.63090268e-17, -1.23853952e-19,
        -9.29126946e-21,  1.68471272e-22, -1.54367074e-24,
         8.34787712e-27, -2.53711116e-29,  3.35685250e-32]]

smax = [ 59.51321215, 113.35577177,  58.70627679,  80.66267064]

ssdot = [[44.15017431,  9.94493352],
       [27.12617116,  7.74348089],
       [11.59016022,  6.5109571 ],
       [18.43266001,  8.3942871 ]]

v0 = 11.17
prob = setup_problem(ssdot, xpoly, ypoly, smax, uc;
N=8, dt=0.2, radius=3., vel_control=true, vel_des = v0)
newton_solve!(prob)

x, y, s, v = trajectory(prob.model, prob.pdtraj.pr)

plt = plot()
plot_trajectory!(plt, prob.model, prob.pdtraj.pr)
xlims!((1000, 1060))
ylims!((970, 1020))
savefig(plt, "./experiments/results/trajectory.png")

plt2 = plot(legend=:bottomright)
for i = 1:prob.model.p
    plot!(range(0.0;length=length(v[i]),step=0.2), v[i],
        label=false)
    scatter!(range(0.0;length=length(v[i]),step=0.2), v[i],
        label="Vehicle $(i)")
end
xlabel!("Time (s)")
ylabel!("Velocity (m/s)")
xlims!((-.1, 1.6))
ylims!((5, 13))

savefig(plt2, "./experiments/results/velocities.png")

max_stars_repo_name: lassepe/DecNashPlanning

max_issues_repo_name: lassepe/DecNashPlanning

max_forks_repo_name: lassepe/DecNashPlanning


text:
CBB
hasprecon(::CBB) = NoPrecon()


text:
module Histogram
using StatsBase

export fit

function fit(values, nbins=15)
  fitted = StatsBase.fit(StatsBase.Histogram, values, nbins=nbins, closed=:left)
  vcat(hcat(fitted.edges[1][1:end-1]...), hcat(fitted.weights...))
end

end


text:
export GeneralizedEuclideanAlgorithm

"""
    abstract type AbstractUnivariateGCDAlgorithm end

Algorithm computing the greatest common divisor of univariate polynomials.
See [`GeneralizedEuclideanAlgorithm`](@ref) and [`SubresultantAlgorithm`](@ref).
"""
abstract type AbstractUnivariateGCDAlgorithm end

"""
    struct GeneralizedEuclideanAlgorithm <: AbstractUnivariateGCDAlgorithm
        primitive_rem::Bool
        skip_last::Bool
    end

Algorithm computing the greatest common divisor of univariate polynomials using
the Euclidean algorithm generalized for polynomials with coefficients over a
a unique factorization domain, see [Knu14, Algorithm E, p. 426-427].

[Knu14] Knuth, D.E., 2014.
*Art of computer programming, volume 2: Seminumerical algorithms.*
Addison-Wesley Professional. Third edition.
"""
struct GeneralizedEuclideanAlgorithm <: AbstractUnivariateGCDAlgorithm
    primitive_rem::Bool
    skip_last::Bool
    function GeneralizedEuclideanAlgorithm(
        primitive_rem::Bool=false,
        skip_last::Bool=false,
    )
        return new(primitive_rem, skip_last)
    end
end

"""
    struct SubresultantAlgorithm <: AbstractUnivariateGCDAlgorithm end

Algorithm computing the greatest common divisor of univariate polynomials using
the Subresultant algorithm, see [Knu14, Algorithm C, p. 428-429].

[Knu14] Knuth, D.E., 2014.
*Art of computer programming, volume 2: Seminumerical algorithms.*
Addison-Wesley Professional. Third edition.
"""
struct SubresultantAlgorithm <: AbstractUnivariateGCDAlgorithm end

Base.lcm(p::APL, q::APL, algo::AbstractUnivariateGCDAlgorithm=GeneralizedEuclideanAlgorithm()) = p * div(q, gcd(p, q, algo))
Base.gcd(α, p::APL, algo::AbstractUnivariateGCDAlgorithm=GeneralizedEuclideanAlgorithm()) = gcd(α, content(p, algo))
Base.gcd(p::APL, α, algo::AbstractUnivariateGCDAlgorithm=GeneralizedEuclideanAlgorithm()) = gcd(content(p, algo), α)

#function MA.promote_operation(::typeof(gcd), P::Type{<:Number}, Q::Type{<:Number})
#    return typeof(gcd(one(P), one(Q)))
#end
function MA.promote_operation(::typeof(gcd), P::Type{<:APL}, Q::Type{<:APL})
    return MA.promote_operation(pseudo_rem, P, Q)
end

"""
    function gcd(p1::AbstractPolynomialLike{T}, p2::AbstractPolynomialLike{S}) where {T, S}

Returns a greatest common divisor of `p1` and `p2`. Note that it does not make
sense, in general, to speak of "the" greatest common divisor of u and v; there
is a set of greatest common divisors, each one being a unit multiple of the
others [Knu14, p. 424].

# Implementation notes

The classical algorithm for computing the `gcd`, commonly referred to as the
Euclidean Algorithm is to use a recursion with the base case `gcd(p, 0) = p`
and the relation `gcd(p1, p2) = gcd(p2, rem(p1, p2))`.
The relation comes from the Euclidean division:
`p1 = q * p2 + r`,
if `g` divides `p1` and `p2` then it divides `r` and
if `g` divides `r` and `p2` then it divides `p1`.

For multivariate polynomials, you may have `rem(p1, p2) = p1` hence this will not
terminate.
To ensure we make progress, we can pick a given variable `xi` and try to find
`q1` and `q2` such that
`q2 * p1 = q1 * p2 + r`
and the degree of `r` in `xi` is strictly smaller than the degree of `p1` in `xi`.
Note that if `g` divides `p1` and `p2` then it divides `r` but
if `g` divides `r` and `p2` then it might divide `q2` and not `p1`.
So what do we do ?
Let `dj` be the degree of `pj` in `xi`.
Suppose we pick `qj` to be the coefficient of `pj` in `xi^dj`.
If `g` divides `q2` then it means that the degree of `g` in `xi` is zero.
Therefore, if it divides `p2` then it also divides the coefficients
of `p2` in `xi^k` for `k = 0, 1, ..., d2`.
This means that if we ensure that these are relatively prime then we won't have
any issue.
So we start by computing a `gcd` `gj` of the coefficients in each degree of
`xi` of `pj`, this is called the [`content`](@ref) of `pj`.
And then we compute `_gcd(p1 / g1, p2 / g2) * gcd(g1, g2)` where we can use the
recursion `_gcd(p1, p2) = _gcd(p2, q2 * p1 - q1 * p2)` where `q1, q2` are as
defined above.
This is the [`GeneralizedEuclideanAlgorithm`](@ref).

[Knu14] Knuth, D.E., 2014.
*Art of computer programming, volume 2: Seminumerical algorithms.*
Addison-Wesley Professional. Third edition.
"""
function Base.gcd(p1::APL{T}, p2::APL{S}, algo::AbstractUnivariateGCDAlgorithm=GeneralizedEuclideanAlgorithm()) where {T, S}
    # If one of these is zero, `shift` should be infinite
    # for this method to work so we exclude these cases.
    if isapproxzero(p1)
        return convert(MA.promote_operation(gcd, typeof(p1), typeof(p2)), p2)
    end
    if isapproxzero(p2)
        return convert(MA.promote_operation(gcd, typeof(p1), typeof(p2)), p1)
    end
    shift1, defl1 = deflation(p1)
    shift2, defl2 = deflation(p2)
    shift = gcd(shift1, shift2)
    defl = mapexponents(gcd, defl1, defl2)
    # We factor out `x.^shift1` from `p1` and
    # `x.^shift2` from `p2`. The `gcd` of these
    # monomials is `x.^shift`.
    # Then, we subsitute `y[i] = x[i]^defl[i]`.
    q1 = deflate(p1, shift1, defl)
    q2 = deflate(p2, shift2, defl)
    g = deflated_gcd(q1, q2, algo)
    return inflate(g, shift, defl)::MA.promote_operation(gcd, typeof(p1), typeof(p2))
end

# Inspired from to `AbstractAlgebra.deflation`
function deflation(p::AbstractPolynomialLike)
    if iszero(p)
        return constantmonomial(p), constantmonomial(p)
    end
    shift = fill(-1, nvariables(p))
    defl = zeros(Int, nvariables(p))
    for mono in monomials(p)
        exps = exponents(mono)
        for i in eachindex(exps)
            exp = exps[i]
            @assert exp >= 0
            if shift[i] == -1
                shift[i] = exp
            elseif exp < shift[i]
                # There are two cases:
                # 1) If `defl[i]` is zero then it means all previous monomials
                #    had degree `shift[i]` so we just set `defl[i]` to
                #    `shift[i] - exp` or equivalently `gcd(0, shift[i] - exp)`.
                # 2) If  `defl[i]` is positive then we have some monomials with
                #    degree `shift[i]` and some with degree
                #    `shift[i] + k * defl[i]` for some `k > 0`. We have
                #    `gcd(shift[i] - exp, shift[i] + k1 * defl[i] - exp, shift[i] + k2 * defl[i] - exp, ...) =`
                #    `gcd(shift[i] - exp, k1 * defl[i], k2 * defl[i], ...)`
                #    Since `gcd(k1, k2, ...) = 1`, this is equal to
                #    `gcd(shift[i] - exp, defl[i])`
                defl[i] = gcd(defl[i], shift[i] - exp)
                shift[i] = exp
            else
                defl[i] = gcd(defl[i], exp - shift[i])
            end
        end
    end
    @assert all(d -> d >= 0, shift)
    @assert all(d -> d >= 0, defl)
    s = prod(variables(p).^shift)::monomialtype(p)
    d = prod(variables(p).^defl)::monomialtype(p)
    return s, d
end

function _zero_to_one_exp(defl::AbstractMonomial)
    # TODO Make it faster by calling something like `mapexponents`.
    return prod(variables(defl).^map(d -> iszero(d) ? one(d) : d, exponents(defl)))
end
function deflate(p::AbstractPolynomialLike, shift, defl)
    if isconstant(shift) && all(d -> isone(d) || iszero(d), exponents(defl))
        return p
    end
    q = MA.operate(deflate, p, shift, _zero_to_one_exp(defl))
    return q
end
function inflate(α, shift, defl)
    return inflate(convert(polynomialtype(shift, typeof(α)), α), shift, defl)
end
function inflate(p::AbstractPolynomialLike, shift, defl)
    if isconstant(shift) && all(d -> isone(d) || iszero(d), exponents(defl))
        return p
    end
    q = MA.operate(inflate, p, shift, _zero_to_one_exp(defl))
    return q
end

function MA.operate(::typeof(deflate), mono::AbstractMonomial, shift, defl)
    mutable_mono = mapexponents(-, mono, shift)
    return mapexponents!(div, mutable_mono, defl)
end
function MA.operate(::typeof(inflate), mono::AbstractMonomial, shift, defl)
    mutable_mono = mapexponents(*, mono, defl)
    return mapexponents!(+, mutable_mono, shift)
end

# Inspired from to `AbstractAlgebra.deflate`
function MA.operate(op::Union{typeof(deflate), typeof(inflate)}, p::AbstractPolynomialLike, shift, defl)
    return polynomial(map(terms(p)) do t
        return term(coefficient(t), MA.operate(op, monomial(t), shift, defl))
    end)
end

function deflated_gcd(p1::APL{T}, p2::APL{S}, algo) where {T, S}
    i1, i2, num_common = _extracted_variable(p1, p2)
    if iszero(i1)
        if iszero(i2)
            return univariate_gcd(p1, p2, algo)
        else
            if isapproxzero(p1)
                return convert(MA.promote_operation(gcd, typeof(p1), typeof(p2)), p2)
            end
            v2 = variables(p2)[i2]
            q2 = isolate_variable(p2, v2)
            g = content(q2, algo)
            return gcd(p1, g, algo)
        end
    else
        if iszero(i2)
            if isapproxzero(p2)
                return convert(MA.promote_operation(gcd, typeof(p1), typeof(p2)), p1)
            end
            v1 = variables(p1)[i1]
            q1 = isolate_variable(p1, v1)
            g = content(q1, algo)
            return gcd(g, p2, algo)
        else
            if num_common > 1
                @assert i1 == i2
                v1 = variables(p1)[i1]
                return multivariate_gcd(p1, p2, v1, algo)
            else
                return univariate_gcd(p1, p2, algo)
            end
        end
    end
end

function Base.gcdx(p1::APL{T}, p2::APL{S}, algo::AbstractUnivariateGCDAlgorithm=GeneralizedEuclideanAlgorithm()) where {T, S}
    i1, i2, num_common = _extracted_variable(p1, p2)
    R = MA.promote_operation(gcd, typeof(p1), typeof(p2))
    if iszero(i1)
        if iszero(i2)
            return univariate_gcdx(p1, p2, algo)
        else
            if isapproxzero(p1)
                return zero(R), one(R), convert(R, p2)
            end
            error("Not implemented yet")
        end
    else
        if iszero(i2)
            if isapproxzero(p2)
                return one(R), zero(R), convert(R, p1)
            end
            error("Not implemented yet")
        else
            if num_common > 1
                @assert i1 == i2
                error("Not implemented yet")
            else
                return univariate_gcdx(p1, p2, algo)
            end
        end
    end
end

# Returns first element in the union of two decreasing vectors
function _extracted_variable(p1, p2)
    v1 = variables(p1)
    v2 = variables(p2)
    i1 = i2 = 1
    best = nothing
    best_var1 = 0
    best_var2 = 0
    num_common = 0
    while i1 <= length(v1) || i2 <= length(v2)
        if i2 > length(v2) || (i1 <= length(v1) && v1[i1] > v2[i2])
            if !iszero(maxdegree(p1, v1[i1]))
                return i1, 0, num_common
            end
            i1 += 1
        elseif i1 > length(v1) || v2[i2] > v1[i1]
            if !iszero(maxdegree(p2, v2[i2]))
                return 0, i2, num_common
            end
            i2 += 1
        else
            @assert v1[i1] == v2[i2]
            v = v1[i1]
            i1 += 1
            i2 += 1
            d1, n1 = deg_num_leading_terms(p1, v)
            d2, n2 = deg_num_leading_terms(p2, v)
            if iszero(d1)
                if iszero(d2)
                    continue
                else
                    return 0, i2-1, num_common
                end
            else
                if iszero(d2)
                    return i1-1, 0, num_common
                end
            end
            if d1 < d2
                d1, d2 = d2, d1
                n1, n2 = n2, n1
            end
            # Heuristic used in `AbstractAlgebra`:
            # https://github.com/Nemocas/AbstractAlgebra.jl/blob/4c6b0a366e550df3db84a665de186111bc3cf8ed/src/generic/MPoly.jl#L4347
            # FIXME what is this based on ? Is there any analysis somewhere comparing different heuristics ?
            cur = max(log(n2) * d1 * d2, log(2) * d2)
            if best === nothing || best > cur
                best = cur
                best_var1 = i1-1
                best_var2 = i2-1
            end
            num_common += 1
        end
    end
    return best_var1, best_var2, num_common
end

function multivariate_gcd(p1::APL, p2::APL, var, algo)
    q = univariate_gcd(isolate_variable(p1, var), isolate_variable(p2, var), algo)
    return sum(coefficient(t) * monomial(t) for t in terms(q))::MA.promote_operation(gcd, typeof(p1), typeof(p2))
end

function isolate_variable(poly::APL, var::AbstractVariable)
    T = typeof(substitute(Subs(), zeroterm(poly), (var,) => (1,)))
    dict = Dict{Int,Vector{T}}()
    for t in terms(poly)
        d = degree(t, var)
        dict_ts = get(dict, d, nothing)
        ts = dict_ts === nothing ? T[] : dict_ts
        push!(ts, subs(t, (var,) => (1,)))
        if dict_ts === nothing
            dict[d] = ts
        end
    end
    return polynomial([
        term(polynomial(ts), var^d) for (d, ts) in dict
    ])
end

"""
    primitive_univariate_gcd(p::APL, q::APL, algo::AbstractUnivariateGCDAlgorithm)

Returns the `gcd` of primitive polynomials `p` and `q` using algorithm `algo`
which is a subtype of [`AbstractUnivariateGCDAlgorithm`](@ref).
"""
function primitive_univariate_gcd end

function not_divided_error(u, v)
    error(
        "Polynomial `$v` of degree `$(maxdegree(v))` and effective",
        " variables `$(effective_variables(v))` does not divide",
        " polynomial `$u` of degree `$(maxdegree(u))` and effective",
        " variables `$(effective_variables(u))`. Did you call",
        " `univariate_gcd` with polynomials with more than one",
        " variable in common ? If yes, call `gcd` instead, otherwise,",
        " please report this.",
    )
end

# If `p` and `q` do not have the same time then the local variables `p` and `q`
# won't be type stable.
function primitive_univariate_gcd(p::APL, q::APL, algo::GeneralizedEuclideanAlgorithm)
    if maxdegree(p) < maxdegree(q)
        return primitive_univariate_gcd(q, p, algo)
    end
    R = MA.promote_operation(gcd, typeof(p), typeof(q))
    u = convert(R, p)
    v = convert(R, q)
    while true
        if isapproxzero(v)
            return u
        elseif isconstant(v)
            # `p` and `q` are primitive so if one of them is constant, it cannot
            # divide the content of the other one.
            return one(R)
        end
        divided, r = pseudo_rem(u, v, algo)
        if !divided
            not_divided_error(u, v)
        end
        if !algo.primitive_rem
            r = primitive_part(r, algo)::R
        end
        u, v = v, r::R
    end
end

function primitive_univariate_gcdx(u0::APL, v0::APL, algo::GeneralizedEuclideanAlgorithm)
    if maxdegree(u0) < maxdegree(v0)
        a, b, g = primitive_univariate_gcdx(v0, u0, algo)
        return b, a, g
    end
    R = MA.promote_operation(gcd, typeof(u0), typeof(v0))
    u = convert(R, u0)
    v = convert(R, v0)
    if isapproxzero(v)
        return one(R), zero(R), u
    elseif isconstant(v)
        # `p` and `q` are primitive so if one of them is constant, it cannot
        # divide the content of the other one.
        return zero(R), one(R), v
    end
    # p * u = q * v + r
    p, q, r = pseudo_divrem(u, v, algo)
    if iszero(q)
        not_divided_error(u, v)
    end
    if iszero(r)
        # Shortcut, does not change the output
        return zero(R), one(R), v
    end
    # TODO
    #if !algo.primitive_rem
    #    r = primitive_part(r, algo)::R
    #end
    # a * v + b * r = g
    # a * v + b * (p * u - q * v) = g
    # b * p * u + (a - b * q) * v = g
    a, b, g = primitive_univariate_gcdx(v, r, algo)
    return p * b, (a - b * q), g
end

function primitive_univariate_gcd(p::APL, q::APL, ::SubresultantAlgorithm)
    error("Not implemented yet")
end

"""
    univariate_gcd(p1::AbstractPolynomialLike, p2::AbstractPolynomialLike, algo::AbstractUnivariateGCDAlgorithm)

Return the *greatest common divisor* of the polynomials `p1` and `p2` that have
at most one variable in common and for which the coefficients are either
`AbstractFloat` or part of a unique factorization domain, e.g., rational numbers,
integers or multivariate polynomials. So `p1` and `p2` should have at most one
variable in common but their coefficients can be multivariate polynomials that
share arbitrarily many variables.

If the coefficients are not `AbstractFloat`, this
1. separates `p1` and `p2` in their [`content`](@ref) and
   [`primitive_part`](@ref) using [`primitive_part_content`](@ref); see
   [Knu14, Algorithm E: E1, p. 426] or [Knu14, Algorithm C: C1, p. 428].
2. Computes the [`gcd`](@ref) of the contents and primitive parts, using
   [`primitive_univariate_gcd`](@ref) for primitive parts.
3. Return the product of these two `gcd`; see
   [Knu14, Algorithm E: E4, p. 427] or [Knu14, Algorithm C: C4, p. 429].

[Knu14] Knuth, D.E., 2014.
*Art of computer programming, volume 2: Seminumerical algorithms.*
Addison-Wesley Professional. Third edition.
"""
function univariate_gcd(p1::APL{S}, p2::APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {S,T}
    return univariate_gcd(algebraic_structure(S, T), p1, p2, algo)
end
function univariate_gcd(::UFD, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm)
    f1, g1 = primitive_part_content(p1, algo)
    f2, g2 = primitive_part_content(p2, algo)
    pp = primitive_univariate_gcd(f1, f2, algo)
    gg = _gcd(g1, g2, algo)#::MA.promote_operation(gcd, typeof(g1), typeof(g2))
    # Multiply each coefficient by the gcd of the contents.
    return mapcoefficientsnz(Base.Fix1(*, gg), pp)
end

function univariate_gcd(::Field, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm)
    return primitive_univariate_gcd(p1, p2, algo)
end

function univariate_gcdx(p1::APL{S}, p2::APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {S,T}
    return univariate_gcdx(algebraic_structure(S, T), p1, p2, algo)
end
function univariate_gcdx(::UFD, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm)
    f1, g1 = primitive_part_content(p1, algo)
    f2, g2 = primitive_part_content(p2, algo)
    a, b, pp = primitive_univariate_gcdx(f1, f2, algo)
    gg = _gcd(g1, g2, algo)#::MA.promote_operation(gcd, typeof(g1), typeof(g2))
    # Multiply each coefficient by the gcd of the contents.
    return g2 * a, g1 * b, g1 * g2 * mapcoefficientsnz(Base.Fix1(*, gg), pp)
end
function univariate_gcdx(::Field, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm)
    return primitive_univariate_gcdx(p1, p2, algo)
end

_gcd(a::APL, b::APL, algo) = gcd(a, b, algo)
_gcd(a, b::APL, algo) = gcd(a, b, algo)
_gcd(a::APL, b, algo) = gcd(a, b, algo)
_gcd(a, b, algo) = gcd(a, b)

_simplifier(a::APL, b::APL, algo) = gcd(a, b, algo)
_simplifier(a, b, algo) = _gcd(a, b, algo)
# Before Julia v1.4, it is not defined.
# After Julia v1.4, it is defined as `gcd of num / lcm of den`.
# We prefer `gcd of den`, otherwise,
# `1/a0 + 1/a1 x + 1/a2 x^2 + ... + 1/an x^n`
# will be transformed into
# `a1*a2*...*an + a0*a2*...*an x + ...`
# which makes the size of the `BigInt`s grow significantly which slows things down.
_simplifier(a::Rational, b::Rational, algo) = gcd(a.num, b.num) // gcd(a.den, b.den)

"""
    content(poly::AbstractPolynomialLike{T}, algo::AbstractUnivariateGCDAlgorithm) where {T}

Return the *content* of the polynomial `poly` over a unique factorization
domain `S` as defined in [Knu14, (3) p. 423].
That is, return the `gcd` of the coefficients of `poly`.
See also [`primitive_part_content`](@ref).

[Knu14] Knuth, D.E., 2014.
*Art of computer programming, volume 2: Seminumerical algorithms.*
Addison-Wesley Professional. Third edition.
"""
function content(poly::APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {T}
    P = MA.promote_operation(gcd, T, T)
    # This is tricky to infer a `content` calls `gcd` which calls `content`, etc...
    # To help Julia break the loop, we annotate the result here.
    return reduce(
        (a, b) -> _simplifier(a, b, algo),
        coefficients(poly),
        init = zero(P),
    )::P
end

"""
    primitive_part_content(poly::AbstractPolynomialLike{T}, algo::AbstractUnivariateGCDAlgorithm) where {T}

Return the *primitive part* of the polynomial `poly` over a unique
factorization domain `S` as defined in [Knu14, (3) p. 423].
That is, return the exact division of `poly` by its [`content`](@ref).
If the content is also needed, call [`primitive_part_content`](@ref)
instead.

[Knu14] Knuth, D.E., 2014.
*Art of computer programming, volume 2: Seminumerical algorithms.*
Addison-Wesley Professional. Third edition.
"""
primitive_part(p::APL, algo::AbstractUnivariateGCDAlgorithm) = primitive_part_content(p, algo)[1]
primitive_part(p::APL{<:AbstractFloat}, algo::AbstractUnivariateGCDAlgorithm) = p

"""
    primitive_part_content(poly::AbstractPolynomialLike{T}, algo::AbstractUnivariateGCDAlgorithm) where {T}

Return the *primitive part* and *content* of the polynomial `poly` over a unique
factorization domain `S` as defined in [Knu14, (3) p. 423]. This is more
efficient to call this function rather than calling [`primitive_part`](@ref) and
[`content`](@ref) separately since computing the primitive part requires
computing the content first and this function avoid computing the content twice.

[Knu14] Knuth, D.E., 2014.
*Art of computer programming, volume 2: Seminumerical algorithms.*
Addison-Wesley Professional. Third edition.
"""
function primitive_part_content(p, algo::AbstractUnivariateGCDAlgorithm)
    g = content(p, algo)
    return mapcoefficientsnz(Base.Fix2(_div, g), p), g
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM.Testing

#= TODO: Fix test.
@testset "calculate mortar matrices and weighted gap vector for 2d model" begin
    nodes = Node[
        [1.0, 1.0],
        [3.0, 1.0],
        [2.0, 2.0],
        [4.0, 4.0]]
    sel = Seg2([1, 2])
    mel = Seg2([3, 4])
    update!([sel, mel], "geometry", nodes)
    calculate_normal_tangential_coordinates!(sel, 0.0)
    sel["master elements"] = [mel]
    prob = MortarProblem("displacement", 2)
    push!(prob, sel)
    a = assemble(prob, 0.0)
    C1 = full(a.C1)
    D = C1[:, [1, 2, 3, 4]]*24
    M = C1[:, [5, 6, 7, 8]]*24
    X = calculate_nodal_vector("geometry", 2, [sel, mel], 0.0)
    @debug begin
        info("nodal vector")
        dump(round(full(X), 3))
    end
    g = -C1*X
    g = full(g)
    @debug begin
        dump(round(D, 3))
        dump(round(M, 3))
        dump(round(g, 3))
    end
    I = eye(2)
    @test isapprox(D, [2*I 4*I; 4*I 14*I])
    @test isapprox(-M, [5*I 1*I; 13*I 5*I])
    @test isapprox(g, 1/6*[0, 2, 0, 7])
end
=#


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

mutable struct VirtualMachineSize <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    numberOfCores::Any # spec type: Union{ Nothing, Int32 } # spec name: numberOfCores
    osDiskSizeInMB::Any # spec type: Union{ Nothing, Int32 } # spec name: osDiskSizeInMB
    resourceDiskSizeInMB::Any # spec type: Union{ Nothing, Int32 } # spec name: resourceDiskSizeInMB
    memoryInMB::Any # spec type: Union{ Nothing, Int32 } # spec name: memoryInMB
    maxDataDiskCount::Any # spec type: Union{ Nothing, Int32 } # spec name: maxDataDiskCount

    function VirtualMachineSize(;name=nothing, numberOfCores=nothing, osDiskSizeInMB=nothing, resourceDiskSizeInMB=nothing, memoryInMB=nothing, maxDataDiskCount=nothing)
        o = new()
        validate_property(VirtualMachineSize, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(VirtualMachineSize, Symbol("numberOfCores"), numberOfCores)
        setfield!(o, Symbol("numberOfCores"), numberOfCores)
        validate_property(VirtualMachineSize, Symbol("osDiskSizeInMB"), osDiskSizeInMB)
        setfield!(o, Symbol("osDiskSizeInMB"), osDiskSizeInMB)
        validate_property(VirtualMachineSize, Symbol("resourceDiskSizeInMB"), resourceDiskSizeInMB)
        setfield!(o, Symbol("resourceDiskSizeInMB"), resourceDiskSizeInMB)
        validate_property(VirtualMachineSize, Symbol("memoryInMB"), memoryInMB)
        setfield!(o, Symbol("memoryInMB"), memoryInMB)
        validate_property(VirtualMachineSize, Symbol("maxDataDiskCount"), maxDataDiskCount)
        setfield!(o, Symbol("maxDataDiskCount"), maxDataDiskCount)
        o
    end
end # type VirtualMachineSize

const _property_map_VirtualMachineSize = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("numberOfCores")=>Symbol("numberOfCores"), Symbol("osDiskSizeInMB")=>Symbol("osDiskSizeInMB"), Symbol("resourceDiskSizeInMB")=>Symbol("resourceDiskSizeInMB"), Symbol("memoryInMB")=>Symbol("memoryInMB"), Symbol("maxDataDiskCount")=>Symbol("maxDataDiskCount"))
const _property_types_VirtualMachineSize = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("numberOfCores")=>"Int32", Symbol("osDiskSizeInMB")=>"Int32", Symbol("resourceDiskSizeInMB")=>"Int32", Symbol("memoryInMB")=>"Int32", Symbol("maxDataDiskCount")=>"Int32")
Base.propertynames(::Type{ VirtualMachineSize }) = collect(keys(_property_map_VirtualMachineSize))
Swagger.property_type(::Type{ VirtualMachineSize }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineSize[name]))}
Swagger.field_name(::Type{ VirtualMachineSize }, property_name::Symbol) =  _property_map_VirtualMachineSize[property_name]

function check_required(o::VirtualMachineSize)
    true
end

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