text:
# Low-level generated functions for working with grids

# _update_buffers => NTuple{N,SMatrix}
# Generate an SArray from the main array and the last SArray
@generated function _update_buffers(
    buffers::Tuple, src::AbstractArray{T}, ::Val{R}, i, j
) where {T,R}
    B = 2R; S = 2R + 1; L = S^2
    newvals = Expr[]
    for n in 0:2B-1
        push!(newvals, :(@inbounds src[i + $n, j + 2R]))
    end
    newbuffers = Expr(:tuple)
    for b in 1:B
        bufvals = Expr(:tuple)
        for n in S+1:L
            push!(bufvals.args, :(@inbounds buffers[$b][$n]))
        end
        for n in b:b+B
            push!(bufvals.args, newvals[n])
        end
        push!(newbuffers.args, :(SArray{Tuple{$S,$S},$T,2,$L}($bufvals)))
    end
    return quote
        return $newbuffers
    end
end

# _initialise_buffers => NTuple{N,SMatrix}
# Generate an SArray from the main array
@generated function _initialise_buffers(src::AbstractArray{T}, ::Val{R}, i, j) where {T,R}
    B = 2R; S = 2R + 1; L = S^2
    columns = []
    zerocol = Expr[]
    for r in 1:2B
        push!(zerocol, :(zero(T)))
    end
    push!(columns, zerocol)
    for c in 0:S-2
        newcol = Expr[]
        for r in 0:2B-1
            push!(newcol, :(@inbounds src[i + $r, j + $c]))
        end
        push!(columns, newcol)
    end
    newbuffers = Expr(:tuple)
    for b in 1:B
        bufvals = Expr(:tuple)
        for c in 1:S, r in b:b+B
            exp = columns[c][r]
            push!(bufvals.args, exp)
        end
        push!(newbuffers.args, :(SArray{Tuple{$S,$S},$T,2,$L}($bufvals)))
    end
    return quote
        return $newbuffers
    end
end

# _getwindow => SMatrix
# Get a single window square from an array, as an SMatrix.
# We use this on GPUs to get a neighborhood window from the main
# grid, which is 10x or more faster than using a view.
# We could possible just use this instead of _update_buffers
# for the sake of simplicity, with some performance loss.
@generated function _getwindow(tile::AbstractArray{T}, ::Neighborhood{R}, i, j) where {T,R}
    S = 2R+1
    L = S^2
    vals = Expr(:tuple)
    for jj in 1:S, ii in 1:S
        push!(vals.args, :(@inbounds tile[$ii + i - 1, $jj + j - 1]))
    end
    return :(SMatrix{$S,$S,$T,$L}($vals))
end

# _readcell
# Returns a single value or NamedTuple of values.
# This occurs for every cell for every rule, so has to be very fast.
@generated function _readcell(data, ::K, I...) where K<:Tuple
    expr = Expr(:tuple)
    keys = map(_unwrap, Tuple(K.parameters))
    for (i, k) in enumerate(keys)
        push!(expr.args, :(@inbounds source(data[$(QuoteNode(k))])[I...]))
    end
    return quote
        keys = $keys
        vals = $expr
        NamedTuple{keys,typeof(vals)}(vals)
    end
end
@inline function _readcell(data::AbstractSimData, ::Val{K}, I...) where K
    @inbounds source(data[K])[I...]
end

# _writecell! => nothing
# Write values to grid/s at index `I`.
# This occurs for every cell for every rule, so has to be very fast.
@generated function _writecell!(data, ::Val{<:CellRule}, wkeys::K, vals::Union{Tuple,NamedTuple}, I...) where K<:Tuple
    expr = Expr(:block)
    keys = map(_unwrap, Tuple(K.parameters))
    for (i, k) in enumerate(keys)
        # MUST write to source(grid) - CellRule doesn't switch grids
        push!(expr.args, :(@inbounds source(data[$(QuoteNode(k))])[I...] = vals[$i]))
    end
    push!(expr.args, :(nothing))
    return expr
end
@inline function _writecell!(data, ::Val{<:CellRule}, wkeys::Val{K}, val, I...) where K
    # MUST write to source(grid) - CellRule doesn't switch grids
    @inbounds source(data[K])[I...] = val
    return nothing
end
@generated function _writecell!(data, ::Val, wkeys::K, vals::Union{Tuple,NamedTuple}, I...) where K<:Tuple
    expr = Expr(:block)
    keys = map(_unwrap, Tuple(K.parameters))
    for (i, k) in enumerate(keys)
        # MUST write to dest(grid) here, not grid K
        # setindex! has overrides for the grid
        push!(expr.args, :(@inbounds dest(data[$(QuoteNode(k))])[I...] = vals[$i]))
    end
    push!(expr.args, :(nothing))
    return expr
end
@inline function _writecell!(data, ::Val, wkeys::Val{K}, val, I...) where K
    # MUST write to dest(grid) here, not grid K
    # setindex! has overrides for the grid
    @inbounds dest(data[K])[I...] = val
    return nothing
end

# _getreadgrids => Union{ReadableGridData,Tuple{ReadableGridData,Vararg}}
# Retrieves `GridData` from a `SimData` object to match the requirements of a `Rule`.
# Returns a `Tuple` holding the key or `Tuple` of keys, and grid or `Tuple` of grids.
@generated function _getreadgrids(::Rule{R,W}, data::AbstractSimData) where {R<:Tuple,W}
    Expr(:tuple,
        Expr(:tuple, (:(Val{$(QuoteNode(key))}()) for key in R.parameters)...),
        Expr(:tuple, (:(data[$(QuoteNode(key))]) for key in R.parameters)...),
    )
end
@generated function _getreadgrids(::Rule{R,W}, data::AbstractSimData) where {R,W}
    :((Val{$(QuoteNode(R))}(), data[$(QuoteNode(R))]))
end

# _getwritegrids => Union{WritableGridData,Tuple{WritableGridData,Vararg}}
# Retrieves `GridData` from a `SimData` object to match the requirements of a `Rule`.
# Returns a `Tuple` holding the key or `Tuple` of keys, and grid or `Tuple` of grids.
@generated function _getwritegrids(::Rule{R,W}, data::AbstractSimData) where {R,W<:Tuple}
    Expr(:tuple,
        Expr(:tuple, (:(Val{$(QuoteNode(key))}()) for key in W.parameters)...),
        Expr(:tuple, (:(WritableGridData(data[$(QuoteNode(key))])) for key in W.parameters)...),
    )
end
@generated function _getwritegrids(::Rule{R,W}, data::AbstractSimData) where {R,W}
    :((Val{$(QuoteNode(W))}(), WritableGridData(data[$(QuoteNode(W))])))
end

# _combinegrids => AbstractSimData
# Combine grids into a new NamedTuple of grids depending
# on the read and write keys required by a rule.
@inline function _combinegrids(data::AbstractSimData, wkeys, wgrids)
    allkeys = map(Val, keys(data))
    allgrids = values(data)
    _combinegrids(data, allkeys, allgrids, wkeys, wgrids)
end
@inline function _combinegrids(data::AbstractSimData, rkeys, rgrids, wkeys, wgrids)
    @set data.grids = _combinegrids(rkeys, rgrids, wkeys, wgrids)
end
@inline function _combinegrids(rkey, rgrids, wkey, wgrids)
    _combinegrids((rkey,), (rgrids,), (wkey,), (wgrids,))
end
@inline function _combinegrids(rkey, rgrids, wkeys::Tuple, wgrids::Tuple)
    _combinegrids((rkey,), (rgrids,), wkeys, wgrids)
end
@inline function _combinegrids(rkeys::Tuple, rgrids::Tuple, wkey, wgrids)
    _combinegrids(rkeys, rgrids, (wkey,), (wgrids,))
end
@generated function _combinegrids(rkeys::Tuple, rgrids::Tuple, wkeys::Tuple, wgrids::Tuple)
    rkeys = _vals2syms(rkeys)
    wkeys = _vals2syms(wkeys)
    keysexp = Expr(:tuple, QuoteNode.(wkeys)...)
    dataexp = Expr(:tuple, :(wgrids...))

    for (i, key) in enumerate(rkeys)
        if !(key in wkeys)
            push!(dataexp.args, :(rgrids[$i]))
            push!(keysexp.args, QuoteNode(key))
        end
    end

    return quote
        keys = $keysexp
        vals = $dataexp
        NamedTuple{keys,typeof(vals)}(vals)
    end
end

# _replacegrids => AbstractSimData
# Replace grids in a NamedTuple with new grids where required.
function _replacegrids(data::AbstractSimData, newkeys, newgrids)
    @set data.grids = _replacegrids(grids(data), newkeys, newgrids)
end
@generated function _replacegrids(allgrids::NamedTuple, newkeys::Tuple, newgrids::Tuple)
    newkeys = map(_unwrap, newkeys.parameters)
    allkeys = allgrids.parameters[1]
    expr = Expr(:tuple)
    for key in allkeys
        if key in newkeys
            i = findfirst(k -> k == key, newkeys)
            push!(expr.args, :(newgrids[$i]))
        else
            push!(expr.args, :(allgrids.$key))
        end
    end

    return quote
        vals = $expr
        NamedTuple{$allkeys,typeof(vals)}(vals)
    end
end
@generated function _replacegrids(allgrids::NamedTuple, newkey::Val, newgrid::GridData)
    newkey = _unwrap(newkey)
    allkeys = allgrids.parameters[1]
    expr = Expr(:tuple)
    for key in allkeys
        if key == newkey
            push!(expr.args, :(newgrid))
        else
            push!(expr.args, :(allgrids.$key))
        end
    end

    return quote
        vals = $expr
        NamedTuple{$allkeys,typeof(vals)}(vals)
    end
end

# _vals2syms => Union{Symbol,Tuple}
# Get symbols from a Val or Tuple type
@inline _vals2syms(x::Type{<:Tuple}) = map(v -> _vals2syms(v), x.parameters)
@inline _vals2syms(::Type{<:Val{X}}) where X = X


text:
using FirstOrderOptimization

lr = LowRankOperator(rand(5), rand(4)')
lr*rand(4)
lr'*rand(5)

svds(lr, nsv=1)

pobs = .5
obs = sprand(5, 4, pobs).>0
iobs = find(obs)
io = IndexingOperator(5, 4, iobs)
io*rand(5,4)
io'*rand(length(iobs))

# svds(io, nsv=1)

ilr = IndexedLowRankOperator(io,rand(length(iobs)))
size(ilr)==(5,4)
ilr*rand(4)
ilr'*rand(5)

svds(ilr, nsv=1)

# test thin update

u,s,v = FirstOrderOptimization.svd(rand(4,5), 2)
A = LowRankOperator(u,spdiagm(s),v, transpose=(:N,:N,:T))
Delta = LowRankOperator(rand(4,1), rand(5,1), transpose=(:N,:T))
Aup = Array(A) + Array(Delta)
thin_update!(A,Delta,1.)
@assert norm(Array(A) - Aup) <= 1e-10

# test dot between sparse and low rank (svd) operator

m,n = size(A)
G = sprand(m,n,.3)
@assert norm(dot(G,A) - dot(full(G),Array(A))) <= 1e-10

ext: jl

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
"""
    polyeval(z, p, N)

Evaluate the polynomial ``\\sum_k c[k] z^{k-1}`` for the coefficients `c[1]`, `c[2]`, ..., `c[N]`;
that is, the coefficients are given in ascending order by power of `z`.  This macro expands
to efficient inline code that uses either Horner's method.

```jldoctest

julia> p = randn(5)
5-element Array{Float64,1}:
  0.106455
  0.0351716
 -0.204764
 -0.87642
  0.260407

julia> polyeval(0.5, p, 5)
-0.020427265391245605

```
"""
function polyeval(x, p::Vector, N)
    return evalpoly(x,@view p[1:N])
end

function polyeval(x, p::Tuple, N)
    return evalpoly(x,p[1:N])
end

"""
    calcz(vm, b0, c0, [EPS, [MAXITER, [relax]]])

Calculates the compressibility factor using a virial equation:

`z = 1 + b₀/z + c₀/vₘ²`

The algorithm transforms the values to a cubic polynomial and solves the roots using
cardan's method.

 * `b0` b₀ parameter of the virial equation
 * `c0` c₀ parameter of the virial equation
"""
function calcz(b0, c0)
    #transformation:
    #z3 = z2 + b₀z + c0
    #z3-z2- b₀z- c0
    _c0,_b0,_1 = promote(c0,b0,1.0)
    poly =(-_c0,-_b0,-_1,_1)
    f0(x) = evalpoly(x,poly)
    xs = _1 #supposing ideal gas
    return Roots.secant_method(f0,xs)

end


text:
module Covers
    import Base.size
    import Base.getindex
    import Base.setindex!

    using DefaultArrays
    export Cover, multicovered, ncovered

    struct Cover{N} <: AbstractArray{Bool,N}
        size::NTuple{N,Int}
        cov::Vector{BitArray{1}}
    end

    Cover(size::Int...) = Cover(size,[falses(s) for s in size])

    size(A::Cover) = A.size

    @inline getindex(A::Cover{1}, i::Int) = _cartesian_getindex(A, i)
    @inline getindex(A::Cover{N}, i::Int) where N = _scalar_getindex(A, i)
    @inline getindex(A::Cover{N}, I::Int...) where N = _cartesian_getindex(A, I...)

    @inline function _scalar_getindex(A::Cover,i::Int)
        I=CartesianIndices(A)[i]
        _cartesian_getindex(A::Cover, Tuple(I)...)
    end

    @inline function _cartesian_getindex(A::Cover{N}, I::Int...) where N
        @boundscheck checkbounds(CartesianIndices(A),I...)
        for i in 1:N
            if A.cov[i][I[i]]
                return true
            end
        end
        return false
    end

    @inline setindex!(A::Cover{1}, v::Bool, i::Int) = _cartesian_setindex!(A, v, i)
    @inline setindex!(A::Cover{N}, v::Bool, i::Int) where N = _scalar_setindex!(A, v, i)
    @inline setindex!(A::Cover{N}, v::Bool, I::Int...) where N = _cartesian_setindex!(A, v, I...)

    @inline function _cartesian_setindex!(A::Cover{N}, v::Bool,I::Int...) where N
        state=A[I...]
        if v && !state
            idx=findmax(size(A))[2]
            A(idx,I[idx],v)
            @goto ret
        end
        if !v && state
            for idx in 1:N
                A(idx,I[idx],v)
            end
            @goto ret
        end
        @label ret
        return v
    end
    @inline function _scalar_setindex!(A::Cover, v::Bool,i::Int)
        I=CartesianIndices(A)[i]
        _cartesian_setindex!(A::Cover,v, Tuple(I)...)
    end

    @inline function (C::Cover)(dim::Int,i::Int,state::Bool)
        C.cov[dim][i]=state
    end

    @inline function multicovered(A::Cover)
        ret=DefaultArray(false,A.size)
        for t in eachindex(A)
            ret[t]=coveredntimes(A, Tuple(t))>1
        end
        ret
    end

    @inline function coveredntimes(A::Cover{N}, I::NTuple{N,Int}) where N
        @boundscheck checkbounds(CartesianIndices(A),I...)
        sum(A.cov[i][I[i]] for i in eachindex(I))
    end

    @inline function ncovered(A::Cover)
        ret=DefaultArray(0,A.size)
        for t in eachindex(A)
            ret[t]=coveredntimes(A, Tuple(t))
        end
        ret
    end

end # module


text:
for quad_degree = 1:20 # Exceeding degree 20 seems unnecessary at this time
    @eval begin
        # Square
        @generated function gauss_quadrature(form::Val{:legendre},
                                  shape::RefSquare,
                                  degree::Val{$quad_degree},
                                  type::Type{T}) where {T}
            # Is there a way to use form and degree here instead?
            line_weights, line_points = gauss_quadrature(Val(:legendre),
                                                         RefLine(),
                                                         Val($quad_degree),
                                                         T)
            weights = Expr[]
            points = Expr[]
            for i = 1:$quad_degree
                for j = 1:$quad_degree
                    push!(weights, :($(line_weights[i]) * $(line_weights[j])))
                    push!(points, :(tuple($(line_points[i][1]), $(line_points[j][1]))))
                end
            end
            return quote
                return (SVector{$(length(line_weights)^2), $T}(tuple($(weights...))),
                        SVector{$(length(line_weights)^2), NTuple{2,$T}}(tuple($(points...))))
            end
        end
        # Cube
        @generated function gauss_quadrature(form::Val{:legendre},
                                  shape::RefCube,
                                  degree::Val{$quad_degree},
                                  type::Type{T}) where {T}
            # Is there a way to use form and degree here instead?
            line_weights, line_points = gauss_quadrature(Val(:legendre),
                                                         RefLine(),
                                                         Val($quad_degree),
                                                         T)
            weights = Expr[]
            points = Expr[]
            for i = 1:$quad_degree
                for j = 1:$quad_degree
                    for k = 1:$quad_degree
                        push!(weights, :($(line_weights[i]) *
                                         $(line_weights[j]) *
                                         $(line_weights[k])
                                        ))
                        push!(points, :(tuple($(line_points[i][1]),
                                              $(line_points[j][1]),
                                              $(line_points[k][1]),
                                             )))
                    end
                end
            end
            return quote
                return (SVector{$(length(line_weights)^3), $T}(tuple($(weights...))),
                        SVector{$(length(line_weights)^3), NTuple{3,$T}}(tuple($(points...))))
            end
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
    Randomwalk()

Recommendation model using random walk with restart on user-item bipartite graph. Implemented algorithm is based on Pixie random walk.

# References
C.  Eksombatchai (2018), [Pixie: A System for Recommending 3+ Billion Items to 200+ Million Users in Real-Time](http://dl.acm.org/citation.cfm?doid=3178876.3186183)
"""
mutable struct Randomwalk <: AbstractRecommender
    adjacency_list::Union{Vector{Int},Nothing}
    offsets::Union{Vector{Int},Nothing}

    user_history::Any
    user2uidx::Union{Dict,Nothing}
    item2iidx::Union{Dict,Nothing}
    iidx2item::Union{Dict,Nothing}
    max_degree::Union{Int,Nothing}

    Randomwalk() = new(nothing, nothing, nothing, nothing, nothing, nothing, nothing)
end

"""
    fit!(model::Randomwalk, table; col_user = :userid, col_item = :itemid)

Build bipartite graph from `table`. One side of the graph collcets user nodes, and the others item nodes. If a user actions an item, an edge is added between them. The graph is undirected, and has no extra weights.
"""
function fit!(model::Randomwalk, table; col_user = :userid, col_item = :itemid, kwargs...)
    model.adjacency_list, model.offsets, model.user2uidx, model.item2iidx =
        build_graph(table, col_user = col_user, col_item = col_item)

    model.iidx2item = Dict(iidx => item for (item, iidx) in model.item2iidx)

    model.user_history = Dict{Int,Vector{Int}}()
    for (userid, history) in
        zip(make_u2i_dataset(table, col_user = col_user, col_item = col_item)...)
        uidx = model.user2uidx[userid]
        history = [model.item2iidx[item] for item in history]
        model.user_history[uidx] = history
    end

    model.max_degree = get_max_degree(model.offsets)
end

@doc raw"""
    predict_u2i(model::Randomwalk, userid::Union{AbstractString,Int}, n::Int64; drop_history = false, terminate_prob = 0.1, total_walk_length = 10000, min_high_visited_candidates = Inf, high_visited_count_threshold = Inf, pixie_walk_length_scaling = false, pixie_multi_hit_boosting = false, aggregate_function = sum)

Make recommendation by random walk with restart. Basic algorithm is as follows:

1. Get items that are already consumed by the user (on the graph, they are connected by one step). We denote them by ``q \in Q``.
2. Starting from each node ``q \in Q``, perform multiple random walks with certain stop probability. Record the visited count of the items on the walk. We denote the counts of item ``p`` on the walk from ``q`` by ``V_q[p]``.
3. Finally aggregate ``V_q[p]`` to ``V[p]``, and recommeds top-scored items. Two mothods for aggregation are provided
- Simple aggregation: Taking sum, ``V[p] = \sum_{q\in Q} V_q[p]``. You can also replace `sum` by, for instance, `maximum`.
- Pixie boosting: ``V[p] = (\sum_{q\in Q} \sqrt{V_q[p]})^2``, putting more importance on the nodes visited by ``q``s.

# Model-specific arguments
- `terminate_prob`: stop probability of one random walk
- `total_walk_length`: total walk length over the multiple walks from ``q``'s.
- `high_visited_count_threshold`: early stopping paramerer. Count up `high_visited_count` when the visited count of certain node reaches this threshold.
- `min_high_visited_candidates`: early stopping parameter. Terminate the walk from some node ``q`` if `high_visited_count` reaches `min_high_visited_candidates`.
- `pixie_walk_length_scaling`: If set to true, the start node ``q`` with more degree will be given more walk length. If false, the walk length is the same over all the nodes ``q \in Q``
- `pixie_multi_hit_boosting`: If true, pixie boosting is adopted for aggregation. If false, simple aggregation is used.
- `aggregate_function`: function used by simple aggregation.
"""
function predict_u2i(
    model::Randomwalk,
    userid::Union{AbstractString,Int},
    n::Int64;
    drop_history = false,
    terminate_prob = 0.1,
    total_walk_length = 10000,
    min_high_visited_candidates = Inf,
    high_visited_count_threshold = Inf,
    pixie_walk_length_scaling = false,
    pixie_multi_hit_boosting = false,
    aggregate_function = sum,
    kwargs...,
)
    if !(userid in keys(model.user2uidx))
        return []
    end

    uidx = model.user2uidx[userid]
    query_nodeids = model.user_history[uidx]

    visited_count = randomwalk_multiple(
        model.adjacency_list,
        model.offsets,
        query_nodeids,
        terminate_prob,
        total_walk_length,
        min_high_visited_candidates,
        high_visited_count_threshold,
        pixie_walk_length_scaling,
        pixie_multi_hit_boosting,
        model.max_degree,
        aggregate_function,
    )

    sorted_idx = sortperm(collect(values(visited_count)), rev = true)
    pred_iidx = collect(keys(visited_count))[sorted_idx]

    # this is very slow
    if drop_history
        filter!(p -> !(p in query_nodeids), pred_iidx)
    end

    pred_items = [model.iidx2item[iidx] for iidx in pred_iidx]

    n = min(n, length(pred_items))
    return pred_items[1:n]
end


text:
"""
    cv_iht(y, x, z; path=1:20, q=5, d=Normal(), l=IdentityLink(), ...)

For each model specified in `path`, performs `q`-fold cross validation and
returns the (averaged) deviance residuals. The purpose of this function is to
find the best sparsity level `k`, obtained from selecting the model with the
minimum out-of-sample error. Note sparsity is enforced on `x` only, unless `zkeep`
keyword is specified.

To check if multithreading is enabled, check output of `Threads.nthreads()`.

# Arguments:
+ `y`: Phenotype vector or matrix. Should be an `Array{T, 1}` (single traits) or
    `Array{T, 2}` (multivariate Gaussian traits). For multivariate traits, each
    column of `y` should be a sample.
+ `x`: Genotype matrix (an `Array{T, 2}` or `SnpLinAlg`). For univariate
    analysis, samples are rows of `x`. For multivariate analysis, samples are
    columns of `x` (i.e. input `Transpose(x)` for `SnpLinAlg`)
+ `z`: Matrix of non-genetic covariates of type `Array{T, 2}` or `Array{T, 1}`.
    For univariate analysis, sample covariates are rows of `z`. For multivariate
    analysis, sample covariates are columns of `z`. If this is not specified, an
    intercept term will be included automatically. If `z` is specified, make sure
    the first column (row) is all 1s to represent the intercept.

# Optional Arguments:
- `path`: Different values of `k` that should be tested. One can input a vector of
    `Int` (e.g. `path=[5, 10, 15, 20]`) or a range (default `path=1:20`).
- `q`: Number of cross validation folds. Larger means more accurate and more computationally
    intensive. Should be larger 2 and smaller than 10. Default `q=5`.
- `d`: Distribution of phenotypes. Specify `Normal()` for quantitative traits,
    `Bernoulli()` for binary traits, `Poisson()` or `NegativeBinomial()` for
    count traits, and `MvNormal()` for multiple quantitative traits.
- `l`: A link function. The recommended link functions are `l=IdentityLink()` for
    quantitative traits, `l=LogitLink()` for binary traits, `l=LogLink()` for Poisson
    distribution, and `l=Loglink()` for NegativeBinomial distribution. For multivariate
    analysis, the choice of link does not matter.
- `zkeep`: BitVector determining whether non-genetic covariates in `z` will be subject
    to sparsity constraint. `zkeep[i] = true` means covariate `i` will NOT be projected.
    Note covariates forced in the model are not subject to sparsity constraints in `path`.
- `est_r`: Symbol (`:MM`, `:Newton` or `:None`) to estimate nuisance parameters for negative binomial regression
- `group`: vector storing group membership for each predictor
- `weight`: vector storing vector of weights containing prior knowledge on each predictor
- `folds`: Vector that separates the sample into `q` disjoint subsets
- `debias`: Boolean indicating whether we should debias at each IHT step. Defaults `false`
- `verbose`: Boolean indicating whether to print progress and mean squared error for
    each `k` in `path`. Defaults `true`
- `max_iter`: is the maximum IHT iteration for a model to converge. Defaults to 100
- `min_iter`: is the minimum IHT iteration before checking for convergence. Defaults to 5.
- `init_beta`: Whether to initialize beta values to univariate regression values.
    Currently only Gaussian traits can be initialized. Default `false`.
- `memory_efficient`: If `true,` it will cause ~1.1 times slow down but one only
    needs to store the genotype matrix (requiring 2np bits for PLINK binary files
    and `8np` bytes for other formats). If `memory_efficient=false`, one may potentially
    store `t` sparse matrices of dimension `8nk` bytes where `k` are the cross validated
    sparsity levels.

# Output
- `mse`: A vector of mean-squared error for each `k` specified in `path`.
"""
function cv_iht(
    y        :: AbstractVecOrMat{T},
    x        :: AbstractMatrix{T},
    z        :: AbstractVecOrMat{T};
    d        :: Distribution = is_multivariate(y) ? MvNormal(T[]) : Normal(),
    l        :: Link = IdentityLink(),
    path     :: AbstractVector{<:Integer} = 1:20,
    q        :: Int64 = 5,
    est_r    :: Symbol = :None,
    group    :: AbstractVector{Int} = Int[],
    weight   :: AbstractVector{T} = T[],
    zkeep    :: BitVector = trues(size(y, 2) > 1 ? size(z, 1) : size(z, 2)),
    folds    :: AbstractVector{Int} = rand(1:q, is_multivariate(y) ? size(x, 2) : size(x, 1)),
    debias   :: Bool = false,
    verbose  :: Bool = true,
    max_iter :: Int = 100,
    min_iter :: Int = 5,
    init_beta :: Bool = false,
    memory_efficient :: Bool = true
    ) where T <: Float

    typeof(x) <: AbstractSnpArray && throw(ArgumentError("x is a SnpArray! Please convert it to a SnpLinAlg first!"))
    check_data_dim(y, x, z)
    verbose && print_iht_signature()

    # preallocated arrays for efficiency
    test_idx  = [falses(length(folds)) for i in 1:Threads.nthreads()]
    train_idx = [falses(length(folds)) for i in 1:Threads.nthreads()]
    V = [initialize(x, z, y, 1, 1, d, l, group, weight, est_r, false, zkeep,
        memory_efficient=memory_efficient) for i in 1:Threads.nthreads()]

    # for displaying cross validation progress
    pmeter = verbose ? Progress(q * length(path), "Cross validating...") : nothing

    # cross validate. TODO: wrap pmap with batch_size keyword to enable distributed CV
    combinations = allocate_fold_and_k(q, path)
    mses = zeros(length(combinations))
    ThreadPools.@qthreads for i in 1:length(combinations)
        fold, sparsity = combinations[i]

        # assign train/test indices
        id = Threads.threadid()
        v = V[id]
        test_idx[id]  .= folds .== fold
        train_idx[id] .= folds .!= fold

        # run IHT on training data with current (fold, sparsity)
        v.k = sparsity
        init_iht_indices!(v, init_beta, train_idx[id], false)
        fit_iht!(v, debias=debias, verbose=false, max_iter=max_iter, min_iter=min_iter)

        # predict on validation data
        v.cv_wts[train_idx[id]] .= zero(T)
        v.cv_wts[test_idx[id]] .= one(T)
        mses[i] = predict!(v)

        # update progres
        verbose && next!(pmeter)
    end

    # weight mses for each fold by their size before averaging
    mse = meanloss(mses, q, folds)

    # find best model size and print cross validation result
    k = path[argmin(mse)] :: Int
    verbose && print_cv_results(mse, path, k)

    return mse
end

# Distributed memory cv_iht: (TODO merge pmap with multithreading)
# function cv_iht(
#     y        :: AbstractVecOrMat{T},
#     x        :: AbstractMatrix{T},
#     z        :: AbstractVecOrMat{T};
#     d        :: Distribution = is_multivariate(y) ? MvNormal(T[]) : Normal(),
#     l        :: Link = IdentityLink(),
#     path     :: AbstractVector{<:Integer} = 1:20,
#     q        :: Int64 = 5,
#     est_r    :: Symbol = :None,
#     group    :: AbstractVector{Int} = Int[],
#     weight   :: AbstractVector{T} = T[],
#     folds    :: AbstractVector{Int} = rand(1:q, is_multivariate(y) ? size(x, 2) : size(x, 1)),
#     debias   :: Bool = false,
#     verbose  :: Bool = true,
#     parallel :: Bool = true,
#     max_iter :: Int = 100,
#     min_iter :: Int = 20,
#     init_beta :: Bool = false
#     ) where T <: Float

#     typeof(x) <: AbstractSnpArray && throw(ArgumentError("x is a SnpArray! Please convert it to a SnpLinAlg first!"))
#     check_data_dim(y, x, z)

#     # preallocated arrays for efficiency
#     test_idx  = [falses(length(folds)) for i in 1:nprocs()]
#     train_idx = [falses(length(folds)) for i in 1:nprocs()]
#     V = [initialize(x, z, y, 1, 1, d, l, group, weight, est_r, init_beta,
#         cv_train_idx=train_idx[i]) for i in 1:nprocs()]

#     # for displaying cross validation progress
#     pmeter = Progress(q * length(path), "Cross validating...")
#     channel = RemoteChannel(()->Channel{Bool}(q * length(path)), 1)
#     @async while take!(channel)
#         next!(pmeter)
#     end

#     # cross validation all (fold, k) combinations
#     combinations = allocate_fold_and_k(q, path)
#     mses = (parallel ? pmap : map)(combinations) do (fold, k)
#         # assign train/test indices
#         id = myid()
#         test_idx[id]  .= folds .== fold
#         train_idx[id] .= folds .!= fold
#         v = V[id]

#         # run IHT on training data with current (fold, k)
#         v.k = k
#         init_iht_indices!(v, init_beta, train_idx[id])
#         fit_iht!(v, debias=debias, verbose=false, max_iter=max_iter, min_iter=min_iter)

#         # predict on validation data
#         v.cv_wts[train_idx[id]] .= zero(T)
#         v.cv_wts[test_idx[id]] .= one(T)
#         mse = predict!(v)

#         # update progres
#         put!(channel, true)

#         return mse
#     end
#     put!(channel, false)

#     #weight mses for each fold by their size before averaging
#     mse = meanloss(mses, q, folds)

#     # find best model size and print cross validation result
#     k = path[argmin(mse)] :: Int
#     verbose && print_cv_results(mse, path, k)

#     return mse
# end

function cv_iht(y::AbstractVecOrMat{T}, x::AbstractMatrix; kwargs...) where T
    z = is_multivariate(y) ? ones(T, 1, size(y, 2)) : ones(T, length(y))
    return cv_iht(y, x, z; kwargs...)
end

"""
    allocate_fold_and_k(q, path)

Returns all combinations of (qᵢ, kᵢ) where `q` is number of cross validation fold
and path is a vector of different `k` to be tested.
"""
function allocate_fold_and_k(q::Int, path::AbstractVector{<:Integer})
    combinations = Tuple{Int, Int}[]
    for fold in 1:q, k in path
        push!(combinations, (fold, k))
    end
    return combinations
end

"""
Runs IHT across many different model sizes specifed in `path` using the full
design matrix. Same as `cv_iht` but **DOES NOT** validate in a holdout set, meaning
that this will definitely induce overfitting as we increase model size.
Use this if you want to quickly estimate a range of feasible model sizes before
engaging in full cross validation.
"""
function iht_run_many_models(
    y        :: AbstractVecOrMat{T},
    x        :: AbstractMatrix,
    z        :: AbstractVecOrMat{T};
    d        :: Distribution = Normal(),
    l        :: Link = canonicallink(d),
    path     :: AbstractVector{Int} = 1:20,
    est_r    :: Symbol = :None,
    group    :: AbstractVector{Int} = Int[],
    weight   :: AbstractVector{T} = Float64[],
    use_maf  :: Bool = false,
    debias   :: Bool = false,
    verbose  :: Bool = true,
    parallel :: Bool = false,
    max_iter :: Int = 100
    ) where T <: Float

    # for each k, fit model and store the loglikelihoods
    results = (parallel ? pmap : map)(path) do k
        if typeof(x) == SnpArray
            xla = SnpLinAlg{T}(x, model=ADDITIVE_MODEL, center=true, scale=true,
                impute=true)
            return fit_iht(y, xla, z, J=1, k=k, d=d, l=l, est_r=est_r, group=group,
                weight=weight, use_maf=use_maf, debias=debias, verbose=false,
                max_iter=max_iter)
        else
            return fit_iht(y, x, z, J=1, k=k, d=d, l=l, est_r=est_r, group=group,
                weight=weight, use_maf=use_maf, debias=debias, verbose=false,
                max_iter=max_iter)
        end
    end

    loglikelihoods = zeros(size(path, 1))
    for i in 1:length(results)
        loglikelihoods[i] = results[i].logl
    end

    #display result and then return
    verbose && print_a_bunch_of_path_results(loglikelihoods, path)
    return loglikelihoods
end

function iht_run_many_models(y::AbstractVecOrMat{T}, x::AbstractMatrix; kwargs...) where T
    z = is_multivariate(y) ? ones(T, 1, size(y, 2)) : ones(T, length(y))
    return iht_run_many_models(y, x, z; kwargs...)
end

function predict!(v::IHTVariable{T, M}) where {T <: Float, M}
    # first update mean μ with estimated (trained) beta (cv weights are handled in deviance)
    update_xb!(v)
    update_μ!(v)

    # Compute deviance residual (MSE for Gaussian response)
    return deviance(v)
end

function predict!(v::mIHTVariable{T, M}) where {T <: Float, M}
    # first update mean μ with estimated (trained) beta
    update_xb!(v)
    update_μ!(v)

    # Compute MSE
    mse = zero(T)
    @inbounds for j in 1:size(v.Y, 2), i in 1:ntraits(v)
        mse += abs2(v.Y[i, j] - v.μ[i, j]) * v.cv_wts[j]
    end
    return mse
end

"""
Scale mean squared errors (deviance residuals) for each fold by its own fold size.
"""
function meanloss(fitloss::AbstractVector, q::Int64, folds::AbstractVector{Int})
    ninfold = zeros(Int, q)
    for fold in folds
        ninfold[fold] += 1
    end

    pathsize = Int(length(fitloss) / q)
    loss = zeros(eltype(fitloss), pathsize)
    for j = 1:q
        wfold = convert(eltype(fitloss), ninfold[j] / length(folds))
        for i = 1:pathsize
            loss[i] += fitloss[i + (j - 1) * pathsize] * wfold
        end
    end

    return loss
end


text:
using KernelAbstractions.Extras.LoopInfo: @unroll

#####
##### Periodic boundary conditions
#####

@inline parent_and_size(c, dim1, dim2) = (parent(c),  size(parent(c))[[dim1, dim2]])

@inline function parent_and_size(c::NTuple, dim1, dim2)
    p = parent.(c)
    p_size = (minimum([size(t, dim1) for t in p]), minimum([size(t, dim2) for t in p]))
    return p, p_size
end

function fill_west_and_east_halo!(c, ::PBCT, ::PBCT, loc, arch, dep, grid, args...; kw...)
    c_parent, yz_size = parent_and_size(c, 2, 3)
    event = launch!(arch, grid, yz_size, fill_periodic_west_and_east_halo!, c_parent, grid.Hx, grid.Nx; dependencies=dep, kw...)
    return event
end

function fill_south_and_north_halo!(c, ::PBCT, ::PBCT, loc, arch, dep, grid, args...; kw...)
    c_parent, xz_size = parent_and_size(c, 1, 3)
    event = launch!(arch, grid, xz_size, fill_periodic_south_and_north_halo!, c_parent, grid.Hy, grid.Ny; dependencies=dep, kw...)
    return event
end

function fill_bottom_and_top_halo!(c, ::PBCT, ::PBCT, loc, arch, dep, grid, args...; kw...)
    c_parent, xy_size = parent_and_size(c, 1, 2)
    event = launch!(arch, grid, xy_size, fill_periodic_bottom_and_top_halo!, c_parent, grid.Hz, grid.Nz; dependencies=dep, kw...)
    return event
end

#####
##### Periodic boundary condition kernels
#####

@kernel function fill_periodic_west_and_east_halo!(c, H::Int, N)
    j, k = @index(Global, NTuple)
    @unroll for i = 1:H
        @inbounds begin
            c[i, j, k]     = c[N+i, j, k] # west
            c[N+H+i, j, k] = c[H+i, j, k] # east
        end
    end
end

@kernel function fill_periodic_south_and_north_halo!(c, H::Int, N)
    i, k = @index(Global, NTuple)
    @unroll for j = 1:H
        @inbounds begin
            c[i, j, k] = c[i, N+j, k]     # south
            c[i, N+H+j, k] = c[i, H+j, k] # north
        end
    end
end

@kernel function fill_periodic_bottom_and_top_halo!(c, H::Int, N)
    i, j = @index(Global, NTuple)
    @unroll for k = 1:H
        @inbounds begin
            c[i, j, k] = c[i, j, N+k]        # top
            c[i, j, N+H+k] = c[i, j, H+k] # bottom
        end
    end
end

####
#### Tupled periodic boundary condition
####

@kernel function fill_periodic_west_and_east_halo!(c::NTuple{M}, H::Int, N) where M
    j, k = @index(Global, NTuple)
    @unroll for n = 1:M
        @unroll for i = 1:H
            @inbounds begin
                  c[n][i, j, k]     = c[n][N+i, j, k] # west
                  c[n][N+H+i, j, k] = c[n][H+i, j, k] # east
            end
        end
    end
end

@kernel function fill_periodic_south_and_north_halo!(c::NTuple{M}, H::Int, N) where M
    i, k = @index(Global, NTuple)
    @unroll for n = 1:M
        @unroll for j = 1:H
            @inbounds begin
                c[n][i, j, k]     = c[n][i, N+j, k] # south
                c[n][i, N+H+j, k] = c[n][i, H+j, k] # north
            end
        end
    end
end

@kernel function fill_periodic_bottom_and_top_halo!(c::NTuple{M}, H::Int, N) where M
    i, j = @index(Global, NTuple)
    @unroll for n = 1:M
        @unroll for k = 1:H
            @inbounds begin
                c[n][i, j, k]     = c[n][i, j, N+k] # top
                c[n][i, j, N+H+k] = c[n][i, j, H+k] # bottom
            end
        end
    end
end


text:
using Dierckx
using Statistics: mean

export BasicBody,Ellipse,Circle,Rectangle,Square,Plate,SplinedBody,NACA4

"""
    BasicBody(x,y[,closuretype=ClosedBody]) <: Body

Construct a body by simply passing in the `x` and `y` coordinate vectors. The last
point will be automatically connected to the first point. The coordinate vectors
are assumed to be expressed in the body-fixed coordinate system. The optional
`closuretype` specifies whether the body is closed (`ClosedBody`) or open (`OpenBody`).
If closed, then the first and last points are assumed joined in operations that
require neighbor points.
"""
mutable struct BasicBody{N,C<:BodyClosureType} <: Body{N,C}
  cent :: Tuple{Float64,Float64}
  α :: Float64

  x̃ :: Vector{Float64}
  ỹ :: Vector{Float64}

  x :: Vector{Float64}
  y :: Vector{Float64}

end

function BasicBody(x::Vector{T},y::Vector{T};closuretype::Type{<:BodyClosureType}=ClosedBody) where {T <: Real}
    @assert length(x) == length(y)
    BasicBody{length(x),closuretype}((0.0,0.0),0.0,x,y,x,y)
end

function Base.show(io::IO, body::BasicBody{N,C}) where {N,C}
    println(io, "Basic pointwise-specified body with $N points")
    println(io, "   Current position: ($(body.cent[1]),$(body.cent[2]))")
    println(io, "   Current angle (rad): $(body.α)")
end

#### Ellipses and circles ####

"""
    Ellipse(a,b,n) <: Body

Construct an elliptical body with semi-major axis `a` and semi-minor axis `b`,
with `n` points distributed on the body perimeter.
"""
mutable struct Ellipse{N} <: Body{N,ClosedBody}
  a :: Float64
  b :: Float64
  cent :: Tuple{Float64,Float64}
  α :: Float64

  x̃ :: Vector{Float64}
  ỹ :: Vector{Float64}

  x :: Vector{Float64}
  y :: Vector{Float64}

end

function Ellipse(a::Real,b::Real,N::Int)
    x̃ = zeros(N)
    ỹ = zeros(N)
    θ = range(0,stop=2π,length=N+1)
    @. x̃ = a*cos(θ[1:N])
    @. ỹ = b*sin(θ[1:N])

    Ellipse{N}(a,b,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ)
end

#Ellipse(a::Real,b::Real,targetsize::Float64;kwargs...) =
#    Ellipse(a,b,_adjustnumber(targetsize,Ellipse,a,b);kwargs...)

"""
    Circle(a,n) <: Body

Construct a circular body with radius `a`
and with `n` points distributed on the body perimeter.
"""
Circle(a::Real,N::Int) = Ellipse(a,a,N)

"""
    Circle(a,targetsize::Float64) <: Body

Construct a circular body with radius `a` with spacing between points set
approximately to `targetsize`.
"""
Circle(a::Real,targetsize::Float64) = Ellipse(a,a,targetsize)

function Base.show(io::IO, body::Ellipse{N}) where {N}
    if body.a == body.b
      println(io, "Circular body with $N points and radius $(body.a)")
    else
      println(io, "Elliptical body with $N points and semi-axes ($(body.a),$(body.b))")
    end
    println(io, "   Current position: ($(body.cent[1]),$(body.cent[2]))")
    println(io, "   Current angle (rad): $(body.α)")
end

#### Rectangles and squares ####

"""
    Rectangle(a,b,na) <: Body

Construct a rectangular body with x̃ side half-length `a` and ỹ side half-length `b`,
with `na` points distributed on the x̃ side (including both corners). The centroid
of the rectangle is placed at the origin (so that the lower left corner is at (-a,-b)).

By default, points are placed at the corners, and the normal vectors are bisectors
between the normals on the adjacent two sides. If the `shifted=true` flag is added,
then the points are shifted counterclockwise by half a segment, so no points
are placed on the corners, and all normals are perpendicular to the sides.
"""
mutable struct Rectangle{N,PS} <: Body{N,ClosedBody}
  a :: Float64
  b :: Float64
  cent :: Tuple{Float64,Float64}
  α :: Float64

  x̃ :: Vector{Float64}
  ỹ :: Vector{Float64}

  x :: Vector{Float64}
  y :: Vector{Float64}

  x̃mid :: Union{Vector{Float64},Nothing}
  ỹmid :: Union{Vector{Float64},Nothing}

  xmid :: Union{Vector{Float64},Nothing}
  ymid :: Union{Vector{Float64},Nothing}

end

Rectangle(a::Real,b::Real,na::Int;shifted=false) = _rectangle(a,b,na,Val(shifted))

function _rectangle(a::Real,b::Real,na::Int,::Val{false})
    x̃, ỹ = _rectangle_points(a::Real,b::Real,na::Int)
    Rectangle{length(x̃),Unshifted}(a,b,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ,nothing,nothing,nothing,nothing)
end

function _rectangle(a::Real,b::Real,na::Int,::Val{true})
    x̃mid, ỹmid = _rectangle_points(a::Real,b::Real,na::Int)
    x̃, ỹ = _midpoints(x̃mid,ỹmid,ClosedBody)
    Rectangle{length(x̃),Shifted}(a,b,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ,x̃mid,ỹmid,x̃mid,ỹmid)
end

function _rectangle_points(a::Real,b::Real,na::Int)
  Δsa = 2a/(na-1)
  nb = ceil(Int,2b/Δsa)+1
  Δsb = 2b/(nb-1)

  N = 2(na-1)+2(nb-1)
  x = zeros(N)
  y = zeros(N)

  @. x[1:na-1] = -a + Δsa*(0:na-2)
  @. y[1:na-1] = -b

  @. x[na:na+nb-2] =  a
  @. y[na:na+nb-2] = -b + Δsb*(0:nb-2)

  @. x[(na+nb-1):(2na+nb-3)] = -a + Δsa*(na-1:-1:1)
  @. y[(na+nb-1):(2na+nb-3)] = b

  @. x[(2na+nb-2):(2na+2nb-4)] = -a
  @. y[(2na+nb-2):(2na+2nb-4)] =  -b + Δsb*(nb-1:-1:1)

  return x, y

end

#Rectangle(a::Real,b::Real,targetsize::Float64;kwargs...) =
#    Rectangle(a,b,_adjustnumber(targetsize,Rectangle,a,b);kwargs...)

centraldiff(b::Rectangle{N,Shifted}) where {N} = _diff(b.xmid,b.ymid,ClosedBody)

function (T::RigidTransform)(b::Rectangle{N,Shifted}) where {N}
  b.xmid, b.ymid = T(b.x̃mid,b.ỹmid)
  b.x, b.y = T(b.x̃,b.ỹ)
  b.α = T.α
  b.cent = T.trans
  return b
end

"""
    Square(a,na) <: Body

Construct a square body with side half-length `a`
and with `na` points distributed on each side (including both corners).
"""
Square(a::Real,na::Int;kwargs...) = Rectangle(a,a,na;kwargs...)

Square(a::Real,targetsize::Float64;kwargs...) = Rectangle(a,a,targetsize;kwargs...)

function Base.show(io::IO, body::Rectangle{N}) where {N}
    if body.a == body.b
      println(io, "Square body with $N points and side half-length $(body.a)")
    else
      println(io, "Rectangular body with $N points and half-lengths ($(body.a),$(body.b))")
    end
    println(io, "   Current position: ($(body.cent[1]),$(body.cent[2]))")
    println(io, "   Current angle (rad): $(body.α)")
end

#### Plates ####

"""
    Plate(length,thick,n,[λ=1.0]) <: Body

Construct a flat plate with length `length` and thickness `thick`,
with `n` points distributed on the body perimeter.

The optional parameter `λ` distributes the points differently. Values between `0.0`
and `1.0` are accepted.

The constructor `Plate(length,n,[λ=1.0])` creates a plate of zero thickness.

Alternatively, either form can be specified with a target spacing in place of `n`.
"""
mutable struct Plate{N,C<:BodyClosureType} <: Body{N,C}
  len :: Float64
  thick :: Float64
  cent :: Tuple{Float64,Float64}
  α :: Float64

  x̃ :: Vector{Float64}
  ỹ :: Vector{Float64}

  x :: Vector{Float64}
  y :: Vector{Float64}

end

function Plate(len::Real,N::Int;λ::Float64=1.0)

    # set up points on plate
    #x = [[len*(-0.5 + 1.0*(i-1)/(N-1)),0.0] for i=1:N]

    Δϕ = π/(N-1)
    Jϕa = [sqrt(sin(ϕ)^2+λ^2*cos(ϕ)^2) for ϕ in range(π-Δϕ/2,stop=Δϕ/2,length=N-1)]
    Jϕ = len*Jϕa/Δϕ/sum(Jϕa)
    x̃ = -0.5*len .+ Δϕ*cumsum([0.0; Jϕ])
    ỹ = zero(x̃)

    Plate{N,OpenBody}(len,0.0,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ)

end

#Plate(a::Real,targetsize::Float64;kwargs...) =
#    Plate(a,_adjustnumber(targetsize,Plate,a);kwargs...)

function Plate(len::Real,thick::Real,N::Int;λ::Float64=1.0)
    # input N is the number of panels on one side only

    # set up points on flat sides
    Δϕ = π/N
    Jϕa = [sqrt(sin(ϕ)^2+λ^2*cos(ϕ)^2) for ϕ in range(π-Δϕ/2,stop=Δϕ/2,length=N)]
    Jϕ = len*Jϕa/Δϕ/sum(Jϕa)
    xtopface = -0.5*len .+ Δϕ*cumsum([0.0; Jϕ])
    xtop = 0.5*(xtopface[1:N] + xtopface[2:N+1])

    Δsₑ = Δϕ*Jϕ[1]
    Nₑ = 2*floor(Int,0.25*π*thick/Δsₑ)
    xedgeface = [0.5*len + 0.5*thick*cos(ϕ) for ϕ in range(π/2,stop=-π/2,length=Nₑ+1)]
    yedgeface = [          0.5*thick*sin(ϕ) for ϕ in range(π/2,stop=-π/2,length=Nₑ+1)]
    xedge = 0.5*(xedgeface[1:Nₑ]+xedgeface[2:Nₑ+1])
    yedge = 0.5*(yedgeface[1:Nₑ]+yedgeface[2:Nₑ+1])

    x̃ = Float64[]
    ỹ = Float64[]
    for xi in xtop
      push!(x̃,xi)
      push!(ỹ,0.5*thick)
    end
    for i = 1:Nₑ
      push!(x̃,xedge[i])
      push!(ỹ,yedge[i])
    end
    for xi in reverse(xtop,dims=1)
      push!(x̃,xi)
      push!(ỹ,-0.5*thick)
    end
    for i = Nₑ:-1:1
      push!(x̃,-xedge[i])
      push!(ỹ,yedge[i])
    end

    Plate{length(x̃),ClosedBody}(len,thick,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ)

end

#Plate(a::Real,b::Real,targetsize::Float64;kwargs...) =
#    Plate(a,b,_adjustnumber(targetsize,Plate,a,b);kwargs...)

function Base.show(io::IO, body::Plate{N}) where {N}
    println(io, "Plate with $N points and length $(body.len) and thickness $(body.thick)")
    println(io, "   Current position: ($(body.cent[1]),$(body.cent[2]))")
    println(io, "   Current angle (rad): $(body.α)")
end

#### Splined body ####

"""
    SplinedBody(X,Δx[,closuretype=ClosedBody]) -> BasicBody

Using control points in `X` (assumed to be N x 2, where N is the number of points), create a set of points
that are uniformly spaced (with spacing `Δx`) on a curve that passes through the control points. A cubic
parametric spline algorithm is used. If the optional parameter `closuretype` is set
to `OpenBody`, then the end points are not joined together.
"""
function SplinedBody(Xpts_raw::Array{Float64,2},Δx::Float64;closuretype::Type{<:BodyClosureType}=ClosedBody)
    # Assume Xpts are in the form N x 2
    Xpts = copy(Xpts_raw)
    if Xpts[1,:] != Xpts[end,:]
        Xpts = vcat(Xpts,Xpts[1,:]')
    end

    spl = (closuretype == ClosedBody) ? ParametricSpline(Xpts',periodic=true) : ParametricSpline(Xpts')
    tfine = range(0,1,length=1001)
    dX = derivative(spl,tfine)

    np = ceil(Int,sqrt(sum(dX.^2)*(tfine[2]-tfine[1]))/Δx)

    tsamp = range(0,1,length=np)
    x = [X[1] for X in spl.(tsamp[1:end-1])]
    y = [X[2] for X in spl.(tsamp[1:end-1])]

    return BasicBody(x,y,closuretype=closuretype)
end

#### NACA 4-digit airfoil ####

"""
    NACA4(cam,pos,thick,np,[len=1.0]) <: Body{N}

Generates points in the shape of a NACA 4-digit airfoil of chord length 1. The
relative camber is specified by `cam`, the position of
maximum camber (as fraction of chord) by `pos`, and the relative thickness
by `thick`. The parameter `np` specifies the number of points on the upper
or lower surface. The optional parameter `len` specifies the chord length,
which defaults to 1.0.

# Example

```jldoctest
julia> w = Bodies.NACA4(0.0,0.0,0.12);
```
"""
mutable struct NACA4{N} <: Body{N,ClosedBody}
  len :: Float64
  camber :: Float64
  pos :: Float64
  thick :: Float64

  cent :: Tuple{Float64,Float64}
  α :: Float64

  x̃ :: Vector{Float64}
  ỹ :: Vector{Float64}

  x :: Vector{Float64}
  y :: Vector{Float64}

end

function NACA4(cam::Real,pos::Real,t::Real,np::Int;len=1.0)

# Here, cam is the fractional camber, pos is the fractional chordwise position
# of max camber, and t is the fractional thickness.

npan = 2*np-2

# Trailing edge bunching
an = 1.5
anp = an+1
x = zeros(np)

θ = zeros(size(x))
yc = zeros(size(x))

for j = 1:np
    frac = Float64((j-1)/(np-1))
    x[j] = 1 - anp*frac*(1-frac)^an-(1-frac)^anp;
    if x[j] < pos
        yc[j] = cam/pos^2*(2*pos*x[j]-x[j]^2)
        if pos > 0
            θ[j] = atan(2*cam/pos*(1-x[j]/pos))
        end
    else
        yc[j] = cam/(1-pos)^2*((1-2*pos)+2*pos*x[j]-x[j]^2)
        if pos > 0
            θ[j] = atan(2*cam*pos/(1-pos)^2*(1-x[j]/pos))
        end
    end
end

xu = zeros(size(x))
yu = xu
xl = xu
yl = yu

yt = t/0.20*(0.29690*sqrt.(x)-0.12600*x-0.35160*x.^2+0.28430*x.^3-0.10150*x.^4)

xu = x-yt.*sin.(θ)
yu = yc+yt.*cos.(θ)

xl = x+yt.*sin.(θ)
yl = yc-yt.*cos.(θ)

rt = 1.1019*t^2;
θ0 = 0
if pos > 0
    θ0 = atan(2*cam/pos)
end
# Center of leading edge radius
xrc = rt*cos(θ0)
yrc = rt*sin(θ0)
θle = collect(0:π/50:2π)
xlec = xrc .+ rt*cos.(θle)
ylec = yrc .+ rt*sin.(θle)

# Assemble data
coords = [xu yu xl yl x yc]
cole = [xlec ylec]

# Close the trailing edge
xpanold = [0.5*(xl[np]+xu[np]); reverse(xl[2:np-1],dims=1); xu[1:np-1]]
ypanold = [0.5*(yl[np]+yu[np]); reverse(yl[2:np-1],dims=1); yu[1:np-1]]

xpan = zeros(npan)
ypan = zeros(npan)
for ipan = 1:npan
    if ipan < npan
        xpan1 = xpanold[ipan]
        ypan1 = ypanold[ipan]
        xpan2 = xpanold[ipan+1]
        ypan2 = ypanold[ipan+1]
    else
        xpan1 = xpanold[npan]
        ypan1 = ypanold[npan]
        xpan2 = xpanold[1]
        ypan2 = ypanold[1]
    end
    xpan[ipan] = 0.5*(xpan1+xpan2)
    ypan[ipan] = 0.5*(ypan1+ypan2)
end
w = ComplexF64[1;reverse(xpan,dims=1)+im*reverse(ypan,dims=1)]*len
w .-= mean(w)

x̃ = real.(w)
ỹ = imag.(w)

NACA4{length(x̃)}(len,cam,pos,t,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ)

end

#NACA4(a::Real,b::Real,c::Real,targetsize::Float64;kwargs...) =
#    NACA4(a,b,c,_adjustnumber(targetsize,NACA4,a,b,c);kwargs...)

function Base.show(io::IO, body::NACA4{N}) where {N}
    println(io, "NACA 4-digit airfoil with $N points and length $(body.len) and thickness $(body.thick)")
    println(io, "   Current position: ($(body.cent[1]),$(body.cent[2]))")
    println(io, "   Current angle (rad): $(body.α)")
end

####

function _adjustnumber(targetsize::Real,shapefcn::Type{T},params...;kwargs...) where {T <: Body}
    ntrial = 501
    return floor(Int,ntrial*mean(dlength(shapefcn(params...,ntrial;kwargs...)))/targetsize)
end

for shape in (:Ellipse,:Rectangle,:Plate,:NACA4)
    @eval RigidBodyTools.$shape(params...;kwargs...) =
        $shape(Base.front(params)...,
          _adjustnumber(Base.last(params),$shape,Base.front(params)...;kwargs...);kwargs...)
end


text:
export print_merge_logs, print_clique_sizes

"""
    AbstractMergeStrategy

An abstract supertype for merge strategies.

A merge strategy determines how the cliques in a clique tree / clique graph are merged to improve computation time of the projection. Each merge strategy should
implement the following functions:

 - initialise!: Initialise the graph / tree from the input clique graph / tree, allocate memory
 - traverse: A method that determines how the clique tree / graph is traversed
 - evaluate: A method that decides whether to merge two cliques
 - update!: A method to update local strategy-related information after a merge
"""
abstract type AbstractMergeStrategy end

"""
    AbstractTreeBasedMerge

An abstract supertype for merge strategies that work with the clique tree.
"""
abstract type AbstractTreeBasedMerge <: AbstractMergeStrategy end

"""
    AbstractGraphBasedMerge

An abstract supertype for merge strategies that work with the clique (intersection) graph.
"""
abstract type AbstractGraphBasedMerge <: AbstractMergeStrategy end

"""
		MergeLog

A struct to analyse the clique merges. Introduced for debugging purposes.
"""
mutable struct MergeLog
	num::Int64 # number of merges
	clique_pairs::Array{Int64, 2} # ordered pair merges
	decisions::Array{Bool, 1} # at what step was merged
	function MergeLog()
		new(0, zeros(0, 2), Array{Bool}(undef, 0))
	end
end

"""
	SuperNodeTree

A structure to represent and analyse the sparsity pattern of the input matrix L.

### Note:
Based on the `merge_strategy` in the constructor, SuperNodeTree might be initialised
as a graph, i.e. seperators `sep` are left empty as well as `snd_child` and `snd_post`.

After cliques have been merged, a valid clique tree will be recomputed from the
consolidated clique graph.
"""
mutable struct SuperNodeTree
	snd::Union{Array{Set{Int64}, 1}, Array{Array{Int64, 1}, 1}} #vertices of supernodes stored in one array (also called residuals)
	snd_par::Array{Int64,1}  # parent of supernode k is supernode j=snd_par[k]
	snd_post::Array{Int64,1} # post order of supernodal elimination tree
	snd_child::Array{Set{Int64},1}
	post::Array{Int64} # post ordering of the vertices in elim tree σ(j) = v
	par::Array{Int64}
	sep::Union{Array{Set{Int64},1}, Array{Array{Int64, 1}, 1}} #vertices of clique seperators
	nBlk::Array{Int64,1} #sizes of submatrizes defined by each clique, sorted by post-ordering, e.g. size of clique with order 3 => nBlk[3]
	num::Int64 # number of supernodes / cliques in tree
	merge_log::MergeLog
	strategy::AbstractMergeStrategy
	function SuperNodeTree(L::SparseMatrixCSC,  merge_strategy::AbstractMergeStrategy)
		par = etree(L)
		child = child_from_par(par)
		post = post_order(par, child)
		# sorting of children from highest order one to lowest make sure that a node always gets
		# added to the supernode with the highest rep vertex
		#sort_children!(child, post)
		# faster algorithm according to Vandenberghe p.308
		degrees = higher_degrees(L)
		snd, snd_par = find_supernodes(par, post, degrees, merge_strategy)

		snd_child = child_from_par(snd_par)
	 	snd_post = post_order(snd_par, snd_child)

		# If the merge strategy is clique graph-based, we give up the tree structure and add the seperators to the supernodes
		# the supernodes then represent the full clique
		# after the clique merging a new clique tree will be computed before psd completion is performed
	 	if typeof(merge_strategy) <: AbstractGraphBasedMerge
			sep = [Set{Int64}() for i = 1:length(snd)]
			add_separators!(L, snd, sep, snd_par, post)
			@. snd_par = -1
			snd_child = [Set{Int64}() for i = 1:length(snd)]
			new(snd, snd_par, snd_post, snd_child, post, par, sep, [0], length(snd_post), MergeLog(), merge_strategy)
	 	# If the merge strategy is tree based keep the supernodes and separators in to different locations
		else
	 		# given the supernodes (clique residuals) find the clique separators
			sep = find_separators(L, snd, snd_par, post, merge_strategy)
			new(snd, snd_par, snd_post, snd_child, post, par, sep, [0], length(snd_post), MergeLog(), merge_strategy)

		end
	end

	function COSMO.SuperNodeTree(cliques::Array{Set{Int64},1}, N::Int64)
	    snd_child = [Set{Int64}() for i = 1:N]
	    snd_par = -1 * ones(Int64, length(cliques))
	    sep = [Set{Int64}() for i = 1:length(cliques)]
	    snd_post = zeros(Int64, length(cliques))
	    new(cliques, snd_par, snd_post, snd_child, collect(1:N), [0], sep, [1], length(cliques), COSMO.MergeLog(), COSMO.CliqueGraphMerge())
	end

	# FIXME: only for debugging purposes
	function SuperNodeTree(snd, par, snd_post, sep, merge_strategy; post::Array{Int64, 1} = [1])
		child = child_from_par(par)
  	new(snd, par, snd_post, child, post, [1], sep, [1], length(snd), MergeLog(), merge_strategy)
	end
end

function sort_children!(child, post)
	for v = 1:length(child)
		child[v] = sort(child[v], by = x -> post[x], rev = true)
	end
end
# -------------------------------------
# FUNCTION DEFINITIONS
# -------------------------------------
# given v=σ^-1(i) it returns i=σ(v)
function invert_order(sigma::Array{Int64,1})
	sigma_inv = zeros(Int64, length(sigma))
	for iii=1:length(sigma)
		sigma_inv[sigma[iii]] = iii
	end
	return sigma_inv
end

# elimination tree algorithm from H.Liu - A Compact Row Storage Scheme for Cholesky Factors Using Elimination Trees
# function etree_liu(g)
# 	N = length(g.adjacency_list)
# 	par = zeros(Int64,N)
# 	ancestor = zeros(Int64,N)

# 	elemSequence = g.reverse_ordering[collect(1:N)]
# 	for iii in elemSequence
# 		for vk in g.adjacency_list[iii]
# 			if g.ordering[vk] < g.ordering[iii]
# 				r = vk
# 				while (ancestor[r] != 0) && (ancestor[r] != iii)
# 					t = ancestor[r]
# 					ancestor[r] = iii
# 					r = t
# 				end
# 				if ancestor[r] == 0
# 					ancestor[r] = iii
# 					par[r] = iii
# 				end
# 			end
# 		end
# 	end

# 	return par
# end

# simplified version of my own elimination tree algorithm with simplified data structure (fastest)
function etree(L)
	N = size(L, 1)
	par = zeros(Int64, N)
	# loop over Vertices of graph
	for i=1:N
		par_ = find_parent_direct(L, i)
		par[i] = par_
	end
	return par
end

# perform a depth-first-search to determine the post order of the tree defined by parent and children vectors
# FIXME: This can be made a lot faster for the case that merges happened, i.e. Nc != length(par)

function post_order(par::Array{Int64,1}, child::Union{Array{Array{Int64, 1}, 1}, Array{Set{Int64}, 1}}, Nc::Int64)

	order = (Nc + 1) * ones(Int64, length(par))
	root = findall(x ->x == 0, par)[1]
	stack = Array{Int64}(undef, 0)
	iii = Nc
	push!(stack, root)
	while !isempty(stack)
		v = pop!(stack)
		order[v] = iii
		iii -= 1
		push!(stack, child[v]...)
	end
	post = collect(1:length(par))
	sort!(post, by = x-> order[x])

	# if merges happened, remove the entries pointing to empty arrays / cliques
	Nc != length(par) && resize!(post, Nc)
	return post
end
post_order(par::Array{Int64,1}, child::Union{Array{Array{Int64, 1}, 1}, Array{Set{Int64},1}}) = post_order(par, child, length(par))

function child_from_par(par::Array{Int64,1})

	child = [Set{Int64}() for i = 1:length(par)]
	for i = 1:length(par)
		par_ = par[i]
		par_ != 0 && push!(child[par_], i)
	end
	return child
end

# ------------------------------------------
# Interfaces function to SuperNodeTree
# ------------------------------------------

function num_cliques(sntree::SuperNodeTree)
	return sntree.num::Int64
end

function get_post_order(sntree::SuperNodeTree)
	return sntree.snd_post
end

function get_post_order(sntree::SuperNodeTree, i::Int64)
	return sntree.snd_post[i]
end
# Using the post order ensures that no empty arrays from the clique merging are returned
function get_snd(sntree::SuperNodeTree, ind::Int64)
		return sntree.snd[sntree.snd_post[ind]]
end

function get_sep(sntree::SuperNodeTree, ind::Int64)
		return sntree.sep[sntree.snd_post[ind]]
end

function get_clique_par(sntree::SuperNodeTree, clique_ind::Int64)
		return sntree.snd_par[sntree.snd_post[clique_ind]]
end
# the block sizes are stored in post order, e.g. if clique 4 (stored in pos 4) has order 2, then nBlk[2] represents the cardinality of clique 4
function get_nBlk(sntree::SuperNodeTree, ind::Int64)
		return sntree.nBlk[ind]::Int64
end

function get_overlap(sntree::SuperNodeTree, ind::Int64)
		return length(sntree.sep[sntree.snd_post[ind]])
end

"Returns the number of rows and the number of overlaps of all the blocks (cliques) represented in the tree after decomposition."
function get_decomposed_dim(sntree::SuperNodeTree, C::DecomposableCones{<: Real})
	dim = 0
	overlaps = 0
	for iii = 1:num_cliques(sntree)
		dim += vec_dim(get_nBlk(sntree, iii), C)
		overlaps += vec_dim(get_overlap(sntree, iii), C)
	end
	return dim::Int64, overlaps::Int64
end

"Given the side dimension of a PSD cone return the number of stored entries."
vec_dim(side_dim::Int64, C::PsdCone{<:Real}) = Base.power_by_squaring(side_dim, 2)
vec_dim(side_dim::Int64, C::PsdConeTriangle{<:Real}) = div(side_dim * (side_dim + 1), 2)

" Return clique with post order `ind` (prevents returning empty arrays due to clique merging)"
get_clique(sntree::SuperNodeTree, ind::Int64) = get_clique(sntree, ind, sntree.strategy)
function get_clique(sntree::SuperNodeTree, ind::Int64, strategy::AbstractTreeBasedMerge)
	c = sntree.snd_post[ind]
	return union(sntree.snd[c], sntree.sep[c])
end

function get_clique(sntree::SuperNodeTree, ind::Int64, strategy::AbstractGraphBasedMerge)
	# if a clique tree has been recomputed, call the method for AbstractTreeBasedMerge types
	if strategy.clique_tree_recomputed
		return get_clique(sntree, ind, COSMO.NoMerge())
	else
		c = sntree.snd_post[ind]
		return sntree.snd[c]
	end
end

function get_clique_by_ind(sntree::SuperNodeTree, ind::Int64)
	return union(sntree.snd[ind], sntree.sep[ind])
end

function print_cliques(sp; reordered = true)
	sntree = sp.sntree
	ordering = sp.ordering
	Nsnd = length(sntree.snd)
	println("Cliques of Graph:")
	println("Reordered = $(reordered)")
	for iii = 1:Nsnd
			if !reordered
				snd = map(x-> ordering[x], sntree.snd[iii])
				sep = map(x-> ordering[x], sntree.sep[iii])
			else
				snd = sntree.snd[iii]
				sep = sntree.sep[iii]
			end
			println("$(iii): res: $(snd), sep: $(sep)")
	end
end
print_cliques(sntree::SuperNodeTree) = print_cliques(sntree, sntree.strategy)
function print_cliques(sntree::SuperNodeTree, strategy::AbstractTreeBasedMerge)
	Nsnd = length(sntree.snd)
	println("Cliques of Graph:")
	for iii = 1:Nsnd
			println("$(iii): res: $(sntree.snd[iii]), sep: $(sntree.sep[iii])")
	end
end

function print_clique_orders(sntree::SuperNodeTree, strategy::AbstractTreeBasedMerge)
	Nsnd = length(sntree.snd)
	println("Cliques of Graph:")
	for iii = 1:Nsnd
		post = sntree.post
		postInv = invert_order(post) #σ-1(v) = j

  	snd_ip = map(x -> postInv[x], sntree.snd[iii])
  	sep_ip = map(x -> postInv[x], sntree.sep[iii])
		println("$(iii): res: $(snd_ip), sep: $(sep_ip)")
	end
end

function print_cliques(sntree::SuperNodeTree, strategy::AbstractGraphBasedMerge)

	if strategy.clique_tree_recomputed
		return print_cliques(sntree, COSMO.NoMerge())
	else
		Nsnd = length(sntree.snd)
		println("Cliques of Graph:")
		for iii = 1:Nsnd
				println("$(iii): clique: $(sntree.snd[iii])")
		end
	end
	return nothing
end

function print_clique_sizes(ws)
  sp_arr = ws.ci.sp_arr
  # print the merge logs for each sparsity pattern
  println(">>> Clique Dimensions:")
  for (iii, sp) in enumerate(sp_arr)
    println("Sparsity Pattern Nr. $(iii), Graph Size: $(length(sp.sntree.par))")

    t = sp.sntree
    Nc = length(t.snd_post)
    # block sizes
    sizes = zeros(Int64, Nc)
    # occurences of block size
    occ = zeros(Int64, Nc)
    for jjj = 1:Nc
      c = t.snd_post[jjj]
      dim = length(t.snd[c]) + length(t.sep[c])
      # try to check if that dimension occured before, if not add, otherwise increment occ
      ind = findfirst(x -> x == dim, sizes)

      if ind != nothing
        occ[ind] = occ[ind] + 1
      else
        sizes[jjj] = dim
        occ[jjj] = 1
      end
    end

    # consolidate
    filter!(x -> x != 0, occ)
    filter!(x -> x != 0, sizes)

    # sort sizes
    p = sortperm(sizes)
    sizes = sizes[p]
    occ = occ[p]

    for (jjj, dim) in enumerate(sizes)
      println("$(occ[jjj])x dim:$(dim)")
    end
  end
end

function check_degree_condition(v::Int64, w::Int64, degrees::Array{Int64,1})
	return degrees[v] > degrees[w] - 1
end

# Algorithm from A. Poten and C. Sun: Compact Clique Tree Data Structures in Sparse Matrix Factorizations (1989)
function pothen_sun(par::Array{Int64,1}, post::Array{Int64,1}, degrees::Array{Int64,1})
	N = length(par)
	sn_ind = -1 * ones(Int64, N) # if snInd[v] < 0 then v is a rep vertex, otherwise v ∈ supernode[snInd[v]]
	supernode_par = -1 * ones(Int64, N)
	children = 	[Array{Int64}(undef, 0) for i = 1:length(par)]

	root_ind = findfirst(x -> x == 0, par)
	# go through parents of vertices in post_order
	for v in post

		if par[v] == 0
			push!(children[root_ind], v)
		else
			push!(children[par[v]], v)
		end

		# parent is not the root
		if par[v] != 0
			if degrees[v] - 1 == degrees[par[v]] && sn_ind[par[v]] == -1
				# Case A: v is a representative vertex
				if sn_ind[v] < 0
					sn_ind[par[v]] = v
					sn_ind[v] -= 1
				# Case B: v is not representative vertex, add to sn_ind[v] instead
				else
					sn_ind[par[v]] = sn_ind[v]
					sn_ind[sn_ind[v]] -= 1
				end
			else
				if sn_ind[v] < 0
					supernode_par[v] = v
				else
					supernode_par[sn_ind[v]] = sn_ind[v]
				end
			end
		end

		# k: rep vertex of the snd that v belongs to
		if sn_ind[v] < 0
			k = v
		else
			k = sn_ind[v]
		end
		# loop over v's children
		v_children = children[v]
		if !isempty(v_children)
			for w in v_children
				if sn_ind[w] < 0
					l = w
				else
					l = sn_ind[w]
				end
				if l != k
					supernode_par[l] = k
				end
			end
		end
	end # loop over vertices

	# representative vertices
	reprv = findall(x-> x < 0, sn_ind)
	# vertices that are the parent of representative vertices
	repr_par = supernode_par[reprv]
	# take into account that all non-representative arrays are removed from the parent structure
	sn_par = zeros(Int64, length(reprv))

	for (iii, rp) in enumerate(repr_par)
		ind = findfirst(x -> x == rp, reprv)
		ind == nothing && (ind = 0)
		sn_par[iii] = ind
	end

	return sn_par, sn_ind
end

function initialise_sets(N::Int64, strategy::AbstractGraphBasedMerge)
	return [Set{Int64}() for i = 1:N]
end

function initialise_sets(N::Int64, strategy::AbstractMergeStrategy)
	return [Int64[] for i = 1:N]
end

function find_supernodes(par::Array{Int64,1}, post::Array{Int64,1}, degrees::Array{Int64,1}, strategy::AbstractMergeStrategy)
	supernode_par, snInd = pothen_sun(par, post, degrees)
	# number of vertices
	N = length(par)
	# number of representative vertices == number of supernodes
	Nrep = length(supernode_par)
	# snode = [Set{Int64}() for i = 1:N]
	snode = initialise_sets(N, strategy)
	for iii = 1:N
		f = snInd[iii]
		if f < 0
			push!(snode[iii], iii)
		else
			push!(snode[f], iii)
		end
	end
	filter!(x -> !isempty(x), snode)
	return snode, supernode_par

end

function find_separators(L, snodes::Union{Array{Set{Int64}, 1}, Array{Array{Int64, 1},1}}, supernode_par::Array{Int64,1}, post::Array{Int64,1}, strategy::AbstractMergeStrategy)
	postInv = invert_order(post)

	Nc = length(supernode_par)
	sep = initialise_sets(Nc, strategy)

	for iii = 1:Nc
		vRep = snodes[iii][1]

		adjPlus = find_higher_order_neighbors(L, vRep)
		deg = length(adjPlus) + 1
		sep[iii] = adjPlus
		setdiff!(sep[iii], snodes[iii])

	end

	return sep

end

function add_separators!(L::SparseMatrixCSC, snodes::Array{Set{Int64},1}, separators::Array{Set{Int64}, 1}, supernode_par::Array{Int64,1}, post::Array{Int64,1})
	postInv = invert_order(post)

	Nc = length(supernode_par)

	for iii = 1:Nc
		snode = snodes[iii]
		sep = separators[iii]
		vRep = minimum(snode)

		adjPlus = find_higher_order_neighbors(L, vRep)
		for neighbor in adjPlus
			if !in(neighbor, snode)
				push!(snode, neighbor)
				push!(sep, neighbor)
			end
		end
	end
	return nothing
end

"""
		reorder_snd_consecutively!(sntree, ordering)

Takes a SuperNodeTree and reorders the vertices in each supernode (and separator) to have consecutive order.

The reordering is needed to achieve equal column structure for the psd completion of the dual variable `Y`. This also modifies `ordering` which maps the vertices in the `sntree` back to the actual location in the not reordered data, i.e.
the primal constraint variable `S` and dual variables `Y`.
"""
function reorder_snd_consecutively!(t::SuperNodeTree, ordering::Array{Int64, 1})

	# determine permutation vector p and permute the vertices in each snd
	p = zeros(Int64,length(t.post))
	snd = t.snd
	sep = t.sep

	k = 1
	for snd_ind in t.snd_post
	  s = snd[snd_ind]
	  l = length(s)
	  p[k:k+l-1] = sort(collect(s))
	  s[:] = collect(k:1:k+l-1)
	  k += l
	end

  # permute the separators as well
  p_inv = invperm(p)
  for sp in t.sep
    sp[:] = map(x -> p_inv[x], sp)
  end

	# # permute the degrees, supernodal parent structure and supernodal post order

	# snd_post_inv = invperm(t.snd_post)
	# root_ind = findfirst(x -> x == 0, t.snd_par)
	# t.snd_par[root_ind] = root_ind
	# t.snd_par = sndpost_inv[t.snd_par[t.snd_post]]
	# t.snd_post = collect(1:1:t.num)

	permute!(ordering, p)
	return nothing
end

function find_higher_order_neighbors(L::SparseMatrixCSC, v::Int64)
	v == size(L, 1) && return 0
	col_ptr = L.colptr
	row_val = L.rowval
	return row_val[col_ptr[v]:col_ptr[v + 1] - 1]
end

function find_parent_direct(L::SparseMatrixCSC, v::Int64)
	v == size(L, 1) && return 0
	col_ptr = L.colptr
	row_val = L.rowval
	return row_val[col_ptr[v]]
end

# findall the cardinality of adj+(v) for all v in V
function higher_degrees(L::SparseMatrixCSC)
	N = size(L, 1)
	degrees = zeros(Int64, N)
	col_ptr = L.colptr
	row_val = L.rowval

	for v = 1:N-1
		degrees[v] = col_ptr[v + 1] - col_ptr[v]
	end
	return degrees
end

# -------------------------------------
# FUNCTION DEFINITIONS
# -------------------------------------

# this assumes a sparse lower triangular matrix L
function connect_graph!(L::SparseMatrixCSC)
	# unconnected blocks don't have any entries below the diagonal in their right-most column
	m = size(L, 1)
	row_val = L.rowval
	col_ptr = L.colptr
	for j = 1:m-1
		connected = false
		for k in col_ptr[j]:col_ptr[j+1]-1
			if row_val[k] > j
				connected  = true
				break
			end
		end
		if !connected
			L[j+1, j] = 1
		end
	end
end

"""
	find_graph!(ci, rows::Array{Int64, 1}, N::Int64, C::AbstractConvexSet)

Given the indices of non-zero rows in `rows`:
- Compute the sparsity pattern and find a chordal extension using `QDLDL` with AMD ordering `F.perm`.
- If unconnected, connect the graph represented by the cholesky factor `L`
"""
function find_graph!(ci, rows::Array{Int64, 1}, N::Int64, C::AbstractConvexSet)
	row_val, col_val = COSMO.row_ind_to_matrix_indices(rows, N, C)
	F = QDLDL.qdldl(sparse(row_val, col_val, ones(length(row_val)), N, N), logical = true)#, perm = collect(1:N))
	nz_ind_map = get_nz_ind_map(rows, N)
	# this takes care of the case that QDLDL returns an unconnected adjacency matrix L
	connect_graph!(F.L)
	ci.L = F.L
	return F.perm, nz_ind_map
end

" A sparse vector that maps the index of svec(i, j) = k to the actual index of where that entry is stored in `A.rowval`."
function get_nz_ind_map(rows::Array{Int64, 1}, N::Int64)
	d = div(N * (N+1), 2)
	nzind = rows
	nzval = collect(1:length(rows))

	return SparseArrays._sparsevector!(nzind, nzval, d)
end

# given an array [rows] that represent the nonzero entries of a vectorized NxN matrix,
# return the rows and columns of the nonzero entries of the original matrix
function row_ind_to_matrix_indices(rows::Array{Int64,1}, N::Int64, ::PsdCone{Float64})
	row_val = zeros(Int64, length(rows))
	col_val = zeros(Int64, length(rows))
	for (ind, r) in enumerate(rows)
		_rem = mod(r, N)
		fl = fld(r, N)
		if _rem == 0
			row_val[ind] = N
			col_val[ind] = fl
		else
			row_val[ind] = _rem
			col_val[ind] = fl + 1
		end
	end
	return row_val, col_val
end

# given an array [rows] that represent the nonzero entries of the vectorized upper triangular part of a NxN matrix,
# return the rows and columns of the nonzero entries of the original matrix
function row_ind_to_matrix_indices(rows::Array{Int64,1}, N::Int64, ::COSMO.PsdConeTriangle{Float64})
	#  allocate conservative here since we don't know how many diagonal entries are contained in row
	row_val = zeros(Int64, 2 * length(rows) + N)
	col_val = zeros(Int64, 2 * length(rows) + N)
	ind = 1
	_step = 1
	for (iii, r) in enumerate(rows)
		i, j = COSMO.svec_to_mat(r)

		row_val[ind] = i
		col_val[ind] = j
		ind += 1

		if i != j
			row_val[ind] = j
			col_val[ind] = i
			ind += 1
		end
	end
	return row_val[1:ind - 1], col_val[1:ind - 1]
end

# Given a linear index find the col of the corresponding upper triangular matrix
function svec_get_col(x::Int64)
	c = (sqrt(8 * x + 1) - 1) / 2
	if c % 1 != 0
		return Int((floor(c) + 1))
	else
		return Int(c)
	end
end

# Given a linear index find the row of the corresponding upper triangular matrix
function svec_get_row(x::Int64)
	c = get_col(x) - 1
	k = (c + 1) * c / 2
	return (x - k)::Int64
end

# Given a linear index find the row and column of the corresponding upper triangular matrix
function svec_to_mat(ind::Int64)
	c = svec_get_col(ind)
	k = (c - 1) * c / 2
	r = Int(ind - k)
	return r::Int64, c::Int64
end


text:
include("constants.jl")
include("bremsstrahlung.jl")

#Plank law
function plank(ν,T)
    b_nu = ( (2 * h * ν.^3) / (c^2)) .* (1 ./ ( exp.( (h .* ν) / (k_b * T) ) - 1))
    return b_nu
end

#Intensity
function intensity(N, ν, T, D, I0)
    I   = [I0]
    Iν  = 0
    tau = []

    for i=2:length(N)

        Δx  = D[i] - D[i-1]
        ttemp = τ([N[i-1] N[i]],[T[i-1] T[i]], ν, Δx)
        push!(tau, ttemp)

        I_temp = I[i-1] * exp(-τ([N[i-1] N[i]],[T[i-1] T[i]], ν, Δx)) +
                 S(ν, T[i]) * (1 - exp(-τ([N[i-1] N[i]],[T[i-1] T[i]], ν, Δx)))

        Iν = I_temp
        push!(I,I_temp)
    end

    #return I, Iν, tau
    return Iν
end

#Opacity
function κ(N,T,ν, opacity=true)
    if opacity == true
        return bremsstrahlung(N, T ,ν)
    else
        return 1.0
    end
end

#Source function
function S(ν, T)
    return plank(ν,T)
end

#Optical depth
function τ(N, T, ν, Δx)

    τ1 = (Δx/2 * 1e5) * ( κ(N[1], T[1], ν) + κ(N[2], T[2], ν) )

    return τ1
end

#Temperature brightness
function tb(I, ν)
    return I * c^2 / (2 * k_b * ν^2)
end

export plank, intensity, κ, S, τ, tb


text:
# TODO: think of a better name for this? and consider using a special container
# type (akin to ModelFrame), since we'd probably like to hold onto the schema,
# too.

################################################################################
# With Formula as a vector of Terms:
#
# entry point is to set_schema on Vector{Term}.  initialize set of terms seen
# already as a Set{Term}.  then map set_schema! on each term.
#
# for individual Term, default method is to do nothing (return the term), and
# add it to the set of terms already seen.
#
# for dealing with categorical terms, need to keep track of what's been seen so
# far, and the context in which the Term we're currently processing occurred (by
# itself, or as part of an interaction, or otherwise).  All that matters for
# determining whether an aliased term is present is the _set_ of variables
# contained in previously encountered terms, so we also need a way of extracting
# that.

function set_schema!(f::Formula, sch::Data.Schema)
    @argcheck !f.schema_set "Schema already set for this formula!"
    already = Set()
    f.term = set_schema(f.term, Set(), sch)
    f.schema_set = true
    f
end

set_schema(terms::AbstractVector{<:Terms.Term}, already::Set, sch::Data.Schema) =
    map(t->set_schema(t, already, sch), terms)

set_schema(::Nothing, ::Set, ::Data.Schema) = nothing

set_schema(t::Terms.FormulaTerm, already::Set, sch::Data.Schema) =
    Terms.FormulaTerm(set_schema(t.lhs, Set(), sch),
                      set_schema(t.rhs, already, sch))

set_schema(t::Terms.Intercept, already::Set, ::Data.Schema) = (push!(already, termsyms(t)); t)

set_schema(t::Terms.FunctionTerm, already::Set, ::Data.Schema) =
    (push!(already, termsyms(t)); t)

function set_schema(t::Terms.Interaction, already::Set, sch::Data.Schema)
    new_t = Terms.Interaction(map(s -> set_schema(s, t, already, sch), t.terms))
    push!(already, termsyms(new_t))
    new_t
end

set_schema(t::Terms.Eval, already::Set, sch::Data.Schema) = set_schema(t, t, already, sch)

function set_schema(t::Terms.Eval, context::T, already::Set, sch::Data.Schema) where T<:Terms.Term
    @debug "$t in context of $context, already seen $already"
    push!(already, termsyms(t))
    if is_categorical(t, sch)
        aliased = aliassyms(t, context)
        @debug "  aliases: $aliased"
        if aliased in already
            # TODO: allow custom contrasts here
            contr = DEFAULT_CONTRASTS()
        else
            # lower-order term that is aliased by full-rank contrasts for this
            # term is NOT present, so use full-rank contrasts and add the
            # aliased term to the set of terms we've seen
            contr = FullDummyCoding()
            push!(already, aliased)
        end
        uniq = sch.metadata[:summaries][name(t)][:unique]
        mat = ContrastsMatrix(contr, uniq)
        invindex = Dict(lev => i for (i,lev) in enumerate(mat.levels))
        N = size(mat, 2)
        Terms.Categorical{N}(name(t), mat, invindex)
    else
        Terms.Continuous(name(t))
    end
end

"""
    termsyms(t::Terms.Term)

Extract the Set of symbols referenced in this term.

This is needed in order to determine when a categorical term should have
standard (reduced rank) or full rank contrasts, based on the context it occurs
in and the other terms that have already been encountered.
"""

termsyms(t::Terms.Term) = Set()
termsyms(t::Terms.Intercept) = Set([1])
termsyms(t::Union{Terms.Eval, Terms.Categorical, Terms.Continuous}) = Set([t.name])
termsyms(t::Terms.Interaction) = mapreduce(termsyms, union, t.terms)
termsyms(t::Terms.FunctionTerm) = Set([t.name])

"""
    aliassyms(t::T, context::S) where {T<:Terms.Term, S<:Terms.Term}

Get the Set of symbols that this term potentially aliases in this context.
"""
aliassyms(::Terms.Eval, context::Terms.Eval) = Set([1])
function aliassyms(t::Terms.Eval, context::Terms.Interaction)
    mapreduce(termsyms, union, c for c in context.terms if t ≠ c)
end


text:
"""
	streamf(filename|io|formatted) -> ::FormattedStream
	streamf(f::Function, filename|io|formatted) -> f(s)

Wrap a resource as a stream of formatted objects, inferring format/coding
automatically if unspecified.

The second signature calls `f` with the opened stream as argument, then closes
the stream. This is intended to be called using `do` block syntax.

# Arguments

	filename::AbstractString
	io::IO
	formatted::Formatted

# Examples

```julia
julia> s = streamf("lysozyme.trr")
[...]

julia> s = streamf(open("lysozyme.trr"))
[...]

julia> s = streamf(specify("lysozyme.dat", "trajectory/x-trr"))
[...]
```
"""
function streamf(f::Formatted, args...; kwargs...)
	format = resolveformat(f)
	coding = resolvecoding(f)
	streamf(f.resource, format, coding, args...; kwargs...)
end

streamf(filename::AbstractString, format::MIME, coding::Union{MIME,Nothing},
		args...; kwargs...) =
		streamf(open(filename), format, coding, args...; kwargs...)

function streamf(io::IO, format::MIME, coding::Union{MIME,Nothing}, args...;
		kwargs...)
	streamer = resolvestreamer(format)
	if coding == nothing
		streamf(streamer, format, io, args...; kwargs...)
	else
		decoder = resolvedecoder(coding)
		io2 = TranscodingStream(decoder(), io)
		streamf(streamer, format, io2, args...; kwargs...)
	end
end

streamf(resource::Union{AbstractString,IO}, args...; kwargs...) =
		streamf(infer(resource), args...; kwargs...)

function streamf(f::Function,
		resource::Union{AbstractString,IO,Formats.FormattedFilename}, args...;
		kwargs...)
	s = streamf(resource, args...; kwargs...)
	results = f(s)
	close(s)
	results
end

"""
	streamf(streamer, mime, io, args...; kwargs...)

Open `io` as a streamf of objects using a specific `streamer` for the `mime`
format.

This version should be specialised for streamers that wish to handle a specific
format. Arguments `args` and `kwargs` make it possible to customize streamer
behavior.

# Arguments

	streamer::FormatHandler
	mime::MIME
	io::IO

# Example

```julia
struct MyIO <: Formats.FormatHandler end

function streamf(::MyIO, ::MIME"image/apng", io::IO)
	# Open io and wrap it to provide a stream of objects
	stream
end
```
"""
streamf

struct FormattedStreamIterator{T<:FormattedStream}
	s::T

	FormattedStreamIterator{T}(s::T) where {T<:FormattedStream} = new(s)
end

function Base.iterate(iter::FormattedStreamIterator)
	seekstart(iter.s)
	iterate(iter, nothing)
end

function Base.iterate(iter::FormattedStreamIterator, ::Nothing)
	if eof(iter.s)
		nothing
	else
		read(iter.s), nothing
	end
end

Base.IteratorSize(::FormattedStreamIterator) = Base.SizeUnknown()

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

"""
	eachval(s::FormattedStream)

Iterate over the values contained in the formatted stream `s`.
"""
eachval(s::T) where {T<:FormattedStream} = FormattedStreamIterator{T}(s)

struct PreallocatedFormattedStreamIterator{T1,T2<:FormattedStream{T1}}
	s::T2
	output::T1

	PreallocatedFormattedStreamIterator{T1,T2}(s::T2, output::T1) where
			{T1,T2<:FormattedStream{T1}} = new(s, output)
end

function Base.iterate(iter::PreallocatedFormattedStreamIterator)
	seekstart(iter.s)
	iterate(iter, nothing)
end

function Base.iterate(iter::PreallocatedFormattedStreamIterator, ::Nothing)
	if eof(iter.s)
		nothing
	else
		read!(iter.s, iter.output), nothing
	end
end

Base.IteratorSize(::PreallocatedFormattedStreamIterator) = Base.SizeUnknown()

Base.eltype(::PreallocatedFormattedStreamIterator{T1,T2}) where
		{T1,T2<:FormattedStream{T1}} = T1

"""
	eachval!(s::FormattedStream, output)

Iterate over the values contained in the formatted stream `s`, reading into
pre-allocated `output`.
"""
eachval!(s::T2, output::T1) where {T1,T2<:FormattedStream{T1}} =
		PreallocatedFormattedStreamIterator{T1,T2}(s, output)


text:
import Base: start, get!, send

"""
    ChatterClient(source::Source) -> ChatterClient

Create a new `ChatterClient` from `source`.

# Arguments
* `source::Source`: Source of communication.
"""
struct ChatterClient
    messages::Vector{Message}  # Ordered oldest to newest.
    source::Source  # Source of sending/receiving.

    ChatterClient(source::Source) = new(Vector(), source)
end

"""
    start(client::ChatterClient) -> Task

Make `client` start listening for new messages.

# Arguments
* `client::ChatterClient`: The client to start.

Returns the listener's `Task`.
"""
function start(client::ChatterClient)
    @async begin
        channel = Channel{Message}(Inf)
        @async listen(client.source, channel)
        while true
            push!(client.messages, take!(channel))
        end
    end
end

"""
    get!(client::ChatterClient; n=1) -> Vector{Message}

Get the `n` oldest unread message(s) from `client`.

# Arguments
* `client::ChatterClient`: Client whose messages are being read.

# Keyword Arguments
* `n::Int=1`: Number of messages to fetch.

Returns messages ordered oldest to newest.
"""
function get!(client::ChatterClient; n::Int=1)
    n = min(n, length(client.messages))
    messages = client.messages[1:n]
    deleteat!(client.messages, 1:n)
    return messages
end

"""
    flush!(client::ChatterClient) -> Vector{Message}

Get all unread messages from `client`.

# Arguments
* `client::ChatterClient`: client whose messages are being read.

Returns messages ordered oldest to newest.
"""
flush!(client::ChatterClient) = get!(client; n=length(client.messages))

"""
    send(client::ChatterClient, msg::AbstractString) -> Int

Send a message to `client`'s source.

# Arguments
* `client::ChatterClient`: Client whose source we're sending to.
* `msg::AbstractString`: Message to be sent.

Returns the number of bytes written.
"""
send(client::ChatterClient, msg::AbstractString) = deliver(client.source, msg)