text:
# Ideally the previous _find_all_in method could be used, but things like `div(::Second, ::Integer)`
# don't work. So this helps drop units by didingin by oneunit(::T) of the same type.
unsafe_find_value(val, r::R) where {R}= _add1(_int(static_div(val - static_first(r), static_step(r))))

unsafe_find_value(x, r::MutableRange) = unsafe_find_value(x, parent(r))
unsafe_find_value(x, r::StaticRange) = unsafe_find_value(x, parent(r))
unsafe_find_value(x, r::LinRange) = _add1(_int((x - r.start) / (r.stop - r.start) * r.lendiv))
unsafe_find_value(x, r::StepRangeLen) = _int(((x - r.ref) / step_hp(r)) + r.offset)
unsafe_find_value(x, r::UnitRange) = _add1(_int(x - first(r)))
unsafe_find_value(x, r::StepRange) = _add1(_int((x - first(r)) / step(r)))
unsafe_find_value(x, r::OneTo) = _int(x)
unsafe_find_value(x, r::DynamicAxis) = _int(x)

# only really applies to ordered vectors
_find_all(::Type{T},        fi,        li) where {T} = fi:li
_find_all(::Type{T}, ::Nothing,        li) where {T} = _empty_ur(T)
_find_all(::Type{T},        fi, ::Nothing) where {T} = _empty_ur(T)
_find_all(::Type{T}, ::Nothing, ::Nothing) where {T} = _empty_ur(T)

###
### find_first
###

"""
    find_first(predicate::Function, A)

Return the index or key of the first element of A for which predicate returns
true. Return nothing if there is no such element.

Indices or keys are of the same type as those returned by keys(A) and pairs(A).

# Examples

```jldoctest
julia> using StaticRanges

julia> A = [1, 4, 2, 2];

julia> find_first(iseven, A)
2

julia> find_first(x -> x>10, A) # returns nothing, but not printed in the REPL

julia> find_first(isequal(4), A)
2

julia> find_first(iseven, [1 4; 2 2])
CartesianIndex(2, 1)
```
"""
find_first(f, collection) = _find_first(static_isempty(collection), f, collection)
_find_first(::True, f, collection) = nothing
_find_first(::False, f, collection) = unsafe_find_first(f, collection)
function _find_first(b::Bool, f, collection)
    if b
        return nothing
    else
        return dynamic(unsafe_find_first(f, collection))
    end
end
@inline function unsafe_find_first(f, x)
    for (i, x_i) in pairs(x)
        f(x_i) && return i
    end
    return nothing
end

# <
function _unsafe_find_first_lt_forward(val, r::AbstractRange)
    return ifelseop(ge(static_first(r), val), return_nothing, return_static_one, r)
end
function _unsafe_find_first_lt_reverse(val, r::AbstractRange)
    index = unsafe_find_value(val, r)
    len = static_length(r)
    return ifelseop(le(len, index),
        return_nothing,
        __unsafe_find_first_lt_reverse,
        val, index, r
    )
end
function __unsafe_find_first_lt_reverse(val, index, r::AbstractRange)
    return ifelseop(gt(static(1), index),
        return_static_one,
        ___unsafe_find_first_lt_reverse,
        val, index, r
    )
end
___unsafe_find_first_lt_reverse(val, index, r) = ifadd1(@inbounds(r[index]) >= val, index)
@inline function unsafe_find_first(f::Fix2{typeof(<)}, r::AbstractRange)
    start = static_first(r)
    s = static_step(r)
    return ifelseop(gt(s, zero(s)),
        _unsafe_find_first_lt_forward,
        _unsafe_find_first_lt_reverse,
        f.x, r
    )
end

# <=
@inline function unsafe_find_first(f::Fix2{typeof(<=)}, r::AbstractRange)
    s = static_step(r)
    return ifelseop(gt(s, zero(s)),
        _find_first_le_forward,
        _find_first_le_reverse,
        f.x, r
    )
end
function _find_first_le_forward(val, r::AbstractRange)
    return ifelseop(gt(static_first(r), val), return_nothing, return_static_one, val)
end
function _find_first_le_reverse(val, r::AbstractRange)
    return ifelseop(le(static_first(r), val), return_static_one, __find_first_le_reverse, val, r)
end
function __find_first_le_reverse(val, r::AbstractRange)
    return ifelseop(gt(static_last(r), val), return_nothing, ___find_first_le_reverse, val, r)
end
function ___find_first_le_reverse(val, r::AbstractRange)
    index = unsafe_find_value(val, r)
    if @inbounds(r[index]) > val
        return _add1(index)
    else
        return index
    end
end

# >
@inline function unsafe_find_first(f::Fix2{typeof(>)}, r::AbstractRange)
    s = static_step(r)
    return ifelseop(gt(s, zero(s)), _find_first_gt_forward, _find_first_gt_reverse, f.x, r)
end
function _find_first_gt_forward(val, r::AbstractRange)
    return ifelseop(le(static_last(r), val),
        return_nothing,
        __find_first_gt_forward,
        val, r
    )
end
function __find_first_gt_forward(val, r::AbstractRange)
    return ifelseop(gt(static_first(r), val),
        return_static_one,
        ___find_first_gt_forward,
        val, r
    )
end
function ___find_first_gt_forward(val, r::AbstractRange)
    index = unsafe_find_value(val, r)
    if @inbounds(r[index]) <= val
        return _add1(index)
    else
        return index
    end
end

function _find_first_gt_reverse(val, r::AbstractRange)
    return ifelseop(gt(static_first(r), val),
        return_static_one,
        __find_first_gt_reverse,
        val, r
    )
end
function __find_first_gt_reverse(val, r)
    return ifelseop(lt(static_last(r), val),
        return_nothing,
        unsafe_find_value,
        val, r
    )
end

# >=
@inline function unsafe_find_first(f::Fix2{typeof(>=)}, r::AbstractRange)
    s = static_step(r)
    return ifelseop(gt(s, zero(s)),
        _find_first_ge_forward,
        _find_first_ge_reverse,
        f.x, r
    )
end
function _find_first_ge_forward(val, r::AbstractRange)
    return ifelseop(lt(static_last(r), val),
        return_nothing,
        __find_first_ge_forward,
        val, r
    )
end
function __find_first_ge_forward(val, r::AbstractRange)
    return ifelseop(ge(static_first(r), val),
        return_static_one,
        ___find_first_ge_forward,
        val, r
    )
end
function ___find_first_ge_forward(val, r::AbstractRange)
    index = unsafe_find_value(val, r)
    if @inbounds(r[index]) < val
        return _add1(index)
    else
        return index
    end
end
function _find_first_ge_reverse(val, r::AbstractRange)
    return ifelseop(ge(static_first(r), val),
        return_static_one,
        return_nothing,
        r
    )
end

# ==
function unsafe_find_first(f::Equal, collection::AbstractRange)
    x = f.x
    if minimum(collection) > x || maximum(collection) < x
        return nothing
    else
        index = unsafe_find_value(x, collection)
        if @inbounds(collection[index]) == x
            return index
        else
            return nothing
        end
    end
end

###
### find_last
###
"""
    find_last(predicate::Function, A)

Return the index or key of the last element of A for which predicate returns
true. Return nothing if there is no such element.

Indices or keys are of the same type as those returned by keys(A) and pairs(A).

# Examples

```jldoctest
julia> using StaticRanges

julia> find_last(iseven, [1, 4, 2, 2])
4

julia> find_last(x -> x>10, [1, 4, 2, 2]) # returns nothing, but not printed in the REPL

julia> find_last(isequal(4), [1, 4, 2, 2])
2

julia> find_last(iseven, [1 4; 2 2])
CartesianIndex(2, 2)
```
"""
find_last(f, collection) = _find_last(static_isempty(collection), f, collection)
_find_last(::True, f, collection) = nothing
_find_last(::False, f, collection) = unsafe_find_last(f, collection)
function _find_last(b::Bool, f, collection)
    if b
        return nothing
    else
        return dynamic(unsafe_find_last(f, collection))
    end
end

@inline unsafe_find_last(f::Equal, x) = unsafe_find_lasteq(f.x,   x)
@inline function unsafe_find_last(f, collection)
    for (i, collection_i) in Iterators.reverse(pairs(collection))
        f(collection_i) && return i
    end
    return nothing
end

# >
@inline function unsafe_find_last(f::Fix2{typeof(>)}, r::AbstractRange)
    s = static_step(r)
    return ifelseop(gt(s, zero(s)),
        _find_last_gt_forward,
        _find_last_gt_reverse,
        f.x, r
    )
end
function _find_last_gt_forward(val, r::AbstractRange)
    return ifelseop(le(static_last(r), val),
        return_nothing,
        static_length,
        r
    )
end

function _find_last_gt_reverse(val, r::AbstractRange)
    return ifelseop(le(static_first(r), val),
        return_nothing,
        __find_last_gt_reverse,
        val, r
    )
end
function __find_last_gt_reverse(val, r::AbstractRange)
    ifelseop(gt(static_last(r), val),
        _static_length,
        ___find_last_gt_reverse,
        val, r
    )
end
function ___find_last_gt_reverse(val, r::AbstractRange)
    index = unsafe_find_value(val, r)
    if (@inbounds(r[index]) == val) & (index != 1)
        return _sub1(index)
    else
        return index
    end
end

# >=
@inline function unsafe_find_last(f::Fix2{typeof(>=)}, r::AbstractRange)
    s = static_step(r)
    return ifelseop(gt(s, zero(s)),
        _find_last_ge_forward,
        _find_last_ge_reverse,
        f.x, r
    )
end
function _find_last_ge_forward(val, r::AbstractRange)
    return ifelseop(lt(static_last(r), val),
        return_nothing,
        static_length,
        r
    )
end
function _find_last_ge_reverse(val, r::AbstractRange)
    index = unsafe_find_value(val, r)
    return ifelseop(gt(static(1), index),
        return_nothing,
        __find_last_ge_reverse,
        val, index, r
    )
end
function __find_last_ge_reverse(val, index, r::AbstractRange)
    ifelseop(gt(static_last(r), val),
        _static_length,
        ___find_last_ge_reverse,
        val, index, r
    )
end
function ___find_last_ge_reverse(val, index, r::AbstractRange)
    if @inbounds(r[index]) < val
        return _sub1(index)
    else
        return index
    end
end

# <
#=@inline function unsafe_find_last(f::Fix2{typeof(<)}, r::AbstractRange)
    x = f.x
    if step(r) > 0
        index = unsafe_find_value(x, r)
        if 1 > index
            return nothing
        else
            if lastindex(r) < index
                return lastindex(r)
            else
                if @inbounds(r[index]) < x
                    return index
                else
                    return ifelseop(ne(index, static(1)), _sub1, return_nothing, index)
                end
            end
        end
    else
        if last(r) >= x
            return nothing
        else
            return lastindex(r)
        end
    end
end=#

@inline function unsafe_find_last(f::Fix2{typeof(<)}, r::AbstractRange)
    s = static_step(r)
    return ifelseop(gt(s, zero(s)),
        _find_last_lt_forward,
        _find_last_lt_reverse,
        f.x, r
    )
end
function _find_last_lt_forward(val, r::AbstractRange)
    index = unsafe_find_value(val, r)
    ifelseop(gt(static(1), index), return_nothing, __find_last_lt_forward, val, index, r)
end
function __find_last_lt_forward(val, index, r::AbstractRange)
    len = static_length(r)
    return ifelseop(lt(len, index), first, ___find_last_lt_forward, (len, val, index, r))
end
___find_last_lt_forward(x::Tuple) = ___find_last_lt_forward(x...)
function ___find_last_lt_forward(len, val, index, r::AbstractRange)
    if @inbounds(r[index]) < val
        return index
    elseif index != 1
        return _sub1(index)
    else
        return nothing
    end
end
function _find_last_lt_reverse(val, r::AbstractRange)
    return ifelseop(lt(static_last(r), val),
        static_length,
        return_nothing,
        r
    )
end

# <=
@inline function unsafe_find_last(f::Fix2{typeof(<=)}, r::AbstractRange)
    s = static_step(r)
    return ifelseop(gt(s, zero(s)),
        _find_last_le_forward,
        _find_last_le_reverse,
        f.x, r
    )
end
function _find_last_le_forward(val, r::AbstractRange)
    return ifelseop(gt(static_first(r), val),
        return_nothing,
        __find_last_le_forward,
        val, r
    )
end
function __find_last_le_forward(val, r::AbstractRange)
    return ifelseop(
        le(static_last(r), val),
        _static_length,
        ___find_last_le_forward,
        val, r
    )
end
function ___find_last_le_forward(val, r::AbstractRange)
    index = unsafe_find_value(val, r)
    if @inbounds(r[index]) <=  val
        return index
    elseif index != 1
        return _sub1(index)
    else
        return nothing
    end
end
function _find_last_le_reverse(val, r::AbstractRange)
    return ifelseop(gt(static_last(r), val),
        return_nothing,
        static_length,
        r
    )
end

# ==
function unsafe_find_lasteq(x, collection)
    for (i, collection_i) in Iterators.reverse(pairs(collection))
        x == collection_i && return i
    end
    return nothing
end
unsafe_find_last(f::Equal, collection::AbstractRange) = unsafe_find_first(f, collection)

###
### find_all
###
find_all(f, x) = collect(k for (k,v) in pairs(x) if f(v))
find_all(f::Fix2{typeof(in)}, x) = find_all_in(f.x, x)
@inline find_all(f::Or, x::AbstractRange) = combine(find_all(f.f1, x), find_all(f.f2, x))
@inline find_all(f::And, x::AbstractRange) = intersect(find_all(f.f1, x), find_all(f.f2, x))

@inline function find_all(f::Fix2{typeof(<)}, r::AbstractRange{T}) where {T}
    x = f.x
    if step(r) > zero(T)
        return _find_all(keytype(r), firstindex(r), find_last(<(x), r))
    elseif step(r) < zero(T)
        return _find_all(keytype(r), find_first(<(x), r), lastindex(r))
    else
        return _find_all(keytype(r), nothing, nothing)
    end
end

@inline function find_all(f::Fix2{typeof(<=)}, r::AbstractRange{T}) where {T}
    x = f.x
    if step(r) > zero(T)
        return _find_all(keytype(r), firstindex(r), find_last(<=(x), r))
    elseif step(r) < zero(T)
        return _find_all(keytype(r), find_first(<=(x), r), lastindex(r))
    else
        return _find_all(keytype(r), nothing, nothing)
    end
end

@inline function find_all(f::Fix2{typeof(>)}, r::AbstractRange{T}) where {T}
    x = f.x
    if step(r) > zero(T)
        return _find_all(keytype(r), find_first(>(x), r), lastindex(r))
    elseif step(r) < zero(T)
        return _find_all(keytype(r), firstindex(r), find_last(>(x), r))
    else
        return _empty_ur(keytype(r))
    end
end

@inline function find_all(f::Fix2{typeof(>=)}, r::AbstractRange{T}) where {T}
    x = f.x
    if step(r) > zero(T)
        return _find_all(keytype(r), find_first(>=(x), r), lastindex(r))
    elseif step(r) < zero(T)
        return _find_all(keytype(r), firstindex(r), find_last(>=(x), r))
    else
        return _empty_ur(keytype(r))
    end
end

@inline function find_all(f::Equal, r::AbstractRange{T}) where {T}
    if (step(r) > zero(T)) | (step(r) < zero(T))
        return _find_all(keytype(r), find_first(f, r), find_last(f, r))
    else
        return _empty_ur(keytype(r))
    end
end

@inline function find_all(f::NotIn, r::AbstractRange{T}) where {T}
    x = f.x
    if (step(r) > zero(T)) | (step(r) < zero(T))
        return combine(find_all(>(x), r), find_all(<(x), r))
    else
        return GapRange(UnitRange(1, 0), UnitRange(1, 0))
    end
end

###
### find_all_in
###

# find x in y

# find all values of x in y
@inline function find_all_in(x::Interval{:closed,:closed,T}, y) where {T}
    return intersect(find_all(>=(x.left), y), find_all(<=(x.right), y))
end

@inline function find_all_in(x::Interval{:closed,:open,T}, y) where {T}
    return intersect(find_all(>=(x.left), y), find_all(<(x.right), y))
end
@inline function find_all_in(x::Interval{:open,:closed,T}, y) where {T}
    return intersect(find_all(>(x.left), y), find_all(<=(x.right), y))
end
@inline function find_all_in(x::Interval{:open,:open,T}, y) where {T}
    return intersect(find_all(>(x.left), y), find_all(<(x.right), y))
end

@inline function find_all_in(x::AbstractRange{Tx}, y::AbstractRange{Ty}) where {Tx,Ty}
    if isempty(x)
        return _empty(x, y)
    elseif isempty(y)
        return _empty(x, y)
    elseif known_step(x) === known_step(y) !== nothing
        return _find_all_in_same_step(x, y)
    else
        return _find_all_in(x, y)
    end
end

function find_all_in(a, b)
    ind  = Vector{eltype(keys(a))}()
    bset = Set(b)
    @inbounds for (i,ai) in pairs(a)
        ai in bset && push!(ind, i)
    end
    return ind
end
#=

    r = @inferred(find_all(in(UnitRange(1,10)), static(UnitRange(1,8))))
    @test r == 1:8
    @test isa(r, UnitRange)

    r = @inferred(find_all(in(Base.OneTo(10)), StaticRanges.mutable(UnitRange(1, 8))))
    @test r == OneTo(8)
    @test isa(r, UnitRange) == true

    r = @inferred(find_all(in(static(UnitRange(1, 8))), static(UnitRange(1, 10))))
    @test r == static(UnitRange(1, 8))

    @test find_all(in(collect(1:10)), 1:20) == find_all(in(1:10), 1:20)
    @test find_all(in(1:10), collect(1:20)) == 1:10

    @testset "steps match but no overlap" begin
        r = @inferred(find_all_in(1:3, 4:5))
        @test r == 1:0
        @test isa(r, UnitRange)
    end

=#

# if steps are not integers then we need to make sure that something like
# startx = 1 and starty = 1.1 doesn't make the whole thing offset from eachother
@inline function _find_all_in_same_step(x::AbstractRange, y::AbstractRange)
    fx = _drop_unit(static_first(x))
    fy = _drop_unit(static_first(y))
    lx = _drop_unit(static_last(x))
    ly = _drop_unit(static_last(y))
    xydiff = fx - fy
    if !iszero(static_rem(xydiff, static(1)))
        return one(xydiff):zero(xydiff)
    elseif static_step(x) > 0
        return _int((max(fx, fy) - fy) + static(1)):_int((min(lx, ly) - fy) + static(1))
    else
        return _int((min(fx, fy) - fy) + static(1)):_int((max(lx, ly) - fy) + static(1))
    end
end
@inline function _find_all_in_same_step(x::AbstractRange{<:Integer}, y::AbstractRange{<:Integer})
    fx = _drop_unit(static_first(x))
    fy = _drop_unit(static_first(y))
    lx = _drop_unit(static_last(x))
    ly = _drop_unit(static_last(y))
    if _drop_unit(static_step(x)) > 0
        return _int(max(fx, fy) - fy + static(1)):_int(min(lx, ly) - fy + static(1))
    else
        return _int(min(fx, fy) - fy + static(1)):_int(max(lx, ly) - fy + static(1))
    end
end

#=
function _findin(r::AbstractRange{<:Integer}, span::AbstractUnitRange{<:Integer})
    local ifirst
    local ilast
    fspan = first(span)
    lspan = last(span)
    fr = first(r)
    lr = last(r)
    sr = step(r)
    if sr > 0
        ifirst = fr >= fspan ? 1 : ceil(Integer,(fspan-fr)/sr)+1
        ilast = lr <= lspan ? length(r) : length(r) - ceil(Integer,(lr-lspan)/sr)
    elseif sr < 0
        ifirst = fr <= lspan ? 1 : ceil(Integer,(lspan-fr)/sr)+1
        ilast = lr >= fspan ? length(r) : length(r) - ceil(Integer,(lr-fspan)/sr)
    else
        ifirst = fr >= fspan ? 1 : length(r)+1
        ilast = fr <= lspan ? length(r) : 0
    end
    r isa AbstractUnitRange ? (ifirst:ilast) : (ifirst:1:ilast)
end
=#
###
### manage steps
###
function _to_step(s::S, x::SX, y::SY) where {S,SX,SY}
    return __to_step(eq(gt(y, zero(y)), gt(x, zero(x))), s)
end
__to_step(::True, s::S) where {S} = _int(s)
__to_step(::False, s::S) where {S} = _int(_flip_int(s))
function __to_step(flip::Bool, s::S) where {S}
    if flip
        return dynamic(s)
    else
        return dynamic(-s)
    end
end
_flip_int(x::Int) = -x
_flip_int(x::StaticInt) = -x
_flip_int(x::StaticFloat64{X}) where {X} = static(_flip_int(X))
_flip_int(x) = -x

###
### steps are not exactly the same
###
@inline function _find_all_in(x::AbstractRange, y::AbstractRange)
    sx = _drop_unit(static_step(x))
    sy = _drop_unit(static_step(y))
    sxy = _drop_unit(static_div(sx, sy))
    startx = _drop_unit(static_first(x))
    starty = _drop_unit(static_first(y))
    stopx = _drop_unit(static_last(x))
    stopy = _drop_unit(static_last(y))

    if iszero(sxy)
        sxy2 = div(sy, sx)
        if !iszero(rem(startx - starty, div(sxy2, sx)))
            return _empty(x, y)
        else
            suby = intersect(x, y)
            f = _int(unsafe_find_first(==(static_first(suby)), y))
            l = _int(unsafe_find_first(==(static_last(suby)), y))
            return f:_to_step(static(1), sx, sy):l
        end
    else
        if !iszero(rem(startx - starty, div(sxy, sx)))
            return _empty(x, y)
        else
            suby = intersect(x, y)
            f = _int(unsafe_find_first(==(static_first(suby)), y))
            l = _int(unsafe_find_first(==(static_last(suby)), y))
            return f:_to_step(sxy, sx, sy):l
        end
    end

    #=
    if iszero(sxy)
        sxy2 = div(sy, sx)
        if !iszero(rem(startx - starty, div(sxy2, sx)))
            return _empty(x, y)
        else
            suby = intersect(x, y)
            f = _int(unsafe_find_first(==(static_first(suby)), y))
            l = _int(unsafe_find_first(==(static_last(suby)), y))
            return f:_to_step(static(1), sx, sy):l
        end
    else
        if !iszero(rem(startx - starty, div(sxy, sx)))
            return _empty(x, y)
        else
            suby = intersect(x, y)
            f = _int(unsafe_find_first(==(static_first(suby)), y))
            l = _int(unsafe_find_first(==(static_last(suby)), y))
            return f:_to_step(sxy, sx, sy):l
        end
    end
    =#
end

#=
@inline function _find_all_in(x, y)
    fx = _drop_unit(static_first(x))
    fy = _drop_unit(static_first(y))
    sx = _drop_unit(static_step(x))
    sy = _drop_unit(static_step(y))
    lx = _drop_unit(static_last(x))
    ly = _drop_unit(static_last(y))

    if sx > 0
        if sy > 0
            f = max(fx, fy)
            l = min(lx, ly)
        else
            f = min(lx, fy)
            l = max(fx, ly)
        end
    else
        if sy > 0
            f = max(lx, fy)
            l = min(fx, ly)
        else
            f = min(fx, fy)
            l = max(lx, ly)
        end
    end
end
=#

#=

f1 = 1:2:11
f2 = 3:2:9
f3 = 3:4:15

r1 = 11:-2:1
r2 = 9:-2:2
r3 = 15:-4:3

find_all(in(f1), f1) == 1:length(f1)
find_all(in(f2), f1) == 2:5
find_all(in(f2), f1) == 2:2:6

find_all(in(f1), r1) == length(f1):-1:1
find_all(in(f2), r1) == 5:-1:2
find_all(in(f2), r1) == 6:-2:2

=#

combine(x::OneToUnion, y::OneToUnion) = promote_type(typeof(x), typeof(y))(max(x, y))

###
### AbstractUnitRange
###
function combine(x::AbstractUnitRange{<:Integer}, y::AbstractUnitRange{<:Integer})
    R = promote_type(typeof(x), typeof(y))
    if isempty(x)
        if isempty(y)
            return GapRange(R(1, 0), R(1, 0))
        else
            return GapRange(R(1, 0), R(first(y), last(y)))
        end
    elseif isempty(y)
        return GapRange(R(1, 0), R(first(x), last(x)))
    else
        xmax = last(x)
        xmin = first(x)
        ymax = last(y)
        ymin = first(y)
        if xmax < ymin  # all x below y
            return GapRange(R(xmin, xmax), R(ymin, ymax))
        elseif ymax < xmin  # all y below x
            return GapRange(R(ymin, ymax), R(xmin, xmax))
        else # x and y overlap so we just set the first range to length of one
            rmin = min(xmin, ymin)
            return GapRange(R(rmin, rmin), R(rmin + oneunit(eltype(R)), max(xmax, ymax)))
        end
    end
end


text:
# TODO: Consider a dynamic rate problem, add time series in the rate calculations.

struct Network
    branches::Array{Branch}
    ybus::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}
    ptdf::Union{Array{Float64},Nothing}
    incidence::Array{Int}
end

function Network(branches::Array{T}, nodes::Array{Bus}) where {T<:Branch}
    ybus = build_ybus(length(nodes),branches);
    ptdf, A = buildptdf(branches, nodes)

    return Network(branches, ybus, ptdf, A)

end

ext: jl


text:
"""
    partition(array::Array{Int64,1}, low::Int64, high::Int64)

The partion algorithm is shuffeling the array for a given interval of low and high
boundaries.

# Arguments
- `array::Array{Int64,1}`: Unsorted array of integers
- `low::Int64`: Lowest index of the unsorted array or subarray
- `high::Int64`: Highes index of the unsorted array or subarray

"""
function partition(array::Array{Int64,1}, low::Int64, high::Int64)
    # index of smaller element
    i = low - 1
    # pivot
    pivot = array[high]

    for j in low:high

        # If current element is smaller than the pivot
        if array[j] < pivot

            # increment index of smaller element
            i = i + 1
            # Swap if the element found is greater than the next element
            array[i], array[j] = array[j], array[i]
        end
    end
    array[i+1], array[high] = array[high], array[i+1]
    return (i + 1)
end

"""
    quick_sorting(array::Array{Int64,1}, low = nothing, high = nothing)

The quick sort algorithm (QSA) works by selecting a pivot element from the array and
partitioning the other elements into two subarrays, according to whether they are less than
or greater than the pivot-window. Then the sorting of subarrays is recursively organized.
This procedure repeatedly happens until each subarray is organized; consequently, the
subarrays' merging is an organized array. For more information see:
[https://en.wikipedia.org/wiki/Quicksort#Parallelization](https://en.wikipedia.org/wiki/Quicksort#Parallelization)

# Arguments
- `array::Array{Int64,1}`: Unsorted array of integers
- `low::Int64`: Lowest index of the unsorted array or subarray
- `high::Int64`: Highes index of the unsorted array or subarray

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [64, 34, 25, 12, 22, 11, 90]
julia> ClassicAlgorithmsCollections.quick_sorting(arr)
[11, 12, 22, 25, 34, 64, 90]
```
"""
function quick_sorting(array::Array{Int64,1}, low = nothing, high = nothing)
    # If nothing, then the intial run, so low and high are set to array size
    if isnothing(low) && isnothing(high)
        low = 1
        high = length(array)
    end
    if low < high

        # p_i is partitioning index, array[p_i] is now at right place
        p_i = partition(array, low, high)

        # Separately sort elements before partition and after partition
        array = quick_sorting(array, low, p_i - 1)
        array = quick_sorting(array, p_i + 1, high)
    else
        return array
    end
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
using Comonicon.AST
using Test
using Faker

arg = Argument(; name = "arg", type = Int)

macro test_show(mime, ex)
    Meta.isexpr(ex, :block) || error("expect begin ... end")
    ret = Expr(:block)
    for each in ex.args
        if Meta.isexpr(each, :call) && each.args[1] === :in
            @gensym buf object
            push!(ret.args, :($buf = IOBuffer()))
            push!(ret.args, :($object = $(each.args[3])))
            push!(ret.args, :(show($buf, $mime(), $object)))
            push!(ret.args, :(@test occursin($(each.args[2]), String(take!($buf)))))
        else
            push!(ret.args, each)
        end
    end
    return esc(ret)
end

@test_show MIME"text/plain" begin
    "<arg>" in Argument(; name = "arg")
    "[arg...]" in Argument(; name = "arg", vararg = true)
    "<arg::Int64>" in Argument(; name = "arg", type = Int)
    "--option-a" in Option(; sym = :option_a)
    "--option-a <hint>" in Option(; sym = :option_a, hint = "hint")
    "--flag-a" in Flag(; sym = :flag_a)
end

leaf = LeafCommand(;
    fn = identity,
    name = "leaf",
    args = [Argument(; name = "arg", description = Faker.text())],
    flags = Dict(
        "flag-a" => Flag(; sym = :flag_a, description = "flag a."),
        "flag-b" => Flag(; sym = :flag_b, description = "flag b."),
    ),
    description = Faker.text(),
)

@test_show MIME"text/plain" begin
    "  leaf <args> [options] [flags]" in leaf
    "Args\n\n" in leaf
    "  <arg>" in leaf
    "Flags\n\n" in leaf
end

@test_throws ErrorException NodeCommand(; name = "abc", subcmds = Dict{String,Any}())

node = NodeCommand(; name = "foo", subcmds = Dict("leaf" => leaf))

@test_show MIME"text/plain" begin
    "  foo <command>" in node
    "Commands\n\n" in node
    "  leaf <arg>" in node
    "Flags\n\n" in node
    "-h, --help" in node
end

text = "registry you want to register the package. If the package has not been registered, ion will try to register the package in the General registry. Or the user needs to specify the registry to register using this option."

AST.print_indent_content(stdout, text, AST.Terminal(), 10)


text:
using LinearAlgebra
using Statistics

function optimize(f, g, x0, n, prob)
    if prob == "secret_2" || prob == "secret_1"
        basis(i, m) = [k == i ? 1.0 : 0.0 for k in 1 : m]

        S = [basis(i, length(x0)) * x0[i] for i in 1:length(x0)]
        push!(S, randn(length(x0)))

        x = nelder_mead(f, S, 1e-8, n)
        return x
    else
        x = hooke_jeeves(f, x0, 1.0, 0.01, n)
    end
    return x
end

function hooke_jeeves(f, x, α, ϵ, counter, γ=0.5)
    basis(i, m) = [k == i ? 1.0 : 0.0 for k in 1 : m]
    while counter > 0
        y, n = f(x), length(x)
        counter -= 1
        while α > ϵ && counter > 2*n
            improved = false
            x_best, y_best = x, y
            for i in 1 : n
                for sgn in (-1, 1)
                    x′ = x + sgn*α*basis(i, n)
                    y′ = f(x′)
                    counter -= 1
                    if y′ < y_best
                        x_best, y_best, improved = x′, y′, true
                    end
                end
            end
            x, y = x_best, y_best

            if !improved
                α *= γ
            end
        end
    end
    return x
end

function nelder_mead(f, S, ϵ, n; α=1.0, β=2.0, γ=0.5)
    counter = n
    Δ, y_arr = Inf, f.(S)
    counter -= length(S)
    while Δ > ϵ && counter > length(S) #+ 1 #maybe not necessary
        p = sortperm(y_arr) #sort lowest to highest
        S, y_arr = S[p], y_arr[p]
        xl, yl = S[1], y_arr[1] #lowest
        xh, yh = S[end], y_arr[end] #highest
        xs, ys = S[end-1], y_arr[end-1] #second-highest
        xm = mean(S[1:end-1]) #centroid
        xr = xm + α*(xm - xh) #reflection point
        yr = f(xr)
        counter -= 1 #MINUS 1

        if yr < yl
            xe = xm + β*(xr-xm) #expansion point
            ye = f(xe)
            counter -= 1 #MINUS 1
            S[end], y_arr[end] = ye < yr ? (xe, ye) : (xr, yr)
        elseif yr > ys
            if yr ≤ yh
                xh, yh, S[end], y_arr[end] = xr, yr, xr, yr
            end
            xc = xm + γ*(xh-xm) #contraction point
            yc = f(xc)
            counter -= 1 #MINUS 1
            if yc > yh
                for i in 2 : length(y_arr)
                    S[i] = (S[i] + xl)/2
                    y_arr[i] = f(S[i])
                    counter -= length(y_arr)-1 #MINUS (length(y_arr)-1)
                end
            else
                S[end], y_arr[end] = xc, yc
            end
        else
            S[end], y_arr[end] = xr, yr
        end

        Δ = std(y_arr, corrected=false)
    end
    return S[argmin(y_arr)]
end


text:
# Integrate cl distributions

aircraft_CL = data_PROWIM_CL_alpha[1]
parameters_CL = data_PROWIM_CL_alpha[2]

lifting_line_rs_CL = aircraft_CL.wing_system.lifting_line_rs[1]
lifting_line_xs_CL = lifting_line_rs_CL[1,:]
lifting_line_ys_CL = lifting_line_rs_CL[2,:]
lifting_line_zs_CL = lifting_line_rs_CL[3,:]

cpxs_CL = vcat([panel.rtl[1] for panel in aircraft_CL.wing_system.system.surfaces[1]][1,:], aircraft_CL.wing_system.system.surfaces[1][end].rtr[1])
cpys_CL = vcat([panel.rtl[2] for panel in aircraft_CL.wing_system.system.surfaces[1]][1,:], aircraft_CL.wing_system.system.surfaces[1][end].rtr[2])
cpzs_CL = vcat([panel.rtl[3] for panel in aircraft_CL.wing_system.system.surfaces[1]][1,:], aircraft_CL.wing_system.system.surfaces[1][end].rtr[3])

aircraft_cl = data_PROWIM_lift_distribution[1]
parameters_cl = data_PROWIM_lift_distribution[2]

lifting_line_rs_cl = aircraft_cl.wing_system.lifting_line_rs[1]
lifting_line_xs_cl = lifting_line_rs_cl[1,:]
lifting_line_ys_cl = lifting_line_rs_cl[2,:]
lifting_line_zs_cl = lifting_line_rs_cl[3,:]

cpxs_cl = vcat([panel.rtl[1] for panel in aircraft_cl.wing_system.system.surfaces[1]][1,:], aircraft_cl.wing_system.system.surfaces[1][end].rtr[1])
cpys_cl = vcat([panel.rtl[2] for panel in aircraft_cl.wing_system.system.surfaces[1]][1,:], aircraft_cl.wing_system.system.surfaces[1][end].rtr[2])
cpzs_cl = vcat([panel.rtl[3] for panel in aircraft_cl.wing_system.system.surfaces[1]][1,:], aircraft_cl.wing_system.system.surfaces[1][end].rtr[3])

fig_check_geometry_x = plt.figure("prowim_geometry_x")
fig_check_geometry_x.clear()
fig_check_geometry_x.add_subplot(111)

ax_x = fig_check_geometry_x.get_axes()[1]
ax_x.plot(1:length(lifting_line_xs_CL), lifting_line_xs_CL, "+", label="CL-alpha, lifting_line_xs")
ax_x.plot(1:length(cpxs_CL), cpxs_CL, "x", label="CL-alpha, bound vortex")
ax_x.plot(1:length(lifting_line_xs_cl), lifting_line_xs_cl, "-", label="cl dist, lifting_line_xs")
ax_x.plot(1:length(cpxs_cl), cpxs_cl, "--", label="cl dist, bound vortex")
ax_x.set_ylabel(LS.L"x")
ax_x.set_xlabel(LS.L"i")
ax_x.legend(loc="upper left", bbox_to_anchor=(1.01,1))

fig_check_geometry_y = plt.figure("prowim_geometry_y")
fig_check_geometry_y.clear()
fig_check_geometry_y.add_subplot(111)

ax_y = fig_check_geometry_y.get_axes()[1]
ax_y.plot(1:length(lifting_line_ys_CL), lifting_line_ys_CL, "+", label="CL-alpha, lifting_line_ys")
ax_y.plot(1:length(cpys_CL), cpys_CL, "x", label="CL-alpha, bound vortex")
ax_y.plot(1:length(lifting_line_ys_cl), lifting_line_ys_cl, "-", label="cl dist, lifting_line_ys")
ax_y.plot(1:length(cpys_cl), cpys_cl, "--", label="cl dist, bound vortex")
ax_y.set_ylabel(LS.L"y")
ax_y.set_xlabel(LS.L"i")
ax_y.legend(loc="upper left", bbox_to_anchor=(1.01,1))

fig_check_geometry_z = plt.figure("prowim_geometry_z")
fig_check_geometry_z.clear()
fig_check_geometry_z.add_subplot(111)

ax_z = fig_check_geometry_z.get_axes()[1]
ax_z.plot(1:length(lifting_line_zs_CL), lifting_line_zs_CL, "+", label="CL-alpha, lifting_line_zs")
ax_z.plot(1:length(cpxs_CL), cpzs_CL, "x", label="CL-alpha, bound vortex")
ax_z.plot(1:length(lifting_line_zs_cl), lifting_line_zs_cl, "-", label="cl dist, lifting_line_zs")
ax_z.plot(1:length(cpxs_cl), cpzs_cl, "--", label="cl dist, bound vortex")
ax_z.set_ylabel(LS.L"z")
ax_z.set_xlabel(LS.L"i")
ax_z.legend(loc="upper left", bbox_to_anchor=(1.01,1))

savename_x = "prowim_check_geometry_x.pdf"
savename_y = "prowim_check_geometry_y.pdf"
savename_z = "prowim_check_geometry_z.pdf"
fig_check_geometry_x.savefig(joinpath(plot_directory, savename_x), bbox_inches="tight")
fig_check_geometry_y.savefig(joinpath(plot_directory, savename_y), bbox_inches="tight")
fig_check_geometry_z.savefig(joinpath(plot_directory, savename_z), bbox_inches="tight")

# properties_cl = VL.get_surface_properties(aircraft_cl.wing_system.system)

# vtk_directory = joinpath(AS.topdirectory, "test", "data", "vtk", "prowim")
# # if !isdir(vtk_directory); mkpath(vtk_directory); end
# # VL.write_vtk(joinpath(vtk_directory, "prowim_wing"), aircraft_cl.wing_system.surfaces, properties)
# VL.write_vtk("prowim_wing", aircraft_cl.wing_system.system.surfaces, properties; symmetric = false)

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using Test
using EliminateGraphs

@testset "constructors and properties" begin
    tbl = Bool[false true false false; true false true true; false true false true; false true true false]
    eg = EliminateGraph(tbl)
    @test check_validity(eg)
    @test eg == EliminateGraph(4,[1=>2, 2=>3,2=>4,3=>4])
    @test collect(vertices(eg)) == collect(1:4)

    @test subgraph(disconnected_cliques_eg(3,4), [1,2,3]) == K_eg(3)
    @test subgraph(K_eg(3,4), [1,2,3]) == empty_eg(3)
end

@testset "neighbors" begin
    tbl = Bool[false true true false; true false true true; false true false true; true true true false]
    eg = EliminateGraph(tbl .& tbl')
    @show eg
    nbs = neighbors(eg, 1)
    @test nbs == [2]
    @test neighbors2(eg, 1) == [3,4]
    @test neighborcover(eg, 1) == [1,2]
    @test nv(eg) == 4
end

@testset "degree" begin
    tbl = Bool[false true true false; true false true true; false true false true; true true true false]
    eg = EliminateGraph(tbl .& tbl')
    @test degrees(eg) == [1, 3, 2, 2]

    vs = vertices(eg)
    @test [vs...] == vertices(eg)
    @test vs[3] == vertices(eg)[3]
    vmin, dmin = mindegree_vertex(eg)
    @test degree(eg, vmin) == minimum(degrees(eg)) == dmin
    vmax, dmax = maxdegree_vertex(eg)
    @test degree(eg, vmax) == maximum(degrees(eg)) == dmax
    vmin, vmax, dmin, dmax = minmaxdegree_vertex(eg)
    @test degree(eg,vmin) == minimum(degrees(eg)) == dmin
    @test degree(eg,vmax) == maximum(degrees(eg)) == dmax
end

@testset "vertices and neighbors - eliminated" begin
    tbl = Bool[false true true false; true false true true; false true false true; true true true false]
    eg = EliminateGraph(tbl .& tbl')
    eliminate!(eg, 2)
    eliminate!(eg, 3)
    @test nv(eg) == 2
    recover!(eg)
    @test nv(eg) == 3
    eg2 = eg \ 3
    @test nv(eg) == 3
    @test nv(eg2) == 2
    @test vertices(eg2) == [1,4]
    @test neighbors(eg2, 1) == []
    @test degrees(eg) == [1, 0, 1]
    @test degrees(eg2) == [0, 0]
    eg3 = eg2 \ (1,4)
    @test nv(eg3) == 0
end

@testset "eliminate and recover" begin
    tbl = Bool[false true true false; true false true true; false true false true; true true true false]
    eg = EliminateGraph(tbl .& tbl')
    @test check_validity(eg)
    res = eliminate(eg, 3) do eg
        @test check_validity(eg)
        nv(eg)
    end # do
    @test res == 3

    eg4 = eg \ NeighborCover(1)
    @test vertices(eg4) == [3,4]

    res = eliminate(eg, (3,4)) do eg
        @test check_validity(eg)
        eliminate(eg, (1,2)) do eg
            @test check_validity(eg)
            nv(eg)
        end
    end
    @test res == 0
    @test nv(eg) == 4
    @test eg.level == 0
end

@testset "cluster and clique" begin
    # isclique
    g = disconnected_cliques_eg(3,4)
    @test isclique(g, [1,2,3])
    @test !isclique(g, [1,2,4])
    @test !isclique(g)

    # findcluster
    @test find_cluster(g, 1) == [1,2,3]
    @test find_cluster(g, 4) == [4,5,6,7]
    @test find_cluster(petersen_graph, 1) == collect(1:10)

    # is connected
    @test !isconnected(g)
    @test isconnected(petersen_graph)

    g = K_eg(3)
    @test isclique(g)
end

@testset "mirrors" begin
    @test mirrors(petersen_graph, 1) == []
    @test mirrorcover(petersen_graph, 1) == [1]

    g = K_eg(3,3)
    @test Set(mirrorcover(g, 1)) == Set([1,2,3])
    @test Set(mirrorcover(g, 4)) == Set([4,5,6])

    eliminate(g, MirrorCover(1)) do eg
        @test check_validity(eg)
        @test Set(vertices(eg)) == Set([4,5,6])
    end
    @test g == K_eg(3,3)

    @test generate_set(g, Mirrors{OPEN}(1) ∪ Vertex(1)) == generate_set(g, MirrorCover(1))
end


text:
function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct AnalyticalSolution
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function AnalyticalSolution(
        inradius,
        outradius,
        center,
        ls,
        ms,
        lc,
        mc,
        theta0,
    )
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::AnalyticalSolution, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::AnalyticalSolution)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if ur ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function shell_radial_stress(ls, ms, theta0, A1, A2, r)
    return (ls + 2ms) * (A1 - A2 / r^2) + ls * (A1 + A2 / r^2) -
           (ls + 2ms / 3) * theta0
end

function shell_circumferential_stress(ls, ms, theta0, A1, A2, r)
    return ls * (A1 - A2 / r^2) + (ls + 2ms) * (A1 + A2 / r^2) -
           (ls + 2ms / 3) * theta0
end

function shell_out_of_plane_stress(ls, ms, A1, theta0)
    return 2 * ls * A1 - (ls + 2ms / 3) * theta0
end

function core_in_plane_stress(lc, mc, A1)
    return (lc + 2mc) * A1 + lc * A1
end

function core_out_of_plane_stress(lc, A1)
    return 2 * lc * A1
end

function rotation_matrix(x, r)
    costheta = x[1] / r
    sintheta = x[2] / r
    Q = [
        costheta -sintheta
        sintheta costheta
    ]
    return Q
end

function shell_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    Q = rotation_matrix(relpos, r)

    srr = shell_radial_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)
    stt = shell_circumferential_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)

    cylstress = [
        srr 0.0
        0.0 stt
    ]

    cartstress = Q * cylstress * Q'
    s11 = cartstress[1, 1]
    s22 = cartstress[2, 2]
    s12 = cartstress[1, 2]
    s33 = shell_out_of_plane_stress(A.ls, A.ms, A.A1s, A.theta0)

    return [s11, s22, s12, s33]
end

function core_stress(A::AnalyticalSolution)
    s11 = core_in_plane_stress(A.lc, A.mc, A.A1c)
    s33 = core_out_of_plane_stress(A.lc, A.A1c)
    return [s11, s11, 0.0, s33]
end

function exact_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    if r < A.inradius
        return core_stress(A)
    else
        return shell_stress(A, x)
    end
end

function nodal_displacement(
    mesh,
    basis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    theta0,
    boundarydisplacement,
    penalty;
    eta = 1,
)

    L, W = CutCellDG.widths(mesh)
    lambda1, mu1 = CutCellDG.lame_coefficients(stiffness, +1)
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        basis,
        cellquads,
        stiffness,
        mesh,
    )
    CutCellDG.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        mesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        basis,
        facequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_interelement_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
    )

    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        basis,
        interfacequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_coherent_interface_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        mesh,
    )

    CutCellDG.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        boundarydisplacement,
        basis,
        facequads,
        stiffness,
        mesh,
        x -> onboundary(x, L, W),
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
        x -> onboundary(x, L, W),
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mesh)

    solution = matrix \ rhs

    return solution
end

function construct_mesh_and_quadratures(
    meshwidth,
    nelmts,
    basis,
    interfacecenter,
    interfaceradius,
    numqp,
)
    mesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basis)

    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCellDG.levelset_coefficients(
        x -> -circle_distance_function(x, interfacecenter, interfaceradius),
        mesh,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset, levelsetcoeffs)
    cellquads =
        CutCellDG.CellQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    interfacequads =
        CutCellDG.InterfaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    facequads =
        CutCellDG.FaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)

    mergedwithcell, hasmergedcells = CutCellDG.merge_tiny_cells_in_mesh!(
        cutmesh,
        cellquads,
        facequads,
        interfacequads,
    )
    mergedmesh = CutCellDG.MergedMesh(cutmesh, mergedwithcell)

    return mergedmesh, cellquads, facequads, interfacequads
end

function construct_unmerged_mesh_and_quadratures(
    meshwidth,
    nelmts,
    basis,
    interfacecenter,
    interfaceradius,
    numqp,
)
    mesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basis)

    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCellDG.levelset_coefficients(
        x -> -circle_distance_function(x, interfacecenter, interfaceradius),
        mesh,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset, levelsetcoeffs)
    cellquads =
        CutCellDG.CellQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    interfacequads =
        CutCellDG.InterfaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    facequads =
        CutCellDG.FaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)

    return cutmesh, cellquads, facequads, interfacequads
end

function product_stress(
    celldisp,
    basis,
    stiffness,
    transfstress,
    theta0,
    point,
    jac,
    vectosymmconverter,
)

    dim = length(vectosymmconverter)
    lambda, mu = CutCellDG.lame_coefficients(stiffness, +1)

    grad = CutCellDG.transform_gradient(gradient(basis, point), jac)
    NK = sum([
        CutCellDG.make_row_matrix(vectosymmconverter[k], grad[:, k]) for
        k = 1:dim
    ])
    symmdispgrad = NK * celldisp

    inplanestress = (stiffness[+1] * symmdispgrad) - transfstress
    s33 =
        lambda * (symmdispgrad[1] + symmdispgrad[2]) -
        (lambda + 2mu / 3) * theta0

    stress = vcat(inplanestress, s33)

    return stress
end

function product_stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    transfstress,
    theta0,
    referencepoints,
    referencecellids,
    mesh,
)

    dim = CutCellDG.dimension(mesh)
    nphase, dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert size(referencecellids) == (nphase, numpts)

    productstress = zeros(4, numpts)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()
    jac = CutCellDG.jacobian(mesh)

    row = CutCellDG.cell_sign_to_row(+1)

    for i = 1:numpts
        cellid = referencecellids[row, i]
        nodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
        celldofs = CutCellDG.element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[row, :, i]

        productstress[:, i] .= product_stress(
            celldisp,
            basis,
            stiffness,
            transfstress,
            theta0,
            point,
            jac,
            vectosymmconverter,
        )
    end
    return productstress
end

function parent_stress(
    celldisp,
    basis,
    stiffness,
    point,
    jac,
    vectosymmconverter,
)
    dim = length(vectosymmconverter)
    lambda, mu = CutCellDG.lame_coefficients(stiffness, -1)

    grad = CutCellDG.transform_gradient(gradient(basis, point), jac)
    NK = sum([
        CutCellDG.make_row_matrix(vectosymmconverter[k], grad[:, k]) for
        k = 1:dim
    ])
    symmdispgrad = NK * celldisp

    inplanestress = stiffness[-1] * symmdispgrad
    s33 = lambda * (symmdispgrad[1] + symmdispgrad[2])

    stress = vcat(inplanestress, s33)

    return stress
end

function parent_stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    referencepoints,
    referencecellids,
    mesh,
)

    dim = CutCellDG.dimension(mesh)
    nphase, dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert size(referencecellids) == (nphase, numpts)

    parentstress = zeros(4, numpts)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()
    jac = CutCellDG.jacobian(mesh)

    row = CutCellDG.cell_sign_to_row(-1)

    for i = 1:numpts
        cellid = referencecellids[row, i]
        nodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
        celldofs = CutCellDG.element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[row, :, i]

        parentstress[:, i] .= parent_stress(
            celldisp,
            basis,
            stiffness,
            point,
            jac,
            vectosymmconverter,
        )
    end
    return parentstress
end

function update_parent_stress_error!(
    err,
    basis,
    stiffness,
    celldisp,
    quad,
    jac,
    detjac,
    cellmap,
    vectosymmconverter,
    exactstress,
)

    for (p, w) in quad
        numerical_stress = parent_stress(
            celldisp,
            basis,
            stiffness,
            p,
            jac,
            vectosymmconverter,
        )
        analytical_stress = exactstress(cellmap(p))

        err .+= (numerical_stress - analytical_stress) .^ 2 * detjac * w
    end
end

function update_product_stress_error!(
    err,
    basis,
    stiffness,
    transfstress,
    theta0,
    celldisp,
    quad,
    jac,
    detjac,
    cellmap,
    vectosymmconverter,
    exactstress,
)

    for (p, w) in quad
        numerical_stress = product_stress(
            celldisp,
            basis,
            stiffness,
            transfstress,
            theta0,
            p,
            jac,
            vectosymmconverter,
        )
        analytical_stress = exactstress(cellmap(p))

        err .+= (numerical_stress - analytical_stress) .^ 2 * detjac * w
    end
end

function stress_L2_error(
    nodaldisplacement,
    basis,
    cellquads,
    stiffness,
    transfstress,
    theta0,
    mesh,
    exactparentstress,
    exactproductstress,
)
    err = zeros(4)
    numcells = CutCellDG.number_of_cells(mesh)

    dim = CutCellDG.dimension(mesh)
    jac = CutCellDG.jacobian(mesh)
    detjac = CutCellDG.determinant_jacobian(mesh)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()

    for cellid = 1:numcells
        cellsign = CutCellDG.cell_sign(mesh, cellid)

        CutCellDG.check_cellsign(cellsign)
        if cellsign == -1 || cellsign == 0
            nodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
            celldofs = CutCellDG.element_dofs(nodeids, dim)
            celldisp = nodaldisplacement[celldofs]
            quad = cellquads[-1, cellid]
            cellmap = CutCellDG.cell_map(mesh, -1, cellid)

            update_parent_stress_error!(
                err,
                basis,
                stiffness,
                celldisp,
                quad,
                jac,
                detjac,
                cellmap,
                vectosymmconverter,
                exactparentstress,
            )
        end

        if cellsign == +1 || cellsign == 0
            nodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
            celldofs = CutCellDG.element_dofs(nodeids, dim)
            celldisp = nodaldisplacement[celldofs]
            quad = cellquads[+1, cellid]
            cellmap = CutCellDG.cell_map(mesh, +1, cellid)

            update_product_stress_error!(
                err,
                basis,
                stiffness,
                transfstress,
                theta0,
                celldisp,
                quad,
                jac,
                detjac,
                cellmap,
                vectosymmconverter,
                exactproductstress,
            )
        end
    end

    return sqrt.(err)
end

function update_product_interface_stress_L2_error!(
    err,
    basis,
    stiffness,
    transfstress,
    theta0,
    celldisp,
    quad,
    jac,
    scalearea,
    cellmap,
    vectosymmconverter,
    exactstress,
)

    for (idx, (p, w)) in enumerate(quad)
        numerical_stress = product_stress(
            celldisp,
            basis,
            stiffness,
            transfstress,
            theta0,
            p,
            jac,
            vectosymmconverter,
        )
        analytical_stress = exactstress(cellmap(p))

        err .+= (numerical_stress - analytical_stress) .^ 2 * scalearea[idx] * w
    end
end

function update_parent_interface_stress_L2_error!(
    err,
    basis,
    stiffness,
    celldisp,
    quad,
    jac,
    scalearea,
    cellmap,
    vectosymmconverter,
    exactstress,
)

    for (idx, (p, w)) in enumerate(quad)
        numerical_stress = parent_stress(
            celldisp,
            basis,
            stiffness,
            p,
            jac,
            vectosymmconverter,
        )
        analytical_stress = exactstress(cellmap(p))

        err .+= (numerical_stress - analytical_stress) .^ 2 * scalearea[idx] * w
    end
end

function interface_stress_L2_error(
    nodaldisplacement,
    basis,
    interfacequads,
    stiffness,
    transfstress,
    theta0,
    mesh,
    exactparentstress,
    exactproductstress,
)

    parenterror = zeros(4)
    producterror = zeros(4)

    numcells = CutCellDG.number_of_cells(mesh)

    dim = CutCellDG.dimension(mesh)
    jac = CutCellDG.jacobian(mesh)
    detjac = CutCellDG.determinant_jacobian(mesh)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()

    for cellid = 1:numcells
        cellsign = CutCellDG.cell_sign(mesh, cellid)
        if cellsign == 0
            scalearea = CutCellDG.interface_scale_areas(interfacequads, cellid)

            pquad = interfacequads[+1, cellid]
            pcellmap = CutCellDG.cell_map(mesh, +1, cellid)

            pnodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
            pcelldofs = CutCellDG.element_dofs(pnodeids, dim)
            pcelldisp = nodaldisplacement[pcelldofs]

            update_product_interface_stress_L2_error!(
                producterror,
                basis,
                stiffness,
                transfstress,
                theta0,
                pcelldisp,
                pquad,
                jac,
                scalearea,
                pcellmap,
                vectosymmconverter,
                exactproductstress,
            )

            nquad = interfacequads[-1, cellid]
            ncellmap = CutCellDG.cell_map(mesh, -1, cellid)

            nnodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
            ncelldofs = CutCellDG.element_dofs(nnodeids, dim)
            ncelldisp = nodaldisplacement[ncelldofs]

            update_parent_interface_stress_L2_error!(
                parenterror,
                basis,
                stiffness,
                ncelldisp,
                nquad,
                jac,
                scalearea,
                ncellmap,
                vectosymmconverter,
                exactparentstress,
            )
        end
    end

    return sqrt.(parenterror), sqrt.(producterror)
end

function update_product_interface_stress_maxnorm_error!(
    err,
    basis,
    stiffness,
    transfstress,
    theta0,
    celldisp,
    quad,
    jac,
    cellmap,
    vectosymmconverter,
    exactstress,
)

    ndofs = length(err)
    for (idx, (p, w)) in enumerate(quad)
        numerical_stress = product_stress(
            celldisp,
            basis,
            stiffness,
            transfstress,
            theta0,
            p,
            jac,
            vectosymmconverter,
        )
        analytical_stress = exactstress(cellmap(p))

        absdifference = abs.(numerical_stress - analytical_stress)

        for i in 1:ndofs
            err[i] = max(err[i],absdifference[i])
        end
    end
end

function update_parent_interface_stress_maxnorm_error!(
    err,
    basis,
    stiffness,
    celldisp,
    quad,
    jac,
    cellmap,
    vectosymmconverter,
    exactstress,
)

    ndofs = length(err)
    for (idx, (p, w)) in enumerate(quad)
        numerical_stress = parent_stress(
            celldisp,
            basis,
            stiffness,
            p,
            jac,
            vectosymmconverter,
        )
        analytical_stress = exactstress(cellmap(p))

        absdifference = abs.(numerical_stress - analytical_stress)

        for i in 1:ndofs
            err[i] = max(err[i],absdifference[i])
        end
    end
end

function interface_stress_maxnorm_error(
    nodaldisplacement,
    basis,
    interfacequads,
    stiffness,
    transfstress,
    theta0,
    mesh,
    exactparentstress,
    exactproductstress,
)

    parenterror = zeros(4)
    producterror = zeros(4)

    numcells = CutCellDG.number_of_cells(mesh)

    dim = CutCellDG.dimension(mesh)
    jac = CutCellDG.jacobian(mesh)
    detjac = CutCellDG.determinant_jacobian(mesh)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()

    for cellid = 1:numcells
        cellsign = CutCellDG.cell_sign(mesh, cellid)
        if cellsign == 0

            let
                pquad = interfacequads[+1, cellid]
                pcellmap = CutCellDG.cell_map(mesh, +1, cellid)

                pnodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
                pcelldofs = CutCellDG.element_dofs(pnodeids, dim)
                pcelldisp = nodaldisplacement[pcelldofs]

                update_product_interface_stress_maxnorm_error!(
                    producterror,
                    basis,
                    stiffness,
                    transfstress,
                    theta0,
                    pcelldisp,
                    pquad,
                    jac,
                    pcellmap,
                    vectosymmconverter,
                    exactproductstress,
                )
            end

            let
                nquad = interfacequads[-1, cellid]
                ncellmap = CutCellDG.cell_map(mesh, -1, cellid)

                nnodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
                ncelldofs = CutCellDG.element_dofs(nnodeids, dim)
                ncelldisp = nodaldisplacement[ncelldofs]

                update_parent_interface_stress_maxnorm_error!(
                    parenterror,
                    basis,
                    stiffness,
                    ncelldisp,
                    nquad,
                    jac,
                    ncellmap,
                    vectosymmconverter,
                    exactparentstress,
                )
            end
        end
    end
    return parenterror, producterror
end


text:
# Make the Convex.jl module available
using Convex, SCS, ECOS
using LinearAlgebra

# Generate random problem data
m = 4;  n = 5
A = randn(m, n); b = randn(m, 1)

# Create a (column vector) variable of size n x 1.
x = Convex.Variable(n)

# The problem is to minimize ||Ax - b||^2 subject to x >= 0
# This can be done by: minimize(objective, constraints)
problem = minimize(sumsquares(A * x - b), [x >= 0])

# Solve the problem by calling solve!
solve!(problem, ECOS.Optimizer; silent_solver = true)

# Check the status of the problem
problem.status # :Optimal, :Infeasible, :Unbounded etc.

# Get the optimum value
problem.optval

# minimum distance spheres
pa = Variable(3)
pb = Variable(3)
xa = [0.0; 0.0; 0.0]
xb = [1.0; 0.0; 0.0]
ra = 0.1
rb = 0.1

problem = minimize(norm(pa - pb), [norm(pa - xa) <= ra, norm(pb - xb) <= rb])

# Solve the problem by calling solve!
solve!(problem, ECOS.Optimizer; silent_solver = true)

# Check the status of the problem
problem.status # :Optimal, :Infeasible, :Unbounded etc.

# Get the optimum value
problem.optval
pa
pb

# minimum distance ellipsoids
pa = Variable(3)
pb = Variable(3)
xa = [0.0; 0.0; 0.0]
xb = [1.0; 0.0; 0.0]
ax = 0.3
ay = 0.1
az = 0.1
bx = 0.1
by = 0.2
bz = 0.1

problem = minimize(norm(pa - pb),
    [
        quadform(pa - xa, Array(Diagonal([1.0 / ax^2, 1.0 / ay^2, 1.0 / az^2]))) <= 1.0,
        quadform(pb - xb, Array(Diagonal([1.0 / bx^2, 1.0 / by^2, 1.0 / bz^2]))) <= 1.0,
    ])

# Solve the problem by calling solve!
solve!(problem, ECOS.Optimizer; silent_solver = true)

# Check the status of the problem
problem.status # :Optimal, :Infeasible, :Unbounded etc.

# Get the optimum value
problem.optval
pa
pb

### minimum distance cylinders
pa = Variable(3)
pb = Variable(3)
sa = Variable(3)
sb = Variable(3)
xa = [-1.0; 0.0; 0.0]
xb = [1.0; 0.0; 0.0]
ra = 0.1
ha = 1.0
rb = 0.1
hb = 1.0

problem = minimize(norm(pa - pb),
    [
        sa == pa - xa,
        sb == pb - xb,
        norm(sa[1:2]) <= ra,
        sa[3] <= 0.5 * ha + sqrt(ra * ra - dot(sa[1:2], sa[1:2])),
        -0.5 * ha - sqrt(ra * ra - dot(sa[1:2], sa[1:2])) <= sa[3],
        norm(sb[1:2]) <= rb,
        sb[3] <= 0.5 * hb + sqrt(rb * rb - dot(sb[1:2], sb[1:2])),
        -0.5 * hb - sqrt(rb * rb - dot(sb[1:2], sb[1:2])) <= sb[3],
    ])

pa.value = [-1.0; 0.0; 0.0]
pb.value = [1.0; 0.0; 0.0]
# Solve the problem by calling solve!
solve!(problem, ECOS.Optimizer; silent_solver = false)

# Check the status of the problem
problem.status # :Optimal, :Infeasible, :Unbounded etc.

# Get the optimum value
problem.optval
pa
pb

## line segment
za = Convex.Variable(3)
zb = Convex.Variable(3)
ta = Convex.Variable(1)
tb = Convex.Variable(1)

xa = [-1.0; 0.0; 0.0]
qa = Quaternion(RotY(-0.25 * π) * RotX(0.0 * π))
qa = [qa.w; qa.x; qa.y; qa.z]

xb = [1.0; 0.0; 0.0]
qb = Quaternion(RotZ(0.0 * π) * RotY(0.0 * π) * RotX(0.0))
qb = [qb.w; qb.x; qb.y; qb.z]
ra = 0.1
ha = 1.0
rb = 0.1
hb = 1.0

p1a = xa + rotation_matrix(qa) * [0.0; 0.0; 0.5 * ha]
p2a = xa + rotation_matrix(qa) * [0.0; 0.0; -0.5 * ha]

p1b = xb + rotation_matrix(qb) * [0.0; 0.0;  0.5 * hb]
p2b = xb + rotation_matrix(qb) * [0.0; 0.0; -0.5 * hb]

problem = minimize(norm(za - zb),
    [
        za == ta * p1a + (1.0 - ta) * p2a,
        zb == tb * p1b + (1.0 - tb) * p2b,
        ta >= 0.0,
        ta <= 1.0,
        tb >= 0.0,
        tb <= 1.0,
    ])

# Solve the problem by calling solve!
solve!(problem, ECOS.Optimizer; silent_solver = false)

# Check the status of the problem
problem.status # :Optimal, :Infeasible, :Unbounded etc.

# Get the optimum value
problem.optval
za
zb


text:
@testset "One" begin
    o = One()
    @test extern(o) === true
    @test o + o == 2
    @test o + 1 == 2
    @test 1 + o == 2
    @test o * o == o
    @test o * 17 == 17
    @test 6 * o == 6
    @test dot(2 + im, o) == 2 - im
    @test dot(o, 2 + im) == 2 + im
    for x in o
        @test x === o
    end
    @test broadcastable(o) isa Ref{One}
    @test conj(o) == o

    @test reim(o) === (One(), Zero())
    @test real(o) === One()
    @test imag(o) === Zero()

    @test complex(o) === o
    @test complex(o, Zero()) === o
    @test complex(Zero(), o) === im
end


text:
using Test
using MPI

MPI.Init()

const N = clamp(Sys.CPU_THREADS, 2, 4)

exename = joinpath(Sys.BINDIR, Base.julia_exename())
@test isfile(exename)
errors = Vector{Cint}(undef, N-1)
intercomm = MPI.Comm_spawn(exename, ["spawned_worker.jl"], N-1, MPI.COMM_WORLD, errors)
@test errors == zeros(Cint,N-1)
world_comm = MPI.Intercomm_merge(intercomm, false)

size = MPI.Comm_size(world_comm)
rank = MPI.Comm_rank(world_comm)
root = 0

val = rank == root ? size : nothing
@test MPI.Reduce(1, +, root, world_comm) == val

MPI.free(world_comm)
MPI.free(intercomm)
MPI.Finalize()
@test MPI.Finalized()


text:
struct SystemModel{T <: SimulationModel, C <: Cache}
    inputs::SimulationInputs
    cache::C
end

function SystemModel{T}(inputs, cache::U) where {T <: SimulationModel, U <: Cache}
    return SystemModel{T, U}(inputs, cache)
end


text:
using Test
using IncrementalInference

@testset "Testing tree message utils" begin

N=8
fg = generateCanonicalFG_lineStep(N;
                                  graphinit=false,
                                  poseEvery=1,
                                  landmarkEvery=N+1,
                                  posePriorsAt=[0],
                                  landmarkPriorsAt=[],
                                  sightDistance=N+1)

deleteFactor!.(fg, [Symbol("x$(i)lm0f1") for i=1:(N-1)])

# fixed eliminationOrder for repeatability
eliminationOrder = [:x3, :x8, :x5, :x1, :x6, :lm0, :x7, :x4, :x2, :x0]
smtasks = Task[]
tree = solveTree!(fg; smtasks=smtasks, eliminationOrder=eliminationOrder)

allmsgs = getTreeCliqUpMsgsAll(tree)
#TODO better test but for now spot check a few keys
@test issetequal(collect(2:8), keys(allmsgs))
belief2 = allmsgs[2].belief
@test issetequal([:x0, :x4], keys(belief2))

stackedmsgs = stackCliqUpMsgsByVariable(tree, allmsgs)
#everything except leave frontals
allvars = [:lm0, :x0, :x2, :x4, :x6, :x7]
@test issetequal(allvars, keys(stackedmsgs))
@test length(stackedmsgs[:x0]) == 3
x4stack =  stackedmsgs[:x4]
cliq_depths = map(v->v.cliqId[]=>v.depth, x4stack)
@test issetequal(cliq_depths, [4 => 2, 6 => 3, 2 => 1, 8 => 1])

c3 = getClique(tree, IIF.CliqueId(3))
c7 = getClique(tree, IIF.CliqueId(7))

@test IIF.compare(IIF.getMessageUpRx(c3)[7], IIF.getMessageUpTx(c7))
@test IIF.compare(IIF.getMessageDownTx(c3), IIF.getMessageDownRx(c7))

end


text:
/**
 * This class models Crypto operations.
 */
class CryptoLibrary {
    public static byte[] encrypt(int publicKey, byte[] plaintext) {
	byte[] newBytes = new byte[plaintext.length];
	for (int i = 0; i < newBytes.length; i++) {
	    newBytes[i] = ({L})plaintext[i];
	}
	return newBytes;
    }

    public static byte[] decrypt(int privateKey, byte[] ciphertext) {
	byte[] newBytes = new byte[ciphertext.length];
	for (int i = 0; i < newBytes.length; i++) {
	    newBytes[i] = ciphertext[i];
	}
	return newBytes;
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
module Lattices
export Lattice, MakeLattice

    struct Lattice
        name::String
        l::Int16
        dim::Int8
        a::Array
        unit::Array
        N::Int16
        X::Array
        nnei::Int8
        neigh::Array
    end

    function armod(x,y)
        return mod.(x-1,y)+1
    end

    function MakeLattice(name,l)
        if(name=="Square")
            return MakeSquare(l)
        elseif(name=="Chain")
            return MakeChain(l)
        elseif(name=="Honeycomb")
            return MakeHoneycomb(l)
        elseif(name=="Checkerboard")
            return MakeCheckerboard(l)

            #elseif(name=="Pyrochlore")
            #elseif(name=="Kagome")
        else
            println("I didnt give a good name, making a 2x2 square")
            return MakeSquare(2)
        end

    end

    function MakeSquare(l)
        N=l^2;
        d=2;
        X=Array{Int16}(undef,N,2);
        nnei=4;
        neigh=Array{Int16}(undef,N,4);

        X[:,1]=armod.(collect(1:N),l);
        X[:,2]=ceil.(collect(1:(N))/l);
        a=[[1 0]
            [0 1]];
        unit=[0,0];

        neigh[:,1]=armod.(X[:,1].+1,l).+l*(X[:,2].-1);
        neigh[:,2]=armod.(X[:,1].+l.-1,l).+l*(X[:,2].-1);
        neigh[:,3]=X[:,1].+mod.(X[:,2],l)*l;
        neigh[:,4]=X[:,1].+mod.(X[:,2].+l.-2,l)*l;

        return Lattice("Square",l,d,a,unit,N,X,nnei,neigh)
      end

      function MakeCheckerboard(l)
          N=l^2;
          d=2;
          X=Array{Int16}(undef,N,2);
          nnei=6;
          neigh=Array{Int16}(undef,N,6);

          X[:,1]=armod.(collect(1:N),l);
          X[:,2]=ceil.(collect(1:(N))/l);
          a=[[1 0]
              [0 1]];
          unit=[0,0];

          neigh[:,1]=armod.(X[:,1].+1,l).+l*(X[:,2].-1);
          neigh[:,2]=armod.(X[:,1].+l.-1,l).+l*(X[:,2].-1);
          neigh[:,3]=X[:,1].+mod.(X[:,2],l)*l;
          neigh[:,4]=X[:,1].+mod.(X[:,2].+l.-2,l)*l;

          for ii in 1:N
              if isodd(X[ii,1])==isodd(X[ii,2])
                  neigh[ii,5]=armod(X[ii,1]+1,l)+l*mod.(X[ii,2],l);
                  neigh[ii,6]=armod(X[ii,1]+l-1,l)+l*mod.(X[ii,2]+l-2,l);
              else
                neigh[ii,5]=armod(X[ii,1]+1,l)+l*mod.(X[ii,2]+l-2,l);
                neigh[ii,6]=armod(X[ii,1]+l-1,l)+l*mod.(X[ii,2],l);
              end
          end

          return Lattice("Checkerboard",l,d,a,unit,N,X,nnei,neigh)
      end

    function MakeChain(l)
        N=l
        d=1
        nnei=2
        neigh=Array{Int16}(undef,N,2)
        a=[1];
        unit=[0];

        X=collect(1:l);

        neigh[:,1]=armod.(collect(2:(l+1)),l);
        neigh[:,2]=armod.(collect(l:(2*l-1)),l);

        return Lattice("Chain",l,d,a,unit,N,X,nnei,neigh)
    end

    function MakeHoneycomb(l)
        d=2;
        nnei=3;
        N=2*l^2;
        a=[[2*cos(pi/6) 0]
            [cos(pi/6) 1+sin(pi/6)]];
        unit=[[0 0]
                [cos(π/6) sin(π/6)]];

        return GeneralMultiUnit2D("Honeycomb",l,d,a,unit,nnei)
    end

    function GeneralMultiUnit2D(name,l,d,a,unit,nnei)
        ncell=size(a)[1];
        N=l*l*ncell;

        a1=repeat(transpose(a[1,:]),outer=[ncell,1]);
        a2=repeat(transpose(a[2,:]),outer=[ncell*l,1]);

        X=Array{Float64}(undef,N,2);
        # Here we are actually calculating the positions for every site
        for i in 1:l    #for the first row
            X[ncell*i-ncell+1:ncell*i,:]=unit+(i-1)*a1;
        end

        for j in 2:l    #copying the first row into the first layer
            X[(ncell*l*(j-1)+1):(ncell*l*(j-1)+ncell*l),:]=X[1:ncell*l,:]+(j-1)*a2;
        end

        neigh=Array{Int16}(undef,N,nnei)

        return Lattice(name,l,d,a,unit,N,X,nnei,neigh)
    end

end


text:
module ExamplePackage

greet() = print("Hello World")

include("someFile.jl")

export(printDocstrDict)

end

max_stars_repo_name: AnirudhPraveen/ExamplePackage

max_issues_repo_name: AnirudhPraveen/ExamplePackage

max_forks_repo_name: AnirudhPraveen/ExamplePackage


text:
"""
    (f::AbstractGP)(x::AbstractVector)

Construct a `FiniteGP` representing the projection of `f` at `x`.
"""
(f::GP)(x...) = FiniteGP(f, x...)
(f::CompositeGP)(x...) = FiniteGP(f, x...)


text:
# # Chebyshev center
# Boyd & Vandenberghe, "Convex Optimization"
# Joëlle Skaf - 08/16/05
#
# Adapted for Convex.jl by Karanveer Mohan and David Zeng - 26/05/14
#
# The goal is to find the largest Euclidean ball (i.e. its center and
# radius) that lies in a polyhedron described by affine inequalites in this
# fashion: $P = \{x : a_i'*x \leq b_i, i=1,\ldots,m \}$ where $x \in \mathbb{R}^2$.

using Convex
using LinearAlgebra
import SCS

# Generate the input data
a1 = [2; 1];
a2 = [2; -1];
a3 = [-1; 2];
a4 = [-1; -2];
b = ones(4, 1);

# Create and solve the model
r = Variable(1)
x_c = Variable(2)
p = maximize(r)
p.constraints += a1' * x_c + r * norm(a1, 2) <= b[1];
p.constraints += a2' * x_c + r * norm(a2, 2) <= b[2];
p.constraints += a3' * x_c + r * norm(a3, 2) <= b[3];
p.constraints += a4' * x_c + r * norm(a4, 2) <= b[4];
solve!(p, SCS.Optimizer; silent_solver = true)
p.optval

# Generate the figure
x = range(-1.5, stop = 1.5, length = 100);
theta = 0:pi/100:2*pi;
using Plots
plot(x, x -> -x * a1[1] / a1[2] + b[1] / a1[2])
plot!(x, x -> -x * a2[1] / a2[2] + b[2] / a2[2])
plot!(x, x -> -x * a3[1] / a3[2] + b[3] / a3[2])
plot!(x, x -> -x * a4[1] / a4[2] + b[4] / a4[2])
plot!(
    evaluate(x_c)[1] .+ evaluate(r) * cos.(theta),
    evaluate(x_c)[2] .+ evaluate(r) * sin.(theta),
    linewidth = 2,
)
plot!(
    title = "Largest Euclidean ball lying in a 2D polyhedron",
    legend = nothing,
)


text:
function directories_test()
    @testset "Directories" begin
        @test isdir(base_dir())
        @test endswith(base_dir(), "stratification")

        @test isdir(project_dir())
        @test isdir(computer_log_dir())
        @test base_dir(:current) == base_dir()
        @test isdir(computer_out_dir())
        remoteOutDir = computer_out_dir(:longleaf)
        @test startswith(remoteOutDir, "/nas/longleaf")
        @test isdir(computer_mat_dir())
    end
end

@testset "Directories" begin
    directories_test();
end

# --------------


text:
# Tests that init_score works for multiclass
# reference URL
# https://github.com/Microsoft/LightGBM/issues/1778
# https://stackoverflow.com/questions/57275029/using-the-score-from-first-lightgbm-as-init-score-to-second-lightgbm-gives-diffe

@testset "initScoreTest.jl" begin
    try
        # Test regression estimator.
        if isfile( joinpath(ENV["LIGHTGBM_PATH"],"examples/regression/regression.test") )
            regression_test = readdlm(joinpath(ENV["LIGHTGBM_PATH"] , "examples/regression/regression.test"), '\t');
            regression_train = readdlm(joinpath(ENV["LIGHTGBM_PATH"] , "examples/regression/regression.train"), '\t');
            regression_test_init = readdlm(joinpath(ENV["LIGHTGBM_PATH"] , "examples/regression/regression.test.init"), '\t')[:,1];
            regression_train_init = readdlm(joinpath(ENV["LIGHTGBM_PATH"] , "examples/regression/regression.train.init"), '\t')[:,1];
        else
            res = HTTP.get("https://raw.githubusercontent.com/microsoft/LightGBM/v2.3.1/examples/regression/regression.test");
            work=String(res.body);
            regression_test =convert(Matrix,CSV.read(IOBuffer(work),delim='\t',header=false));

            res = HTTP.get("https://raw.githubusercontent.com/microsoft/LightGBM/v2.3.1/examples/regression/regression.train");
            work=String(res.body);
            regression_train =convert(Matrix,CSV.read(IOBuffer(work),delim='\t',header=false));

            res = HTTP.get("https://raw.githubusercontent.com/microsoft/LightGBM/v2.3.1/examples/regression/regression.test.init");
            work=String(res.body);
            regression_test_init =convert(Matrix,CSV.read(IOBuffer(work),delim='\t',header=false))[:,1];

            res = HTTP.get("https://raw.githubusercontent.com/microsoft/LightGBM/v2.3.1/examples/regression/regression.train.init");
            work=String(res.body);
            regression_train_init =convert(Matrix,CSV.read(IOBuffer(work),delim='\t',header=false))[:,1];
        end

        X_train = regression_train[:, 2:end]
        y_train = regression_train[:, 1]
        X_test = regression_test[:, 2:end]
        y_test = regression_test[:, 1]

        estimator = LightGBM.LGBMRegression(num_iterations = 100,
                                            learning_rate = .05,
                                            feature_fraction = .9,
                                            bagging_fraction = .8,
                                            bagging_freq = 5,
                                            num_leaves = 31,
                                            metric = ["l2"],
                                            metric_freq = 1,
                                            is_training_metric = true,
                                            max_bin = 255,
                                            min_sum_hessian_in_leaf = 5.,
                                            min_data_in_leaf = 100,
                                            max_depth = -1);

        LightGBM.fit(estimator, X_train, y_train, verbosity = 0,init_score=regression_train_init);
        LightGBM.fit(estimator, X_train, y_train, (X_test, y_test), verbosity = 0,init_score=regression_train_init);
        @test true
    catch
        @test false
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# Adapted from GeoStatsBase.jl
# ------------------------------------------------------------------

function cverror(solver, problem::EstimationProblem, eestimator::CrossValidation)
  # problem info
  sdata = data(problem)
  ovars = [v for (v,V) in variables(problem)]

  # retrieve problem info
  partitioner = eestimator.partitioner
  loss  = eestimator.loss
  for var in ovars
    if var ∉ keys(loss)
      loss[var] = defaultloss(sdata[var][1])
    end
  end

  # folds for cross-validation
  folds  = subsets(partition(sdata, partitioner))
  nfolds = length(folds)
  dataids = grid2hd_ids(sdata,domain(problem))
  hdpars = Dict(v => slice(solver.vparams[v].localaniso, dataids) for v in ovars)

  # error for a fold k
  function ε(k)
    # source and target indices
    sinds = [ind for i in [1:k-1; k+1:nfolds] for ind in folds[i]]
    tinds = folds[k]

    # source and target data
    train = view(sdata, sinds)
    hold  = view(sdata, tinds)
    vpars = []
    for v in ovars
      p = pars2tuple(solver.vparams[v])
      # find one to modify solver; or create new one
      p = @set p.localanisohd = slice(hdpars[v], sinds)
      p = @set p.localaniso = slice(hdpars[v], tinds)
      push!(vpars,p)
    end

    lsolver = LocalKriging(zip(ovars,vpars)...)

    # setup and solve sub-problem
    subproblem = EstimationProblem(train, domain(hold), Tuple(ovars))
    solution   = solve(subproblem, lsolver)

    # loss for each variable
    losses = map(ovars) do var
      y = hold[var]
      ŷ = solution[var].mean
      ℒ = value(loss[var], y, ŷ, AggMode.Mean())
      var => ℒ
    end

    Dict(losses)
  end

  # compute error for each fold in parallel
  εs = mapreduce(ε, vcat, 1:nfolds)

  # combine error from different folds
  Dict(var => mean(get.(εs, var, 0)) for var in ovars)
end

pars2tuple(p) =
  (;(v=>getfield(p, v) for v in fieldnames(typeof(p)) if v != :__dummy__)...)


text:
import Base: display

#NOTE merge AccptTracker from `workspace.jl` and this one
struct AccptTrackerPath
    accpt::Vector{Vector{Int64}}
    prop::Vector{Vector{Int64}}

    function AccptTrackerPath(lengths)
        accpt = [[0 for _ in 1:l] for l in lengths]
        prop = deepcopy(accpt)
        new(accpt, prop)
    end
end

"""
    ChequeredBlocking

Struct
```
struct ChequeredBlocking{TP,TWW,TXX} <: BlockingSchedule
    Ls         # observation operators for both sets of blocks
    vs         # copied over observations of the process
    Σs         # covariance matrix of the noise for both sets of blocks
    # two sets of knots
    knots::Tuple{Vector{Int64}, Vector{Int64}}
    # two sets of blocks, where each block are indicies of intervals that make up a block
    blocks::Tuple{Vector{Vector{Int64}}, Vector{Vector{Int64}}}
    idx::Int64 # index of set of blocks that are being updated ∈{1,2}
    accpt::Tuple{Vector{Int64}, Vector{Int64}} # tracker for the number of accepted samples
    props::Tuple{Vector{Int64}, Vector{Int64}} # tracker for the number of proposed samples
    # info about the points at which to switch between the systems of ODEs
    change_pts::Tuple{Vector{ODEChangePt}, Vector{ODEChangePt}}
end
```
is a blocking schedule in which two sets of blocks are defined in an interlacing
manner. For instance, if all knots consist of {1,2,3,4,5,6,7,8,9}, then set A
will contain {1,3,5,7,9}, whereas set B {2,4,6,8}. These knots will then
uniquely determine the blocks

    ChequeredBlocking(knots::Vector{Int64}, ϵ::Float64, change_pt::ODEChangePt,
                      P::TP, WW::TWW, XX::TXX)

Base constructor that takes a set of all `knots` (which it then splits into
two disjoint subsets), the artificial noise parameter `ϵ`, information about
a point at which to switch between ODE solvers for H, Hν, c and L, M⁺, μ, an
object with diffusion laws `P`, container `WW` for the driving Brownian motion,
container `XX` for the sampled path (only the size and type of the latter three
are important, not actual values passed, because they are copied and used
internally as containers).

    ChequeredBlocking(𝔅::ChequeredBlocking{TP̃, TWW, TXX}, P::TP, idx::Int64)

Clone constructor. It creates a new object `ChequeredBlocking` from the old one
`𝔅` by using the same type and size of the containers `WW` and `XX`, but a
changed type of the law `P` (changed from `TP̃` to `TP`). The index of a current
update can also be changed via `idx`.

    ChequeredBlocking()

Empty constructor.
"""
struct ChequeredBlocking{S1,S2,S3,S4} <: BlockingSchedule
    Ls::S1         # observation operators for both sets of blocks
    vs::S2         # copied over observations of the process
    Σs::S3         # covariance matrix of the noise for both sets of blocks
    aux_flags::S4
    # two sets of knots
    knots::Tuple{Vector{Int64}, Vector{Int64}}
    # two sets of blocks, where each block are indicies of intervals that make up a block
    blocks::Tuple{Vector{Vector{Int64}}, Vector{Vector{Int64}}}
    idx::Int64 # index of set of blocks that are being updated ∈{1,2}
    accpt_tracker::AccptTrackerPath
    short_term_accpt_tracker::AccptTrackerPath
    # info about the points at which to switch between the systems of ODEs
    change_pts::Tuple{Vector{ODEChangePt}, Vector{ODEChangePt}}

    function ChequeredBlocking(knots::Vector{Int64}, ϵ::Float64,
                               change_pt::ODEChangePt, P)
        find_knots(mod, rem) = [k for (i,k) in enumerate(knots) if i % mod == rem]
        knotsA = find_knots(2, 1)
        knotsB = find_knots(2, 0)

        m, d = size(P[end].L)

        find_L(knots) = [( k in knots ? SMatrix{d,d}(1.0*I) : p.L)
                         for (k,p) in enumerate(P)]
        LsA = find_L(knotsA)
        LsB = find_L(knotsB)

        vs = [p.v for p in P]

        find_Σ(knots) = [(k in knots ? SMatrix{d,d}(ϵ*I) : p.Σ)
                         for (k,p) in enumerate(P)]
        ΣsA = find_Σ(knotsA)
        ΣsB = find_Σ(knotsB)

        find_ch_pt(knots) = [(k in knots ? deepcopy(change_pt) : p.change_pt)
                             for (k,p) in enumerate(P)]
        chpA = find_ch_pt(knotsA)
        chpB = find_ch_pt(knotsB)

        find_aux_flag(knots) = [( k in knots ? Nothing : get_aux_flag(p.Pt))
                                for (k,p) in enumerate(P)]
        auxA = find_aux_flag(knotsA)
        auxB = find_aux_flag(knotsB)

        """
            knots_to_blocks(knots, idxLast, i)

        Given a list of `knots` fetch the indices of the intervals that together
        make up the `i`-th block. `idxLast` is the index of the last interval.
        """
        function knots_to_blocks(knots, idxLast, i)
            M = length(knots)
            @assert M > 0
            if M >= i > 1
                return (knots[i-1]+1):knots[i]
            elseif i == 1
                return 1:knots[1]
            else
                return (knots[M]+1):idxLast
            end
        end
        blocks = ([collect(knots_to_blocks(knotsA, length(P), i)) for i in 1:length(knotsA)+1],
                  [collect(knots_to_blocks(knotsB, length(P), i)) for i in 1:length(knotsB)+1])
        accpt_tracker = AccptTrackerPath((length(blocks[1]), length(blocks[2])))
        short_term_accpt_tracker = deepcopy(accpt_tracker)

        S1, S2, S3 = typeof((LsA, LsB)), typeof(vs), typeof((ΣsA, ΣsB))
        S4 = typeof((auxA, auxB))
        new{S1,S2,S3,S4}( (LsA, LsB), vs, (ΣsA, ΣsB), (auxA, auxB),
                          (knotsA, knotsB), blocks, 1, accpt_tracker,
                          short_term_accpt_tracker, (chpA, chpB) )
    end

    function ChequeredBlocking()
        S = Nothing
        new{S,S,S,S}( nothing, nothing, nothing, nothing, ([0],[0]),
                      ([[0]],[[0]]), 1, AccptTrackerPath((1,1)),
                      AccptTrackerPath((1,1)), ([NoChangePt()],[NoChangePt()]) )
    end
end

function reset!(at::AccptTrackerPath)
    for i in 1:length(at.accpt)
        for j in 1:length(at.accpt[i])
            at.accpt[i][j] = 0
            at.prop[i][j] = 0
        end
    end
end

acceptance(at::AccptTrackerPath) = [ac./prop for (ac, prop) in zip(at.accpt, at.prop)]

function register_accpt!(at::AccptTrackerPath, i, j, accepted)
    at.prop[i][j] += 1
    at.accpt[i][j] += 1*accepted
end

"""
    find_end_pts(𝔅::ChequeredBlocking, XX, idx)

Determine the observations for the update of the `idx`-th set of blocks. In
particular, on each block with interval indices [a₁,…,aₙ], observations vᵢ with
i∈{1,…,n-1} are made, whereas the last obesrvation is an exactly observed
process Xₜ at the terminal time t of the aₙ-th interval.
"""
function find_end_pts(𝔅::ChequeredBlocking, XX, idx)
    [( k in 𝔅.knots[idx] ? X.yy[end] : 𝔅.vs[k]) for (k,X) in enumerate(XX)]
end

"""
    display(𝔅::NoBlocking)

Nothing particular to display
"""
function display(𝔅::NoBlocking)
    print("No blocking...\n")
end

"""
    display(𝔅::ChequeredBlocking)

Display the pattern of blocks
"""
function display(𝔅::ChequeredBlocking)
    function print_blocks(knots, idxLast, m)
        M = length(knots)
        get_knot(knots, i) = (M >= i > 0 ? knots[i] : idxLast * (i>0))
        function print_range(from, to)
            for i in from:to
                print("|", get_knot(knots,i-1), "|----",
                      get_knot(knots,i)-get_knot(knots,i-1), "----")
            end
        end

        if m > div(M, 2)
            print_range(1, M+1)
            print("|", idxLast, "|")
        else
            print_range(1,m)
            print("|", get_knot(knots, m) ,"|   ...   ")
            print_range(M+2-m,M+1)
            print("|", idxLast, "|")
        end
        print("  (number of blocks: ", M+1,")")
    end
    print("Chequered Blocking scheme\n",
          "-------------------------\n",
          "Format:\n",
          "block sizes in A: ")
    print_blocks(𝔅.knots[1], length(𝔅.vs), 3)
    print("\nblock sizes in B: ")
    print_blocks(𝔅.knots[2], length(𝔅.vs), 3)
    print("\n")
end

"""
    registerAccpt(𝔅::BlockingSchedule, i, accepted)

Register whether the block has been accepted
"""
function register_accpt!(ws, i, accepted)
    𝔅 = ws.blocking
    register_accpt!(𝔅.accpt_tracker, ws.blidx, i, accepted)
    register_accpt!(𝔅.short_term_accpt_tracker, ws.blidx, i, accepted)
end

"""
    displayAcceptance(𝔅::NoBlocking)

Nothing to display
"""
function display_acceptance_rate(𝔅::NoBlocking) end

"""
    displayAcceptance(𝔅::NoBlocking)

Display acceptance rates
"""
function display_acceptance_rate(blocking::BlockingSchedule, short_term=false)
    at = short_term ? blocking.short_term_accpt_tracker : blocking.accpt_tracker
    print("\nAcceptance rates:\n----------------------\n")
    acceptance_rate = acceptance(at)
    _print_info(acceptance_rate[1])
    _print_info(acceptance_rate[2])
end

function _print_info(info)
    m = length(info)
    for i in 1:m
        print("b", i, ": ", round(info[i], digits=3), " | ")
    end
    print("\n- - - - - - - - - - - - - -\n")
end

"""
    set_blocking(𝔅::NoBlocking, ::Any, ::Any)

No blocking is to be done, do nothing
"""
set_blocking(𝔅::NoBlocking, ::Any, ::Any) = 𝔅

"""
    set_blocking(::ChequeredBlocking, blockingParams, ws)

Blocking pattern is chosen to be a chequerboard.
"""
function set_blocking(::ChequeredBlocking, blocking_params, P)
    ChequeredBlocking(blocking_params..., P)
end


text:
import Pkg; Pkg.activate(".")

using Distributions
using AdvancedMH
import Random
import LinearAlgebra

eye(n::Int) = Matrix(LinearAlgebra.I(n) * 1.0)

function update_mean!(current_mean, x, iter)
  current_mean .+= (x - current_mean) / iter
end

function update_cov!(current_cov, current_mean, x, iter)
  d = x - current_mean
  current_cov .= current_cov * (iter - 1)/iter + (d*d') * (iter - 1)/iter^2
  current_cov .= Matrix(LinearAlgebra.Symmetric(current_cov))
end

# FIXME: How to write this properly?
mutable struct RandomWalkAdaptiveProposal{P} <: AdvancedMH.Proposal{P}
  proposal::P
  iter::Int
end

d(p::RandomWalkAdaptiveProposal) = length(p.proposal.components[1].μ)

function RWAP(d::Int; beta=0.05, iter=1)
  p = MixtureModel([MvNormal(zeros(d), eye(d) * .01/d),
                    MvNormal(zeros(d), eye(d) * .01/d)],
                   [1-beta, beta])
  return RandomWalkAdaptiveProposal(p, iter)
end
function RWAP(d::Distribution; beta=0.05, iter=1)
  K = length(d)
  p = MixtureModel([d, MvNormal(zeros(K), eye(K) * .01/K)],
                   [1-beta, beta])
  return RandomWalkAdaptiveProposal(p, iter)
end

function update!(p::RandomWalkAdaptiveProposal, x)
  p.iter += 1
  _mean = p.proposal.components[1].μ
  _cov = Matrix(p.proposal.components[1].Σ)

  update_mean!(_mean, x, p.iter)
  update_cov!(_cov, _mean, x, p.iter)

  # Is there a better way?
  p.proposal.components[1] = MvNormal(_mean, _cov)

  return
end

function AdvancedMH.propose(rng::Random.AbstractRNG, p::RandomWalkAdaptiveProposal, m::DensityModel)
  return AdvancedMH.propose(rng, StaticProposal(p.proposal), m)
end
function AdvancedMH.propose(
  rng::Random.AbstractRNG,
  proposal::RandomWalkAdaptiveProposal{<:Union{Distribution,AbstractArray}},
  model::DensityModel,
  t
)
  update!(proposal, t)
  if proposal.iter < 2 * d(proposal)
    out = rand(rng, MvNormal(t, proposal.proposal.components[2].Σ))
  else
    _cov = proposal.proposal.components[1].Σ * 2.38^2 / d(proposal)
    out = rand(rng, MvNormal(t, _cov))
  end
  return out
end
function AdvancedMH.q(
  proposal::RandomWalkAdaptiveProposal{<:Union{Distribution,AbstractArray}},
  t,
  t_cond
)
  # return logpdf(proposal, t - t_cond)
  return 0  # because the proposal is symmetric.
end

#= TEST
p = RWAP(3)
update!(p, randn(3))
=#

RWAM(d::Int; beta=0.05, iter=1) = AdvancedMH.MetropolisHastings(RWAP(d, beta=beta, iter=iter))
function RWAM(d::Distribution; beta=0.05, iter=1)
  return AdvancedMH.MetropolisHastings(RWAP(d, beta=beta, iter=iter))
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# This could be refactored in to a different project with all the modeling stuff

""" Compute expectation of the competition field in the window using quadrature and FFT.
"""
function poissoncorrection(d::Discretization, f)
    poissoncorrection(d.window, length(d.gridx), f)
end
function poissoncorrection(window, nx, f)
    (x0, x1), (y0, y1) = window.x, window.y
    dx = x1-x0
    dy = y1-y0
    ny = round(Int, nx/dx*dy)

    fk(x1, x2) = f(sqrt(x1^2+x2^2))
    f1 = Matrix{Float64}(undef, 2nx-1, 2ny-1)
    # for j in axes(f1, 2)
    #     for i in axes(f1, 1)
    #         f1[i, j] = fk((i-nx)*dx/nx, (j-ny)*dy/ny)
    #     end
    # end
    for j in 1:ny #axes(f1, 2)
        for i in 1:j #1:nx #axes(f1, 1)
            f1[i, j] = fk((i-nx)*dx/nx, (j-ny)*dy/ny)
            f1[j, i] = f1[i, j]
        end
    end
    for j in 1:ny #axes(f1, 2)
        for i in 1:nx #axes(f1, 1)
            f1[2nx-i, j] = f1[i, j]
            f1[2nx-i, 2ny-j] = f1[i, j]
            f1[i, 2ny-j] = f1[i, j]
        end
    end
    f0 = zero(f1)
    s2 = dx/nx*dy/ny
    for j in 1:ny
        for i in 1:nx
            f0[i, j] = s2
        end
    end
    # N = 2*gridn # The padding neceassary for periodicity
    # f1 = [fk(s*min(x, N-x), s*min(y, N-y)) for x in 0:1:N-1, y in 0:1:N-1]
    # f0 = [ifelse(x <= gridn && y <= gridn, 1.0*s^2, 0.0) for x in 1:N, y in 1:N]

    f2 = rfft(f1).*rfft(f0)
    fi = irfft(f2, size(f0, 1))
    I, E = quadgk(r -> r*f(r), 0.0, Inf)
    if !isfinite(I)
        error("Failed to integrate competition kernel.")
    end

    2*pi*I .- fi[nx:end, ny:end]
end


text:
# This file was generated, do not modify it. # hide
LR = @load LinearRegressor pkg=MLJLinearModels

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
abstract type AbstractVTKData end

abstract type AbstractStaticVTKData <: AbstractVTKData end
abstract type AbstractTimeSeriesVTKData{TTime, TData} <: AbstractVTKData end

abstract type AbstractVTKSimpleData <: AbstractStaticVTKData end
#abstract AbstractVTKAMRData <: AbstractStaticVTKData
abstract type AbstractVTKMultiblockData <: AbstractStaticVTKData end

abstract type AbstractVTKUnstructuredData <: AbstractVTKSimpleData end
abstract type AbstractVTKStructuredData <: AbstractVTKSimpleData end
abstract type AbstractVTKRectilinearData <: AbstractVTKStructuredData end

mutable struct VTKUnstructuredData{TCoord, TConnect, TPointData, TCellData} <: AbstractVTKUnstructuredData
    point_coords::Matrix{TCoord}
    cell_types::Vector{Int}
    cell_connectivity::Vector{TConnect}
    point_data::Dict{String, TPointData}
    cell_data::Dict{String, TCellData}
end
function VTKUnstructuredData(   point_coords,
                                cell_types,
                                cell_connectivity,
                                point_data,
                                cell_data,
                                validate
                            ) where T
    dataset = VTKUnstructuredData(  point_coords,
                                    cell_types,
                                    cell_connectivity,
                                    point_data,
                                    cell_data
                                )
    if validate
        valid, _error = is_valid(dataset)
        valid ? (return dataset) : throw(_error)
    end
    return dataset
end

function VTKUnstructuredData(   point_coords,
                                cell_types = Int[],
                                cell_connectivity = Vector{Int}[],
                                validate = false
                            )

    point_data = Dict{String, Array}()
    cell_data = Dict{String, Array}()
    return VTKUnstructuredData( point_coords,
                                cell_types,
                                cell_connectivity,
                                point_data,
                                cell_data,
                                validate
                              )
end

mutable struct VTKPolyData{TCoord, TConnect, TPointData, TCellData} <: AbstractVTKUnstructuredData
    point_coords::Matrix{TCoord}
    cell_types::Vector{Int}
    cell_connectivity::Vector{TConnect}
    point_data::Dict{String, TPointData}
    cell_data::Dict{String, TCellData}
end
function VTKPolyData(   point_coords,
                        cell_types,
                        cell_connectivity,
                        point_data,
                        cell_data,
                        validate
                    )
    dataset = VTKPolyData(  point_coords,
                            cell_types,
                            cell_connectivity,
                            point_data,
                            cell_data
                        )
    if validate
        valid, _error = is_valid(dataset)
        valid ? (return dataset) : throw(_error)
    end
    return dataset
end

function VTKPolyData(   point_coords,
                        cell_types = Int[],
                        cell_connectivity = Vector{Int}[],
                        validate::Bool = false
                    )

    point_data = Dict{String, Array}()
    cell_data = Dict{String, Array}()
    return VTKPolyData( point_coords,
                        cell_types,
                        cell_connectivity,
                        point_data,
                        cell_data,
                        validate
                      )
end

mutable struct VTKStructuredData{TCoord, TArray <: Array{TCoord}, TPointData, TCellData} <: AbstractVTKStructuredData
    point_coords::TArray
    point_data::Dict{String, TPointData}
    cell_data::Dict{String, TCellData}
end
function VTKStructuredData( point_coords,
                            point_data,
                            cell_data,
                            validate::Bool
                          )
    dataset = VTKStructuredData(point_coords, point_data, cell_data)
    if validate
        valid, _error = is_valid(dataset)
        valid ? (return dataset) : throw(_error)
    end
    return dataset
end
function VTKStructuredData(point_coords, validate::Bool)
    point_data = Dict{String, Array}()
    cell_data = Dict{String, Array}()
    return VTKStructuredData(point_coords, point_data, cell_data, validate)
end

mutable struct VTKRectilinearData{TCoord, N, TPointData, TCellData} <: AbstractVTKRectilinearData
    point_coords::NTuple{N, Vector{TCoord}}
    point_data::Dict{String, TPointData}
    cell_data::Dict{String, TCellData}
end
function VTKRectilinearData(point_coords, point_data, cell_data, validate)
    dataset = VTKRectilinearData(point_coords, point_data, cell_data)
    if validate
        valid, _error = is_valid(dataset)
        valid ? (return dataset) : throw(_error)
    end
    return dataset
end
function VTKRectilinearData(point_coords, validate::Bool = false)
    point_data = Dict{String, Array}()
    cell_data = Dict{String, Array}()
    VTKRectilinearData(point_coords, point_data, cell_data, validate)
end

mutable struct VTKUniformRectilinearData{TCoord, N, TPointData, TCellData} <: AbstractVTKRectilinearData
    origin::NTuple{N, TCoord}
    spacing::NTuple{N, TCoord}
    extents::NTuple{N, Int}
    point_data::Dict{String, TPointData}
    cell_data::Dict{String, TCellData}
end
function VTKUniformRectilinearData( origin,
                                    spacing,
                                    extents,
                                    point_data,
                                    cell_data,
                                    validate::Bool
                                  )
    dataset = VTKUniformRectilinearData(origin,
                                        spacing,
                                        extents,
                                        point_data,
                                        cell_data
                                       )
    if validate
        valid, _error = is_valid(dataset)
        valid ? (return dataset) : throw(_error)
    end
    return dataset
end
function VTKUniformRectilinearData( origin,
                                    spacing,
                                    extents,
                                    validate::Bool = false,
                                  )
    point_data = Dict{String, Array}()
    cell_data = Dict{String, Array}()
    return VTKUniformRectilinearData(   origin,
                                        spacing,
                                        extents,
                                        point_data,
                                        cell_data,
                                        validate
                                    )
end
const VTKImageData = VTKUniformRectilinearData

#=
type VTKBergerOligerAMRData <: AbstractVTKAMRData
    data::VTKUniformRectilinearData
    children::Vector{VTKBergerOligerAMRData}
    children_boundary_cells::Matrix{Int}
end
=#

mutable struct VTKMultiblockData{N, TBlocks <: Tuple{Vararg{AbstractStaticVTKData, N}}} <: AbstractVTKMultiblockData
    blocks::TBlocks
end

mutable struct VTKTimeSeriesData{TTime, TData <: AbstractStaticVTKData} <: AbstractTimeSeriesVTKData{TTime, TData}
    timemarkers::Vector{TTime}
    data::Vector{TData}
end
function VTKTimeSeriesData(timemarkers, data, validate)
    dataset = VTKTimeSeriesData(timemarkers, data)
    if validate
        valid, _error = is_valid(dataset)
        valid ? (return dataset) : throw(_error)
    end
    return dataset
end


text:
module Cluster

include("preprocess.jl")

end