text:
abstract type Basis{T} <: LazyQuasiMatrix{T} end
abstract type Weight{T} <: LazyQuasiVector{T} end

struct WeightLayout <: AbstractQuasiLazyLayout end
abstract type AbstractBasisLayout <: AbstractQuasiLazyLayout end
abstract type AbstractWeightedBasisLayout <: AbstractBasisLayout end
struct BasisLayout <: AbstractBasisLayout end
struct SubBasisLayout <: AbstractBasisLayout end
struct MappedBasisLayout <: AbstractBasisLayout end
struct WeightedBasisLayout{Basis} <: AbstractWeightedBasisLayout end
const SubWeightedBasisLayout = WeightedBasisLayout{SubBasisLayout}
const MappedWeightedBasisLayout = WeightedBasisLayout{MappedBasisLayout}

SubBasisLayouts = Union{SubBasisLayout,SubWeightedBasisLayout}
WeightedBasisLayouts = Union{WeightedBasisLayout,SubWeightedBasisLayout,MappedWeightedBasisLayout}
MappedBasisLayouts = Union{MappedBasisLayout,MappedWeightedBasisLayout}

struct AdjointBasisLayout{Basis} <: AbstractQuasiLazyLayout end
const AdjointSubBasisLayout = AdjointBasisLayout{SubBasisLayout}
const AdjointMappedBasisLayout = AdjointBasisLayout{MappedBasisLayout}

MemoryLayout(::Type{<:Basis}) = BasisLayout()
MemoryLayout(::Type{<:Weight}) = WeightLayout()

adjointlayout(::Type, ::Basis) where Basis<:AbstractBasisLayout = AdjointBasisLayout{Basis}()
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::Basis) where Basis<:AbstractBasisLayout = WeightedBasisLayout{Basis}()

# A sub of a weight is still a weight
sublayout(::WeightLayout, _) = WeightLayout()
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{Map,AbstractVector}}) = MappedBasisLayout()

## Weighted basis interface
unweighted(P::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{AbstractQuasiVector,AbstractQuasiMatrix}}) = last(P.args)
unweighted(V::SubQuasiArray) = view(unweighted(parent(V)), parentindices(V)...)
weight(P::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{AbstractQuasiVector,AbstractQuasiMatrix}}) = first(P.args)
weight(V::SubQuasiArray) = weight(parent(V))[parentindices(V)[1]]

unweighted(a::AbstractQuasiArray) = unweighted(MemoryLayout(a), a)
# Default is lazy
ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
pinv(J::Basis) = apply(pinv,J)

function _equals(::AbstractBasisLayout, ::AbstractBasisLayout, A, B)
    axes(A) == axes(B) && throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
    false
end

_equals(::SubBasisLayouts, ::SubBasisLayouts, A::SubQuasiArray, B::SubQuasiArray) = parentindices(A) == parentindices(B) && parent(A) == parent(B)
_equals(::MappedBasisLayouts, ::MappedBasisLayouts, A::SubQuasiArray, B::SubQuasiArray) = parentindices(A) == parentindices(B) && demap(A) == demap(B)
_equals(::AbstractWeightedBasisLayout, ::AbstractWeightedBasisLayout, A, B) = weight(A) == weight(B) && unweighted(A) == unweighted(B)

@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)
for op in (:+, :-)
    @eval @inline copy(L::Ldiv{Lay,BroadcastLayout{typeof($op)},<:Any,<:AbstractQuasiVector}) where Lay<:MappedBasisLayouts =
        copy(Ldiv{Lay,LazyLayout}(L.A,L.B))
end

@inline function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)}})
    a,b = arguments(L.B)
    (L.A\a)-(L.A\b)
end

@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)},<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

@inline function copy(P::Ldiv{<:AbstractBasisLayout,<:AbstractBasisLayout})
    A, B = P.A, P.B
    A == B || throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
    SquareEye{eltype(eltype(P))}((axes(A,2),)) # use double eltype for array-valued
end
@inline function copy(P::Ldiv{<:SubBasisLayouts,<:SubBasisLayouts})
    A, B = P.A, P.B
    parent(A) == parent(B) ||
        throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
    Eye{eltype(eltype(P))}((axes(A,2),axes(B,2)))
end

@inline function copy(P::Ldiv{<:MappedBasisLayouts,<:MappedBasisLayouts})
    A, B = P.A, P.B
    demap(A)\demap(B)
end

function transform_ldiv_if_columns(P::Ldiv{<:MappedBasisLayouts,<:Any,<:Any,<:AbstractQuasiVector}, ::OneTo)
    A,B = P.A, P.B
    demap(A) \ B[invmap(basismap(A))]
end

function transform_ldiv_if_columns(P::Ldiv{<:MappedBasisLayouts,<:Any,<:Any,<:AbstractQuasiMatrix}, ::OneTo)
    A,B = P.A, P.B
    demap(A) \ B[invmap(basismap(A)),:]
end

copy(L::Ldiv{<:MappedBasisLayouts,ApplyLayout{typeof(*)}}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{<:MappedBasisLayouts,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = transform_ldiv(L.A, L.B)

@inline copy(L::Ldiv{<:AbstractBasisLayout,<:SubBasisLayouts}) = apply(\, L.A, ApplyQuasiArray(L.B))
@inline function copy(L::Ldiv{<:SubBasisLayouts,<:AbstractBasisLayout})
    P = parent(L.A)
    kr, jr = parentindices(L.A)
    layout_getindex(apply(\, P, L.B), jr, :) # avoid sparse arrays
end

# default to transform for expanding weights
copy(L::Ldiv{<:AbstractBasisLayout,WeightLayout}) = transform_ldiv(L.A, L.B)

# multiplication operators, reexpand in basis A
@inline function _broadcast_mul_ldiv(::Tuple{Any,AbstractBasisLayout}, A, B)
    a,b = arguments(B)
    @assert a isa AbstractQuasiVector # Only works for vec .* mat
    ab = (A * (A \ a)) .* b # broadcasted should be overloaded
    MemoryLayout(ab) isa BroadcastLayout && return transform_ldiv(A, ab)
    A \ ab
end

@inline function _broadcast_mul_ldiv(::Tuple{Any,ApplyLayout{typeof(*)}}, A, B)
    a,b = arguments(B)
    @assert a isa AbstractQuasiVector # Only works for vec .* mat
    args = arguments(*, b)
    *(A \ (a .* first(args)), tail(args)...)
end

function _broadcast_mul_ldiv(::Tuple{ScalarLayout,Any}, A, B)
    a,b = arguments(B)
    a * (A \ b)
end

function _broadcast_mul_ldiv(::Tuple{ScalarLayout,ApplyLayout{typeof(*)}}, A, B)
    a,b = arguments(B)
    a * (A \ b)
end

_broadcast_mul_ldiv(::Tuple{ScalarLayout,AbstractBasisLayout}, A, B) =
    _broadcast_mul_ldiv((ScalarLayout(),UnknownLayout()), A, B)
_broadcast_mul_ldiv(_, A, B) = copy(Ldiv{typeof(MemoryLayout(A)),UnknownLayout}(A,B))

copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)}}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)
copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)

# ambiguity
copy(L::Ldiv{<:MappedBasisLayouts,BroadcastLayout{typeof(*)}}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)
copy(L::Ldiv{<:MappedBasisLayouts,BroadcastLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)

# expansion
_grid(_, P) = error("Overload Grid")

_grid(::MappedBasisLayout, P) = invmap(parentindices(P)[1])[grid(demap(P))]
_grid(::SubBasisLayout, P) = grid(parent(P))
_grid(::WeightedBasisLayouts, P) = grid(unweighted(P))
grid(P) = _grid(MemoryLayout(P), P)

struct TransformFactorization{T,Grid,Plan,IPlan} <: Factorization{T}
    grid::Grid
    plan::Plan
    iplan::IPlan
end

TransformFactorization{T}(grid, plan) where T = TransformFactorization{T,typeof(grid),typeof(plan),Nothing}(grid, plan, nothing)

"""
    TransformFactorization(grid, plan)

associates a planned transform with a grid. That is, if `F` is a `TransformFactorization`, then
`F \\ f` is equivalent to `F.plan * f[F.grid]`.
"""
TransformFactorization(grid, plan) = TransformFactorization{promote_type(eltype(eltype(grid)),eltype(plan))}(grid, plan)

TransformFactorization{T}(grid, ::Nothing, iplan) where T = TransformFactorization{T,typeof(grid),Nothing,typeof(iplan)}(grid, nothing, iplan)

"""
    TransformFactorization(grid, nothing, iplan)

associates a planned inverse transform with a grid. That is, if `F` is a `TransformFactorization`, then
`F \\ f` is equivalent to `F.iplan \\ f[F.grid]`.
"""
TransformFactorization(grid, ::Nothing, iplan) = TransformFactorization{promote_type(eltype(eltype(grid)),eltype(iplan))}(grid, nothing, iplan)

grid(T::TransformFactorization) = T.grid
function size(T::TransformFactorization, k)
    @assert k == 2 # TODO: make consistent
    size(T.plan,1)
end

function size(T::TransformFactorization{<:Any,<:Any,Nothing}, k)
    @assert k == 2 # TODO: make consistent
    size(T.iplan,2)
end

\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiVector{T}) where T = a.iplan \  convert(Array{T}, b[a.grid])
\(a::TransformFactorization, b::AbstractQuasiVector) = a.plan * convert(Array, b[a.grid])
\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractVector) = a.iplan \  b
\(a::TransformFactorization, b::AbstractVector) = a.plan * b
ldiv!(ret::AbstractVecOrMat, a::TransformFactorization, b::AbstractVecOrMat) = mul!(ret, a.plan, b)
ldiv!(ret::AbstractVecOrMat, a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractVecOrMat) = ldiv!(ret, a.iplan, b)

\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiMatrix{T}) where T = a \  convert(Array{T}, b[a.grid,:])
\(a::TransformFactorization, b::AbstractQuasiMatrix) = a \ convert(Array, b[a.grid,:])
\(a::TransformFactorization, b::AbstractMatrix) = ldiv!(Array{promote_type(eltype(a),eltype(b))}(undef,size(a,2),size(b,2)), a, b)

function _factorize(::AbstractBasisLayout, L, dims...; kws...)
    p = grid(L)
    TransformFactorization(p, nothing, factorize(L[p,:]))
end

"""
    ProjectionFactorization(F, inds)

projects a factorization to a subset of coefficients. That is, if `P` is a `ProjectionFactorization`
then `P \\ f` is equivalent to `(F \\ f)[inds]`
"""
struct ProjectionFactorization{T, FAC<:Factorization{T}, INDS} <: Factorization{T}
    F::FAC
    inds::INDS
end

\(a::ProjectionFactorization, b::AbstractQuasiVector) = (a.F \ b)[a.inds]
\(a::ProjectionFactorization, b::AbstractQuasiMatrix) = (a.F \ b)[a.inds,:]
\(a::ProjectionFactorization, b::AbstractVector) = (a.F \ b)[a.inds]

_factorize(::SubBasisLayout, L, dims...; kws...) = ProjectionFactorization(factorize(parent(L), dims...; kws...), parentindices(L)[2])

"""
    MappedFactorization(F, map)

remaps a factorization to a different domain. That is, if `M` is a `MappedFactorization`
then `M \\ f` is equivalent to `F \\ f[map]`
"""
struct MappedFactorization{T, FAC<:Factorization{T}, MAP} <: Factorization{T}
    F::FAC
    map::MAP
end

\(a::MappedFactorization, b::AbstractQuasiVector) = a.F \ view(b, a.map)
\(a::MappedFactorization, b::AbstractVector) = a.F \ b
\(a::MappedFactorization, b::AbstractQuasiMatrix) = a.F \ view(b, a.map, :)

function invmap end

function _factorize(::MappedBasisLayout, L, dims...; kws...)
    kr, jr = parentindices(L)
    P = parent(L)
    MappedFactorization(factorize(view(P,:,jr), dims...; kws...), invmap(parentindices(L)[1]))
end

plan_ldiv(A, B::AbstractQuasiVector) = factorize(A)
plan_ldiv(A, B::AbstractQuasiMatrix) = factorize(A, size(B,2))

transform_ldiv(A::AbstractQuasiArray{T}, B::AbstractQuasiArray{V}, _) where {T,V} = plan_ldiv(A, B) \ B
transform_ldiv(A, B) = transform_ldiv(A, B, size(A))

copy(L::Ldiv{<:AbstractBasisLayout}) = transform_ldiv(L.A, L.B)
# TODO: redesign to use simplifiable(\, A, B)
copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = transform_ldiv(L.A, L.B)
copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
# A BroadcastLayout of unknown function is only knowable pointwise
transform_ldiv_if_columns(L, _) = ApplyQuasiArray(\, L.A, L.B)
transform_ldiv_if_columns(L, ::OneTo) = transform_ldiv(L.A,L.B)
transform_ldiv_if_columns(L) = transform_ldiv_if_columns(L, axes(L.B,2))
copy(L::Ldiv{<:AbstractBasisLayout,<:BroadcastLayout}) = transform_ldiv_if_columns(L)
# Inclusion are QuasiArrayLayout
copy(L::Ldiv{<:AbstractBasisLayout,QuasiArrayLayout}) = transform_ldiv(L.A, L.B)
# Otherwise keep lazy to support, e.g., U\D*T
copy(L::Ldiv{<:AbstractBasisLayout,<:AbstractLazyLayout}) = transform_ldiv_if_columns(L)
copy(L::Ldiv{<:AbstractBasisLayout,ZerosLayout}) = Zeros{eltype(L)}(axes(L)...)

transform_ldiv_if_columns(L::Ldiv{<:Any,<:ApplyLayout{typeof(hcat)}}, ::OneTo) = transform_ldiv(L.A, L.B)
transform_ldiv_if_columns(L::Ldiv{<:Any,<:ApplyLayout{typeof(hcat)}}, _) = hcat((Ref(L.A) .\ arguments(hcat, L.B))...)

"""
    WeightedFactorization(w, F)

weights a factorization `F` by `w`.
"""
struct WeightedFactorization{T, WW, FAC<:Factorization{T}} <: Factorization{T}
    w::WW
    F::FAC
end

_factorize(::WeightedBasisLayouts, wS, dims...; kws...) = WeightedFactorization(weight(wS), factorize(unweighted(wS), dims...; kws...))

\(F::WeightedFactorization, b::AbstractQuasiVector) = F.F \ (b ./ F.w)

##
# Algebra
##

struct ExpansionLayout{Lay} <: AbstractLazyLayout end
applylayout(::Type{typeof(*)}, ::Lay, ::Union{PaddedLayout,AbstractStridedLayout}) where Lay <: AbstractBasisLayout = ExpansionLayout{Lay}()

basis(v::ApplyQuasiArray{<:Any,N,typeof(*)}) where N = v.args[1]
coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any}}) where N = v.args[2]

function unweighted(lay::ExpansionLayout, a)
    wP,c = arguments(lay, a)
    unweighted(wP) * c
end

LazyArrays._mul_arguments(::ExpansionLayout, A) = LazyArrays._mul_arguments(ApplyLayout{typeof(*)}(), A)
copy(L::Ldiv{Bas,<:ExpansionLayout}) where Bas<:AbstractBasisLayout = copy(Ldiv{Bas,ApplyLayout{typeof(*)}}(L.A, L.B))
copy(L::Mul{<:ExpansionLayout,Lay}) where Lay = copy(Mul{ApplyLayout{typeof(*)},Lay}(L.A, L.B))

function _broadcastbasis(::typeof(+), _, _, a, b)
    try
        a ≠ b && error("Overload broadcastbasis(::typeof(+), ::$(typeof(a)), ::$(typeof(b)))")
    catch
        error("Overload broadcastbasis(::typeof(+), ::$(typeof(a)), ::$(typeof(b)))")
    end
    a
end

_broadcastbasis(::typeof(+), ::MappedBasisLayouts, ::MappedBasisLayouts, a, b) = broadcastbasis(+, demap(a), demap(b))[basismap(a), :]
function _broadcastbasis(::typeof(+), ::SubBasisLayout, ::SubBasisLayout, a, b)
    kr_a,jr_a = parentindices(a)
    kr_b,jr_b = parentindices(b)
    @assert kr_a == kr_b # frist axes must match
    view(broadcastbasis(+, parent(a), parent(b)), kr_a, union(jr_a,jr_b))
end
_broadcastbasis(::typeof(+), ::SubBasisLayout, _, a, b) = broadcastbasis(+, parent(a), b)
_broadcastbasis(::typeof(+), _, ::SubBasisLayout, a, b) = broadcastbasis(+, a, parent(b))

broadcastbasis(::typeof(+), a, b) = _broadcastbasis(+, MemoryLayout(a), MemoryLayout(b), a, b)
broadcastbasis(::typeof(+), a, b, c...) = broadcastbasis(+, broadcastbasis(+, a, b), c...)

broadcastbasis(::typeof(-), a, b) = broadcastbasis(+, a, b)

@eval function layout_broadcasted(::NTuple{2,ExpansionLayout}, ::typeof(-), f, g)
    S,c = arguments(f)
    T,d = arguments(g)
    ST = broadcastbasis(-, S, T)
    ST * ((ST \ S) * c - (ST \ T) * d)
end

_plus_P_ldiv_Ps_cs(P, ::Tuple{}, ::Tuple{}) = ()
_plus_P_ldiv_Ps_cs(P, Q::Tuple, cs::Tuple) = tuple((P \ first(Q)) * first(cs), _plus_P_ldiv_Ps_cs(P, tail(Q), tail(cs))...)
function layout_broadcasted(::Tuple{Vararg{ExpansionLayout}}, ::typeof(+), fs...)
    Ps = first.(arguments.(fs))
    cs = last.(arguments.(fs))
    P = broadcastbasis(+, Ps...)
    P * +(_plus_P_ldiv_Ps_cs(P, Ps, cs)...)  # +((Ref(P) .\ Ps .* cs)...)
end

function layout_broadcasted(::NTuple{2,ExpansionLayout}, ::typeof(*), a, f)
    axes(a,1) == axes(f,1) || throw(DimensionMismatch())
    P,c = arguments(f)
    (a .* P) * c
end

function layout_broadcasted(::Tuple{ExpansionLayout{<:AbstractWeightedBasisLayout},AbstractBasisLayout}, ::typeof(*), a, P)
    axes(a,1) == axes(P,1) || throw(DimensionMismatch())
    wQ,c = arguments(a)
    w,Q = weight(wQ),unweighted(wQ)
    uP = unweighted(a) .* P
    w .* uP
end

# function layout_broadcasted(::Tuple{PolynomialLayout,MappedBasisLayouts}, ::typeof(*), x, P)
#     axes(a,1) == axes(P,1) || throw(DimensionMismatch())
#     wQ,c = arguments(a)
#     w,Q = arguments(wQ)
#     uP = unweighted(a) .* P
#     w .* uP
# end

function _equals(::ExpansionLayout, ::ExpansionLayout, f, g)
    S,c = arguments(f)
    T,d = arguments(g)
    ST = broadcastbasis(+, S, T)
    (ST \ S) * c == (ST \ T) * d
end

function QuasiArrays._mul_summary(::ExpansionLayout, io::IO, f)
    P,c = arguments(f)
    show(io, P)
    print(io, " * ")
    show(io, c)
end

##
# Multiplication for mapped and subviews x .* view(P,...)
##

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    T = promote_type(eltype(x), eltype(C))
    x == axes(C,1) || throw(DimensionMismatch())
    P = parent(C)
    kr,jr = parentindices(C)
    y = axes(P,1)
    Y = P \ (y .* P)
    X = kr.A \ (Y     - kr.b * I)
    P[kr, :] * view(X,:,jr)
end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
    T = promote_type(eltype(x), eltype(C))
    x == axes(C,1) || throw(DimensionMismatch())
    P = parent(C)
    kr,_ = parentindices(C)
    y = axes(P,1)
    Y = P \ (y .* P)
    X = kr.A \ (Y     - kr.b * I)
    P[kr, :] * X
end

# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::AbstractQuasiVector, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
#     T = promote_type(eltype(f), eltype(C))
#     axes(f,1) == axes(C,1) || throw(DimensionMismatch())
#     P = parent(C)
#     kr,jr = parentindices(C)
#     (f[invmap(kr)] .* P)[kr,jr]
# end

# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::AbstractQuasiVector, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
#     T = promote_type(eltype(f), eltype(C))
#     axes(f,1) == axes(C,1) || throw(DimensionMismatch())
#     P = parent(C)
#     kr,jr = parentindices(C)
#     (f[invmap(kr)] .* P)[kr,jr]
# end

broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::Broadcasted, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) =
    broadcast(*, materialize(f), C)

function layout_broadcasted(::Tuple{PolynomialLayout,WeightedBasisLayout}, ::typeof(*), x, C)
    x == axes(C,1) || throw(DimensionMismatch())
    weight(C) .* (x .* unweighted(C))
end

## materialize views

# materialize(S::SubQuasiArray{<:Any,2,<:ApplyQuasiArray{<:Any,2,typeof(*),<:Tuple{<:Basis,<:Any}}}) =
#     *(arguments(S)...)

# mass matrix
# y = p(x), dy = p'(x) * dx
# \int_a^b f(y) g(y) dy = \int_{-1}^1 f(p(x))*g(p(x)) * p'(x) dx

_sub_getindex(A, kr, jr) = A[kr, jr]
_sub_getindex(A, ::Slice, ::Slice) = A

@simplify function *(Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}},
             B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    A = Ac'
    PA,PB = parent(A),parent(B)
    kr,jr = parentindices(B)
    _sub_getindex((PA'PB)/kr.A,parentindices(A)[2],jr)
end

# Differentiation of sub-arrays

# avoid stack overflow from unmaterialize Derivative() * parent()
_der_sub(DP, inds...) = DP[inds...]
_der_sub(DP::ApplyQuasiMatrix{T,typeof(*),<:Tuple{Derivative,Any}}, kr, jr) where T = ApplyQuasiMatrix{T}(*, DP.args[1], view(DP.args[2], kr, jr))

# need to customise simplifiable so can't use @simplify
simplifiable(::typeof(*), A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}})= simplifiable(*, Derivative(axes(parent(B),1)), parent(B))
simplifiable(::typeof(*), Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = simplifiable(*, Bc', Ac')
function mul(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}})
    axes(A,2) == axes(B,1) || throw(DimensionMismatch())
    P = parent(B)
    _der_sub(Derivative(axes(P,1))*P, parentindices(B)...)
end
mul(Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = mul(Bc', Ac')'

simplifiable(::typeof(*), A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) = simplifiable(*, Derivative(axes(parent(B),1)), parent(B))
simplifiable(::typeof(*), Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = simplifiable(*, Bc', Ac')
function mul(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    axes(A,2) == axes(B,1) || throw(DimensionMismatch())
    P = parent(B)
    kr,jr = parentindices(B)
    (Derivative(axes(P,1))*P*kr.A)[kr,jr]
end
mul(Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = mul(Bc', Ac')'

# we represent as a Mul with a banded matrix
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractVector}}) = SubBasisLayout()
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractVector}}) = MappedBasisLayout()
sublayout(::WeightedBasisLayouts, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractVector}}) = MappedWeightedBasisLayout()
sublayout(::WeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractVector}}) = SubWeightedBasisLayout()
sublayout(::MappedWeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractVector}}) = MappedWeightedBasisLayout()

@inline sub_materialize(::AbstractBasisLayout, V::AbstractQuasiArray) = V
@inline sub_materialize(::AbstractBasisLayout, V::AbstractArray) = V

demap(x) = x
demap(x::BroadcastQuasiArray) = BroadcastQuasiArray(x.f, map(demap, arguments(x))...)
demap(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Any,Slice}}) = parent(V)
demap(V::SubQuasiArray{<:Any,1}) = parent(V)
function demap(V::SubQuasiArray{<:Any,2})
    kr, jr = parentindices(V)
    demap(parent(V)[kr,:])[:,jr]
end
function demap(wB::ApplyQuasiArray{<:Any,N,typeof(*)}) where N
    a = arguments(wB)
    *(demap(first(a)), tail(a)...)
end

basismap(x::SubQuasiArray) = parentindices(x)[1]
basismap(x::BroadcastQuasiArray) = basismap(x.args[1])

##
# SubLayout behaves like ApplyLayout{typeof(*)}

combine_mul_styles(::SubBasisLayouts) = combine_mul_styles(ApplyLayout{typeof(*)}())
_mul_arguments(::SubBasisLayouts, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
arguments(::SubBasisLayouts, A) = arguments(ApplyLayout{typeof(*)}(), A)
call(::SubBasisLayouts, ::SubQuasiArray) = *

combine_mul_styles(::AdjointSubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
_mul_arguments(::AdjointSubBasisLayout, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
arguments(::AdjointSubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
call(::AdjointSubBasisLayout, ::SubQuasiArray) = *

copy(M::Mul{AdjointSubBasisLayout,<:SubBasisLayouts}) = apply(*, arguments(M.A)..., arguments(M.B)...)

function arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
    A = parent(V)
    _,jr = parentindices(V)
    first(jr) ≥ 1 || throw(BoundsError())
    P = _BandedMatrix(Ones{Int}((Base.OneTo(1),axes(jr,1))), axes(A,2), first(jr)-1,1-first(jr))
    A,P
end

function arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractVector}})
    A = parent(V)
    _,jr = parentindices(V)
    first(jr) ≥ 1 || throw(BoundsError())
    P = Eye{Int}((axes(A,2),))[:,jr]
    A,P
end

####
# sum
####

function __sum(::SubBasisLayout, Vm, dims)
    @assert dims == 1
    sum(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
end

__sum(::AdjointBasisLayout, Vm::AbstractQuasiMatrix, dims) = permutedims(sum(Vm'; dims=(isone(dims) ? 2 : 1)))

function __sum(::MappedBasisLayouts, V, dims)
    kr = basismap(V)
    @assert kr isa AbstractAffineQuasiVector
    sum(demap(V); dims=dims)/kr.A
end

__sum(::ExpansionLayout, A, dims) = __sum(ApplyLayout{typeof(*)}(), A, dims)
__cumsum(::ExpansionLayout, A, dims) = __cumsum(ApplyLayout{typeof(*)}(), A, dims)

include("splines.jl")


text:
using DFControl, Test

Client.rm_pseudoset!("test", "localhost_test")

@test_throws ErrorException Client.pseudos("localhost_test", "test", [:Ni])

ext: jl


text:
include("OpCodes.jl")
using .OpCodes

function run1(data)
    function move!(ic, nsteps, pos, origin)
        pos ∈ positions && return
        push!(positions, pos)
        dirs = setdiff(1:4, origin)
        for dir in dirs
            put!(inp, dir)
            res = take!(outp)
            field[(pos .+ directions[dir])] = res
            if res > 0
                res == 2 && (shortest_path_length = nsteps + 1)
                orig = back[dir]
                move!(ic, nsteps + 1, pos .+ directions[dir], orig)
                put!(inp, orig)
                throwaway = take!(outp)
            end
        end
    end
    inp, outp = Channel{Int}(1), Channel{Int}(1)
    ic = IntComputer(data, inp, outp)
    nsteps = 0
    positions = Set{Tuple{Int, Int}}()
    directions = ((1,0), (-1,0), (0,-1), (0,1))
    back = (2, 1, 4, 3)
    shortest_path_length = 99999
    field = Dict((0,0)=>1)
    @async compute!(ic)
    move!(ic, 0, (0,0), 5)
    shortest_path_length, field
end
result1, field = run1(readline("data/Dec15_data.txt"))

function run2(field)
    neighbors(pos) = (pos .+ a for a in ((1,0), (-1,0), (0,1), (0,-1)))
    current = Set([k for (k,v) in field if v == 2])
    space = Set([k for (k,v) in field if v > 0])
    pop!(space, first(current))
    counter = 0
    while true
        counter += 1
        next = Set([v for v in Iterators.flatten(neighbors.(current)) if v ∈ space])
        @show counter, length(space), length(next), next
        [pop!(space, n) for n in next]
        length(space) == 0 && return counter
        current = next
    end
end
run2(field)


text:
using Pkg

function addandpin(spec)
    Pkg.add(PackageSpec(; spec...))
    Pkg.pin(spec.name)
end

# packages without external dependencies
# that have to be manually instlled

pkg1 = [(name="StatsBase", version="0.31.0"),
        (name="TimeZones", version="0.9.2"),
        (name="BSON", version="0.2.3"),
        (name="Revise", version="2.1.8"),
        (name="Distributions", version="0.21.0"),
        (name="Clp", version="0.6.2"),
        (name="HTTP", version="0.8.5"),
        (name="Gumbo", version="0.5.1"),
        (name="StringEncodings", version="0.3.1"),
        (name="ZMQ", version="1.0.0"),
        (name="CodecZlib", version="0.6.0"),
        (name="JSON", version="0.21.0"),
        (name="BenchmarkTools", version="0.4.2"),
        (name="JuliaWebAPI", version="0.6.0"),
        (name="FileIO", version="1.0.7"),
        (name="ProfileView", version="0.4.1"),
        (name="StaticArrays", version="0.11.0"),
        (name="ForwardDiff", version="0.10.3"),
        (name="Optim", version="0.19.1"),
        (name="JuMP", version="0.19.2"),
        (name="JLD2", version="0.1.2"),
        (name="XLSX", version="0.5.4"),
        (name="Cbc", version="0.6.2"),
        (name="DataFrames", version="0.19.2"),
        (name="CSV", version="0.5.9"),
        (name="DataFramesMeta", version="0.5.0"),
        (name="Feather", version="0.5.3"),
        (name="FreqTables", version="0.3.1"),
        (name="OnlineStats", version="0.27.0"),
        (name="MySQL", version="0.7.0"),
        (name="Cascadia", version="0.4.0"),
        (name="UnicodePlots", version="1.1.0"),
        (name="ParallelDataTransfer", version="0.5.0")]

foreach(addandpin, pkg1)

# packages that depend on Python Anaconda
# whose manual installation is optional

pkg2 = [(name="Conda", version="1.3.0"),
        (name="PyCall", version="1.91.2"),
        (name="PyPlot", version="2.8.1"),
        (name="Plots", version="0.26.0"),
        (name="StatsPlots", version="0.10.2")]

foreach(addandpin, pkg2)

# packages that have external dependencies
# that have to be manually installed

pkg3 = [(name="RCall", version="0.13.3"),
        (name="JDBC", version="0.4.1"),
        (name="LibPQ", version="0.11.0"),
        (name="Gurobi", version="0.6.0")]

foreach(addandpin, pkg3)


text:
module QPDAS

global const DEBUG = false

export QuadraticProgram, BoxConstrainedQP, solve!, update!, ldiv2!

using LinearAlgebra, SparseArrays

abstract type AbstractCholeskySpecial{T,MT} <: Factorization{T} end

# Special type that allowes for solving M\b with some rows/columns "deleted"
include("choleskySpecial.jl")
include("choleskySpecialShifted.jl")
#include("LDLTSpecial.jl")
# The type for representing the dual (semidefinite) problem
include("boxConstrainedQP.jl")
# Teh main type for solving positive definite quadratic programs
include("quadraticProgram.jl")

end # module

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
# Token types
abstract type Token end

mutable struct TextToken <: Token
_type::String
value::String
end

mutable struct TagToken <: Token
_type::String
value::String
ltag::String
rtag::String
indent::String
TagToken(_type, value, ltag, rtag, indent="") = new(_type, value, ltag, rtag, indent)
end

mutable struct SectionToken <: Token
_type::String
value::String
ltag::String
rtag::String
collector::Vector
SectionToken(_type, value, ltag, rtag) = new(_type, value, ltag, rtag, Any[])
end

mutable struct MustacheTokens
tokens::Vector{Token}
end
MustacheTokens() = MustacheTokens(Token[])

Base.length(tokens::MustacheTokens) = length(tokens.tokens)
Base.lastindex(tokens::MustacheTokens) = lastindex(tokens.tokens)
Base.getindex(tokens::MustacheTokens, ind) = getindex(tokens.tokens, ind)
Base.pop!(tokens::MustacheTokens) = pop!(tokens.tokens)

function Base.push!(tokens::MustacheTokens, token::Token)
    # squash if possible
    if length(tokens) == 0
        push!(tokens.tokens, token)
        return
    end

    lastToken = tokens[end]
    if !falsy(lastToken) && lastToken._type == "text" ==  token._type
        lastToken.value *= token.value
    else
        push!(tokens.tokens, token)
    end
end

mutable struct AnIndex
    ind::Int
    value::String
end
Base.string(ind::AnIndex) = string(ind.value)

### We copy some functions from Tokenize
function peekchar(io::Base.GenericIOBuffer)
    if !io.readable || io.ptr > io.size
        return '∅'
    end
    ch, _ = readutf(io)
    return ch
end

# this implementation is copied from Tokenize,
@inline function utf8_trailing(i)
    if i < 193
        return 0
    elseif i < 225
        return 1
    elseif i < 241
        return 2
    elseif i < 249
        return 3
    elseif i < 253
        return 4
    else
        return 5
    end
end

const utf8_offset = [0x00000000
                    0x00003080
                    0x000e2080
                    0x03c82080
                    0xfa082080
                    0x82082080]

function readutf(io, offset = 0)
    ch = convert(UInt8, io.data[io.ptr + offset])
    if ch < 0x80
        return convert(Char, ch), 0
    end
    trailing = utf8_trailing(ch + 1)
    c::UInt32 = 0
    for j = 1:trailing
        c += ch
        c <<= 6
        ch = convert(UInt8, io.data[io.ptr + j + offset])
    end
    c += ch
    c -= utf8_offset[trailing + 1]
    return convert(Char, c), trailing
end

####

## this is modified from dpeekchar in Tokenize
function peekaheadmatch(io::IOBuffer, m=['{','{'])
    if !io.readable || io.ptr > io.size
        return false
    end
    i,N = 1, length(m)
    ch1, trailing = readutf(io)
    ch1 != m[1] && return false
    offset = 0
    while true
        i == N && return true
        i += 1

        offset += trailing + 1
        if io.ptr + offset > io.size
            return false
        end
        chn, trailing = readutf(io, offset)
        chn != m[i] && return false
    end
    return false
end

# ...{{.... -> (..., {{....)
function scan_until!(io, tags, firstline)
    out = IOBuffer()
    while !end_of_road(io)
        if peekaheadmatch(io, tags)
            val = String(take!(out))
            close(out)
            return val, firstline
        else
            c = read(io, Char)
            firstline = firstline && !(c== '\n')
            print(out, c)
        end
    end
    val = String(take!(out))
    close(out)
    return val, firstline
end

# skip over tags
function scan_past!(io, tags = ('{','{'))
    for i in 1:length(tags)
        t = read(io, Char)
    end
end

# end of io stream
function end_of_road(io)
    !io.readable && return true
    io.ptr > io.size && return true
    return false
end

# stip whitespace at head of io stream
WHITESPACE = (' ', '\t', '\r')
function popfirst!_whitespace(io)

    while !end_of_road(io)
        c::Char = peekchar(io)
        !(c in WHITESPACE) && break
        read(io, Char)
    end

end

# is tag possibly standalone
# check the left side
# XXX This could be part of scan_until!
# then we could avoid regular expressions here
function is_l_standalone(txt, firstline)
    if firstline
        occursin(r"^ *$", txt)
    else
        out = occursin(r"\n *$", txt)
        out
    end
end

# check the last side
function is_r_standalone(io)
    end_of_road(io) && return true
    offset = 0
    ch::Char, trailing::Int = readutf(io)
    while true
        ch == '\n' && return true
        !(ch in WHITESPACE) && return false
        offset += trailing + 1
        io.ptr + offset > io.size && return true
        ch, trailing = readutf(io, offset)
    end
    return false
end

## Make the intial set of tokens before nesting
function make_tokens(template, tags)

    ltag, rtag = tags
    ltags = collect(ltag)
    rtags = collect(rtag)

    sections = MustacheTokens()
    tokens = MustacheTokens()

    firstline = true
    io = IOBuffer(template)
    while !end_of_road(io)
        # we have
        # ...text... {{ tag }} ...rest...
        # each lap makes token of ...text... and {{tag}} leaving ...rest...

        b0 = 0
        text_value, firstline = scan_until!(io, ltags, firstline)
        b1 = position(io)

        if end_of_road(io)
            token = TextToken("text", text_value)
            push!(tokens, token)
            return tokens
        end
        scan_past!(io, ltags)
        text_token = TextToken("text", text_value)

        # grab tag token
        t0 = position(io)
        token_value, firstline = scan_until!(io, rtags, firstline)
        t1 = position(io)

        if end_of_road(io)
            throw(ArgumentError("tag not closed: $token_value $ltag $rtag"))
        end
        scan_past!(io, rtags)

        # what kinda tag did we get?
        token_value = stripWhitespace(token_value)
        _type = token_value[1:1]
        if _type in  ("#", "^", "/", ">", "<", "!", "|", "=", "{", "&")
            # type is first, we peel it off, also strip trailing = and },
            # as necessary

            token_value = stripWhitespace(token_value[2:(end-(_type == "="))])

            if _type == "{"
                # strip "}" if present in io
                c = peekchar(io)
                c == '}' && read(io, Char)
                # strip "}" if  present in token
                if token_value[end]=='}'
                    token_value = token_value[1:end-1]
                end
            end

            if _type == "="

                tag_token = TagToken("=", token_value, ltag, rtag)
                ts = String.(split(token_value, r"[ \t]+"))
                length(ts) != 2 && throw(ArgumentError("Change of delimiter must be a pair. Identified: $ts"))
                ltag, rtag = ts
                ltags, rtags = collect(ltag), collect(rtag)

            elseif _type == ">"

                # get indentation
                if firstline && occursin(r"^\h*$", text_value)
                    m = match(r"^(\h*)$", text_value)
                else
                    m = match(r"\n([\s\t\h]*)$", text_value)
                end
                indent = m == nothing ? "" : m.captures[1]
                tag_token = TagToken(_type, token_value, ltag, rtag, indent)

            elseif _type in ("#", "^", "|")

                tag_token = SectionToken(_type, token_value, ltag, rtag)

            else

                tag_token = TagToken(_type, token_value, ltag, rtag)

            end
        else

            token_value = stripWhitespace(token_value)
            tag_token = TagToken("name", token_value, ltag, rtag)

        end

        if tag_token._type == "name"
            firstline = false
        end

        # are we standalone?
        # yes if firstline and all whitespace
        # yes if !firstline and match \nwhitespace
        # AND
        # match rest with space[EOF,\n]

        standalone = is_r_standalone(io) && is_l_standalone(text_value, firstline)
        firstline = false

        if standalone && _type in ("!", "^", "/", "#", "<", ">", "|", "=")

            # we strip text_value back to \n or beginning
            text_token.value = replace(text_token.value, r" *$" => "")
            # strip whitspace at outset of io and "\n"
            popfirst!_whitespace(io)

            if !end_of_road(io)
                read(io, Char) # will be \n
                firstline = true
            end

        end
        push!(tokens, text_token)
        push!(tokens, tag_token)

        # account for matching/nested sections
        if _type == "#" || _type == "^" || _type ==  "|"
            push!(sections, tag_token)
        elseif _type == "/"
            ## section nestinng
            if length(sections) == 0
                throw(ArgumentError("Unopened section $token_value at $t0"))
            end

            openSection = pop!(sections)
            if openSection.value != token_value
                throw(ArgumentError("Unclosed section: " * openSection.value * " at $t0"))
            end
        end

    end

    if length(sections) > 0
        openSection = pop!(sections)
        throw(ArgumentError("Unclosed section " * string(openSection.value)))
    end

    return(tokens)
end

## Forms the given array of `tokens` into a nested tree structure where
## tokens that represent a section have two additional items: 1) an array of
## all tokens that appear in that section and 2) the index in the original
## template that represents the end of that section.
function nestTokens(tokens)
    tree = Array{Any}(undef, 0)
    collector = tree
    sections = MustacheTokens() #Array{Any}(undef, 0)

    for i in 1:length(tokens)
        token = tokens[i]
        ## a {{#name}}...{{/name}} will iterate over name
        ## a {{^name}}...{{/name}} does ... if we have no name
        ## start nesting
        if token._type == "^" || token._type == "#" || token._type == "|"
            push!(sections, token)
            push!(collector, token)
            token.collector = Array{Any}(undef, 0)
            collector = token.collector
        elseif token._type == "/"
            section = pop!(sections)
            if length(sections) > 0
                push!(sections[end].collector, token)
                collector = sections[end].collector
            else
                collector = tree
            end
        else
            push!(collector, token)
        end
    end

    return(tree)
end

## In lambdas with section this is used to go from the tokens to an unevaluated string
## This might have issues with space being trimmed
function toString(tokens)
    io = IOBuffer()
    for token in tokens
        write(io, _toString(Val{Symbol(token._type)}(), token))
    end
    out = String(take!(io))
    close(io)
    out
end

_toString(::Val{:name}, t) = t.ltag * t.value * t.rtag
_toString(::Val{:text}, t) = t.value
_toString(::Val{Symbol("^")}, t) = t.ltag * "^" * t.value * t.rtag
_toString(::Val{Symbol("|")}, t) = t.ltag * "|" * t.value * t.rtag
_toString(::Val{Symbol("/")}, t) = t.ltag * "/" * t.value * t.rtag
_toString(::Val{Symbol(">")}, t) = t.ltag * ">" * t.value * t.rtag
_toString(::Val{Symbol("<")}, t) = t.ltag * "<" * t.value * t.rtag
_toString(::Val{Symbol("&")}, t) = t.ltag * "&" * t.value * t.rtag
_toString(::Val{Symbol("{")}, t) = t.ltag * "{" * t.value * t.rtag
_toString(::Val{Symbol("=")}, t) = ""
function _toString(::Val{Symbol("#")}, t)
    out = t.ltag * "#" * t.value * t.rtag
    if !isempty(t.collector)
        out *= toString(t.collector)
    end
    out
end

# render tokens with values given in context
function renderTokensByValue(value, io, token, writer, context, template, args...)
    inverted = token._type == "^"
    if (inverted && falsy(value))
        _renderTokensByValue(value, io, token, writer, context, template, args...)
    elseif Tables.istable(value)
        if isempty(Tables.rows(value))
            return nothing
        else
            rows = Tables.rows(value)
            sch = Tables.schema(rows)
            rD = Dict()
            for row in rows
                if sch == nothing
                    _renderTokensByValue(getfield(row,1), io, token, writer, context, template, args...)
                else
                    # create a dictionary, though perhaps getfield is general enough
                    Tables.eachcolumn(sch, row) do val, col, name
                        rD[name] = val
                    end
                    renderTokens(io, token.collector, writer, ctx_push(context, rD), template, args...)
                end
            end
        end
    elseif is_dataframe(value) # XXX remove once istable(x::DataFrame) == true works
        for i in 1:size(value)[1]
            renderTokens(io, token.collector, writer, ctx_push(context, value[i,:]), template, args...)
        end
    elseif !falsy(value)
        _renderTokensByValue(value, io, token, writer, context, template, args...)
    else
        nothing
        # value == nothing
    end
end

# function renderTokensByValue(value, io, token, writer, context, template, args...)
#     if Tables.istable(value)
#         if isempty(Tables.rows(value))
#             @show "no length", token._type

#             return nothing
#         else
#             for row in Tables.rows(value)
#                 _renderTokensByValue(getfield(row,1), io, token, writer, context, template, args...)
#             end
#         end
#     elseif is_dataframe(value) # XXX remove once istable(x::DataFrame) == true works
#         for i in 1:size(value)[1]
#             renderTokens(io, token.collector, writer, ctx_push(context, value[i,:]), template, args...)
#         end
#     else
#         inverted = token._type == "^"
#         if (inverted && falsy(value)) || !falsy(value)
#             _renderTokensByValue(value, io, token, writer, context, template, args...)
#         end
#     end
# end

## Helper function for dispatch based on value in renderTokens
function _renderTokensByValue(value::AbstractDict, io, token, writer, context, template, args...)
        renderTokens(io, token.collector, writer, ctx_push(context, value), template, args...)
end

function _renderTokensByValue(value::AbstractArray, io, token, writer, context, template, args...)
   inverted = token._type == "^"
   if (inverted && falsy(value))
       renderTokens(io, token.collector, writer, ctx_push(context, ""), template, args...)
   else
       n = length(value)
       for (i,v) in enumerate(value)
           renderTokens(io, token.collector, writer, ctx_push(context, v), template, (i,n))
       end
    end
end

## ## DataFrames
## function renderTokensByValue(value::DataFrames.DataFrame, io, token, writer, context, template)
##     ## iterate along row, Call one for each row
##     for i in 1:size(value)[1]
##         renderTokens(io, token.collector, writer, ctx_push(context, value[i,:]), template)
##     end
## end

## what to do with an index value `.[ind]`?
## We have `.[ind]` being of a leaf type (values are not pushed onto a Context) so of simple usage
function _renderTokensByValue(value::AnIndex, io, token, writer, context, template, idx=(0,0))
    idx_match = (value.ind == idx[1]) || (value.ind==-1 && idx[1] == idx[2])
    if token._type == "#" || token._type == "|"
        # print if match
        if idx_match
            renderTokens(io, token.collector, writer, context, template, idx)
        end
    elseif token._type == "^"
        # print if *not* a match
        if !idx_match
            renderTokens(io, token.collector, writer, context, template, idx)
        end
    else
        renderTokens(io, token.collector, writer, ctx_push(context, value.value), template)
    end
end

function _renderTokensByValue(value::Function, io, token, writer, context, template, args...)
    ## function get's passed
    # When the value is a callable
    # object, such as a function or lambda, the object will
    # be invoked and passed the block of text. The text
    # passed is the literal block, unrendered. {{tags}} will
    # not have been expanded - the lambda should do that on
    # its own. In this way you can implement filters or
    # caching.

    #    out = (value())(token.collector, render)
    if token._type == "name"
        out = render(value(), context.view)
    elseif token._type == "|"
        # pass evaluated values
        view = context.parent.view
        sec_value = render(MustacheTokens(token.collector), view)
        out = render(value(sec_value), view)
    else
        ## How to get raw section value?
        ## desc: Lambdas used for sections should receive the raw section string.
        ## Lambdas used for sections should parse with the current delimiters.
        sec_value = toString(token.collector)
        view = context.parent.view
        tpl = value(sec_value)
        out = render(parse(tpl, (token.ltag, token.rtag)),  view)

    end
    write(io, out)

end

function _renderTokensByValue(value::Any, io, token, writer, context, template, args...)
    inverted = token._type == "^"
    if (inverted && falsy(value)) || !falsy(value)
        renderTokens(io, token.collector, writer, context, template)
    end
end

## Low-level function that renders the given `tokens` using the given `writer`
## and `context`. The `template` string is only needed for templates that use
## higher-order sections to extract the portion of the original template that
## was contained in that section.
function renderTokens(io, tokens, writer, context, template, idx=(0,0))
    for i in 1:length(tokens)
        token = tokens[i]
        tokenValue = token.value
        if token._type == "#" || token._type == "|"
            ## iterate over value if Dict, Array or DataFrame,
            ## or display conditionally
            value = lookup(context, tokenValue)
            ctx = isa(value, AnIndex) ? context : Context(value, context)
            renderTokensByValue(value, io, token, writer, ctx, template, idx)

            # if !isa(value, AnIndex)
            #    context = Context(value, context)
            # end
            # renderTokensByValue(value, io, token, writer, context, template, idx)
            # if !isa(value, AnIndex)
            #     context = ctx_pop(context)
            # end

        elseif token._type == "^"

            ## display if falsy, unlike #
            value = lookup(context, tokenValue)
            if !isa(value, AnIndex)
                context = Context(value, context)

                if falsy(value)
                    renderTokensByValue(value, io, token, writer, context, template, idx)
                end
            else
                # for indices falsy test is performed in
                # _renderTokensByValue(value::AnIndex,...)
                renderTokensByValue(value, io, token, writer, context, template, idx)
            end

        elseif token._type == ">"
            ## partials: desc: Each line of the partial should be indented before rendering.
            fname = stripWhitespace(tokenValue)
            value = lookup(context, tokenValue)

            tpl_string = ""

            if isfile(fname)
                found_partial = true
                indent = token.indent
                buf = IOBuffer()
                for (rowno, l) in enumerate(eachline(fname, keep=true))
                    # we don't strip indent from first line, so we don't indent that
                    print(buf, rowno > 0 ? indent : "", l)
                end
                tpl_string = String(take!(buf))
            else
                value = lookup(context, fname)
                if !falsy(value)
                    indent = token.indent
                    slashn = ""
                    # don't indent if last \n
                    if occursin(r"\n$", value)
                        value = chomp(value)
                        slashn = "\n"
                    end
                    buf = IOBuffer()
                    l = split(value, r"[\n]")
                    print(buf, indent*join(l, "\n"*indent))
                    tpl_string = String(take!(buf)) * slashn
                    close(buf)
                end
            end

            if tpl_string != ""
                ## Delimiters set in a parent template should not affect a partial.
                ## Delimiters set in a partial should not affect the parent template.
                tpl = parse(tpl_string, ("{{", "}}"))
                ## Issue 82, what context to pass along to the partial
                ## this makes "scope" of partial just the immediate value, not
                ## the parent
                renderTokens(io, tpl, writer, Context(context.view,nothing), template,idx)
            else
                # can't find partial. Issue #83 requests signal, but this doesn't
                # seem to match spec. Signal would go here.
            end

        elseif token._type == "<"
            ## partials without parse
            fname = stripWhitespace(tokenValue)
            if isfile(fname)
                print(io, open(x -> read(x, String), fname))
            else
                @warn("File $fname not found")
            end

        elseif token._type == "&"
            value = lookup(context, tokenValue)
            if !falsy(value)
                ## desc: A lambda's return value should parse with the default delimiters.
                ##       parse(value()) ensures that
                val = isa(value, Function) ? render(parse(value()), context.view) : value
                print(io, val)
            end

        elseif token._type == "{"
            value = lookup(context, tokenValue)
            if !falsy(value)
                val = isa(value, Function) ? render(parse(value()), context.view) : value
                print(io, val)
            end

        elseif token._type == "name"
            value = lookup(context, tokenValue)
            if !falsy(value)
                val = isa(value, Function) ? render(parse(value()), context.view) : value
                print(io, escape_html(val))
            end

        elseif token._type == "text"
            print(io, string(tokenValue))
        end

    end

end


text:
struct ChebyshevUltrasphericalConstants{D,T}
    λ::T
    M::Int
    N::Int
    nM₀::Int
    αN::T
    K::Int
end

function ChebyshevUltrasphericalConstants(c::AbstractVector{T},λ::T;M::Int=7,D::Bool=FORWARD) where T
    N = length(c)-1
    if modλ(λ) == 0 return ChebyshevUltrasphericalConstants{D,T}(λ,M,N,0,zero(T),0) end
    λm = modλ(λ)
    if D == FORWARD
        nM₀ = min(floor(Int,(eps(T)*sqrtpi*2^(M-λm-1)*factorial(M)/pochhammer(λm,M)/pochhammer(1-λm,M))^(-1/(M+1/2))),N)
        αN = min(one(T)/log(N/nM₀),half(T))
        K = N > 0 ? ceil(Int,log(N/nM₀)/log(1/αN)) : 0
    else#if D == BACKWARD
        nM₀ = min(floor(Int,(eps(T)*sqrtpi*2^(M-λm-1)*factorial(M)/pochhammer(λm,M)/pochhammer(1-λm,M))^(-1/(M+1/2))),2N)
        αN = min(one(T)/log(2N/nM₀),half(T))
        K = N > 0 ? ceil(Int,log(2N/nM₀)/log(1/αN)) : 0
    end
    ChebyshevUltrasphericalConstants{D,T}(λ,M,N,nM₀,αN,K)
end

struct ChebyshevUltrasphericalIndices
    i₁::Vector{Int}
    i₂::Vector{Int}
    j₁::Vector{Int}
    j₂::Vector{Int}
end

function ChebyshevUltrasphericalIndices(λ::T,CUC::ChebyshevUltrasphericalConstants{D,T},tempmindices::Vector{T},tempsin::Vector{T},tempsinλ::Vector{T}) where {D,T}
    M,N,nM₀,αN,K = getconstants(CUC)

    i₁,i₂ = zeros(Int,K+1),zeros(Int,K+1)
    j₁,j₂ = zeros(Int,K+1),zeros(Int,K+1)

    if D == FORWARD
        i₁[1],i₂[1],j₂[1] = 1,N+1,N+1
        for k=1:K j₁[k] = ceil(Int,αN^k*N) end
    else#if D == BACKWARD
        i₁[1],i₂[1],j₂[1] = 1,2N+1,2N+1
        for k=1:K j₁[k] = ceil(Int,αN^k*2N) end
    end

    for k=1:K
        if j₁[k] < nM₀ break end
        i₁[k+1],i₂[k+1] = findmindices!(tempmindices,λ,j₁[k],M,tempsin,tempsinλ)
    end
    for k=1:K j₂[k+1] = j₁[k]-1 end

    ChebyshevUltrasphericalIndices(i₁,i₂,j₁,j₂)
end

mutable struct ChebyshevUltrasphericalPlan{D,T,DCT,DST,SA} <: FastTransformPlan{D,T}
    CUC::ChebyshevUltrasphericalConstants{D,T}
    CUI::ChebyshevUltrasphericalIndices
    p₁::DCT
    p₂::DST
    rp::RecurrencePlan{T}
    c₁::Vector{T}
    c₂::SA
    um::Vector{T}
    vm::Vector{T}
    θ::Vector{T}
    tempsin::Vector{T}
    tempsin2::Vector{T}
    tempsinλ::Vector{T}
    tempsinλm::Vector{T}
    tempmindices::Vector{T}
    cnλ::Vector{T}
    cnmλ::Vector{T}
    w::Vector{T}
    anλ::Vector{T}
    c_cheb2::Vector{T}
    pr::Vector{T}
    function ChebyshevUltrasphericalPlan{D,T,DCT,DST,SA}(CUC::ChebyshevUltrasphericalConstants{D,T},CUI::ChebyshevUltrasphericalIndices,p₁::DCT,p₂::DST,rp::RecurrencePlan{T},c₁::Vector{T},c₂::SA,um::Vector{T},vm::Vector{T},θ::Vector{T},tempsin::Vector{T},tempsin2::Vector{T},tempsinλ::Vector{T},tempsinλm::Vector{T},tempmindices::Vector{T},cnλ::Vector{T},cnmλ::Vector{T}) where {D,T,DCT,DST,SA}
        P = new{D,T,DCT,DST,SA}()
        P.CUC = CUC
        P.CUI = CUI
        P.p₁ = p₁
        P.p₂ = p₂
        P.rp = rp
        P.c₁ = c₁
        P.c₂ = c₂
        P.um = um
        P.vm = vm
        P.θ = θ
        P.tempsin = tempsin
        P.tempsin2 = tempsin2
        P.tempsinλ = tempsinλ
        P.tempsinλm = tempsinλm
        P.tempmindices = tempmindices
        P.cnλ = cnλ
        P.cnmλ = cnmλ
        P
    end
    function ChebyshevUltrasphericalPlan{D,T,DCT,DST,SA}(CUC::ChebyshevUltrasphericalConstants{D,T},CUI::ChebyshevUltrasphericalIndices,p₁::DCT,p₂::DST,rp::RecurrencePlan{T},c₁::Vector{T},c₂::SA,um::Vector{T},vm::Vector{T},θ::Vector{T},tempsin::Vector{T},tempsin2::Vector{T},tempsinλ::Vector{T},tempsinλm::Vector{T},tempmindices::Vector{T},cnλ::Vector{T},cnmλ::Vector{T},w::Vector{T},anλ::Vector{T},c_cheb2::Vector{T},pr::Vector{T}) where {D,T,DCT,DST,SA}
        P = ChebyshevUltrasphericalPlan{D,T,DCT,DST,SA}(CUC,CUI,p₁,p₂,rp,c₁,c₂,um,vm,θ,tempsin,tempsin2,tempsinλ,tempsinλm,tempmindices,cnλ,cnmλ)
        P.w = w
        P.anλ = anλ
        P.c_cheb2 = c_cheb2
        P.pr = pr
        P
    end
    function ChebyshevUltrasphericalPlan{D,T,DCT,DST,SA}(CUC::ChebyshevUltrasphericalConstants{D,T}) where {D,T,DCT,DST,SA}
        P = new{D,T,DCT,DST,SA}()
        P.CUC = CUC
        P
    end
end

function ForwardChebyshevUltrasphericalPlan(c_ultra::AbstractVector{T},λ::T,M::Int) where T
    # Initialize constants
    CUC = ChebyshevUltrasphericalConstants(c_ultra,λ;M=M,D=FORWARD)
    if modλ(λ) == 0 return ChebyshevUltrasphericalPlan{FORWARD,T,Any,Any,Any}(CUC) end
    M,N,nM₀,αN,K = getconstants(CUC)
    λ = modλ(λ)

    # Initialize recurrence coefficients
    rp = RecurrencePlan(λ-half(λ),λ-half(λ),N+1)

    # Initialize temporary arrays
    c₁,c₂,um,vm = zero(c_ultra),view(zero(c_ultra),2:N),zero(c_ultra),zero(c_ultra)

    # Initialize DCT-I and DST-I plans
    p₁,p₂ = applyTN_plan(c₁),applyUN_plan(c₂)

    # Clenshaw-Curtis points
    θ = N > 0 ? T[k/N for k=zero(T):N] : T[0]

    # Initialize sines and cosines
    tempsin = sinpi.(θ)
    @inbounds for i = 1:N÷2 tempsin[N+2-i] = tempsin[i] end
    tempsin2 = sinpi.(θ./2)
    tempsinλ,tempmindices = zero(c_ultra),zero(c_ultra)
    @inbounds for i=1:N+1 tempsinλ[i] = tempsin[i]^λ end
    tempsinλm = similar(tempsinλ)

    # Initialize normalizing constant
    cnλ = Cnλ.(0:N,λ)
    cnmλ = similar(cnλ)

    # Get indices
    CUI = ChebyshevUltrasphericalIndices(λ,CUC,tempmindices,tempsin,tempsinλ)

    ChebyshevUltrasphericalPlan{FORWARD,T,typeof(p₁),typeof(p₂),typeof(c₂)}(CUC,CUI,p₁,p₂,rp,c₁,c₂,um,vm,θ,tempsin,tempsin2,tempsinλ,tempsinλm,tempmindices,cnλ,cnmλ)
end

function BackwardChebyshevUltrasphericalPlan(c_ultra::AbstractVector{T},λ::T,M::Int) where T
    # Initialize constants
    CUC = ChebyshevUltrasphericalConstants(c_ultra,λ;M=M,D=BACKWARD)
    if modλ(λ) == 0 return ChebyshevUltrasphericalPlan{BACKWARD,T,Any,Any,Any}(CUC) end
    M,N,nM₀,αN,K = getconstants(CUC)
    λ = modλ(λ)

    # Array of almost double the size of the coefficients
    c_cheb2 = N > 0 ? zeros(T,2N+1) : T[0]

    # Initialize recurrence coefficients
    rp,pr = RecurrencePlan(λ-half(λ),λ-half(λ),2N+1),zero(c_cheb2)

    # Initialize temporary arrays
    c₁,c₂,um,vm = zero(c_cheb2),view(zero(c_cheb2),2:2N),zero(c_cheb2),zero(c_cheb2)

    # Initialize DCT-I and DST-I plans
    p₁,p₂ = applyTN_plan(c₁),applyUN_plan(c₂)

    # Clenshaw-Curtis nodes and weights
    θ = N > 0 ? T[k/2N for k=zero(T):2N] : T[0]
    w = N > 0 ? clenshawcurtisweights!(chebyshevjacobimoments1(T, 2N+1, λ-half(λ), λ-half(λ)), p₁) : T[0]

    # Initialize sines and cosines
    tempsin = sinpi.(θ)
    @inbounds for i = 1:N tempsin[2N+2-i] = tempsin[i] end
    tempsin2 = sinpi.(θ/2)
    tempsinλ,tempmindices = zero(c_cheb2),zero(c_cheb2)
    @inbounds for i=1:2N+1 tempsinλ[i] = tempsin[i]^λ end
    tempsinλm = similar(tempsinλ)

    # Initialize normalizing constant
    cnλ = Cnλ.(0:2N,λ)
    cnmλ = similar(cnλ)

    # Initialize orthonormality constants
    anλ = Anαβ.(0:N,λ-half(λ),λ-half(λ))

    # Get indices
    CUI = ChebyshevUltrasphericalIndices(λ,CUC,tempmindices,tempsin,tempsinλ)

    ChebyshevUltrasphericalPlan{BACKWARD,T,typeof(p₁),typeof(p₂),typeof(c₂)}(CUC,CUI,p₁,p₂,rp,c₁,c₂,um,vm,θ,tempsin,tempsin2,tempsinλ,tempsinλm,tempmindices,cnλ,cnmλ,w,anλ,c_cheb2,pr)
end

getplanλ(CUC::ChebyshevUltrasphericalConstants) = CUC.λ
getplanλ(plan::ChebyshevUltrasphericalPlan) = getplanλ(plan.CUC)

getconstants(CUC::ChebyshevUltrasphericalConstants) = CUC.M,CUC.N,CUC.nM₀,CUC.αN,CUC.K
getconstants(plan::ChebyshevUltrasphericalPlan) = getconstants(plan.CUC)

getindices(CUI::ChebyshevUltrasphericalIndices) = CUI.i₁,CUI.i₂,CUI.j₁,CUI.j₂
getindices(plan::ChebyshevUltrasphericalPlan) = getindices(plan.CUI)

getplan(plan::ChebyshevUltrasphericalPlan{FORWARD})  = plan.p₁,plan.p₂,plan.rp,plan.c₁,plan.c₂,plan.um,plan.vm,plan.θ,plan.tempsin,plan.tempsin2,plan.tempsinλ,plan.tempsinλm,plan.tempmindices,plan.cnλ,plan.cnmλ
getplan(plan::ChebyshevUltrasphericalPlan{BACKWARD}) = plan.p₁,plan.p₂,plan.rp,plan.c₁,plan.c₂,plan.um,plan.vm,plan.θ,plan.tempsin,plan.tempsin2,plan.tempsinλ,plan.tempsinλm,plan.tempmindices,plan.cnλ,plan.cnmλ,plan.w,plan.anλ,plan.c_cheb2,plan.pr


text:
import Pkg
Pkg.add("Documenter")
using Documenter, BosonSampling
push!(LOAD_PATH, "./src")

DocMeta.setdocmeta!(BosonSampling, :DocTestSetup, :(using MyPackage); recursive=true)

makedocs(
    source = "./src/",
    sitename = "BosonSampling.jl",
    modules = [BosonSampling],
    authors = "Benoit Seron, Antoine Restivo",
    format = Documenter.HTML(prettyurls=false, sidebar_sitename=false),
    doctest = false,
    pages = [
        "About" => "about.md",
        "Tutorials" => Any[
                "Installation" => "tutorial/installation.md",
                "Basic usage" => "tutorial/basic_usage.md",
                "Implementing new models" => "tutorial/user_defined_models.md",
                "Samplers" => "tutorial/boson_samplers.md",
                "Partitions" => "tutorial/partitions.md",
                "Bunching" => "tutorial/bunching.md",
                "Certification" => "tutorial/certification.md",
                "Optimization" => "tutorial/optimization.md",
                "Counting statistics" => "tutorial/compute_distr.md",
                "Circuits" => "tutorial/circuits.md",
                "Permanent conjectures" => "tutorial/permanent_conjectures.md"
                        ],
        "Benchmarks" => "benchmarks/bench.md",
        "API" => Any[
            "Types" => Any[
                "Inputs" => "types/input.md",
                "Events" => "types/events.md",
                "Interferometers" => "types/interferometers.md",
                "Measurements" => "types/measurements.md",
                "Partitions" => "types/partitions.md",
                "Type utilities" => "types/type_functions.md"],
            "Functions" => Any[
                "Certification" => "functions/bayesian.md",
                "Bunching" => "functions/bunching.md",
                "Distributions" => "functions/distributions.md",
                "Partitions" => "functions/partitions.md",
                "Permanent conjectures" => "functions/permanent_conjectures.md",
                "Tools" => "functions/proba_tools.md",
                "Samplers" => "functions/samplers.md",
                "Scattering" => "functions/scattering.md",
                "Special_matrices" => "functions/special_matrices.md",
                "Visualization" => "functions/visualize.md"]
                ]
    ]
)

deploydocs(
    repo = "github.com/benoitseron/BosonSampling.jl.git",
    target = "build",
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#"],
)


text:
using DSPopt
using Test
using SparseArrays

const dsp = DSPopt.dspenv

@testset "Initializing DSPopt" begin
    @test dsp.p != C_NULL
    @test length(dsp.numRows) == 0
    @test length(dsp.numCols) == 0
    @test isnan(dsp.primVal)
    @test isnan(dsp.dualVal)
    @test length(dsp.colVal) == 0
    @test length(dsp.rowVal) == 0
    @test dsp.nblocks == -1
    @test dsp.block_ids == []
    @test dsp.is_stochastic == false
    @test dsp.solve_type == DSPopt.DW
    @test isnothing(dsp.comm)
    @test dsp.comm_size == 1
    @test dsp.comm_rank == 0
end

@testset "Check DSP versions" begin
    major = DSPopt.getVersionMajor(dsp)
    minor = DSPopt.getVersionMinor(dsp)
    patch = DSPopt.getVersionPatch(dsp)
    if major > 1
        @test true
    elseif major == 1
        if minor > 4 || (minor == 4 && patch >= 1)
            @test true
        else
            @test false
        end
    else
        @test false
    end
    @test DSPopt.getVersion(dsp) == "$(major).$(minor).$(patch)"
end

@testset "Setting options" begin
    @testset "param:" begin
        DSPopt.setoptions!(Dict(:param => "params.txt"))
    end
    @testset "is_stochastic: $i" for i in [true, false]
        DSPopt.setoptions!(Dict(:is_stochastic => i))
        @test dsp.is_stochastic == i
    end
    @testset "solve_type: $t" for t in instances(DSPopt.Methods)
        DSPopt.setoptions!(Dict(:solve_type => t))
        @test dsp.solve_type == t
    end
end

@testset "Farmer example: stochastic form" begin
    include("farmer_stoc.jl")
end

@testset "Farmer example: block form" begin
    include("farmer_block.jl")
end

@testset "dcap" begin
    include("dcap.jl")
end

if DSPopt.getVersionMajor(dsp) >= 2
    @testset "Farmer example: stochastic quadratic form" begin
        include("farmer_qcp.jl")
    end

    @testset "Farmer example2: stochastic quadratic form" begin
        include("farmer_qcp2.jl")
    end

    @testset "Farmer example3: stochastic quadratic form" begin
        include("farmer_qcp3.jl")
    end
end

@testset "Distributionally robust extension" begin
    include("Dro.jl")
end

@testset "Freeing DSPopt" begin
    DSPopt.freeEnv(dsp)
    @test dsp.p == C_NULL
    @test length(dsp.numRows) == 0
    @test length(dsp.numCols) == 0
    @test isnan(dsp.primVal)
    @test isnan(dsp.dualVal)
    @test length(dsp.colVal) == 0
    @test length(dsp.rowVal) == 0
    @test dsp.nblocks == -1
    @test dsp.block_ids == []
    @test dsp.is_stochastic == false
    @test dsp.solve_type == DSPopt.Dual
    @test isnothing(dsp.comm)
    @test dsp.comm_size == 1
    @test dsp.comm_rank == 0
end


text:
function GaussIntegrator(arg0::Pair)
    return GaussIntegrator((Pair,), arg0)
end

function GaussIntegrator(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return GaussIntegrator((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function get_number_of_points(obj::GaussIntegrator)
    return jcall(obj, "getNumberOfPoints", jint, ())
end

function get_point(obj::GaussIntegrator, arg0::jint)
    return jcall(obj, "getPoint", jdouble, (jint,), arg0)
end

function get_weight(obj::GaussIntegrator, arg0::jint)
    return jcall(obj, "getWeight", jdouble, (jint,), arg0)
end

function integrate(obj::GaussIntegrator, arg0::UnivariateFunction)
    return jcall(obj, "integrate", jdouble, (UnivariateFunction,), arg0)
end


text:
using DelayEmbeddings, StatsBase, LinearAlgebra
export PseudoPeriodic, noiseradius

"""
    PseudoPeriodic(d, τ, ρ, shift = true)

Create surrogates suitable for pseudo-periodic signals. They retain the periodic structure
of the signal, while inter-cycle dynamics that are either deterministic or correlated
noise are destroyed (for appropriate `ρ` choice).
Therefore these surrogates are suitable to test the null hypothesis
that the signal is a periodic orbit with uncorrelated noise[^Small2001].

Arguments `d, τ, ρ` are as in the paper, the embedding dimension, delay time and
noise radius. The method works by performing a delay coordinates ambedding via the
library [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystems.jl/dev/embedding/reconstruction/).
See its documentation for choosing appropriate values for `d, τ`. For `ρ`, we have implemented
the method proposed in the paper in the function [`noiseradius`](@ref).

The argument `shift` is not discussed in the paper. If `shift=false`
we adjust the algorithm so that there is little phase shift between the
periodic component of the original and surrogate data.

See also [`CycleShuffle`](@ref).

[^Small2001]: Small et al., Surrogate test for pseudoperiodic time series data, [Physical Review Letters, 87(18)](https://doi.org/10.1103/PhysRevLett.87.188101)
"""
struct PseudoPeriodic{T<:Real} <: Surrogate
    d::Int
    τ::Int
    ρ::T
    shift::Bool
end
PseudoPeriodic(d, t, r) = PseudoPeriodic(d, t, r, true)

function surrogenerator(x::AbstractVector, pp::PseudoPeriodic, rng = Random.default_rng())
    # in the following symbol `y` stands for `s` of the paper
    d, τ = getfield.(Ref(pp), (:d, :τ))
    N = length(x)
    z = embed(x, d, τ)
    Ñ = length(z)
    w = zeros(eltype(z), Ñ-1) # weights vector
    y = Dataset([z[1] for i in 1:N])
    s = similar(x)
    init = (y = y, w = w, z = z)
    return SurrogateGenerator(pp, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:PseudoPeriodic})()
    y, z, w = getfield.(Ref(sg.init), (:y, :z, :w))
    ρ, shift = getfield.(Ref(sg.method), (:ρ, :shift))
    pseudoperiodic!(sg.s, y, sg.x, z, w, ρ, shift, sg.rng)
end
# Low-level method, also used in `noiseradius`
function pseudoperiodic!(s, y, x, z, w, ρ, shift, rng)
    N, Ñ = length.((x, z))
    y[1] = shift ? rand(rng, z.data) : z[1]
    @inbounds for i in 1:N-1
        w .= (exp(-norm(z[t] - y[i])/ρ) for t in 1:Ñ-1)
        j = sample(rng, 1:Ñ-1, pweights(w))
        y[i+1] = z[j+1]
    end
    for (i, p) in enumerate(y); s[i] = p[1]; end
    return s
end

"""
    noiseradius(x::AbstractVector, d::Int, τ, ρs, n = 1) → ρ

Use the proposed* algorithm of[^Small2001] to estimate optimal `ρ` value for
[`PseudoPeriodic`](@ref) surrogates, where `ρs` is a vector of possible `ρ` values.
*The paper is ambiguous about exactly what to calculate. Here we count how many times
we have pairs of length-2 that are identical in `x` and its surrogate, but **are not**
also part of pairs of length-3.

This function directly returns the arg-maximum of the evaluated distribution of these counts
versus `ρ`, use `TimeseriesSurrogates._noiseradius` with same arguments to get the actual
distribution. `n` means to repeat τhe evaluation `n` times, which increases accuracy.

[^Small2001]: Small et al., Surrogate test for pseudoperiodic time series data, [Physical Review Letters, 87(18)](https://doi.org/10.1103/PhysRevLett.87.188101)
"""
function noiseradius(x::AbstractVector, d::Int, τ, ρs, n = 1)
    l2n = _noiseradius(x, d, τ, ρs, n)
    return ρs[argmax(l2n)]
end

function _noiseradius(x::AbstractVector, d::Int, τ, ρs, n = 1)
    l2n = noiseradius(surrogenerator(x, PseudoPeriodic(d, τ, ρs[1])), ρs, n)
end

function noiseradius(sg::SurrogateGenerator{<:PseudoPeriodic}, ρs, n = 1)
    l2n = zero(ρs) # length-2 number of points existing in both timeseries
    y, z, w = getfield.(Ref(sg.init), (:y, :z, :w))
    x, s = sg.x, sg.s
    N = length(sg.x)
    @inbounds for _ in 1:n
        for (ℓ, ρ) in enumerate(ρs)
            s̄ = pseudoperiodic!(s, y, x, z, w, ρ, sg.method.shift, sg.rng)
            for i in 1:N-1
                # TODO: This can be optimized heavily: checking x[i+2] already tells us
                # that we shouldn't check x[i+1] on the next iteration.
                l2n[ℓ] += count(j -> s̄[j]==x[i] && s̄[j+1]==x[i+1] && s̄[j+2]!=x[i+2], 1:N-2)
            end
        end
    end
    # TODO: here we can directly find maximum and return it as a single number
    return l2n
end


text:
function quicksort!(A, lo, hi)
    i = lo
    j = hi
    while i < hi
        pivot = A[(lo + hi) >>> 1]
        while i <= j
            while A[i] < pivot
                i += 1
            end

            while A[j] > pivot
                j -= 1
            end

            if i <= j
                A[i], A[j] = A[j], A[i]
                i, j = i+1, j-1
            end
        end

        if lo < j
            quicksort!(A, lo, j)
        end

        lo = i
        j = hi
    end
end

max_stars_repo_name: mjaglarz/BA-thesis

max_issues_repo_name: mjaglarz/BA-thesis

max_forks_repo_name: mjaglarz/BA-thesis


text:
# %% imports
using Statistics
using Plots
gr()

# %% based on a spike history of the simulation layout

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function spike_histogram(arch::Array{Int64,1}, spike_history::Array{Array{Int64}})
    p = Any[]
    for layer in eachindex(arch)
        count_layer = zeros(arch[layer])
        for i in 1:length(spike_history)
            count_layer .+= count_spikes(arch, spike_history[i], layer)
        end
        p_tmp = bar(count_layer,
                    xlabel = "neuron",
                    ylabel = "spiking distribution",
                    legend = false)
        push!(p,p_tmp)

    end
    display(plot(p..., layout = length(arch)))
end

function layer_spike_hist(arch::Array{Int64}, spike_history::Array{Array{Int64}}, layer::Int64)
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    gap = 10 - (length(arch) -1) + 1
    neurons = Array{Int64}[]
    for neuron in interval[1]:interval[2]
        spikebin = [(neuron in spike_history[i]) for i in 1:length(spike_history)]
        neuronspike = [sum(spikebin[i:i+gap-1]) for i in 1:gap:length(spike_history)-gap]
        push!(neurons, neuronspike)
    end
    display(plot(neurons,
                 xlabel = "Time step",
                 ylabel = "Number spikes per neuron",
                 legend = false))
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
@testset "417.pacific-atlantic-water-flow.jl" begin
    @test pacific_atlantic([1 2 2 3 5
                            3 2 3 4 4
                            2 4 5 3 1
                            6 7 1 4 5
                            5 1 1 2 4]) ==
          Set([(2, 4), (2, 5), (4, 2), (5, 1), (4, 1), (3, 3), (1, 5)])
end


text:
module Mods

import Base: (==), (+), (-), (*), (inv), (/), (//), (^), hash, show
import Base: rand, conj, iszero

export Mod, modulus, value, AbstractMod
export is_invertible
export CRT

abstract type AbstractMod <: Number end

"""
`Mod{m}(v)` creates a modular number in mod `m` with value `mod(v,m)`.
"""
struct Mod{N,T} <: AbstractMod
    val::T
end

# safe constructors (slower)
function Mod{N}(x::T) where {T<:Union{Integer,Complex{<:Integer}},N}
    @assert N isa Integer && N > 1 "modulus must be at least 2"
    Mod{N,T}(x)
end

# type casting
Mod{N,T}(x::Mod{N,T2}) where {T,N,T2} = Mod{N,T}(T(x.val))
Mod{N,T}(x::Mod{N,T}) where {T,N} = x

show(io::IO, z::Mod{N}) where {N} = print(io, "Mod{$N}($(value(z)))")
show(io::IO, ::MIME"text/plain", z::Mod{N}) where {N} = show(io, z)

"""
`modulus(a::Mod)` returns the modulus of this `Mod` number.
```
julia> a = Mod{13}(11);

julia> modulus(a)
13
```
"""
modulus(a::Mod{N}) where {N} = N

"""
`value(a::Mod)` returns the value of this `Mod` number.
```
julia> a = Mod{13}(11);

julia> value(a)
11
```
"""
value(a::Mod{N}) where {N} = mod(a.val, N)

Base.abs(a::Mod{N,<:Real} where {N}) = abs(value(a))

function hash(x::Mod, h::UInt64 = UInt64(0))
    v = value(x)
    m = modulus(x)
    return hash(v, hash(m, h))
end

# Test for equality
iszero(x::Mod{N,T}) where {N,T} = iszero(mod(x.val, N))
==(x::Mod{N,T1}, y::Mod{M,T2}) where {M,N,T1,T2} = false
==(x::Mod{N,T1}, y::Mod{N,T2}) where {N,T1,T2} = iszero(value(x - y))

# Easy arithmetic
@inline function +(x::Mod{N,T}, y::Mod{N,T}) where {N,T}
    s, flag = Base.add_with_overflow(x.val, y.val)
    if !flag
        return Mod{N,T}(s)
    end
    t = widen(x.val) + widen(y.val)    # add with added precision
    return Mod{N,T}(mod(t, N))
end

function -(x::Mod{M,T}) where {M,T<:Signed}
    if (x.val isa BigInt || x.val == typemin(T))
        return Mod{M,T}(-mod(x.val, M))
    else
        return Mod{M,T}(-x.val)
    end
end

function -(x::Mod{M,T}) where {M,T<:Unsigned}
    return Mod{M,T}(M - value(x))
end

-(x::Mod, y::Mod) = x + (-y)

@inline function *(x::Mod{N,T}, y::Mod{N,T}) where {N,T}
    p, flag = Base.mul_with_overflow(x.val, y.val)
    if !flag
        return Mod{N,T}(p)
    else
        q = widemul(x.val, y.val)         # multipy with added precision
        return Mod{N,T}(mod(q, N)) # return with proper type
    end
end

# Division stuff
"""
`is_invertible(x::Mod)` determines if `x` is invertible.
"""
function is_invertible(x::Mod{M})::Bool where {M}
    return gcd(x.val, M) == 1
end

"""
`inv(x::Mod)` gives the multiplicative inverse of `x`.
"""
@inline function inv(x::Mod{M,T}) where {M,T}
    Mod{M,T}(_invmod(x.val, M))
end
_invmod(x::Unsigned, m::Unsigned) = invmod(x, m)
# faster version of `Base.invmod`, only works for for signed types
@inline function _invmod(x::Signed, m::Signed)
    (g, v, _) = gcdx(x, m)
    if g != 1
        error("$x (mod $m) is not invertible")
    end
    return v
end

function /(x::Mod{N,T}, y::Mod{N,T}) where {N,T}
    return x * inv(y)
end

(//)(x::Mod, y::Mod) = x / y
(//)(x::Number, y::Mod{N}) where {N} = x / y
(//)(x::Mod{N}, y::Number) where {N} = x / y

Base.promote_rule(::Type{Mod{M,T1}}, ::Type{Mod{N,T2}}) where {M,N,T1,T2<:Number} =
    error("can not promote types `Mod{$M,$T1}`` and `Mod{$N,$T2}`")
Base.promote_rule(::Type{Mod{M,T1}}, ::Type{Mod{M,T2}}) where {M,T1,T2<:Number} =
    Mod{M,promote_type(T1, T2)}
Base.promote_rule(::Type{Mod{M,T1}}, ::Type{T2}) where {M,T1,T2<:Number} =
    Mod{M,promote_type(T1, T2)}
Base.promote_rule(::Type{Mod{M,T1}}, ::Type{Rational{T2}}) where {M,T1,T2} =
    Mod{M,promote_type(T1, T2)}

# Operations with rational numbers
Mod{N}(k::Rational) where {N} = Mod{N}(numerator(k)) / Mod{N}(denominator(k))
Mod{N,T}(k::Rational{T2}) where {N,T,T2} = Mod{N,T}(numerator(k)) / Mod{N,T}(denominator(k))

# Random
rand(::Type{Mod{N}}, args::Integer...) where {N} = rand(Mod{N,Int}, args...)
rand(::Type{Mod{N,T}}) where {N,T} = Mod{N}(rand(T))
rand(::Type{Mod{N,T}}, dims::Integer...) where {N,T} = Mod{N}.(rand(T, dims...))

# Chinese remainder theorem functions
"""
    `CRT([T=BigInt, ]m1, m2,...)`

Chinese Remainder Theorem.

```
julia> CRT(Int, Mod{11}(4), Mod{14}(814))
92

julia> 92%11
4

julia> 92%14
8

julia> CRT(Mod{9223372036854775783}(9223372036854775782), Mod{9223372036854775643}(9223372036854775642))
85070591730234614113402964855534653468
```

!!! note

    `CRT` uses `BigInt` by default to prevent potential integer overflow.
    If you are confident that numbers do not overflow in your application,
    please specify an optional type parameter as the first argument.
"""
function CRT(::Type{T}, remainders, primes) where {T}
    length(remainders) == length(primes) || error("size mismatch")
    isempty(remainders) && return zero(T)
    primes = convert.(T, primes)
    M = prod(primes)
    Ms = M .÷ primes
    ti = _invmod.(Ms, primes)
    mod(sum(convert.(T, remainders) .* ti .* Ms), M)
end

function CRT(::Type{T}, rs::Mod...) where {T}
    CRT(T, value.(rs), modulus.(rs))
end
CRT(rs::Mod...) = CRT(BigInt, rs...)

include("GaussMods.jl")
include("iterate.jl")

end # end of module Mods


text:
using Ogre, PyPlot

type GraphValues
    P
    ρ
end
function GraphValues(filename::AbstractString)
    graph = readcsv(filename)
    P = graph[:, 1]
    ρ = graph[:, 2]
    GraphValues(P, ρ)
end

function main()
    graph_h2o    = GraphValues("data/eos/graphs/seager/h2o.csv")
    graph_mgsio3 = GraphValues("data/eos/graphs/seager/mgsio3.csv")
    graph_fe     = GraphValues("data/eos/graphs/seager/fe.csv")

    pressures    = logspace(5, 14)

    h2o          = map(Ogre.h2o_seager, pressures)
    mgsio3       = map(Ogre.mgsio3_seager, pressures)
    fe           = map(Ogre.fe_seager, pressures)
    my_h2o_300   = map(P -> Ogre.my_h2o_full(P, 300), pressures)
    my_h2o_500   = map(P -> Ogre.my_h2o_full(P, 500), pressures)
    my_h2o_800   = map(P -> Ogre.my_h2o_full(P, 800), pressures)
    my_h2o_1200  = map(P -> Ogre.my_h2o_full(P, 1200), pressures)

    fig = figure()
    ax = subplot(111)

    plot(pressures, h2o, "k--", label="BME/TFD (isothermal ice VII)")
    plot(pressures, my_h2o_300, label="Realistic water EOS (300K)")
    plot(pressures, my_h2o_500, label="(500K)")
    plot(pressures, my_h2o_800, label="(800K)")
    plot(pressures, my_h2o_1200, label="(1200K)")
    plot(graph_h2o.P, graph_h2o.ρ, label="Seager H2O")
    plot(graph_mgsio3.P, graph_mgsio3.ρ, label="Seager MgSiO3")
    plot(graph_fe.P, graph_fe.ρ, label="Seager Fe")

    xlabel("Pressure (P) / Pa")
    ylabel(L"Density ($\rho$) / kg m$^{-3}$")
    xlim(1e5, 1e14)
    # ylim(1e1, 1e5)
    xscale("log")
    yscale("log")

    legend(loc="best")
    tight_layout()
end

main()


text:
using NPZ

abstract type Dataset end

length(datasett::Dataset) = throw("unimplemented")

batchify(x, n) = [x[i:min(i + n - 1, length(x))] for i in 1:n:length(x)]

mutable struct FashionMNIST <: Dataset
    xs::AbstractMatrix{<:Real}
    ys::AbstractVector{<:Integer}
end

length(dataset::FashionMNIST) = size(dataset.xs, 1)

function FashionMNIST(dirpath::String, set::String)
    xs = npzread("$(dirpath)/fashion-$(set)-imgs.npz")
    ys = npzread("$(dirpath)/fashion-$(set)-labels.npz")
    FashionMNIST(reshape(xs, :, size(xs, 3))', ys)
end

Base.getindex(dataset::FashionMNIST, i::Integer) = dataset.xs[i, :], dataset.ys[i]
Base.getindex(dataset::FashionMNIST, r::UnitRange{Int64}) = dataset.xs[r, :], dataset.ys[r]

function shuffle!(dataset::FashionMNIST)
    ord = randperm(length(dataset))
    dataset.xs = dataset.xs[ord, :]
    dataset.ys = dataset.ys[ord]
    return dataset
end

max_stars_repo_name: taraspiotr/DeepJulia

max_issues_repo_name: taraspiotr/DeepJulia

max_forks_repo_name: taraspiotr/DeepJulia


text:
# class constructor
QQ(::Type{T}) where T<:Integer = QQ{T}
QQ(::Type{ZZ{T}}) where T = QQ{T}

category_trait(::Type{<:QQ}) = FieldTrait

# construction
basetype(::Type{<:QQ{T}}) where T = ZZ{T}
depth(::Type{<:QQ}) = 1
issimpler(a::T, b::T) where T<:QQ = QQ(abs(a.num),a.den) < QQ(abs(b.num),b.den)
copy(a::QQ) = typeof(a)(a.num,a.den)

QQ{T}(a::QQ) where T = QQ{T}(T(a.num), T(a.den), NOCHECK)
QQ{T}(a::ZZ) where T = QQ{T}(T(a.val), one(T), NOCHECK)
QQ{T}(a::Integer) where T = QQ{T}(T(a), one(T), NOCHECK)
QQ{T}(a::Rational) where T = QQ{T}(T(a.num), T(a.den), NOCHECK)

QQ(a::QQ{T}) where T = a
QQ(a::ZZ{T}) where T = QQ{T}(a)
QQ(a::T) where T = QQ{T}(a)
QQ(a::Rational{T}) where T = QQ{T}(a)

function QQ{T}(num::Integer, den::Integer) where T<:Integer
    iszero(num) && iszero(den) && throw(ArgumentError("invalid rational: zero($T)//zero($T)"))
    num == den && return one(QQ{T})
    R = promote_type(T, typeof(signed(num)), typeof(signed(den)))
    num1, den1 = Base.divgcd(num, den)
    num2, den2 = T(num1), T(den1)
    if signbit(den2)
        den2 = -den2
        signbit(den2) && throw(ArgumentError("invalid rational: denominator $den"))
        num2 = -num2
    end
    QQ{T}(num2, den2, NOCHECK)
end
function QQ(a::S, b::T) where {S<:Integer,T<:Integer}
    R = promote_type(typeof(signed(a)),typeof(signed(b)))
    QQ{R}(R(a), R(b))
end
QQ(num::T, den::T) where T = QQ{T}(num, den)
Rational(a::QQ{T}) where T = Rational(a.num, a.den)
//(a::ZZ{T}, b::ZZ{T}) where T = QQ(Rational(a.val, b.val))

# promotion and conversion
_promote_rule(::Type{QQ{T}}, ::Type{QQ{S}}) where {S,T} = QQ{promote_type(S,T)}
_promote_rule(::Type{QQ{T}}, ::Type{ZZ{S}}) where {S,T} = QQ{promote_type(S,T)}
promote_rule(::Type{QQ{T}}, ::Type{S}) where {S<:Integer,T} = QQ{promote_type(S,T)}
promote_rule(::Type{QQ{T}}, ::Type{Rational{S}}) where {S<:Integer,T} = QQ{promote_type(S,T)}

# induced homomorphism
function (h::Hom{F,R,S})(p::QQ{<:R}) where {F,R,S}
    QQ(h.f(a.num), h.f(a.den))
end

# operations for QQ
for op in (:+, :- , :*)
    @eval begin
        ($op)(a::QQ{T}, b::QQ{T}) where T = QQ(($op)(Rational(a), Rational(b)))
    end
end
==(a::T, b::T) where T<:QQ = Rational(a) == Rational(b)
/(a::T, b::T) where T<:QQ = iszero(b) ? throw(DivideError()) : QQ(Rational(a) / Rational(b))
-(a::QQ{T}) where T = QQ{T}(checked_neg(a.num),a.den)
divrem(a::T, b::T) where T<:QQ = (a / b, zero(T))
div(a::T, b::T) where T<:QQ = a / b
rem(a::T, b::T) where T<:QQ = zero(T)
isless(a::T, b::T) where T<:QQ = isless(Rational(a), Rational(b))

isunit(a::QQ) = !iszero(a.num)
isone(a::QQ) = a.num == a.den
iszero(a::QQ) = iszero(a.num)
zero(::Type{QQ{T}}) where T = QQ{T}(zero(T), one(T), NOCHECK)
one(::Type{QQ{T}}) where T = QQ{T}(one(T), one(T), NOCHECK)
hash(a::QQ, h::UInt) = hash(Rational(a), h)

function show(io::IO, a::QQ)
    if isone(a.den)
        show(io, a.num)
    else
        print(io, a.num, '/', a.den)
    end
end


text:
#################################################################
# This file defines the abstract distribution type
# AbstractDistribution: the abstract super type for the transition and observation distributions
# DiscreteDistribution: discrete distributions support state indexing and length functions
#################################################################

@pomdp_func create_transition_distribution(pomdp::POMDP)
@pomdp_func create_observation_distribution(pomdp::POMDP)
@pomdp_func rand(rng::AbstractRNG, d::AbstractDistribution, sample::Any)
@pomdp_func pdf(d::AbstractDistribution, x::Any)
@pomdp_func iterator(d::AbstractDistribution)
@pomdp_func domain(d::AbstractDistribution)

abstract DiscreteDistribution <: AbstractDistribution

@pomdp_func Base.length(d::DiscreteDistribution)
@pomdp_func weight(d::DiscreteDistribution, i::Int)
@pomdp_func index(pomdp::POMDP, d::DiscreteDistribution, i::Int)

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
include("util/buffer.jl")
include("util/sumtree.jl")

import Random
# import Base.getindex, Base.size
import DataStructures

abstract type AbstractReplay end

function sample
end

function push!
end

abstract type AbstractWeightedReplay <: AbstractReplay end

mutable struct ExperienceReplay{CB} <: AbstractReplay
    buffer::CB
end

ExperienceReplay(size, types, shapes, column_names) = begin
    cb = CircularBuffer(size, types, shapes, column_names)
    ExperienceReplay(cb)
end

ExperienceReplayDef(size, obs_size, obs_type) =
    ExperienceReplay(size,
                     (obs_type, Int, obs_type, Float32, Bool),
                     (obs_size, 1, obs_size, 1, 1),
                     (:s, :a, :sp, :r, :t))

Base.length(er::ExperienceReplay) = length(er.buffer)
Base.getindex(er::ExperienceReplay, idx) = er.buffer[idx]

Base.push!(er::ExperienceReplay, experience) = Base.push!(er.buffer, experience)

sample(er::ExperienceReplay, batch_size) =
    sample(Random.GLOBAL_RNG, er, batch_size)

function sample(rng::Random.AbstractRNG, er::ExperienceReplay, batch_size)
    idx = rand(rng, 1:length(er), batch_size)
    return er[idx]
end

mutable struct OnlineReplay{CB<:DataStructures.CircularBuffer, T<:Tuple} <: AbstractReplay
    buffer::CB
    column_names::T
end

OnlineReplay(size, types, column_names) =
    OnlineReplay(DataStructures.CircularBuffer{Tuple{types...}}(size), tuple(column_names...))

# size(er::OnlineReplay) = size(er.buffer)
@forward OnlineReplay.buffer Base.lastindex
function getindex(er::OnlineReplay, idx)
    data = er.buffer[idx]
    NamedTuple{er.column_names}((getindex.(data, i) for i in 1:length(er.column_names)))
end
isfull(er::OnlineReplay) = DataStructures.isfull(er.buffer)
Base.push!(er::OnlineReplay, experience) = push!(er.buffer, experience)

function sample(er::OnlineReplay, batch_size; rng=Random.GLOBAL_RNG)
    @assert batch_size <= size(er.buffer)[1]
    return er[(end-batch_size+1):end]
end

warmup(::OnlineReplay, x) = x

mutable struct WeightedExperienceReplay{CB<:CircularBuffer} <: AbstractWeightedReplay
    buffer::CB
    sumtree::SumTree
end

WeightedExperienceReplay(size, types, column_names) =
    WeightedExperienceReplay(
        CircularBuffer(size, types, column_names),
        SumTree{Int}(size))

@forward WeightedExperienceReplay.buffer getindex

function Base.push!(er::WeightedExperienceReplay, experience, weight)
    idx = push!(er.buffer, experience)
    add!(er.sumtree, weight, idx)
    return
end

function sample(er::WeightedExperienceReplay, batch_size; rng=Random.GLOBAL_RNG)
    batch_idx, batch_priorities, idx = sample(er.sumtree, batch_size; rng=rng)
    return er[idx]
end


text:
# =========================================================================== #
## Signal
# =========================================================================== #

abstract type AbstractSignal end

mutable struct Signal{OutType, InType} <: AbstractSignal
    start_pin::OutType
    end_pin::InType
    right_sibling_signal::Signal

    function Signal{OutType, InType}(start_pin::OutType, end_pin::InType) where {InType <: AbstractInputPin, OutType <: AbstractOutputPin}
        return new(start_pin, end_pin)
    end
end

function Signal(start_pin::OutType, end_pin::InType) where {InType <: AbstractInputPin, OutType <: AbstractOutputPin}
    return Signal{OutType, InType}(start_pin, end_pin)
end

function has_right_sibling(signal::S) where {S <: AbstractSignal}
    return isdefined(signal, :right_sibling_signal)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia