text:
struct Boat_ode_params
    bh2o_0::Float64
    angle_ranges::Boat_angle_ranges # in radians
    timing::Boat_timing
    x1va::Float64
    x1air::Float64
end

Boat_ode_params(θ0, bs::Boat_settings) =
    Boat_ode_params(bh2o0(θ0), deg2rad(bs.angle_ranges), bs.timing, bs.x1va, bs.x1air)

function Boat_sim_pars(z, θ, u, w, p::Boat_ode_params, t)
    sθ, cθ = sincos(θ)
    R = SA[cθ -sθ; sθ cθ]
    U = SA[u, w]
    Boat_sim_pars(
        get_bh2o(z, θ, p),
        Boat_angles(p.angle_ranges, p.timing, t),
        η(t, p.timing),
        p.x1va,
        p.x1air,
        (R * U)...,
    )
end

function boatode!(dy, y, p::Boat_ode_params, t)
    # y contains int(u,t), int(w, t), the - position level coordinates
    # as well as u, w, the1  - velocity level coordinates
    _, z, θ, u, w, θ′ = y

    bsp = Boat_sim_pars(z, θ, u, w, p, t)

    M = system_lhs(bsp.boat_angles)
    b = system_rhs(u, w, θ′, bsp.boat_angles)
    F = Q(θ, bsp)

    dy[1] = bsp.x′
    dy[2] = bsp.z′
    dy[3] = θ′
    dy[4:6] = M \ (b + F)
end

function static_equilibrium(θ, settings::Boat_settings)
    (Q_g(
        θ,
        deg2rad(settings.angle_ranges.γ_OA[1]),
        deg2rad(settings.angle_ranges.θ_k[1]),
        deg2rad(settings.angle_ranges.θ_t[1]),
    )+Q_BUOY(θ, bh2o0(θ)))[3]
end

function find_the0(themin, themax, settings::Boat_settings)
    bisectionMethod((the) -> static_equilibrium(the, settings), themin, themax, 1e-2)
end


text:
###
# padded
####

struct PaddedLayout{L} <: MemoryLayout end
for op in (:hcat, :vcat)
    @eval begin
        applylayout(::Type{typeof($op)}, ::A, ::ZerosLayout) where A = PaddedLayout{A}()
        applylayout(::Type{typeof($op)}, ::ScalarLayout, ::ScalarLayout, ::ZerosLayout) = PaddedLayout{ApplyLayout{typeof($op)}}()
        applylayout(::Type{typeof($op)}, ::A, ::PaddedLayout) where A = PaddedLayout{ApplyLayout{typeof($op)}}()
        applylayout(::Type{typeof($op)}, ::ScalarLayout, ::ScalarLayout, ::PaddedLayout) = PaddedLayout{ApplyLayout{typeof($op)}}()
    end
end
applylayout(::Type{typeof(hvcat)}, _, ::A, ::ZerosLayout...) where A = PaddedLayout{A}()
cachedlayout(::A, ::ZerosLayout) where A = PaddedLayout{A}()
sublayout(::PaddedLayout{L}, ::Type{I}) where {L,I<:Tuple{AbstractUnitRange}} =
    PaddedLayout{typeof(sublayout(L(), I))}()
sublayout(::PaddedLayout{L}, ::Type{I}) where {L,I<:Tuple{AbstractUnitRange,AbstractUnitRange}} =
    PaddedLayout{typeof(sublayout(L(), I))}()
sublayout(::PaddedLayout{Lay}, sl::Type{<:Tuple{AbstractUnitRange,Integer}}) where Lay =
    PaddedLayout{typeof(sublayout(Lay(), sl))}()
sublayout(::PaddedLayout{Lay}, sl::Type{<:Tuple{Integer,AbstractUnitRange}}) where Lay =
    PaddedLayout{typeof(sublayout(Lay(), sl))}()
transposelayout(::PaddedLayout{L}) where L = PaddedLayout{typeof(transposelayout(L))}()

paddeddata(A::CachedArray{<:Any,N,<:Any,<:Zeros}) where N = cacheddata(A)
_vcat_paddeddata(A, B::Zeros) = A
_vcat_paddeddata(A, B) = Vcat(A, paddeddata(B))
_vcat_paddeddata(A, B, C...) = Vcat(A, _vcat_paddeddata(B, C...))
paddeddata(A::Vcat) = _vcat_paddeddata(A.args...)

_hcat_paddeddata(A, B::Zeros) = A
_hcat_paddeddata(A, B) = Hcat(A, paddeddata(B))
_hcat_paddeddata(A, B, C...) = Hcat(A, _hcat_paddeddata(B, C...))
paddeddata(A::Hcat) = _hcat_paddeddata(A.args...)

paddeddata(A::Transpose) = transpose(paddeddata(parent(A)))
paddeddata(A::Adjoint) = paddeddata(parent(A))'

_hvcat_paddeddata(N, A, B::Zeros...) = A
paddeddata(A::ApplyMatrix{<:Any,typeof(hvcat)}) = _hvcat_paddeddata(A.args...)

const DualOrPaddedLayout{Lay} = Union{PaddedLayout{Lay},DualLayout{PaddedLayout{Lay}}}

function colsupport(lay::DualOrPaddedLayout{Lay}, A, j) where Lay
    P = paddeddata(A)
    MemoryLayout(P) == lay && return colsupport(UnknownLayout, P, j)
    j̃ = j ∩ axes(P,2)
    cs = colsupport(P,j̃)
    isempty(j̃) ? convert(typeof(cs), Base.OneTo(0)) : cs
end
function rowsupport(lay::DualOrPaddedLayout{Lay}, A, k) where Lay
    P = paddeddata(A)
    MemoryLayout(P) == lay && return rowsupport(UnknownLayout, P, k)
    k̃ = k ∩ axes(P,1)
    rs = rowsupport(P,k̃)
    isempty(k̃) ? convert(typeof(rs), Base.OneTo(0)) : rs
end

function _vcat_resizedata!(::DualOrPaddedLayout{Lay}, B, m...) where Lay
    any(iszero,m) || Base.checkbounds(paddeddata(B), m...)
    B
end

resizedata!(B::Vcat, m...) = _vcat_resizedata!(MemoryLayout(B), B, m...)

function ==(A::CachedVector{<:Any,<:Any,<:Zeros}, B::CachedVector{<:Any,<:Any,<:Zeros})
    length(A) == length(B) || return false
    n = max(A.datasize[1], B.datasize[1])
    resizedata!(A,n); resizedata!(B,n)
    view(A.data,OneTo(n)) == view(B.data,OneTo(n))
end

function _copyto!(::PaddedLayout, ::PaddedLayout, dest::AbstractVector, src::AbstractVector)
    length(src) ≤ length(dest)  || throw(BoundsError())
    src_data = paddeddata(src)
    n = length(src_data)
    resizedata!(dest, n) # if resizeable, otherwise this is a no-op
    dest_data = paddeddata(dest)
    copyto!(view(dest_data,OneTo(n)), src_data)
    zero!(view(dest_data,n+1:length(dest_data)))
    dest
end

function _copyto!(::PaddedLayout, ::ZerosLayout, dest::AbstractVector, src::AbstractVector)
    zero!(paddeddata(dest))
    dest
end

function zero!(::PaddedLayout, A)
    zero!(paddeddata(A))
    A
end

function ArrayLayouts._norm(::PaddedLayout, A, p)
    dat = paddeddata(A)
    if MemoryLayout(dat) isa PaddedLayout
        Base.invoke(norm, Tuple{Any,Real}, dat, p)
    else
        norm(dat, p)
    end
end

# special case handle broadcasting with padded and cached arrays
function _paddedpadded_broadcasted(op, A::AbstractVector, B::AbstractVector)
    a,b = paddeddata(A),paddeddata(B)
    n,m = length(a),length(b)
    dat = if n ≤ m
        [broadcast(op, a, _view_vcat(b,1:n)); broadcast(op, zero(eltype(A)), _view_vcat(b,n+1:m))]
    else
        [broadcast(op, _view_vcat(a,1:m), b); broadcast(op, _view_vcat(a,m+1:n), zero(eltype(B)))]
    end
    CachedArray(dat, broadcast(op, Zeros{eltype(A)}(length(A)), Zeros{eltype(B)}(length(B))))
end

function _paddedpadded_broadcasted(op, A::AbstractMatrix{T}, B::AbstractMatrix{V}) where {T,V}
    a,b = paddeddata(A),paddeddata(B)
    (a_m,a_n) = size(a)
    (b_m,b_n) = size(b)
    m,n = min(a_m,b_m),min(a_n,b_n)
    dat = if a_m ≤ b_m && a_n ≤ b_n
        [broadcast(op, a, _view_vcat(b,1:a_m,1:a_n)) broadcast(op, zero(T), _view_vcat(b,1:a_m,a_n+1:b_n));
         broadcast(op, zero(T), _view_vcat(b,a_m+1:b_m,1:b_n))]
    elseif a_m ≤ b_m
        [broadcast(op, _view_vcat(a,1:a_m,1:b_n), _view_vcat(b,1:a_m,1:b_n)) broadcast(op, _view_vcat(a,1:a_m,b_n+1:a_n), zero(V));
         broadcast(op, zero(T), _view_vcat(b,a_m+1:b_m,1:b_n)) broadcast(op, Zeros{T}(b_m-a_m, a_n-b_n), Zeros{V}(b_m-a_m, a_n-b_n))]
    elseif b_n ≤ a_n # && b_m < a_m
        [broadcast(op, view(a,1:b_m,1:b_n), _view_vcat(b,1:b_m,1:b_n)) broadcast(op, _view_vcat(a,1:b_m,b_n+1:a_n), zero(V));
         broadcast(op, _view_vcat(a,b_m+1:a_m,1:a_n), zero(V))]
    else  # b_m < a_m && a_n < b_n
        [broadcast(op, _view_vcat(a,1:b_m,1:a_n), _view_vcat(b,1:b_m,1:a_n)) broadcast(op, zero(T),  _view_vcat(b,1:b_m,a_n+1:b_n));
         broadcast(op, _view_vcat(a,b_m+1:a_m,1:a_n), zero(V)) broadcast(op, Zeros{T}(a_m-b_m, b_n-a_n), Zeros{V}(a_m-b_m, b_n-a_n))]
    end
    PaddedArray(dat, max.(size(A), size(B))...)
end

layout_broadcasted(::PaddedLayout, ::PaddedLayout, op, A::AbstractVector, B::AbstractVector) =
    _paddedpadded_broadcasted(op, A, B)
layout_broadcasted(::PaddedLayout, ::PaddedLayout, op, A::AbstractMatrix, B::AbstractMatrix) =
    _paddedpadded_broadcasted(op, A, B)
layout_broadcasted(::PaddedLayout, ::PaddedLayout, ::typeof(*), A::Vcat{<:Any,1}, B::AbstractVector) =
    _paddedpadded_broadcasted(*, A, B)
layout_broadcasted(::PaddedLayout, ::PaddedLayout, ::typeof(*), A::AbstractVector, B::Vcat{<:Any,1}) =
    _paddedpadded_broadcasted(*, A, B)

function layout_broadcasted(_, ::PaddedLayout, op, A::AbstractVector, B::AbstractVector)
    b = paddeddata(B)
    m = length(b)
    zB = Zeros{eltype(B)}(size(B)...)
    CachedArray(convert(Array,broadcast(op, view(A,1:m), b)), broadcast(op, A, zB))
end

function layout_broadcasted(::PaddedLayout, _, op, A::AbstractVector, B::AbstractVector)
    a = paddeddata(A)
    n = length(a)
    zA = Zeros{eltype(A)}(size(A)...)
    CachedArray(convert(Array,broadcast(op, a, view(B,1:n))), broadcast(op, zA, B))
end

function layout_broadcasted(::PaddedLayout, ::CachedLayout, op, A::AbstractVector, B::AbstractVector)
    a = paddeddata(A)
    n = length(a)
    resizedata!(B,n)
    Bdata = cacheddata(B)
    b = view(Bdata,1:n)
    zA1 = Zeros{eltype(A)}(size(Bdata,1)-n)
    zA = Zeros{eltype(A)}(size(A)...)
    CachedArray([broadcast(op, a, b); broadcast(op, zA1, @view(Bdata[n+1:end]))], broadcast(op, zA, B.array))
end

function layout_broadcasted(::CachedLayout, ::PaddedLayout, op, A::AbstractVector, B::AbstractVector)
    b = paddeddata(B)
    n = length(b)
    resizedata!(A,n)
    Adata = cacheddata(A)
    a = view(Adata,1:n)
    zB1 = Zeros{eltype(B)}(size(Adata,1)-n)
    zB = Zeros{eltype(B)}(size(B)...)
    CachedArray([broadcast(op, a, b); broadcast(op, @view(Adata[n+1:end]), zB1)], broadcast(op, A.array, zB))
end

layout_broadcasted(lay::ApplyLayout{typeof(vcat)}, ::PaddedLayout, op, A::AbstractVector, B::AbstractVector) =
    layout_broadcasted(lay, lay, op, A, B)
layout_broadcasted(::PaddedLayout, lay::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector) =
    layout_broadcasted(lay, lay, op, A, B)

# special case for * to preserve Vcat Structure

layout_broadcasted(::ApplyLayout{typeof(vcat)}, lay::PaddedLayout, ::typeof(*), A::AbstractVector, B::AbstractVector) =
    layout_broadcasted(UnknownLayout(), lay, *, A, B)
layout_broadcasted(lay::PaddedLayout, ::ApplyLayout{typeof(vcat)}, ::typeof(*), A::AbstractVector, B::AbstractVector) =
    layout_broadcasted(lay, UnknownLayout(), *, A, B)

for op in (:+, :-)
    @eval function layout_broadcasted(::PaddedLayout, ::PaddedLayout, ::typeof($op), A::Vcat{T,1}, B::Vcat{V,1}) where {T,V}
        a,b = paddeddata(A),paddeddata(B)
        if a isa Number && b isa Number
            Vcat($op(a, b), Zeros{promote_type(T,V)}(max(length(A),length(B))-1))
        elseif a isa Number
            m = length(b)
            dat = [broadcast($op, a, b[1]); broadcast($op,view(b,2:m))]
            Vcat(convert(Array,dat), Zeros{promote_type(T,V)}(max(length(A),length(B))-length(dat)))
        elseif b isa Number
            n = length(a)
            dat = [broadcast($op, a[1], b); view(a,2:n)]
            Vcat(convert(Array,dat), Zeros{promote_type(T,V)}(max(length(A),length(B))-length(dat)))
        else
            n,m = length(a),length(b)
            dat = if n == m
                broadcast($op, a, b)
            elseif n > m
                [broadcast($op, view(a,1:m), b); view(a,m+1:n)]
            else # n < m
                [broadcast($op, a, view(b,1:n)); broadcast($op,view(b,n+1:m))]
            end
            Vcat(convert(Array,dat), Zeros{promote_type(T,V)}(max(length(A),length(B))-length(dat)))
        end
    end

    @eval function layout_broadcasted(::PaddedLayout, ::PaddedLayout, ::typeof($op), A::Vcat{T,2}, B::Vcat{V,2}) where {T,V}
        a,b = paddeddata(A),paddeddata(B)
        n,m = size(a,1),size(b,1)
        dat = if n == m
            broadcast($op, a, b)
        elseif n > m
            [broadcast($op, view(a,1:m,:), b); view(a,m+1:n,:)]
        else # n < m
            [broadcast($op, a, view(b,1:n,:)); broadcast($op,view(b,n+1:m,:))]
        end
        Vcat(convert(Array,dat), Zeros{promote_type(T,V)}(max(size(A,1),size(B,1))-size(dat,1),size(a,2)))
    end
end

function layout_broadcasted(::PaddedLayout, ::PaddedLayout, ::typeof(*), A::Vcat{T,1}, B::Vcat{V,1}) where {T,V}
    a,b = paddeddata(A),paddeddata(B)
    if a isa Number || b isa Number
        Vcat(a[1]*b[1], Zeros{promote_type(T,V)}(max(length(A),length(B))-1))
    else
        n = min(length(a),length(b))
        dat = broadcast(*, view(a,1:n), view(b,1:n))
        Vcat(convert(Array,dat), Zeros{promote_type(T,V)}(max(length(A),length(B))-n))
    end
end

function layout_broadcasted(_, ::PaddedLayout, ::typeof(*), A::AbstractVector{T}, B::Vcat{V,1}) where {T,V}
    b = paddeddata(B)
    m = length(b)
    dat = broadcast(*, view(A,1:m), b)
    Vcat(convert(Array,dat), Zeros{promote_type(T,V)}(max(length(A),length(B))-m))
end

function layout_broadcasted(::PaddedLayout, _, ::typeof(*), A::Vcat{T,1}, B::AbstractVector{V}) where {T,V}
    a = paddeddata(A)
    n = length(a)
    dat = broadcast(*, a, view(B,1:n))
    Vcat(convert(Array,dat), Zeros{promote_type(T,V)}(max(length(A),length(B))-n))
end

layout_broadcasted(::PaddedLayout, lay::ApplyLayout{typeof(vcat)}, ::typeof(*), A::Vcat{<:Any,1}, B::AbstractVector) = layout_broadcasted(lay, lay, *, A, B)
layout_broadcasted(lay::ApplyLayout{typeof(vcat)}, ::PaddedLayout, ::typeof(*), A::AbstractVector, B::Vcat{<:Any,1}) = layout_broadcasted(lay, lay, *, A, B)

layout_broadcasted(_, _, op, A, B) = Base.Broadcast.Broadcasted{typeof(Base.BroadcastStyle(Base.BroadcastStyle(typeof(A)),Base.BroadcastStyle(typeof(B))))}(op, (A, B))

###
# Dot/Axpy
###

struct Axpy{StyleX,StyleY,T,XTyp,YTyp}
    α::T
    X::XTyp
    Y::YTyp
end

Axpy(α::T, X::XTyp, Y::YTyp) where {T,XTyp,YTyp} = Axpy{typeof(MemoryLayout(XTyp)), typeof(MemoryLayout(YTyp)), T, XTyp, YTyp}(α, X, Y)
materialize!(d::Axpy{<:Any,<:Any,<:Number,<:AbstractArray,<:AbstractArray}) = Base.invoke(BLAS.axpy!, Tuple{Number,AbstractArray,AbstractArray}, d.α, d.X, d.Y)
function materialize!(d::Axpy{<:PaddedLayout,<:PaddedLayout,U,<:AbstractVector{T},<:AbstractVector{V}}) where {U,T,V}
    x = paddeddata(d.X)
    resizedata!(d.Y, length(x))
    y = paddeddata(d.Y)
    BLAS.axpy!(d.α, x, view(y,1:length(x)))
    y
end
axpy!(α, X, Y) = materialize!(Axpy(α,X,Y))
BLAS.axpy!(α, X::LazyArray, Y::AbstractArray) = axpy!(α, X, Y)
BLAS.axpy!(α, X::SubArray{<:Any,N,<:LazyArray}, Y::AbstractArray) where N = axpy!(α, X, Y)

###
# l/rmul!
###

function materialize!(M::Lmul{ScalarLayout,<:PaddedLayout})
    lmul!(M.A, paddeddata(M.B))
    M.B
end

function materialize!(M::Rmul{<:PaddedLayout,ScalarLayout})
    rmul!(paddeddata(M.A), M.B)
    M.A
end

_norm2(::PaddedLayout, a) = norm(paddeddata(a),2)
_norm1(::PaddedLayout, a) = norm(paddeddata(a),1)
_normInf(::PaddedLayout, a) = norm(paddeddata(a),Inf)
_normp(::PaddedLayout, a, p) = norm(paddeddata(a),p)

function copy(D::Dot{layA, layB}) where {layA<:PaddedLayout,layB<:PaddedLayout}
    a = paddeddata(D.A)
    b = paddeddata(D.B)
    T = eltype(D)
    if MemoryLayout(a) isa layA && MemoryLayout(b) isa layB
        return convert(T, dot(Array(a),Array(b)))
    end
    length(a) == length(b) && return convert(T, dot(a,b))
    # following handles scalars
    ((length(a) == 1) || (length(b) == 1)) && return convert(T, a[1] * b[1])
    m = min(length(a), length(b))
    convert(T, dot(view(a, 1:m), view(b, 1:m)))
end

function copy(D::Dot{<:PaddedLayout})
    a = paddeddata(D.A)
    m = length(a)
    v = view(D.B, 1:m)
    if MemoryLayout(a) isa PaddedLayout
        convert(eltype(D), dot(Array(a), v))
    else
        convert(eltype(D), dot(a, v))
    end

end

function copy(D::Dot{<:Any, <:PaddedLayout})
    b = paddeddata(D.B)
    m = length(b)
    v = view(D.A, 1:m)
    if MemoryLayout(b) isa PaddedLayout
        convert(eltype(D), dot(v, Array(b)))
    else
        convert(eltype(D), dot(v, b))
    end
end

_vcat_sub_arguments(::PaddedLayout, A, V) = _vcat_sub_arguments(ApplyLayout{typeof(vcat)}(), A, V)

_lazy_getindex(dat, kr...) = view(dat, kr...)
_lazy_getindex(dat::Number, _...) = dat

function sub_paddeddata(_, S::SubArray{<:Any,1,<:AbstractVector})
    dat = paddeddata(parent(S))
    (kr,) = parentindices(S)
    kr2 = kr ∩ axes(dat,1)
    _lazy_getindex(dat, kr2)
end

function sub_paddeddata(_, S::SubArray{<:Any,1,<:Any,<:Tuple{Integer,Any}})
    P = parent(S)
    (k,jr) = parentindices(S)
    # need to resize in case dat is empty... not clear how to take a vector view of a 0-dimensional matrix in a type-stable way
    resizedata!(P, k, 1); dat = paddeddata(P)
    if k in axes(dat,1)
        _lazy_getindex(dat, k, jr ∩ axes(dat,2))
    else
        _lazy_getindex(dat, first(axes(dat,1)), jr ∩ Base.OneTo(0))
    end
end
function sub_paddeddata(_, S::SubArray{<:Any,1,<:Any,<:Tuple{Any,Integer}})
    P = parent(S)
    (kr,j) = parentindices(S)
    # need to resize in case dat is empty... not clear how to take a vector view of a 0-dimensional matrix in a type-stable way
    resizedata!(P, 1, j); dat = paddeddata(P)
    if j in axes(dat,2)
        _lazy_getindex(dat, kr ∩ axes(dat,1), j)
    else
        _lazy_getindex(dat, kr ∩ Base.OneTo(0), first(axes(dat,2)))
    end
end

function sub_paddeddata(_, S::SubArray{<:Any,2})
    P = parent(S)
    (kr,jr) = parentindices(S)
    dat = paddeddata(P)
    kr2 = kr ∩ axes(dat,1)
    jr2 = jr ∩ axes(dat,2)
    _lazy_getindex(dat, kr2, jr2)
end

paddeddata(S::SubArray) = sub_paddeddata(MemoryLayout(parent(S)), S)

function _padded_sub_materialize(v::AbstractVector{T}) where T
    dat = paddeddata(v)
    if MemoryLayout(dat) isa PaddedLayout
        Vcat(dat, Zeros{T}(length(v) - length(dat)))
    else
        Vcat(sub_materialize(dat), Zeros{T}(length(v) - length(dat)))
    end
end

sub_materialize(::PaddedLayout, v::AbstractVector, _) = _padded_sub_materialize(v)

function sub_materialize(l::PaddedLayout, v::AbstractMatrix, _)
    dat = paddeddata(v)
    PaddedArray(MemoryLayout(dat) isa PaddedLayout ? dat : sub_materialize(dat), size(v)...)
end

function layout_replace_in_print_matrix(::PaddedLayout{Lay}, f::AbstractVecOrMat, k, j, s) where Lay
    # avoid infinite-loop
    f isa SubArray && return Base.replace_in_print_matrix(parent(f), Base.reindex(f.indices, (k,j))..., s)
    data = paddeddata(f)
    MemoryLayout(data) isa PaddedLayout{Lay} && return layout_replace_in_print_matrix(UnknownLayout(), f, k, j, s)
    k in axes(data,1) && j in axes(data,2) && return _replace_in_print_matrix(data, k, j, s)
    Base.replace_with_centered_mark(s)
end

# avoid ambiguity in LazyBandedMatrices
copy(M::Mul{<:DiagonalLayout,<:PaddedLayout}) = copy(Lmul(M))

# Triangular columns

sublayout(::TriangularLayout{'U','N', ML}, INDS::Type{<:Tuple{KR,Integer}}) where {KR,ML} =
    sublayout(PaddedLayout{typeof(sublayout(ML(), INDS))}(), Tuple{KR})

sublayout(::TriangularLayout{'L','N', ML}, INDS::Type{<:Tuple{Integer,JR}}) where {JR,ML} =
    sublayout(PaddedLayout{typeof(sublayout(ML(), INDS))}(), Tuple{JR})

function sub_paddeddata(::TriangularLayout{'U','N'}, S::SubArray{<:Any,1,<:AbstractMatrix,<:Tuple{Any,Integer}})
    P = parent(S)
    (kr,j) = parentindices(S)
    view(triangulardata(P), kr ∩ (1:j), j)
end

function sub_paddeddata(::TriangularLayout{'L','N'}, S::SubArray{<:Any,1,<:AbstractMatrix,<:Tuple{Integer,Any}})
    P = parent(S)
    (k,jr) = parentindices(S)
    view(triangulardata(P), k, jr ∩ (1:k))
end

###
# setindex
###

@inline ndims(A::Applied{<:Any,typeof(setindex)}) = ndims(A.args[1])
@inline eltype(A::Applied{<:Any,typeof(setindex)}) = eltype(A.args[1])
axes(A::ApplyArray{<:Any,N,typeof(setindex)}) where N = axes(A.args[1])

function getindex(A::ApplyVector{T,typeof(setindex)}, k::Integer) where T
    P,v,kr = A.args
    convert(T, k in kr ? v[something(findlast(isequal(k),kr))] : P[k])::T
end

function getindex(A::ApplyMatrix{T,typeof(setindex)}, k::Integer, j::Integer) where T
    P,v,kr,jr = A.args
    convert(T, k in kr && j in jr ? v[something(findlast(isequal(k),kr)),something(findlast(isequal(j),jr))] : P[k,j])::T
end

const PaddedArray{T,N,M} = ApplyArray{T,N,typeof(setindex),<:Tuple{Zeros,M,Vararg{Any,N}}}
const PaddedVector{T,M} = PaddedArray{T,1,M}
const PaddedMatrix{T,M} = PaddedArray{T,2,M}

MemoryLayout(::Type{<:PaddedArray{T,N,M}}) where {T,N,M} = PaddedLayout{typeof(MemoryLayout(M))}()
paddeddata(A::PaddedArray) = A.args[2]

PaddedArray(A::AbstractArray{T,N}, n::Vararg{Integer,N}) where {T,N} = PaddedArray(A, map(oneto,n))
PaddedArray(A::AbstractArray{T,N}, ax::NTuple{N,Any}) where {T,N} = ApplyArray{T,N}(setindex, Zeros{T,N}(ax), A, axes(A)...)
PaddedArray(A::T, n::Vararg{Integer,N}) where {T<:Number,N} = ApplyArray{T,N}(setindex, Zeros{T,N}(n...), A, ntuple(_ -> OneTo(1),N)...)
(PaddedArray{T,N} where T)(A, n::Vararg{Integer,N}) where N = PaddedArray(A, n...)

BroadcastStyle(::Type{<:PaddedArray{<:Any,N}}) where N = LazyArrayStyle{N}()

ext: jl


text:
# shock capturning functions

"""
  Main entry point for shock capturing.  Applies the specified shock capturing
  scheme to every element

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * data: an [`AbstractShockCaputring`](@ref) object
"""
function applyShockCapturing(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts,
                             capture::AbstractVolumeShockCapturing)

  # unlike AbstractFaceShockCapturing, nothing to do here

  sensor = getShockSensor(capture)
  calcShockCapturing(mesh, sbp, eqn, opts, sensor, capture)

  return nothing
end

function applyShockCapturing(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts,
                             capture::AbstractFaceShockCapturing)

  sensor = getShockSensor(capture)
  shockmesh = eqn.params.shockmesh
  setupShockCapturing(mesh, sbp, eqn, opts, sensor, capture, shockmesh)

  # call shock capturing scheme
  calcShockCapturing(mesh, sbp, eqn, opts, capture, shockmesh)

  return nothing
end

"""
  This function sets up the `ShockedElements` object and calls `allocateArrays`
  on the `AbstractFaceShockCapturing` object afterwards.  This completes
  the preparations required before the shock capturing term can be calculated

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * sensor
   * capture
"""
@noinline function setupShockCapturing(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts,
                             sensor::AbstractShockSensor,
                             capture::AbstractFaceShockCapturing,
                             shockmesh::ShockedElements)

  if mesh.commsize > 1 &&
    getParallelData(eqn.shared_data) != PARALLEL_DATA_ELEMENT

    error("shock capturing requires PARALLEL_DATA_ELEMENT")
  end

  assertReceivesWaited(eqn.shared_data)  # parallel communication for the regular
                                     # face integrals should already have
                                     # finished parallel communication

  # re-using the shockmesh from one iteration to the next makes allows
  # the algorithm to re-use existing arrays that depend on the number of
  # elements with the shock in it.
  reset(shockmesh)
  numEl_shock = 0
  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i)
    coords_i = ro_sview(mesh.coords, :, :, i)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i)
    jac_i = ro_sview(mesh.jac, :, i)

    val = isShockElement(eqn.params, sbp, sensor, q_i, i, coords_i, dxidx_i,
                            jac_i)
    if val
      # push to shockmesh
      #println("element ", i, " viscoscity = ", ee)
      numEl_shock += 1
      push!(shockmesh, i)
    end
  end

  #println("counted ", numEl_shock, " elements with nonzero viscoscity")
  completeShockElements(mesh, shockmesh)
  @assert shockmesh.numShock == numEl_shock

  allocateArrays(capture, mesh, shockmesh)

  return nothing
end

#------------------------------------------------------------------------------
# Debugging/testing

using PumiInterface
using apf

function writeShockSensorField(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts,
                               sensor::AbstractShockSensor) where {Tsol, Tres}

  fname = "shock_sensor"

  f_ptr = apf.findField(mesh.m_ptr, fname)
  if f_ptr == C_NULL
    fshape = apf.getConstantShapePtr(mesh.dim)
    f_ptr = apf.createPackedField(mesh, fname, 2, fshape)
  end

  vals = zeros(Float64, 2)
  numEl_shock = 0

  Se = zeros(Tres, mesh.dim, mesh.numNodesPerElement)
  ee = zeros(Tres, mesh.dim, mesh.numNodesPerElement)
  for i=1:mesh.numEl
    q_i = sview(eqn.q, :, :, i)
    coords_i = ro_sview(mesh.coords, :, :, i)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i)
    jac_i = sview(mesh.jac, :, i)

    is_nonzero = getShockSensor(eqn.params, sbp, sensor, q_i, i, coords_i, dxidx_i,
                            jac_i, Se, ee)

    #val1 = computeL2Norm(eqn.params, sbp, jac_i, Se)
    #val2 = computeL2Norm(eqn.params, sbp, jac_i, ee)
    val1 = zero(Tres); val2 = zero(Tres)
    for j=1:mesh.numNodesPerElement
      for d=1:mesh.dim
        val1 += Se[d, j]
        val2 += ee[d, j]
      end
    end

    val1 /= mesh.dim*mesh.numNodesPerElement
    val2 /= mesh.dim*mesh.numNodesPerElement

    vals[1] = val1
    vals[2] = val2
    apf.setComponents(f_ptr, mesh.elements[i], 0, vals)
    if is_nonzero
      numEl_shock += 1
    end
  end

  numEl_shock = MPI.Allreduce(numEl_shock, MPI.SUM, eqn.comm)

  if mesh.myrank == 0
    println("wrote shock field with ", numEl_shock, " elements with non-zero viscoscity")
  end

  return nothing
end


text:
function add_chan!(S::GphysData, C::GphysChannel, strict::Bool)
  if isempty(S)
    push!(S, C)
    return 1
  end
  j = strict ? channel_match(S, C) : findid(S, C.id)
  if j > 0
    T1 = t_extend(S.t[j], C.t[1,2], length(C.x), S.fs[j])
    if T1 != nothing
      S.t[j] = T1
    end
    append!(S.x[j], C.x)
  else
    push!(S, C)
    j = S.n
  end
  return j
end


text:
module ImageModifierProject

using Sockets
using HTTP
using LibPQ

const editorHtml = """
<!DOCTYPE html>

  <meta charset=utf-8>
  <meta name="viewport" content="width=device-width, initial-scale=1">

<canvas id=main></canvas>

"""

mutable struct Image
  width::Int
  height::Int
  data::Array{UInt8,2}
end
Image(w,h) = Image(w, h, zeros(UInt8, w, h))

function imgToStr(img::Image, x::Int=0, y::Int=0, w::Int=img.width, h::Int=img.height; onlyData = false)
  data = join(vec(string.(img.data[1+x : x+w, 1+y : y+h], base=16))) # Column-major.
  if onlyData
    data
  else
    string(x, " ", y, " ", w, " ", h, " ", data)
  end
end
function strToImg(w::Int, h::Int, data::String)
  if w*h != length(data)
    Image(w, h)
  else
  end
  data = reshape(parse.(UInt8, collect(data), base=16), (w, h))
  Image(w, h, data)
end

conn = nothing
scratchpad = nothing
had_updates = false
function createDB(conn)
  execute(conn, """
  CREATE TABLE IF NOT EXISTS image (
    id integer PRIMARY KEY,
    width integer NOT NULL,
    height integer NOT NULL,
    data text NOT NULL
  )
  """)
end
function saveToDB(conn, img)
  data = imgToStr(img; onlyData = true)
  result = execute(conn, """
  INSERT INTO image (id, width, height, data)
  VALUES (0, \$1, \$2, \$3)
  ON CONFLICT (id) DO UPDATE
    SET width = excluded.width, height = excluded.height, data = excluded.data;
  """, [img.width, img.height, data])
  # Sure hope there are no exceptions, in this no-try-catch code.
end
function loadFromDB(conn, width, height)
  try
    result = execute(conn, """
    SELECT width, height, data FROM image WHERE id=0
    """)
    w, h, data = result[1,1], result[1,2], result[1,3]
    strToImg(width, height, data)
  catch err
    Image(width, height)
  end
end
function periodicallySaveImageToDB(t)
  global had_updates
  if !had_updates return end
  had_updates = false
  saveToDB(conn, scratchpad)
end

function __init__() # Ah yes, Julia, execute non-__init__ code during pre-compilation. Makes sense.
  # Every hour, save the image to the DB if it changed.
  global conn, scratchpad, had_updates
  conn = LibPQ.Connection(ENV["DATABASE_URL"])
  createDB(conn)
  scratchpad = loadFromDB(conn, 500, 500)
  Timer(periodicallySaveImageToDB, 600; interval = 600)

  # Yes, ARGS is totally a compile-time constant, and pre-compilation should really complain about ARGS[1].
  ar = ARGS
  push!(ar, "8081")
  if size(ar)[1] < 1
    throw("Too few args; expected the port")
  end
  port = parse(Int64, ar[1])
  server = Sockets.listen(Sockets.InetAddr("0.0.0.0", port))

  connections = []
  println("Serving on port ", port, "…")
  HTTP.serve(server=server; stream=true) do stream::HTTP.Stream
    if stream.message.target == "/"
      write(stream, editorHtml)
      stream.message.response.status = 200
    elseif stream.message.target == "/image"
      # Set up a WebSocket: send image & updates, receive updates.
      stream.message.response.status = 200
      HTTP.WebSockets.upgrade(stream) do ws
        # Send the whole image, then receive one-pixel updates unless they are too fast.
        global had_updates
        img = scratchpad
        push!(connections, ws)
        write(ws, imgToStr(img) * "\n\n")
        last_update = time()
        try
          while !eof(ws)
            data = String(readavailable(ws)) # Still don't know if it's one-read-per-message, or shitty. No info on the Web; not about to do stress-tests to know.
            had_updates = true
            values = map(x -> tryparse(Int64, x), split(strip(data), " "))
            if size(values)[1] != 3
              continue
            end
            x, y, value = values
            if x < 0 || y < 0 || value < 0 || x >= img.width || y >= img.height || value >= 16
              continue
            end
            new_update = time()
            if new_update - last_update > .5
              # Update and broadcast, at most once per second.
              img.data[x+1, y+1] = value
              update_msg = string(imgToStr(img, x, y, 1, 1))
              for conn in connections
                write(conn, update_msg)
              end
            else
              # Fail.
              write(ws, string(imgToStr(img, x, y, 1, 1)))
            end
            last_update = new_update
          end
        catch
          # Socket went away (disconnected), most likely. Ignore it.
        end
        filter!(x -> x !== ws, connections)
      end
    else
      write(stream, "404 not found (this is a good error page, what are you talking about)")
      stream.message.response.status = 404
    end
  end
end

end # module

max_stars_repo_name: Antipurity/image-modifier-project

max_issues_repo_name: Antipurity/image-modifier-project

max_forks_repo_name: Antipurity/image-modifier-project


text:
@testset "Simulator: codegen" begin

    robots = [hopper, biped, quadruped]

    for robot in robots
        # kinematics
        contact_kinematics = RoboDojo.eval(Symbol(String(RoboDojo.name(robot)) * "_contact_kinematics"))
        contact_kinematics_jacobians = RoboDojo.eval(Symbol(String(RoboDojo.name(robot)) * "_contact_kinematics_jacobians"))

        # codegen
        mass_matrix, dynamics_bias = RoboDojo.codegen_dynamics(robot)
        r_model, rz_model, rθ_model = RoboDojo.codegen_residual(robot, mass_matrix, dynamics_bias, contact_kinematics, contact_kinematics_jacobians)

        # evaluate methods
        nz = num_var(robot)
        nθ = num_data(robot, nf=length(friction_coefficients(robot)))
        r0 = zeros(nz)
        z0 = rand(nz)
        θ0 = rand(nθ)
        μ0 = ones(1)
        rz0 = rand(nz, nz)
        rθ0 = rand(nz, nθ)

        eval(r_model)(r0, z0, θ0, μ0)
        eval(rz_model)(rz0, z0, θ0)
        eval(rθ_model)(rθ0, z0, θ0)

        # test methods
        @test sum(r0) != 0.0
        @test sum(rz0) != 0.0
        @test sum(rθ0) != 0.0
    end
end


text:
using Pkg
pkg"precompile"

using Test, MPI

nprocs = clamp(Sys.CPU_THREADS, 2, 4)

@info "Testing: test_cman_julia.jl"
run(`$(Base.julia_cmd()) $(joinpath(@__DIR__, "test_cman_julia.jl")) $nprocs`)

if VERSION >= v"1.5"
    @info "Testing: test_cman_mpi.jl"
    mpiexec() do cmd
        run(`$cmd -n $nprocs $(Base.julia_cmd()) $(joinpath(@__DIR__, "test_cman_mpi.jl"))`)
    end
else
    @warn "MPI_TRANSPORT_ALL broken on Julia versions < 1.5\nSee https://github.com/JuliaParallel/MPIClusterManagers.jl/issues/9"
end

@info "Testing: test_cman_tcp.jl"
mpiexec() do cmd
    run(`$cmd -n $nprocs $(Base.julia_cmd()) $(joinpath(@__DIR__, "test_cman_tcp.jl"))`)
end


text:
import LightQuery
import Documenter: doctest

doctest(LightQuery)


text:
"""
Constructs a DGMesh composite type, a subtype of the abstract type Mesh.

    DGMesh(
        model::SFFM.Model,
        Nodes::Array{Float64,1},
        NBases::Int;
        Fil::IndexDict=IndexDict(),
        Basis::String = "legendre",
    )

# Arguments
- `model`: a Model object
- `Nodes::Array{Float64,1}`: (K+1)×1 array, specifying the edges of the cells
- `NBases::Int`: specifying the number of bases within each cell (same for all
    cells)
- `Fil::IndexDict`: (optional) A dictionary of the sets Fᵢᵐ, they
    keys are Strings specifying i and m, i.e. `"2+"`, the values are BitArrays of
    boolean values which specify which cells of the stencil correspond to Fᵢᵐ. If no
    value specified then `Fil` is generated automatically evaluating ``r_i(x)`` at
    the modpoint of each cell.
- `Basis::String`: a string specifying whether to use the `"lagrange"` basis or
    the `"legendre"` basis

# Output
- a Mesh object with fieldnames:
    - `:NBases`: the number of bases in each cell
    - `:CellNodes`: Array of nodal points (cell edges + GLL points)
    - `:Fil`: As described in the arguments
    - `:Δ`:A vector of mesh widths, Δ[k] = x_{k+1} - x_k
    - `:NIntervals`: The number of cells
    - `:Nodes`: the cell edges
    - `:TotalNBases`: `NIntervals*NBases`
    - `:Basis`: a string specifying whether the

# Examples
TBC

#
A blank initialiser for a DGMesh.

    DGMesh()

Used for initialising a blank plot only. There is no reason to call this, ever.
"""
struct DGMesh <: Mesh
    Nodes::Array{<:Real,1}
    NBases::Int
    Fil::IndexDict
    Basis::String
end
# Convenience constructors
function DGMesh(
    model::SFFM.Model,
    Nodes::Array{<:Real,1},
    NBases::Int;
    Fil::SFFM.IndexDict=SFFM.IndexDict(),
    Basis::String = "lagrange",
    v::Bool = false,
)
    if isempty(Fil)
        Fil = MakeFil(model, Nodes)
    end

    mesh = DGMesh(Nodes, NBases, Fil, Basis)

    v && println("UPDATE: DGMesh object created with fields ", fieldnames(SFFM.DGMesh))
    return mesh
end
function DGMesh()
    DGMesh(
        [0.0],
        0,
        IndexDict(),
        "",
    )
end

"""

    NBases(mesh::DGMesh)

Number of bases in a cell
"""
NBases(mesh::DGMesh) = mesh.NBases

"""

    NIntervals(mesh::Mesh)

Total number of cells for a mesh
"""
NIntervals(mesh::Mesh) = length(mesh.Nodes) - 1

"""

    Δ(mesh::Mesh)

The width of each cell
"""
Δ(mesh::Mesh) = mesh.Nodes[2:end] - mesh.Nodes[1:end-1]

"""

    TotalNBases(mesh::DGMesh)

Total number of bases in the stencil
"""
TotalNBases(mesh::Mesh) = NBases(mesh) * NIntervals(mesh)

"""

    CellNodes(mesh::DGMesh)

The positions of the GLJ nodes within each cell
"""
function CellNodes(mesh::DGMesh)
    meshNBases = NBases(mesh)
    meshNIntervals = NIntervals(mesh)
    cellNodes = zeros(Float64, NBases(mesh), NIntervals(mesh))
    if meshNBases > 1
        z = Jacobi.zglj(meshNBases, 0, 0) # the LGL nodes
    elseif meshNBases == 1
        z = 0.0
    end
    for i = 1:meshNIntervals
        # Map the LGL nodes on [-1,1] to each cell
        cellNodes[:, i] .= (mesh.Nodes[i+1] + mesh.Nodes[i]) / 2 .+ (mesh.Nodes[i+1] - mesh.Nodes[i]) / 2 * z
    end
    cellNodes[1,:] .+= sqrt(eps())
    if meshNBases>1
        cellNodes[end,:] .-= sqrt(eps())
    end
    return cellNodes
end

"""

    Basis(mesh::DGMesh)

Returns mesh.Basis; either "lagrange" or "legendre"
"""
Basis(mesh::DGMesh) = mesh.Basis

"""
Construct a generalised vandermonde matrix.

    vandermonde( nBases::Int)

Note: requires Jacobi package Pkg.add("Jacobi")

# Arguments
- `nBases::Int`: the degree of the basis

# Output
- a tuple with keys
    - `:V::Array{Float64,2}`: where `:V[:,i]` contains the values of the `i`th
        legendre polynomial evaluate at the GLL nodes.
    - `:inv`: the inverse of :V
    - `:D::Array{Float64,2}`: where `V.D[:,i]` contains the values of the derivative
        of the `i`th legendre polynomial evaluate at the GLL nodes.
"""
function vandermonde(nBases::Int)
    if nBases > 1
        z = Jacobi.zglj(nBases, 0, 0) # the LGL nodes
    elseif nBases == 1
        z = 0.0
    end
    V = zeros(Float64, nBases, nBases)
    DV = zeros(Float64, nBases, nBases)
    if nBases > 1
        for j = 1:nBases
            # compute the polynomials at gauss-labotto quadrature points
            V[:, j] = Jacobi.legendre.(z, j - 1) .* sqrt((2 * (j - 1) + 1) / 2)
            DV[:, j] = Jacobi.dlegendre.(z, j - 1) .* sqrt((2 * (j - 1) + 1) / 2)
        end
        # Compute the Gauss-Lobatto weights for numerical quadrature
        w =
            2.0 ./ (
                nBases *
                (nBases - 1) *
                Jacobi.legendre.(Jacobi.zglj(nBases, 0, 0), nBases - 1) .^ 2
            )
    elseif nBases == 1
        V .= [1/sqrt(2)]
        DV .= [0]
        w = [2]
    end
    return (V = V, inv = inv(V), D = DV, w = w)
end

function local_dg_operators(
    mesh::DGMesh;
    v::Bool = false,
)
    ## Construct local blocks
    V = vandermonde(NBases(mesh))
    if Basis(mesh) == "legendre"
        Dw = (
            DwInv = LinearAlgebra.diagm(0 => ones(Float64, NBases(mesh))),
            Dw = LinearAlgebra.diagm(0 => ones(Float64, NBases(mesh))),
        ) # function weights are not available for legendre basis as this is
        # in density land
        MLocal = Matrix{Float64}(LinearAlgebra.I(NBases(mesh)))
        GLocal = V.inv * V.D
        MInvLocal = Matrix{Float64}(LinearAlgebra.I(NBases(mesh)))
        Phi = V.V[[1; end], :]
    elseif Basis(mesh) == "lagrange"
        Dw = (
            DwInv = LinearAlgebra.diagm(0 => 1.0 ./ V.w),
            Dw = LinearAlgebra.diagm(0 => V.w),
        )# function weights so that we can work in probability land as
        # opposed to density land

        MLocal = Dw.DwInv * V.inv' * V.inv * Dw.Dw
        GLocal = Dw.DwInv * V.inv' * V.inv * (V.D * V.inv) * Dw.Dw
        MInvLocal = Dw.DwInv * V.V * V.V' * Dw.Dw
        Phi = (V.inv*V.V)[[1; end], :]
    end

    PosDiagBlock = -Dw.DwInv * Phi[end, :] * Phi[end, :]' * Dw.Dw
    NegDiagBlock = Dw.DwInv * Phi[1, :] * Phi[1, :]' * Dw.Dw
    UpDiagBlock = Dw.DwInv * Phi[end, :] * Phi[1, :]' * Dw.Dw
    LowDiagBlock = Dw.DwInv * Phi[1, :] * Phi[end, :]' * Dw.Dw

    out = (
        G = GLocal, M = MLocal, MInv = MInvLocal, V = V, Phi = Phi, Dw = Dw,
        PosDiagBlock = PosDiagBlock, NegDiagBlock = NegDiagBlock,
        UpDiagBlock = UpDiagBlock, LowDiagBlock = LowDiagBlock,
        )
    v && println("UPDATE: Local operators created; ", keys(out))
    return out
end

"""
Creates the DG approximation to the generator `B`.

    MakeLazyGenerator(
        model::SFFM.Model,
        mesh::DGMesh,
        Matrices::NamedTuple;
    )

# Arguments
- `model`: A Model object
- `mesh`: A Mesh object
- `Matrices`: A Matrices tuple from `MakeMatrices`

# Output
- A Generator with fields `:BDict, :B, :QBDidx`
    - `:BDict::Dict{String,Array{Float64,2}}`: a dictionary storing Bᵢⱼˡᵐ with
        keys string(i,j,ℓ,m), and values Bᵢⱼˡᵐ, i.e. `B.BDict["12+-"]` = B₁₂⁺⁻
    - `:B::SparseArrays.SparseMatrixCSC{Float64,Int64}`:
        `NPhases(model)*TotalNBases(mesh)×NPhases(model)*TotalNBases(mesh)`, the
        global approximation to `B`
    - `:QBDidx::Array{Int64,1}`: `NPhases(model)*TotalNBases(mesh)×1` vector of
        integers such such that `:B[QBDidx,QBDidx]` puts all the blocks relating
        to cell `k` next to each other
"""
function MakeLazyGenerator(
    model::SFFM.Model,
    mesh::DGMesh;
    v::Bool = false,
)

    m = local_dg_operators(mesh; v=v)
    blocks = (m.LowDiagBlock*m.MInv*2, (m.G+m.PosDiagBlock)*m.MInv*2,
        -(m.G+m.NegDiagBlock)*m.MInv*2, m.UpDiagBlock*m.MInv*2)

    boundary_flux = (
        upper = (in = (m.Dw.DwInv * m.Phi[end, :]*2)[:],
                out = (m.Phi[end, :]' * m.Dw.Dw * m.MInv)[:]),
        lower = (in = (m.Dw.DwInv * m.Phi[1, :]*2)[:],
                out = (m.Phi[1, :]' * m.Dw.Dw * m.MInv)[:])
    )

    T = model.T
    C = model.C
    delta = Δ(mesh)
    D = LinearAlgebra.I(SFFM.NBases(mesh))
    pmidx = falses(size(T))

    out = SFFM.LazyGenerator(blocks,boundary_flux,T,C,delta,D,pmidx,mesh.Fil)
    v && println("UPDATE: LazyGenerator object created with keys ", keys(out))
    return out
end

"""
Uses Eulers method to integrate the matrix DE ``f'(x) = f(x)D`` to
approxiamte ``f(y)``.

    EulerDG(
        D::Union{Array{<:Real,2},SparseArrays.SparseMatrixCSC{Float64,Int64}},
        y::Real,
        x0::Array{<:Real};
        h::Float64 = 0.0001,
    )

# Arguments
- `D::Union{Array{<:Real,2},SparseArrays.SparseMatrixCSC{Float64,Int64}}`:
    the matrix ``D`` in the system of ODEs ``f'(x) = f(x)D``.
- `y::Real`: the value where we want to evaluate ``f(y)``.
- `x0::Array{<:Real}`: a row-vector initial condition.
- `h::Float64`: a stepsize for theEuler scheme.

# Output
- `f(y)::Array`: a row-vector approximation to ``f(y)``
"""
function EulerDG(
    D::AbstractArray{<:Real,2},
    y::Real,
    x0::AbstractArray{<:Real,2};
    h::Float64 = 0.0001,
)
    x = x0
    for t = h:h:y
        dx = h * (x * D)
        x = x + dx
    end
    return x
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
using GeodesicLearning
using Test

@testset "GeodesicLearning.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: chelate/GeodesicLearning

max_issues_repo_name: chelate/GeodesicLearning

max_forks_repo_name: chelate/GeodesicLearning


text:
export merge_map

import Base: show
import DataStructures: Deque

"""
    merge_map(::Type{R}, mappingFn::F = identity; concurrent::Int = typemax(Int)) where { R, F <: Function }

Creates a `merge_map` operator, which returns an Observable that emits the result of applying the projection
function `mappingFn` to each item emitted by the source Observable and merging the results of the Observables
obtained from this transformation.

# Arguments
- `::Type{R}`: the type of data of output Observables after projection with `mappingFn`
- `mappingFn::F`: projection function with `(data) -> Observable{R}` signature
- `concurrent::Int`: optional, default is `typemax(Int)`, maximum number of input Observables being subscribed to concurrently

# Producing

Stream of type `<: Subscribable{R}`

# Examples
```julia
using Rocket

source = from([ 0, 0 ]) |> merge_map(Int, d -> from([ 1, 2, 3 ], scheduler = AsyncScheduler(0)))
subscribe!(source, logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Data: 3
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`RightTypedOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
merge_map(::Type{R}, mappingFn::F = identity; concurrent::Int = typemax(Int)) where { R, F <: Function } = MergeMapOperator{R, F}(mappingFn, concurrent)

struct MergeMapOperator{R, F} <: RightTypedOperator{R}
    mappingFn  :: F
    concurrent :: Int
end

function on_call!(::Type{L}, ::Type{R}, operator::MergeMapOperator{R, F}, source) where { L, R, F }
    return proxy(R, source, MergeMapProxy{L, R, F}(operator.mappingFn, operator.concurrent))
end

struct MergeMapProxy{L, R, F} <: ActorSourceProxy
    mappingFn  :: F
    concurrent :: Int
end

actor_proxy!(::Type{R}, proxy::MergeMapProxy{L, R, F}, actor::A) where { L, R, F, A } = MergeMapActor{L, R, F, A}(proxy.mappingFn, proxy.concurrent, actor)

# m - main
mutable struct MergeMapActorProps{L}
    msubscription :: Teardown
    ismcompleted  :: Bool
    isdisposed    :: Bool

    active_listeners :: Vector{Any}
    pending_data     :: Deque{L}

    MergeMapActorProps{L}() where L = new(voidTeardown, false, false, Vector{Any}(), Deque{L}())
end

struct MergeMapActor{L, R, F, A} <: Actor{L}
    mappingFn  :: F
    concurrent :: Int
    actor      :: A
    props      :: MergeMapActorProps{L}

    MergeMapActor{L, R, F, A}(mappingFn::F, concurrent::Int, actor::A) where { L, R, F, A } = begin
        return new(mappingFn, concurrent, actor, MergeMapActorProps{L}())
    end
end

getactive(actor::MergeMapActor)  = actor.props.active_listeners
getpending(actor::MergeMapActor) = actor.props.pending_data

pushactive!(actor::MergeMapActor, active)   = push!(getactive(actor), active)
pushpending!(actor::MergeMapActor, pending) = push!(getpending(actor), pending)

isdisposed(actor::MergeMapActor)   = actor.props.isdisposed
ismcompleted(actor::MergeMapActor) = actor.props.ismcompleted
isicompleted(actor::MergeMapActor) = length(getactive(actor)) + length(getpending(actor)) === 0

setmcompleted!(actor::MergeMapActor, value::Bool) = actor.props.ismcompleted = value

function seticompleted!(actor::MergeMapActor, inner)
    filter!(listener -> listener !== inner, getactive(actor))
    if length(getpending(actor)) !== 0
        attach_source_with_map!(actor, popfirst!(getpending(actor)))
    elseif length(getactive(actor)) === 0 && ismcompleted(actor)
        complete!(actor)
    end
end

function attach_source_with_map!(actor::M, data::L) where { L, R, M <: MergeMapActor{L, R} }
    if length(getactive(actor)) < actor.concurrent
        inner = MergeMapInnerActor{R, M}(actor, MergeMapInnerActorProps())
        pushactive!(actor, inner)
        setsubscription!(inner, subscribe!(actor.mappingFn(data), inner))
    else
        pushpending!(actor, data)
    end
end

mutable struct MergeMapInnerActorProps
    subscription :: Teardown

    MergeMapInnerActorProps() = new(voidTeardown)
end

struct MergeMapInnerActor{R, M} <: Actor{R}
    main  :: M
    props :: MergeMapInnerActorProps
end

getsubscription(actor::MergeMapInnerActor)                = actor.props.subscription
setsubscription!(actor::MergeMapInnerActor, subscription) = actor.props.subscription = subscription

on_next!(actor::MergeMapInnerActor{R}, data::R) where R = next!(actor.main.actor, data)
on_error!(actor::MergeMapInnerActor, err)               = error!(actor.main, err)
on_complete!(actor::MergeMapInnerActor)                 = begin
    seticompleted!(actor.main, actor)
end

function on_next!(actor::M, data::L) where { L, R, M <: MergeMapActor{L, R} }
    if !isdisposed(actor)
        attach_source_with_map!(actor, data)
    end
end

function on_error!(actor::MergeMapActor, err)
    if !isdisposed(actor)
        dispose!(actor)
        error!(actor.actor, err)
    end
end

function on_complete!(actor::MergeMapActor)
    setmcompleted!(actor, true)
    if !isdisposed(actor) && isicompleted(actor)
        dispose!(actor)
        complete!(actor.actor)
    end
end

function dispose!(actor::MergeMapActor)
    actor.props.isdisposed = true
    unsubscribe!(actor.props.msubscription)
    for listener in getactive(actor)
        unsubscribe!(getsubscription(listener))
    end
    empty!(getactive(actor))
    empty!(getpending(actor))
end

struct MergeMapSource{L, S} <: Subscribable{L}
    source :: S
end

source_proxy!(::Type{R}, proxy::MergeMapProxy{L, R, F}, source::S) where { L, R, F, S } = MergeMapSource{L, S}(source)

function on_subscribe!(source::MergeMapSource, actor::MergeMapActor)
    actor.props.msubscription = subscribe!(source.source, actor)
    return MergeMapSubscription(actor)
end

struct MergeMapSubscription{A} <: Teardown
    actor :: A
end

as_teardown(::Type{ <: MergeMapSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::MergeMapSubscription)
    dispose!(subscription.actor)
    return nothing
end

Base.show(io::IO, ::MergeMapOperator{R})   where {    R } = print(io, "MergeMapOperator($R)")
Base.show(io::IO, ::MergeMapProxy{L, R})   where { L, R } = print(io, "MergeMapProxy($L, $R)")
Base.show(io::IO, ::MergeMapActor{L, R})   where { L, R } = print(io, "MergeMapActor($L -> $R)")
Base.show(io::IO, ::MergeMapInnerActor{R}) where {    R } = print(io, "MergeMapInnerActor($R)")
Base.show(io::IO, ::MergeMapSource{S})     where S        = print(io, "MergeMapSource($S)")
Base.show(io::IO, ::MergeMapSubscription)                 = print(io, "MergeMapSubscription()")


text:
using DataStructures, Test
import DataStructures: IntSet

@testset "IntSet" begin
    @testset "Construction, collect" begin
        data_in = (1,5,100)
        s = IntSet(data_in)
        data_out = collect(s)
        @test all(map(d->in(d,data_out), data_in))
        @test length(data_out) == length(data_in)
    end

    @testset "eltype, empty" begin
        @test eltype(IntSet()) == Int
        @test eltype(typeof(IntSet())) == Int
        @test isequal(empty(IntSet([1,2,3])), IntSet())
    end

    @testset "show" begin
        @test sprint(show, IntSet()) == "IntSet([])"
        @test sprint(show, IntSet([1,2,3])) == "IntSet([1, 2, 3])"
        @test occursin("...,", sprint(show, complement(IntSet())))
    end

    @testset "Core Functionality" begin
        s = IntSet([0,1,10,20,200,300,1000,10000,10002])
        @test last(s) == 10002
        @test first(s) == 0
        @test length(s) == 9
        @test pop!(s) == 10002
        @test length(s) == 8
        @test popfirst!(s) == 0
        @test length(s) == 7
        @test !in(0,s)
        @test !in(10002,s)
        @test in(10000,s)
        @test_throws ArgumentError first(IntSet())
        @test_throws ArgumentError last(IntSet())
        t = copy(s)
        sizehint!(t, 20000) #check that hash does not depend on size of internal Array{UInt32, 1}
        @test hash(s) == hash(t)
        @test hash(complement(s)) == hash(complement(t))
    end

    @testset "setdiff / symdiff" begin
        @test setdiff(IntSet([1, 2, 3, 4]), IntSet([2, 4, 5, 6])) == IntSet([1, 3])
        @test symdiff(IntSet([1, 2, 3, 4]), IntSet([2, 4, 5, 6])) == IntSet([1, 3, 5, 6])
    end

    @testset "setdiff!" begin
        s2 = IntSet([1, 2, 3, 4])
        setdiff!(s2, IntSet([2, 4, 5, 6]))

        @test s2 == IntSet([1, 3])
    end

    @testset "issue #7851" begin
        @test_throws ArgumentError IntSet(-1)
        @test !(-1 in IntSet(0:10))
    end

    #@testset "issue #8570" begin
    #    # This requires 2^29 bytes of storage, which is too much for a simple test
    #    # s = IntSet(2^32)
    #    @test length(s) == 1
    #    for b in s; b; end
    #end

    @testset "Copy, copy!, empty" begin
        s1 = IntSet([1,2,3])
        s2 = empty(s1)
        copy!(s2, s1)
        s3 = copy(s2)
        @test s3 == s2 == s1
        @test collect(s3) == collect(s2) == [1,2,3]
    end

    @testset "complement! / pop!" begin
        c1 = complement!(IntSet())
        s1 = IntSet([1,2,3])

        pop!(c1, 1)
        pop!(c1, 2)
        pop!(c1, 3)
        c2 = empty(c1)
        copy!(c2, c1)
        c3 = copy(c2)
        c4 = complement(s1)
        @test c1 == c2 == c3 == c4
        @test c4 === sizehint!(c4, 100)
        @test c1 == c4
        @test last(c1) == typemax(Int)-1
        @test last(complement(IntSet())) == typemax(Int)-1
    end

    @testset "Push, union" begin
        # Push, union
        s1 = IntSet()
        @test_throws ArgumentError push!(s1, -1)
        push!(s1, 1, 10, 100, 1000)
        @test collect(s1) == [1, 10, 100, 1000]
        push!(s1, 606)
        @test collect(s1) == [1, 10, 100, 606, 1000]
        s2 = IntSet()
        @test s2 === union!(s2, s1)
        s3 = IntSet([1, 10, 100])
        union!(s3, [1, 606, 1000])
        s4 = union(IntSet([1, 100, 1000]), IntSet([10, 100, 606]))
        @test s1 == s2 == s3 == s4

        c1 = complement(s1)
        @test !(1 in c1)
        push!(c1, 1)
        @test 1 in c1
        push!(c1, 10, 100, 10)
        @test collect(complement(c1)) == [606, 1000]
        c2 = complement(IntSet([606, 1000]))
        @test c2 === union!(c2, c1)
        c3 = union!(complement(IntSet([10, 606, 1000])), complement(IntSet([1, 606, 1000, 2000])))
        @test c3 == union(complement(IntSet([1, 606, 1000, 2000])), complement(IntSet([10, 606, 1000])))
        c4 = union!(complement(IntSet([1, 10, 606, 1000, 1001])), complement(IntSet([606, 1000])))
        @test c4 == union(complement(IntSet([606, 1000])), complement(IntSet([1, 10, 606, 1000, 1001])))
        c5 = union!(complement(IntSet([10, 606, 1000])), IntSet([0, 10, 20, 30]))
        @test c5 == union(IntSet([0, 10, 20, 30]), complement(IntSet([10, 606, 1000])))
        c6 = union!(complement(IntSet([10, 606, 1000])), IntSet([10, 4000, 3]))
        @test c6 == union(IntSet([10, 4000, 3]), complement(IntSet([10, 606, 1000])))
        @test c1 == c2 == c3 == c4 == c5 == c6
    end

    @testset "pop!, delete!" begin
        s = IntSet(1:2:10)
        @test pop!(s, 1) == 1
        @test !(1 in s)
        @test_throws KeyError pop!(s, 1)
        @test_throws ArgumentError pop!(s, -1)
        @test_throws ArgumentError pop!(s, -1, 1)
        @test_throws ArgumentError pop!(()->throw(ErrorException()), s, -1)
        @test pop!(s, 1, 0) == 0
        @test s === delete!(s, 1)
        for i in s; pop!(s, i); end
        @test isempty(s)
        x = 0
        @test 1 == pop!(()->(x+=1), s, 100)
        @test x == 1
        push!(s, 100)
        @test pop!(()->throw(ErrorException()), s, 100) == 100
        push!(s, 1:2:10...)
        @test pop!(s) == 9
        @test pop!(s) == 7
        @test popfirst!(s) == 1
        @test popfirst!(s) == 3
        @test collect(s) == [5]
        empty!(s)
        @test isempty(s)

        c = complement(IntSet())
        @test pop!(c, 1) == 1
        @test !(1 in c)
        @test_throws KeyError pop!(c, 1)
        @test_throws ArgumentError pop!(c, -1)
        @test_throws ArgumentError pop!(c, -1, 1)
        @test_throws ArgumentError pop!(()->throw(ErrorException()), c, -1)
        @test pop!(c, 1, 0) == 0
        @test c === delete!(c, 1)
        @test popfirst!(c) == 0
        @test popfirst!(c) == 2
        @test_throws ArgumentError pop!(c)
        @test collect(complement(c)) == [0,1,2]
        @test empty!(c) == IntSet()
    end

    @testset "Intersect" begin
        @test isempty(intersect(IntSet()))
        @test isempty(intersect(IntSet(1:10), IntSet()))
        @test isempty(intersect(IntSet(), IntSet(1:10)))
        @test isempty(intersect(IntSet(), complement(IntSet())))
        @test isempty(intersect(IntSet(), complement(IntSet(1:10))))
        @test isempty(intersect(complement(IntSet()), IntSet()))
        @test isempty(intersect(complement(IntSet(1:10)), IntSet()))

        @test intersect(IntSet([1,2,3])) == IntSet([1,2,3])
        @test intersect(complement!(IntSet()), IntSet(1)) ==
              intersect(IntSet(1), complement!(IntSet())) == IntSet(1)

        @test intersect(IntSet(0:7), IntSet(3:10)) ==
              intersect(IntSet(3:10), IntSet(0:7)) == IntSet(3:7)
        @test intersect(complement(IntSet([0:2; 11:16])), IntSet(0:7)) ==
              intersect(IntSet(0:7), complement(IntSet([0:2; 11:16]))) == IntSet(3:7)

        @test intersect(complement(IntSet(5:12)), complement(IntSet(7:10))) ==
              intersect(complement(IntSet(7:10)), complement(IntSet(5:12))) == complement(IntSet(5:12))

        @test intersect(IntSet(0:10), IntSet(1:4), 0:5, [1,2,10]) == IntSet(1:2)
    end

    @testset "Setdiff" begin
        s1 = IntSet(1:100)
        setdiff!(s1, IntSet(1:2:100))
        s2 = setdiff(IntSet(1:100), IntSet(1:2:100))
        @test s1 == s2 == IntSet(2:2:100)
        @test collect(s1) == collect(2:2:100)

        s1 = IntSet(1:10)
        s2 = complement(IntSet(3:5))
        @test setdiff(s1, s2) == setdiff(s1, [0:2; 6:100]) == IntSet(3:5)
        @test isempty(setdiff(complement(IntSet()), complement(IntSet())))
        @test setdiff(complement(IntSet()), complement(IntSet(3:5))) == IntSet(3:5)
        @test setdiff(complement(IntSet(1:5)), complement(IntSet(3:10))) == IntSet([6, 7, 8, 9, 10])
        @test setdiff(complement(IntSet(2:2:10)), IntSet(1:5)) == complement(IntSet([1:5; 6:2:10]))
        @test setdiff!(complement(IntSet()), complement(IntSet())) == IntSet()
        @test setdiff!(complement(IntSet(0:2:10)), complement(IntSet(0:10))) == IntSet(1:2:9)
    end

    @testset "Symdiff" begin
        @test symdiff(IntSet([1, 2, 3, 4]), IntSet([2, 4, 5, 6])) ==
              symdiff(IntSet([2, 4, 5, 6]), IntSet([1, 2, 3, 4])) ==
              symdiff(IntSet([1, 2, 3, 4]), [2, 4, 5, 6]) ==
              symdiff(IntSet([2, 4, 5, 6]), [1, 2, 3, 4]) == IntSet([1, 3, 5, 6])
        @test symdiff(complement(IntSet()), IntSet(2:3)) ==
              symdiff(IntSet(2:3), complement(IntSet())) ==
              symdiff(complement(IntSet()), 2:3) ==complement(IntSet(2:3))

        @test symdiff(complement(IntSet(2:7)), IntSet(5:10)) ==
              symdiff(IntSet(5:10), complement(IntSet(2:7))) == complement(IntSet([2:4; 8:10]))
        @test symdiff(complement(IntSet(3:7)), complement(IntSet(5:10))) == IntSet([3:4; 8:10])
    end

    @testset "Subsets, equality" begin
        @test IntSet(2:2:10) < IntSet(1:10)
        @test !(IntSet(2:2:10) < IntSet(2:2:10))
        @test IntSet(2:2:10) <= IntSet(2:10)
        @test IntSet(2:2:10) <= IntSet(2:2:10)
        @test IntSet(1) < complement!(IntSet())
        @test IntSet(1) <= complement!(IntSet())
        @test !(IntSet(1) < complement!(IntSet(1)))
    end

    @testset "Test logic against Set" begin
        p = IntSet([0,1,4,5])
        resize!(p.bits, 6)
        q = IntSet([0,2,4,6])
        resize!(q.bits, 8)
        p′ = complement(p)
        q′ = complement(q)
        function collect10(itr)
            r = eltype(itr)[]
            for i in itr
                i > 10 && break
                push!(r, i)
            end
            r
        end
        a = Set(p)
        b = Set(q)
        a′ = Set(collect10(p′))
        b′ = Set(collect10(q′))
        for f in (union, intersect, setdiff, symdiff)
            @test collect(f(p, p)) == sort(collect(f(a, a)))
            @test collect(f(q, q)) == sort(collect(f(b, b)))
            @test collect(f(p, q)) == sort(collect(f(a, b)))
            @test collect(f(q, p)) == sort(collect(f(b, a)))

            @test collect10(f(p′, p)) == sort(collect(f(a′, a)))
            @test collect10(f(q′, q)) == sort(collect(f(b′, b)))
            @test collect10(f(p′, q)) == sort(collect(f(a′, b)))
            @test collect10(f(q′, p)) == sort(collect(f(b′, a)))

            @test collect10(f(p, p′)) == sort(collect(f(a, a′)))
            @test collect10(f(q, q′)) == sort(collect(f(b, b′)))
            @test collect10(f(p, q′)) == sort(collect(f(a, b′)))
            @test collect10(f(q, p′)) == sort(collect(f(b, a′)))

            @test collect10(f(p′, p′)) == sort(collect(f(a′, a′)))
            @test collect10(f(q′, q′)) == sort(collect(f(b′, b′)))
            @test collect10(f(p′, q′)) == sort(collect(f(a′, b′)))
            @test collect10(f(q′, p′)) == sort(collect(f(b′, a′)))
        end
    end

    @testset "other 1" begin
        s = IntSet()
        push!(s, 0, 2, 100)
        @test 0 in s
        @test !(1 in s)
        @test 2 in s
        @test 100 in s
        @test !(101 in s)
        @test !(1000 in s)
        @test first(s) == 0
        @test last(s) == 100
        @test s == IntSet([0, 2, 100])
        push!(s, 1000)
        @test [i for i in s] == [0, 2, 100, 1000]
        @test pop!(s) == 1000
        @test s == IntSet([0, 2, 100])
        @test hash(s) == hash(IntSet([0, 2, 100]))
    end

    @testset "other 2" begin
        b = 0:1000
        s = IntSet(b)
        @test collect(s) == collect(b)
        @test length(s) == length(b)
        @test pop!(s, 100) == 100
        @test collect(s) == [0:99; 101:1000]
        @test_throws KeyError pop!(s, 100)
        @test pop!(s, 100, 0) == 0
        @test pop!(s, 99, 0) == 99
        @test pop!(()->1, s, 99) == 1
        @test pop!(()->1, s, 98) == 98
    end

    @testset "show" begin
        show(IOBuffer(), IntSet())
        show(IOBuffer(), complement(IntSet()))
    end

end # @testset IntSet

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    BinaryMaximum(x::AbstractIntVar, y::AbstractIntVar, z::AbstractIntVar, trailer)

BinaryMaximum constraint, states that `x == max(y, z)`
"""
struct BinaryMaximumBC <: Constraint
    x::AbstractIntVar
    y::AbstractIntVar
    z::AbstractIntVar
    active::StateObject{Bool}

    function BinaryMaximumBC(x::AbstractIntVar, y::AbstractIntVar, z::AbstractIntVar, trailer)
        constraint = new(x, y, z, StateObject{Bool}(true, trailer))
        for xi in [x, y, z]
            addOnDomainChange!(xi, constraint)
        end
        return constraint
    end
end

"""
    propagate!(constraint::BinaryMaximumBC, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`BinaryMaximumBC` propagation function.
"""
function propagate!(constraint::BinaryMaximumBC, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    min_x, min_y, min_z = minimum(constraint.x.domain), minimum(constraint.y.domain), minimum(constraint.z.domain)
    max_x, max_y, max_z = maximum(constraint.x.domain), maximum(constraint.y.domain), maximum(constraint.z.domain)
    max_of_min = max(min_y, min_z)
    max_of_max = max(max_y, max_z)
    min_of_min = min(min_y, min_z)

    # feasibility
    if min_x > max_of_max || max_x < min_of_min
        return false
    end

    # y & z pruning x
    if min_x < max_of_min || max_x > max_of_max
        prunedX = vcat(removeBelow!(constraint.x.domain, max_of_min), removeAbove!(constraint.x.domain, max_of_max))
        addToPrunedDomains!(prunedDomains, constraint.x, prunedX)
        triggerDomainChange!(toPropagate, constraint.x)
    end

    # x pruning y & z
    if min_x > max_z # if x can only be associated to y
        # x trying to prune y
        prunedY = vcat(removeBelow!(constraint.y.domain, min_x), removeAbove!(constraint.y.domain, max_x))
        if !isempty(prunedY)
            addToPrunedDomains!(prunedDomains, constraint.y, prunedY)
            triggerDomainChange!(toPropagate, constraint.y)
        end
    elseif min_x > max_y # if x can only be associated to z
        # x trying to prune z
        prunedZ = vcat(removeBelow!(constraint.z.domain, min_x), removeAbove!(constraint.z.domain, max_x))
        if !isempty(prunedZ)
            addToPrunedDomains!(prunedDomains, constraint.z, prunedZ)
            triggerDomainChange!(toPropagate, constraint.z)
        end
    end

    # deactivation
    if isbound(constraint.x) + isbound(constraint.y) + isbound(constraint. z) >= 2
        setValue!(constraint.active, false)
    end

    return !isempty(constraint.x.domain) & !isempty(constraint.y.domain) & !isempty(constraint.z.domain)
end

variablesArray(constraint::BinaryMaximumBC) = [constraint.x, constraint.y, constraint.z]

function Base.show(io::IO, ::MIME"text/plain", con::BinaryMaximumBC)
    println(io, string(typeof(con)), ": ", con.x.id, " == max(", con.y.id, con.z.id, "), ", "active = ", con.active)
    println(io, "   ", con.x)
    println(io, "   ", con.y)
    print(io, "   ", con.z)
end

function Base.show(io::IO, con::BinaryMaximumBC)
    print(io, typeof(con), ": ", con.x.id, " == max(", con.y.id, ", ", con.z.id, ")")
end


text:
function Base.foreach(f, d::AbstractDictionary, d2::AbstractDictionary, ds::AbstractDictionary...)
    if sharetokens(d, d2, ds...)
        @inbounds for t in tokens(d)
            f(gettokenvalue(d, t), gettokenvalue(d2, t), map(x -> @inbounds(gettokenvalue(x, t)), ds)...)
        end
    else
        @inbounds for i in keys(d)
            f(d[i], d2[i], map(x -> @inbounds(x[i]), ds)...)
        end
    end
    return nothing
end

# (avoid an _apply for the two-input case)
function Base.foreach(f, d::AbstractDictionary, d2::AbstractDictionary)
    if sharetokens(d, d2)
        @inbounds for t in tokens(d)
            f(gettokenvalue(d, t), gettokenvalue(d2, t))
        end
    else
        @inbounds for i in keys(d)
            f(d[i], d2[i])
        end
    end
    return nothing
end

function Base.foreach(f, d::AbstractDictionary)
    @inbounds for x in d
        f(x)
    end
    return nothing
end


text:
@testset "get_prenex_normal_form" begin
    @customtest String(get_prenex_normal_form(Formula("(~?{y}*{g}P(f(y))) & (?{x}P(x) | *{y}P(y))"))) ==
        "*{##5}?{##6}?{##9}*{##10}(~P(f(##5)) & (P(##9) | P(##10)))"
    @customtest String(get_prenex_normal_form(Formula("*{x}?{y}P(x, g(y, f(x))) | ~Q(z) | ~*{x}R(x,y)"))) ==
        "*{##1}?{##2}?{##7}(P(##1, g(##2, f(##1))) | (~Q(z) | ~R(##7, y)))"
end


text:
using Gym, Test

include("env.jl")


text:
module DelegationPattern

# account
export Account
export account_number, balance, date_opened
export deposit!, withdraw!, transfer!

# savings account
export SavingsAccount, SavingsAccount2
export interest_rate, accrue_daily_interest!

# package dependencies
using Dates: Date
using Lazy: @forward

include("Account.jl")
include("SavingsAccount.jl")
include("SavingsAccount2.jl")

end # module


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

include("utils/TikzUtils.jl")

using DataStructures

using TikzPictures
import PGFPlots: Plots, Axis, GroupPlot

using Base.Test

function make_fieldmap(num_aircraft::Int)

  M = OrderedDict{ASCIIString, Tuple{Int64, ASCIIString}}()

  for i = 1:num_aircraft

    i1 = i - 1

    M["State.x.$i1"] = (i, "wm.x")
    M["State.y.$i1"] = (i, "wm.y")
    M["State.h.$i1"] = (i, "wm.z")
    M["State.dx.$i1"] = (i, "wm.vx")
    M["State.dy.$i1"] = (i, "wm.vy")
    M["State.dh.$i1"] = (i, "wm.vz")
    M["Input.own.dz.$i1"] = (i, "ra_detailed.ownInput.dz")
    M["Input.own.z.$i1"] = (i, "ra_detailed.ownInput.z")
    M["Input.own.psi.$i1"] = (i, "ra_detailed.ownInput.psi")
    M["Input.own.h.$i1"] = (i, "ra_detailed.ownInput.h")
    M["Input.own.modes.$i1"] = (i, "ra_detailed.ownInput.modes")

    for j = 1:(num_aircraft - 1)

      j1 = j - 1 #LL uses 0 indexing

      #start from 0, count to n-1, skip self
      lid = getll_id(i1, j1)

      M["Input.intruder.modes.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].modes")
      M["Input.intruder.sr.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].sr")
      M["Input.intruder.chi.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].chi")
      M["Input.intruder.z.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].z")
      M["Input.intruder.cvc.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].cvc")
      M["Input.intruder.vrc.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].vrc")
      M["Input.intruder.vsb.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].vsb")
      M["Input.intruder.equipage.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].equipage")
      M["Input.intruder.protection\_mode.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].protection\_mode")
      M["Input.intruder.quant.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].quant")
      M["Input.intruder.sensitivity\_index.$i1.$lid"] = (i, "ra_detailed.intruderInput[$j].sensitivity\_index")
    end

    M["Output.cc.$i1"] = (i, "ra_detailed.ownOutput.cc")
    M["Output.vc.$i1"] = (i, "ra_detailed.ownOutput.vc")
    M["Output.ua.$i1"] = (i, "ra_detailed.ownOutput.ua")
    M["Output.da.$i1"] = (i, "ra_detailed.ownOutput.da")
    M["Output.target\_rate.$i1"] = (i, "ra_detailed.ownOutput.target\_rate")
    M["Output.turn\_off\_aurals.$i1"] = (i, "ra_detailed.ownOutput.turn\_off\_aurals")
    M["Output.crossing.$i1"] = (i, "ra_detailed.ownOutput.crossing")
    M["Output.alarm.$i1"] = (i, "ra_detailed.ownOutput.alarm")
    M["Output.alert.$i1"] = (i, "ra_detailed.ownOutput.alert")
    M["Output.dh\_min.$i1"] = (i, "ra_detailed.ownOutput.dh\_min")
    M["Output.dh\_max.$i1"] = (i, "ra_detailed.ownOutput.dh\_max")
    M["Output.sensitivity\_index.$i1"] = (i, "ra_detailed.ownOutput.sensitivity\_index")
    M["Output.ddh.$i1"] = (i, "ra_detailed.ownOutput.ddh")

    for j = 1:(num_aircraft - 1)

      j1 = j - 1 #LL uses 0 indexing

      #start from 0, skip self
      lid = getll_id(i1, j1)

      M["Output.intruder.cvc.$i1.$lid"] = (i, "ra_detailed.intruderOutput[$j].cvc")
      M["Output.intruder.vrc.$i1.$lid"] = (i, "ra_detailed.intruderOutput[$j].vrc")
      M["Output.intruder.vsb.$i1.$lid"] = (i, "ra_detailed.intruderOutput[$j].vsb")
      M["Output.intruder.tds.$i1.$lid"] = (i, "ra_detailed.intruderOutput[$j].tds")
      M["Output.intruder.code.$i1.$lid"] = (i, "ra_detailed.intruderOutput[$j].code")
    end
  end

  return M
end

#assumes that both ids start counts at 0
getll_id(listowner_id::Int, item::Int) = item >= listowner_id ? item + 1 : item

function pgfplotField{T <: Real}(llvec::Vector{T}, slvec::Vector{T}, var_name::AbstractString = "")

  g = GroupPlot(2, 1, groupStyle = "horizontal sep = 2.2cm, vertical sep = 2.2cm")

  #plot overlay
  plotArray = Plots.Plot[]
  push!(plotArray, Plots.Linear([1:length(llvec)], llvec,
                               style="mark options={color=blue}", mark="*"))
  push!(plotArray, Plots.Linear([1:length(slvec)], slvec,
                               style="mark options={color=red}", mark="*"))
  ax = Axis(plotArray,
            xlabel = "Timestep",
            ylabel = var_name,
            title = "Overlay")
  push!(g, ax)

  #plot difference
  plotArray = Plots.Plot[]
  t_end = min(length(llvec), length(slvec))
  push!(plotArray,Plots.Linear([1:t_end], abs(llvec[1:t_end] - slvec[1:t_end]),
                               style="mark options={color=blue}", mark="*"))
  ax = Axis(plotArray,
            xlabel = "Timestep",
            ylabel = var_name,
            title = "Absolute Difference")
  push!(g, ax)

  tp = PGFPlots.plot(g)
  use_geometry_package!(tp, landscape = false) # portrait

  return tp
end

function getvec(data::Array{Any, 2}, field::AbstractString)

  v = []
  headers = data[1, :] # 1st row is headers
  i = findfirst(x -> x == field, headers)

  if i != 0
    try
      v = demote_vectype(data[2:end, i]) #convert to lowest common type
      v = convert(Vector{Float64}, v) #data starts at row 2
    catch
      warn("Cannot convert $field: $(data[2:end, i])")
      v = []
    end
  end

  return v
end

function demote_vectype(v::Vector)
  # Converts vector type to lowest common type amongst elements
  # Fixes the overly general container type problem, i.e., Vector{Any}

  T = promote_type(map(typeof, v)...) #determine the lowest common type

  return convert(Vector{T}, v)
end

import Base.convert
function convert{T<:AbstractString}(::Type{Vector{Float64}}, v::Vector{T})

  vv = map(lowercase, v)
  vv = map(parse, vv)
  vv = map(float64, vv)

  return vv
end

function namegen(fileroot::AbstractString, indices::Vector{Int64}, extension::AbstractString)
  map(i -> string(fileroot, i, extension), indices)
end

function reorder2ll{T<:AbstractString}(sislesfiles::Vector{T}, llfile::AbstractString;
                                  outfile::AbstractString = "reorderedcsv.csv")

  num_aircraft = length(sislesfiles)
  fieldmap = make_fieldmap(num_aircraft)

  lldat = readdlm(llfile)
  llheaders = squeeze(lldat[1, :], 1)

  sldats = map(readcsv, sislesfiles)

  # skip units row
  map!(d -> d[[1,3:end], :], sldats)

  outcsv = Array(Any, size(sldats[1], 1), 0)

  for llfield in llheaders

    if llfield == "time" #skip...
        continue
    end

    if !haskey(fieldmap, llfield)
      error("reorder2ll::No such key! $llfield")
    end

    slfile, slfield = fieldmap[llfield]

    slvec = getvec(sldats[slfile], slfield)

    if isempty(slvec)
      error("$slfield not found in SISLES csv file.")
    end

    outcsv = hcat(outcsv, vcat(slfield, slvec)) #TODO: remove the dynamic growing...
  end

  f = open(outfile, "w")
  writecsv(f, outcsv)
  close(f)

  return outcsv
end

function reorderedcsvdiff(slcsv::AbstractString, llcsv::AbstractString; thresh::Float64=1.0, use_ll_names::Bool=false)

  fileroot, fileext = splitext(slcsv)
  outfile = string(fileroot, "_emph", fileext)

  sldat = readcsv(slcsv)
  lldat = readdlm(llcsv)
  lldat = lldat[:, 2:end] #skip time col

  if use_ll_names
    sldat[1, :] = lldat[1, :]
  end

  #TODO: avoid the loops
  for c = 1:size(lldat, 2)
    for r = 2:size(lldat, 1) #skip header

      a = abs(lldat[r, c]-sldat[r, c])
      sldat[r, c] = a >= thresh ? a : 0

    end
  end

  f = open(outfile, "w")
  writecsv(f, sldat)
  close(f)

  return sldat
end

function ll_csv_diff{T<:AbstractString}(sislesfiles::Vector{T}, llfile::AbstractString;
                                  slskipfirst::Bool = false,
                                  outfile::AbstractString = "csvdiff.pdf")

  num_aircraft = length(sislesfiles)
  fieldmap = make_fieldmap(num_aircraft)

  lldat = readdlm(llfile)
  sldats = map(readcsv, sislesfiles)

  firstdiff = Tuple{ASCIIString,Int64}[]

  # skip units row
  map!(d -> d[[1,3:end], :], sldats)

  if slskipfirst
    map!(d -> d[[1,3:end], :], sldats)
  end

  td = TikzDocument()

  for (llfield, slfieldtup) in fieldmap

    slfile, slfield = slfieldtup

    llvec = getvec(lldat, llfield)
    slvec = getvec(sldats[slfile], slfield)

    if isempty(llvec) || isempty(slvec)
      llnotfound = isempty(llvec) ? "$llfield not found in LL file." : ""
      slnotfound = isempty(slvec) ? "$slfield not found in SISLES csv file." : ""

      warn("$llnotfound $slnotfound")

      continue
    end

    llfield = replace(llfield, "\_", "\\_")
    slfield = replace(slfield, "\_", "\\_")
    tp = pgfplotField(llvec, slvec, llfield)

    cap = "Variable: $llfield vs. $slfield"

    add_to_document!(td, tp, cap)

    ind = findfirst(x -> x >= 1.0, abs(llvec-slvec))
    push!(firstdiff, (llfield, ind))
  end

  TikzPictures.save(PDF(outfile), td)

  sort!(firstdiff, by = x -> x[2])

  return firstdiff
end

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((3, true, false, true),),
    ((3, false, true, true),),
    ((3, false, true, false),),
    ]
insts["fast"] = [
    ((10, true, false, true),),
    ((10, false, true, true),),
    ((10, false, true, false),),
    ((50, true, false, true),),
    ((50, false, true, true),),
    ((50, false, true, false),),
    ((200, true, false, true),),
    ((200, false, true, true),),
    ((200, false, true, false),),
    ]
return (PortfolioNative, insts)


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "GLVertexBuffer.jl" begin

   # function GLVertexBuffer()
   @testset "GLVertexBuffer" begin
      @test
      @test
      @test
      @test
   end

   # function GLVertexBuffer(vector::Vector{Float32})
   @testset "GLVertexBuffer" begin
      @test
      @test
      @test
      @test
   end

   # function releaseGpuResources(buffer::GLVertexBuffer)
   @testset "releaseGpuResources" begin
      @test
      @test
      @test
      @test
   end

   # function enableAttribute(location::Int32,buffer::GLVertexBuffer,num_components::Int64)
   @testset "enableAttribute" begin
      @test
      @test
      @test
      @test
   end

   # function disableAttribute(location::Int32,buffer::GLVertexBuffer)
   @testset "disableAttribute" begin
      @test
      @test
      @test
      @test
   end

end


text:
isdefined(Base, :__precompile__) && __precompile__()

module MAVs
using NLOptControl
using VehicleModels

include("CaseModule.jl")
using .CaseModule

include("AutonomousControl.jl")
#include("SharedControl.jl")
#include("SimpleModel.jl")
#include("PathFollowing.jl")

using .AutonomousControl
#using .SharedControl
#using .SimpleModel
#using .PathFollowing

export

      # CaseModule.jl
      Case,
      defineCase,
      setWeights!,
      setMisc!,
      setupResults,
      case2dfs,  # currently not working
      dataSet,
      Obs,
      defineObs,

      # AutonomousControl.jl
      initializeAutonomousControl,
      autonomousControl!,
      updateAutoParams!

      # SharedControl.jl
    #  initializeSharedControl!,
    #  checkFeasibility!,
    #  sharedControl!,
    #  getPlantData!,
    #  sendOptData!,
    #  ExternalModel,

      # SimpleModel.jl
  #    initializeSimpleModel,
    #  updateSimpleModel,
  #    runSimpleModel,
  #    compareSimpleModels,
#      resetDesignVariables,

      # PathFollowing.jl
    #  initializePathFollowing,
    #  updatePathParams!
end


text:
struct RecordCOVID
    # {
    #     "dateRep" : "26/10/2020",
    #     "day" : "26",
    #     "month" : "10",
    #     "year" : "2020",
    #     "cases" : 65,
    #     "deaths" : 3,
    #     "countriesAndTerritories" : "Afghanistan",
    #     "geoId" : "AF",
    #     "countryterritoryCode" : "AFG",
    #     "popData2019" : 38041757,
    #     "continentExp" : "Asia",
    #     "Cumulative_number_for_14_days_of_COVID-19_cases_per_100000" : "2.71806583"
    # }
    dateRep::Date
    day::String

end


text:
using Test, ArgParse
using SignalingDimension

try print_test_results catch err
    include("../utilities.jl")
    using .ScriptUtilities: print_test_results
end

# parsing command line arguments
arg_table = ArgParseSettings()
@add_arg_table! arg_table begin
    "-N"
        help = "max number of inputs/outputs to consider"
        arg_type = Int64
        default = 45
    "--dir"
        help = "directory to print results .txt to"
        arg_type = String
        default = "./script/facet_verifications/results/"
    "--stdout"
        help = "Print output to STDOUT. If not flagged, output is printed to `.txt` file."
        action = :store_true
end
args = parse_args(ARGS, arg_table)

# verification script: verifies all maximum likelihood facets with N ∈ [3, N_max]
function verify_maximum_likelihood_facet(N_max :: Int64)
    @testset "Testing all maximum likelihood facets of size `($N_max, $N_max)` or less." begin
        @testset "scanning over all `d ∈ [2, $N - 1]`" for N in 3:N_max
            @testset "Verifying N = $N and d = $d" for d in 2:N-1
                vertices = aff_ind_maximum_likelihood_vertices(N,d)
                facet = maximum_likelihood_facet(N,d)

                @test verify_facet(d, facet, vertices)
            end
        end
    end
end

# outputting result to file or stdout
print_test_results(
    verify_maximum_likelihood_facet,
    params=[args["N"]],
    stdout=args["stdout"],
    dir=args["dir"]
)


text:
using XGPaint
using Healpix

## Load halos from HDF5 files, establish a CIB model and cosmology
halo_pos, halo_mass = read_halo_catalog_hdf5(
    ENV["SCRATCH"] * "/websky_halos-light.hdf5")
cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0)
model = CIB_Planck2013{Float32}(nside=8192)

## Write one chunk to disk

function write_chunk(
    output_dir, chunk_index, model, cosmo,
    pos, mass, freqs)
    # Allocate some arrays and fill them up for centrals and satellites
    @time sources = generate_sources(model, cosmo, pos, mass);

    # Deposit the sources into maps
    fluxes_cen = Array{Float32, 1}(undef, sources.N_cen)
    fluxes_sat = Array{Float32, 1}(undef, sources.N_sat)
    m = Map{Float64,RingOrder}(model.nside)

    # loop over all frequencies and paint sources to appropriate freq map
    @time begin
        for freq in freqs

            XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, model, sources,
                fluxes_cen, fluxes_sat)

            filename = "$(output_dir)/cib_$(freq).fits"

            if chunk_index > 1
                m0 = Healpix.readMapFromFITS(filename, 1, Float32)
                m.pixels = m.pixels + m0.pixels
            end
            Healpix.saveToFITS(m, "!$(filename)", typechar="D")
        end
    end

end

##
function run_all_chunks(output_dir, halo_pos, halo_mass, freqs; N_chunks=2)
    # provide views into halo positions and masses for chunks of the halos
    N_halos = size(halo_mass, 1)
    chunksize = trunc(Integer, N_halos / N_chunks + 1)
    chunks = chunk(N_halos, chunksize)
    for chunk_index in 1:length(chunks)
    left_ind, right_ind = chunks[chunk_index]
        println("Chunk ", chunk_index, "/", length(chunks),
            " ", left_ind, " ", right_ind)
        pos = @view halo_pos[:, left_ind:right_ind]
        mass = @view halo_mass[left_ind:right_ind]
        write_chunk(output_dir, chunk_index, model, cosmo,
            pos, mass, freqs)
    end
end
## compute on all chunks, on all halos

# freqs = [
#     "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4",
#     "63.9", "67.8", "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143",
#     "153", "164", "189", "210", "217", "232", "256", "275", "294", "306", "314",
#     "340", "353", "375", "409", "467", "525", "545", "584", "643", "729", "817",
#     "857", "906", "994", "1080"
# ]
freqs = ["143"]

scratch_dir = ENV["SCRATCH"]
println("SCRATCH: ", scratch_dir)
run_all_chunks(scratch_dir, halo_pos, halo_mass, freqs)

##