text:
## CachedOperator

abstract type AbstractCachedArray{T,N} <: LazyArray{T,N} end
const AbstractCachedVector{T} = AbstractCachedArray{T,1}
const AbstractCachedMatrix{T} = AbstractCachedArray{T,2}

mutable struct CachedArray{T,N,DM<:AbstractArray{T,N},M<:AbstractArray{T,N}} <: AbstractCachedArray{T,N}
    data::DM
    array::M
    datasize::NTuple{N,Int}
    function CachedArray{T,N,DM,M}(data::DM, array::M, datasize::NTuple{N,Int}) where {T,N,DM<:AbstractArray{T,N},M<:AbstractArray{T,N}}
        for d in datasize
            d < 0 && throw(ArgumentError("Datasize must be 0 or more"))
        end
        new{T,N,DM,M}(data, array, datasize)
    end
end

CachedArray(data::AbstractArray{T,N}, array::AbstractArray{T,N}, datasize::NTuple{N,Int}) where {T,N} =
    CachedArray{T,N,typeof(data),typeof(array)}(data, array, datasize)

function CachedArray(data::AbstractArray{T,N}, array::AbstractArray{V,N}, datasize::NTuple{N,Int}) where {T,V,N}
    TV = promote_type(T,V)
    CachedArray(convert(AbstractArray{TV,N}, data), convert(AbstractArray{TV,N}, array), datasize)
end

const CachedVector{T,DM<:AbstractVector{T},M<:AbstractVector{T}} = CachedArray{T,1,DM,M}
const CachedMatrix{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} = CachedArray{T,2,DM,M}

# CachedArray(data::AbstractArray{T,N}, array::AbstractArray{T,N}, sz::NTuple{N,Int}) where {T,N} =
#     CachedArray{T,N,typeof(data),typeof(array)}(data, array, sz)
CachedArray(data::AbstractArray, array::AbstractArray) = CachedArray(data, array, size(data))

# function CachedArray(::Type{Diagonal}, array::AbstractMatrix{T}) where T
#     axes(array,1) == axes(array,2) || throw(DimensionMismatch("Matrix must be square to cache as diagonal"))
#     CachedArray(Diagonal(Vector{T}(undef, size(array,1))), array)
# end

CachedArray(::Type{Array}, array::AbstractArray{T,N}) where {T,N} = CachedArray(Array{T,N}(undef, ntuple(zero,N)), array)
CachedArray(array::AbstractArray{T,N}) where {T,N} = CachedArray(similar(array, ntuple(zero,N)), array)

"""
    cache(array::AbstractArray)

Caches the entries of an array.
"""
cache(::Type{MT}, O::AbstractArray) where {MT<:AbstractArray} = CachedArray(MT,O)
cache(A::AbstractArray) = _cache(MemoryLayout(typeof(A)), A)
_cache(_, O::AbstractArray) = CachedArray(O)
_cache(_, O::CachedArray) = CachedArray(copy(O.data), O.array, O.datasize)
_cache(::AbstractStridedLayout, O::AbstractArray) = copy(O)

cacheddata(A::AbstractCachedArray) = view(A.data,OneTo.(A.datasize)...)

convert(::Type{AbstractArray{T}}, S::CachedArray{T}) where T = S
convert(::Type{AbstractArray{T,N}}, S::CachedArray{T,N}) where {T,N} = S
convert(::Type{AbstractArray{T}}, S::CachedArray{<:Any,N}) where {T,N} = convert(AbstractArray{T,N}, S)
convert(::Type{AbstractArray{T,N}}, S::CachedArray{<:Any,N}) where {T,N} =
    CachedArray(convert(AbstractArray{T}, S.data), convert(AbstractArray{T}, S.array), S.datasize)

axes(A::CachedArray) = axes(A.array)
size(A::CachedArray) = size(A.array)
length(A::CachedArray) = length(A.array)

@propagate_inbounds function Base.getindex(B::AbstractCachedArray{T,N}, kj::Vararg{Integer,N}) where {T,N}
    @boundscheck checkbounds(B, kj...)
    resizedata!(B, kj...)
    B.data[kj...]
end

@propagate_inbounds function Base.getindex(B::AbstractCachedArray{T,1}, k::Integer) where T
    @boundscheck checkbounds(B, k)
    resizedata!(B, k)
    B.data[k]
end

@propagate_inbounds function Base.setindex!(B::AbstractCachedArray{T,N}, v, kj::Vararg{Integer,N}) where {T,N}
    @boundscheck checkbounds(B, kj...)
    resizedata!(B,kj...)
    @inbounds B.data[kj...] = v
    v
end

_maximum(ax, I::AbstractUnitRange{Int}) = maximum(I)
_maximum(ax, I) = maximum(ax[I])
_maximum(ax, ::Colon) = maximum(ax)
function getindex(A::AbstractCachedArray, I...)
    @boundscheck checkbounds(A, I...)
    resizedata!(A, _maximum.(axes(A), I)...)
    A.data[I...]
end

getindex(A::AbstractCachedMatrix, I::Integer) = A[Base._to_subscript_indices(A, I)...]

@inline getindex(A::AbstractCachedMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange) = layout_getindex(A, kr, jr)
@inline getindex(A::AbstractCachedMatrix, kr::AbstractVector, jr::AbstractVector) = layout_getindex(A, kr, jr)
@inline getindex(A::AbstractCachedMatrix, k::Integer, jr::AbstractVector) = layout_getindex(A, k, jr)
@inline getindex(A::AbstractCachedMatrix, k::Integer, ::Colon) = layout_getindex(A, k, :)
@inline getindex(A::AbstractCachedMatrix, kr::AbstractVector, ::Colon) = layout_getindex(A, kr, :)
@inline getindex(A::AbstractCachedMatrix, kr::AbstractUnitRange, ::Colon) = layout_getindex(A, kr, :)
@inline getindex(A::AbstractCachedMatrix, ::Colon, j::Integer) = layout_getindex(A, :, j)
@inline getindex(A::AbstractCachedMatrix, ::Colon, jr::AbstractVector) = layout_getindex(A, :, jr)
@inline getindex(A::AbstractCachedMatrix, ::Colon, jr::AbstractUnitRange) = layout_getindex(A, :, jr)
@inline getindex(A::AbstractCachedMatrix, ::Colon, ::Colon) = layout_getindex(A, :, :)

getindex(A::AbstractCachedVector, ::Colon) = copy(A)
getindex(A::AbstractCachedVector, ::Slice) = copy(A)

function cache_getindex(_, A::AbstractVector, I, J...)
    n = _maximum(axes(A,1), I)
    @boundscheck checkbounds(A, I, J...)
    isempty(I) || resizedata!(A, n)
    A.data[I]
end
# allow dispatch on resize length for infinite arrays
cache_getindex(A::AbstractVector, I, J...) = cache_getindex(length(I), A, I, J...)

getindex(A::AbstractCachedVector, I, J...) = cache_getindex(A, I, J...)
getindex(A::AbstractCachedVector, I::AbstractVector) = cache_getindex(A, I)

function getindex(A::AbstractCachedVector, I::CartesianIndex)
    resizedata!(A, Tuple(I)...)
    A.data[I]
end

function getindex(A::AbstractCachedArray, I::CartesianIndex)
    resizedata!(A, Tuple(I)...)
    A.data[I]
end

## Array caching

resizedata!(B::CachedArray, mn...) = resizedata!(MemoryLayout(B.data), MemoryLayout(B.array), B, mn...)
resizedata!(B::AbstractCachedArray, mn...) = resizedata!(MemoryLayout(B.data), UnknownLayout(), B, mn...)
resizedata!(A::AbstractArray, mn...) = A # don't do anything

function cache_filldata!(B, inds...)
    B.data[inds...] .= view(B.array,inds...)
end

function _vec_resizedata!(B::AbstractVector, n)
    n ≤ 0 && return B
    @boundscheck checkbounds(Bool, B, n) || throw(ArgumentError("Cannot resize beyound size of operator"))

    # increase size of array if necessary
    olddata = cacheddata(B)
    ν, = B.datasize
    n = max(ν,n)
    if n > length(B.data) # double memory to avoid O(n^2) growing
        B.data = similar(B.data, min(2n,length(B)))
        B.data[axes(olddata,1)] = olddata
    end
    if n > ν
        cache_filldata!(B, ν+1:n)
        B.datasize = (n,)
    end

    B
end

resizedata!(_, _, B::AbstractVector, n) = _vec_resizedata!(B, n)
resizedata!(_, _, B::AbstractVector, n::Integer) = _vec_resizedata!(B, n)

function resizedata!(_, _, B::AbstractArray{<:Any,N}, nm::Vararg{Integer,N}) where N
    @boundscheck checkbounds(Bool, B, nm...) || throw(ArgumentError("Cannot resize beyound size of operator"))

    # increase size of array if necessary
    olddata = cacheddata(B)
    νμ = size(olddata)
    nm = max.(νμ,nm)
    if νμ ≠ nm
        B.data = similar(B.data, nm...)
        B.data[axes(olddata)...] = olddata
    end

    for k in 1:N-1
        inds = tuple(axes(B.data)[1:k-1]...,νμ[k]+1:nm[k],Base.OneTo.(B.datasize[k+1:end])...)
        cache_filldata!(B, inds...)
    end
    let k = N
        inds = tuple(axes(B.data)[1:k-1]...,νμ[k]+1:nm[k])
        cache_filldata!(B, inds...)
    end
    B.datasize = nm

    B
end

convexunion(a::AbstractVector) = a

function convexunion(a::AbstractVector, b::AbstractVector)
    isempty(a) && return b
    isempty(b) && return a
    min(minimum(a),minimum(b)):max(maximum(a),maximum(b))
end

convexunion(a::AbstractVector, b::AbstractVector, c::AbstractVector...) =
    convexunion(convexunion(a,b), c...)

convexunion(a::Number, b::Number, c...) = convexunion(a:a, b:b, c...)
convexunion(a::Number, b, c...) = convexunion(a:a, b, c...)
convexunion(a, b::Number, c...) = convexunion(a, b:b, c...)

_tounitrange(a) = first(a):last(a)

function colsupport(A::CachedMatrix, i)
    isempty(i) && return 1:0
    _tounitrange(minimum(i) ≤ A.datasize[2] ? convexunion(colsupport(A.array, i),colsupport(A.data,i) ∩ Base.OneTo(A.datasize[1])) : colsupport(A.array, i))
end
colsupport(A::CachedVector, i) =
    convexunion(colsupport(A.array, i),colsupport(A.data,i) ∩ Base.OneTo(A.datasize[1]))
function rowsupport(A::CachedMatrix, i)
    isempty(i) && return 1:0
    minimum(i) ≤ A.datasize[1] ? convexunion(rowsupport(A.array, i),rowsupport(A.data,i) ∩ Base.OneTo(A.datasize[2])) : rowsupport(A.array, i)
end

###
# special for zero cache
###

function zero!(A::CachedArray{<:Any,N,<:Any,<:Zeros}) where N
    zero!(A.data)
    A
end
function cache_getindex(_, A::CachedVector{<:Any,<:AbstractVector,<:AbstractFill{<:Any,1}}, I::AbstractVector)
    @boundscheck checkbounds(A, I)
    CachedArray(A.data[I ∩ OneTo(A.datasize[1])], A.array[oneto(length(I))])
end

###
# MemoryLayout
####

struct CachedLayout{Data,Array} <: MemoryLayout end

cachedlayout(::Data, ::Array) where {Data,Array} = CachedLayout{Data,Array}()
MemoryLayout(C::Type{CachedArray{T,N,DAT,ARR}}) where {T,N,DAT,ARR} = cachedlayout(MemoryLayout(DAT), MemoryLayout(ARR))

#####
# broadcasting
#
# We want broadcasting for numbers with concaenations to pass through
# to take advantage of special implementations of the sub-components
######

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

broadcasted(::LazyArrayStyle, op, A::CachedArray) = CachedArray(broadcast(op, cacheddata(A)), broadcast(op, A.array))
broadcasted(::LazyArrayStyle, op, A::CachedArray, c::Number) = CachedArray(broadcast(op, cacheddata(A), c), broadcast(op, A.array, c))
broadcasted(::LazyArrayStyle, op, c::Number, A::CachedArray) = CachedArray(broadcast(op, c, cacheddata(A)), broadcast(op, c, A.array))
broadcasted(::LazyArrayStyle, op, A::CachedArray, c::Ref) = CachedArray(broadcast(op, cacheddata(A), c), broadcast(op, A.array, c))
broadcasted(::LazyArrayStyle, op, c::Ref, A::CachedArray) = CachedArray(broadcast(op, c, cacheddata(A)), broadcast(op, c, A.array))

function layout_broadcasted(::CachedLayout, _, op, A::AbstractVector, B::AbstractVector)
    dat = cacheddata(A)
    n = length(dat)
    m = length(B)
    CachedArray(broadcast(op, dat, view(B,1:n)), broadcast(op, A.array, B))
end

function layout_broadcasted(_, ::CachedLayout, op, A::AbstractVector, B::AbstractVector)
    dat = cacheddata(B)
    n = length(dat)
    m = length(A)
    CachedArray(broadcast(op, view(A,1:n), dat), broadcast(op, A, B.array))
end

function layout_broadcasted(::CachedLayout, ::CachedLayout, op, A::AbstractVector, B::AbstractVector)
    n = max(A.datasize[1],B.datasize[1])
    resizedata!(A,n)
    resizedata!(B,n)
    Adat = view(cacheddata(A),1:n)
    Bdat = view(cacheddata(B),1:n)
    CachedArray(broadcast(op, Adat, Bdat), broadcast(op, A.array, B.array))
end

function layout_broadcasted(op, A, B)
    if length(A) ≠ length(B)
        (length(A) == 1 || length(B) == 1) && error("Internal error: Scalar-like broadcasting not yet supported.")
        throw(DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths $(length(A)) and $(length(B))"))
    end
    layout_broadcasted(MemoryLayout(A), MemoryLayout(B), op, A, B)
end

broadcasted(::LazyArrayStyle, op, A::CachedVector, B::AbstractVector) = layout_broadcasted(op, A, B)
broadcasted(::LazyArrayStyle, op, A::AbstractVector, B::CachedVector) = layout_broadcasted(op, A, B)
broadcasted(::LazyArrayStyle, op, A::CachedVector, B::CachedVector) = layout_broadcasted(op, A, B)
broadcasted(::LazyArrayStyle, op, A::Broadcasted, B::CachedVector) = broadcast(op, materialize(A), B)

broadcasted(::LazyArrayStyle, op, A::SubArray{<:Any,1,<:CachedMatrix}, B::CachedVector) = layout_broadcasted(op, A, B)
broadcasted(::LazyArrayStyle, op, A::SubArray{<:Any,1,<:CachedMatrix}, B::AbstractVector) = layout_broadcasted(op, A, B)
broadcasted(::LazyArrayStyle, op, A::CachedVector, B::SubArray{<:Any,1,<:CachedMatrix}) = layout_broadcasted(op, A, B)
broadcasted(::LazyArrayStyle, op, A::AbstractVector, B::SubArray{<:Any,1,<:CachedMatrix}) = layout_broadcasted(op, A, B)

broadcasted(::LazyArrayStyle{1}, op, a::CachedVector, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), op, a, b)
broadcasted(::LazyArrayStyle{1}, op, a::Zeros{<:Any,1}, b::CachedVector) = broadcast(DefaultArrayStyle{1}(), op, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::CachedVector, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), *, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::Zeros{<:Any,1}, b::CachedVector) = broadcast(DefaultArrayStyle{1}(), *, a, b)

###
# norm
###

# allow overloading for special backends, e.g., padded
_norm2(_, a) = sqrt(norm(cacheddata(a),2)^2 + norm(@view(a.array[a.datasize[1]+1:end]),2)^2)
_norm1(_, a) = norm(cacheddata(a),1) + norm(@view(a.array[a.datasize[1]+1:end]),1)
_normInf(_, a) = max(norm(cacheddata(a),Inf), norm(@view(a.array[a.datasize[1]+1:end]),Inf))
_normp(_, a, p) = (norm(cacheddata(a),p)^p + norm(@view(a.array[a.datasize[1]+1:end]),p)^p)^inv(p)

norm1(a::CachedVector) = _norm1(MemoryLayout(a), a)
norm2(a::CachedVector) = _norm2(MemoryLayout(a), a)
normInf(a::CachedVector) = _normInf(MemoryLayout(a), a)
normp(a::CachedVector, p) = _normp(MemoryLayout(a), a, p)

###
# fill!/lmul!/rmul!
###

function fill!(a::CachedArray, x)
    fill!(a.data, x)
    fill!(a.array, x)
    a
end

function rmul!(a::CachedArray, x::Number)
    rmul!(a.data, x)
    rmul!(a.array, x)
    a
end

function lmul!(x::Number, a::CachedArray)
    lmul!(x, a.data)
    lmul!(x, a.array)
    a
end

lmul!(x::Number, a::SubArray{<:Any,N,<:CachedArray}) where N = ArrayLayouts.lmul!(x, a)
rmul!(a::SubArray{<:Any,N,<:CachedArray}, x::Number) where N = ArrayLayouts.rmul!(a, x)

###
# copy
###

# need to copy data to prevent mutation. `a.array` is never changed so does not need to be
# copied
copy(a::CachedArray) = CachedArray(copy(a.data), a.array, a.datasize)
copy(a::Adjoint{<:Any,<:CachedArray}) = copy(parent(a))'
copy(a::Transpose{<:Any,<:CachedArray}) = transpose(copy(parent(a)))

###
# special maximum/minimum/sum
# TODO: when view returns a Fill this can be generalised
###

function maximum(a::CachedVector{<:Any,<:Any,<:AbstractFill})
    data = cacheddata(a)
    if length(data) < length(a)
        max(maximum(data), getindex_value(a.array))
    else
        maximum(data)
    end
end

function minimum(a::CachedVector{<:Any,<:Any,<:AbstractFill})
    data = cacheddata(a)
    if length(data) < length(a)
        min(minimum(data), getindex_value(a.array))
    else
        minimum(data)
    end
end

function sum(a::CachedVector{<:Any,<:Any,<:AbstractFill})
    data = cacheddata(a)
    sum(data) + sum(a.array[length(data)+1:end])
end

permutedims(a::CachedMatrix) = CachedArray(permutedims(a.data), permutedims(a.array), reverse(a.datasize))
permutedims(a::CachedVector) = CachedArray(permutedims(a.data), permutedims(a.array), (1,a.datasize[1]))

##
# hide type of AbstactArray
##
const CachedAbstractArray{T,N} = CachedArray{T,N,Array{T,N},AbstractArray{T,N}}
const CachedAbstractVector{T} = CachedAbstractArray{T,1}
const CachedAbstractMatrix{T} = CachedAbstractArray{T,2}

CachedAbstractArray{T,N}(data::AbstractArray{T,N}, array::AbstractArray{T,N}) where {T,N} = CachedAbstractArray{T,N}(data, array, size(data))
CachedAbstractArray{T,N}(array::AbstractArray{T,N}) where {T,N} = CachedAbstractArray{T,N}(Array{T,N}(undef, ntuple(zero,N)), array)
CachedAbstractArray(array::AbstractArray{T,N}) where {T,N} = CachedAbstractArray{T,N}(array)
CachedAbstractArray(data::AbstractArray{T,N}, array::AbstractArray{T,N}) where {T,N} = CachedAbstractArray{T,N}(data, array)

# TODO: uncomment when Julia v1.6 is required
# (CachedAbstractArray{V,N} where V)(array::AbstractArray{T,N}) where {T,N} = CachedAbstractArray{T,N}(array)
CachedAbstractVector(array::AbstractVector{T}) where T = CachedAbstractVector{T}(array)
CachedAbstractMatrix(array::AbstractMatrix{T}) where T = CachedAbstractMatrix{T}(array)

broadcasted(::LazyArrayStyle, op, A::CachedAbstractArray) = CachedAbstractArray(broadcast(op, cacheddata(A)), broadcast(op, A.array))
function broadcasted(::LazyArrayStyle, op, A::CachedAbstractVector, B::CachedAbstractVector)
    n = max(A.datasize[1],B.datasize[1])
    resizedata!(A,n)
    resizedata!(B,n)
    Adat = view(cacheddata(A),1:n)
    Bdat = view(cacheddata(B),1:n)
    CachedAbstractArray(broadcast(op, Adat, Bdat), broadcast(op, A.array, B.array))
end
broadcasted(::LazyArrayStyle, op, A::CachedAbstractArray, c::Number) = CachedAbstractArray(broadcast(op, cacheddata(A), c), broadcast(op, A.array, c))
broadcasted(::LazyArrayStyle, op, c::Number, A::CachedAbstractArray) = CachedAbstractArray(broadcast(op, c, cacheddata(A)), broadcast(op, c, A.array))
broadcasted(::LazyArrayStyle, op, A::CachedAbstractArray, c::Ref) = CachedAbstractArray(broadcast(op, cacheddata(A), c), broadcast(op, A.array, c))
broadcasted(::LazyArrayStyle, op, c::Ref, A::CachedAbstractArray) = CachedAbstractArray(broadcast(op, c, cacheddata(A)), broadcast(op, c, A.array))

###
# copyto!
###

function _copyto!(_, ::CachedLayout, dest::AbstractArray{T,N}, src::AbstractArray{V,N}) where {T,V,N}
    resizedata!(src, size(dest)...)
    copyto!(dest, view(cacheddata(src), axes(dest)...))
end

###
# SubArray
###

sublayout(::CachedLayout{MLAY,ALAY}, ::Type{I}) where {MLAY,ALAY,I} =
    cachedlayout(sublayout(MLAY(),I), sublayout(ALAY,I))

function resizedata!(V::SubArray, n::Integer...)
    resizedata!(parent(V), getindex.(parentindices(V), n)...)
    V
end

function resizedata!(v::SubArray{<:Any,1,<:AbstractMatrix,<:Tuple{AbstractVector,Integer}}, m::Integer)
    kr,j = parentindices(v)
    resizedata!(parent(v), maximum(view(kr,1:m)), j)
    v
end

function resizedata!(v::SubArray{<:Any,1,<:AbstractMatrix,<:Tuple{Integer,AbstractVector}}, m::Integer)
    k,jr = parentindices(v)
    resizedata!(parent(v), k, maximum(view(jr,1:m)))
    v
end

function cacheddata(V::SubArray)
    P = parent(V)
    data = cacheddata(P)
    view(data, intersect.(axes(data), parentindices(V))...)
end


text:
# changed node numbering on 5/28 when readSubTree will start internal
# node numbers at -2 to avoid error in undirectedOtherNetworks

if !(@isdefined doalltests) doalltests = false; end

@testset "test: auxiliary" begin
global net
net = readTopology("((((B:102.3456789)#H1)#H2,((D:0.00123456789,C,#H2:::0.123456789)S1,(#H1,A_coolname)S2)S3)S4);")
# using PhyloPlots; plot(net, :R, showEdgeNumber=true, showNodeNumber=true);
s = IOBuffer()
@test_logs printEdges(s, net)
@test String(take!(s)) == """
edge parent child  length  hybrid isMajor gamma   containRoot inCycle istIdentitiable
1    2      1      102.346 false  true    1       false       -1      false
2    3      2              true   true            false       -1      true
3    10     3              true   true    0.8765  true        -1      true
4    6      4      0.001   false  true    1       true        -1      false
5    6      5              false  true    1       true        -1      false
6    6      3              true   false   0.1235  true        -1      true
7    9      6              false  true    1       true        -1      true
8    8      2              true   false           true        -1      true
9    8      7              false  true    1       true        -1      false
10   9      8              false  true    1       true        -1      true
11   10     9              false  true    1       true        -1      true
"""
close(s); s = IOBuffer()
@test_logs printNodes(s, net)
@test String(take!(s)) == """
node leaf  hybrid hasHybEdge name       inCycle edges'numbers
1    true  false  false      B          -1      1
2    false true   true       H1         -1      1    2    8
3    false true   true       H2         -1      2    3    6
4    true  false  false      D          -1      4
5    true  false  false      C          -1      5
6    false false  true       S1         -1      4    5    6    7
7    true  false  false      A_coolname -1      9
8    false false  true       S2         -1      8    9    10
9    false false  false      S3         -1      7    10   11
10   false false  true       S4         -1      3    11
"""
@test_throws ErrorException PhyloNetworks.getMinorParent(net.node[1])
@test PhyloNetworks.getMajorParent(net.node[1]).number == 2
@test PhyloNetworks.getMinorParent(net.node[3]).number == 6
@test PhyloNetworks.getMajorParent(net.node[3]).number == 10
@test PhyloNetworks.getMajorParentEdge(net.node[6]).number == 7
@test_throws ErrorException PhyloNetworks.getMinorParentEdge(net.node[6])
@test PhyloNetworks.getMajorParentEdge(net.node[2]).number == 2
@test PhyloNetworks.getMinorParentEdge(net.node[2]).number == 8
@test [n.number for n in PhyloNetworks.getChildren(net.node[4])] == [] # leaf
@test [n.number for n in PhyloNetworks.getChildren(net.node[2])] == [1] # hybrid node
@test [n.number for n in PhyloNetworks.getChildren(net.node[9])] == [6,8] # tree node
@test [n.number for n in PhyloNetworks.getChildren(net.node[10])] == [3,9] # at root
@test [n.number for n in PhyloNetworks.getChildren(net.node[6])] == [4,5,3] # polytomy
@test PhyloNetworks.getParent(net.edge[8]).number == 8
@test [n.number for n in PhyloNetworks.getParents(net.node[3])] == [10, 6]
@test [n.number for n in PhyloNetworks.getParents(net.node[6])] == [9]
@test_throws ErrorException deleteleaf!(net, net.node[9])
n = deepcopy(net)
@test_logs deleteleaf!(n, n.node[7])
@test n.numNodes == 8; @test n.numEdges == 9;
@test_logs deleteleaf!(net, net.node[7], simplify=false)
deleteleaf!(net, 4, simplify=false); deleteleaf!(net, 5, simplify=false)
@test net.numNodes == 5; @test net.numEdges == 6;
end

@testset "testing directEdges! and re-rootings" begin
# on a tree, then on a network with h=2"

if doalltests
tre = readTopology("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
tre.edge[1].isChild1=false; tre.edge[17].isChild1=false
directEdges!(tre)
tre.edge[1].isChild1  || error("directEdges! didn't correct the direction of 1st edge")
tre.edge[17].isChild1 || error("directEdges! didn't correct the direction of 17th edge")
for i=1:18 tre.edge[i].containRoot=false; end;
tre.root = 9;
directEdges!(tre);
!tre.edge[9].isChild1 || error("directEdges! didn't correct the direction of 9th edge")
for i=1:18
 tre.edge[i].containRoot || error("directEdges! didn't correct containRoot of $(i)th edge.")
end
tre = readTopology("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
rootatnode!(tre, -9); ## clau: previously -8
end

global net
net = readTopology("(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,(((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);");
# 5th node = node number -7 (clau: previously -6).
net.root = 5
@test_logs directEdges!(net);
@test !net.edge[12].isChild1 # or: directEdges! didn't correct the direction of 12th edge
@test !net.edge[23].isChild1 # or: directEdges! didn't correct the direction of 23th edge"
@test [!net.edge[i].containRoot for i in [8;collect(13:17)]] == [true for i in 1:6]
# or: "directEdges! didn't correct containRoot below a hyb node, $(i)th edge."
@test [net.edge[i].containRoot for i in [9,5,18,2]] == [true for i in 1:4]
# or: "directEdges! didn't correct containRoot of hyb edges."
# plot(net, showNodeNumber=true, showEdgeLength=false, showEdgeNumber=true)
@test_logs rootatnode!(net, -10); # or: rootatnode! complained, node -10
@test_throws PhyloNetworks.RootMismatch rootatnode!(net, "M"; verbose=false);
# println("the rootmismatch about node 5 is good and expected.")
@test_logs rootonedge!(net, 9); # or: rootonedge! complained, edge 9
@test_logs PhyloNetworks.fuseedgesat!(27, net);
# earlier warning: """node 1 is a leaf. Will create a new node if needed, to set taxon "Ag" as outgroup."""
@test_logs rootatnode!(net, "Ag"); # need for new node
# earlier: """node 1 is a leaf. Will create a new node if needed, to set taxon "Ag" as outgroup."""
@test_logs rootatnode!(net, "Ag"); # no need this time
@test length(net.node) == 27 # or: wrong # of nodes after rootatnode! twice on same outgroup
# earlier: """node 10 is a leaf. Will create a new node if needed, to set taxon "Ap" as outgroup."""
@test_logs rootatnode!(net, "Ap");
@test length(net.node) == 27 # or: wrong # of nodes, after 3rd rooting with outgroup
@test_logs rootonedge!(net, 5);

# example with one hybridization below another

if doalltests
net = readTopology("((((((((1,2),3),4),(5)#H1),(#H1,(6,7))))#H2,(8,9)),(#H2,10));");
# sum([!e.containRoot for e in net.edge]) # only 4.
directEdges!(net); # or error("directEdges! says that the root position is incompatible with hybrids")
sum([!e.containRoot for e in net.edge]) == 16 ||
 error("directEdges! wrong on net with 2 stacked hybrids");
plot(net, showEdgeNumber=true, showEdgeLength=false, showNodeNumber=true);
net = readTopology("((((((((1,2),3),4),(5)#H1),(#H1,(6,7))))#H2,(8,9)),(#H2,10));");
net.root=19; # node number -13 (clau: previously -12)
directEdges!(net); # or error("directEdges! says that the root position is incompatible with hybrids");
end

net = readTopology("((((((((1,2),3),4),(5)#H1),(#H1,(6,7))))#H2,(8,9)),(#H2,10));");
net.root=15; # node number -5 (clau: previously -4)
@test_throws PhyloNetworks.RootMismatch directEdges!(net);
# occursin(r"non-leaf node 9 had 0 children",e.msg))
@test_logs rootatnode!(net, -13); # or: rootatnode complained...
@test_throws PhyloNetworks.RootMismatch rootatnode!(net, -5); # verbose = true this time
# occursin(r"non-leaf node 9 had 0 children", e.msg))
@test_throws PhyloNetworks.RootMismatch rootatnode!(net,"H2"; verbose=false); #try rethrow();
# occursin(r"hybrid edge 17 conflicts", e.msg))
# earlier: """node 12 is a leaf. Will create a new node if needed, to set taxon "10" as outgroup."""
@test_logs rootatnode!(net,"10");

end # of testset for directEdges! and re-rootings

@testset "testing preorder!" begin
# on a tree, then on a network with h=2
global net
tre = readTopology("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
net = readTopology("(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,(((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);");

if doalltests
preorder!(tre)
## clau: previously [-1,10,-2,-9,9,8,-3,-8,7,6,-4,5,-5,4,-6,3,-7,2,1];
nodeN = [-2,10,-3,-10,9,8,-4,-9,7,6,-5,5,-6,4,-7,3,-8,2,1];
for i=1:length(tre.node)
  tre.nodes_changed[i].number==nodeN[i] ||
    error("node pre-ordered $i is node number $(tre.nodes_changed[i].number) instead of $(nodeN[i])")
end
end

@test_logs preorder!(net)
## clau previously: [-1,14,-14,13,12,-2,-9,11,-10,10,-3,-4,-5,-6,-7,5,6,4,3,2,-12,9,-13,8,7,1];
nodeN = [-2,14,-15,13,12,-3,-10,11,-11,10,-4,-5,-6,-7,-8,5,6,4,3,2,-13,9,-14,8,7,1];
@test [n.number for n in net.nodes_changed] == nodeN

cui3str = "(Xmayae,((Xhellerii,(((Xclemenciae_F2,Xmonticolus):1.458,(((((Xmontezumae,(Xnezahuacoyotl)#H26:0.247::0.804):0.375,((Xbirchmanni_GARC,Xmalinche_CHIC2):0.997,Xcortezi):0.455):0.63,(#H26:0.0::0.196,((Xcontinens,Xpygmaeus):1.932,(Xnigrensis,Xmultilineatus):1.401):0.042):2.439):2.0)#H7:0.787::0.835,(Xmaculatus,(Xandersi,(Xmilleri,((Xxiphidium,#H7:9.563::0.165):1.409,(Xevelynae,(Xvariatus,(Xcouchianus,(Xgordoni,Xmeyeri):0.263):3.532):0.642):0.411):0.295):0.468):0.654):1.022):0.788):1.917)#H27:0.149::0.572):0.668,Xalvarezi):0.257,(Xsignum,#H27:1.381::0.428):4.669);"
net3  = readTopology(cui3str);
deleteleaf!(net3,"Xhellerii"); deleteleaf!(net3,"Xsignum");
deleteleaf!(net3,"Xmayae", simplify=false);
# now: net3 has a 2-cycle
directEdges!(net3)
preorder!(net3)
@test [n.number for n in net3.nodes_changed] == [-3,25,-6,-8,-20,-21,-22,-23,-25,-26,-27,-28,24,23,22,21,20,-24,15,-10,-16,-17,-19,14,13,-18,12,11,-11,-14,10,-15,9,8,-12,7,6,5,19,18,17,16,-7,4,3,26]

end # of testset for preorder!

@testset "testing cladewiseorder!" begin # on a tree then network with h=2
# cladewiseorder! is used for plotting: to avoid crossing edges in main tree

if doalltests
cladewiseorder!(tre)
nodeN = collect(19:-1:1);
for i=1:length(tre.node)
  tre.cladewiseorder_nodeIndex[i]==nodeN[i] ||
    error("node clade-wise ordered $i is $(tre.cladewiseorder_nodeIndex[i])th node instead of $(nodeN[i])th")
end
end

@test_logs cladewiseorder!(net)
nodeN = collect(26:-1:1);
@test net.cladewiseorder_nodeIndex == nodeN
end # of testset for cladewiseorder!

@testset "testing rotate!" begin
# to change the order of children edges at a given node

if doalltests
net = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
rotate!(net, -5) ## clau: previously -4
[e.number for e in net.node[13].edge] == [14,12,15] || error("rotate didn't work at node -5"); ## clau: previously -4
plot(net); # just to check no error.
end

global net
net=readTopology("(4,((1,(2)#H7:::0.864):2.069,(6,5):3.423):0.265,(3,#H7:::0.136):10.0);");
@test_logs rotate!(net, -2, orderedEdgeNum=[1,12,9])
@test [e.number for e in net.node[12].edge] == [1,12,9] # or: rotate didn't work at node -2

end # of testset for rotate

ext: jl


text:
macro template(arr...)
    feat = filter(item->item.head == :tuple, arr[1].args)
    nfeat = length(feat)
    res = :(Int[])
    for i = 1:nfeat
        val = :( tohash(len, $(i), $(feat[i].args...)); )
        push!(res.args, val)
    end
    res
end

function tohash(len::Int, arr::Int...)
    res = 1
    for i in arr
        res += i
        res += (res << 10)
        res $= (res >> 6)
    end
    res += (res << 3)
    res $= (res >> 11)
    res += (res << 15)
    abs(res) % len + 1
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
@testset "glossary" begin
    @testset "isolate_gloss" begin
        gloss = "like"
        isolike(x) = isolate_gloss(x, gloss)
        @test isolike("") == [""]
        @test isolike("word") == ["word"]
        @test isolike("like") == ["like"]
        @test isolike("likeword") == ["like", "word"]
        @test isolike("wordlikeword") == ["word", "like", "word"]
        @test isolike("likelike") == ["like", "like"]
        @test isolike("wordlikewordlike") == ["word", "like", "word", "like"]
    end

    @testset "multiple gloss" begin
        glossaries = ["like", "Manuel", "USA"]
        @test isolate_gloss("wordlikeUSAwordManuelManuelwordUSA" , glossaries) == ["word", "like", "USA", "word", "Manuel", "Manuel", "word", "USA"]
    end

    @testset "regex gloss" begin
        glossaries = [r"<country>\w*</country>", r"<name>\w*</name>", r"\d+"]
        @test isolate_gloss("wordlike<country>USA</country>word10001word<name>Manuel</name>word<country>USA</country>", glossaries) == ["wordlike", "<country>USA</country>", "word", "10001", "word", "<name>Manuel</name>", "word", "<country>USA</country>"]
    end

    @testset "bpe with gloss" begin
        glossaries = ["like", "Manuel", "USA"]

        #mocking bpe
        bpefile = joinpath(dirname(@__FILE__), "data/bpe.out")
        bpe = Bpe(bpefile; sepsym="@@", endsym="", glossaries = glossaries)
        import BytePairEncoding: Bpe
        function (b::Bpe)(x::String)::Tuple
            if x ∈ b.glossaries
                tuple(x)
            else
                l = length(x)
                tuple(x[1:div(l,2)], x[div(l,2)+1:end])
            end
        end

        @test join(segment(bpe, "wordlikeword likeManuelword"), " ") == "wo@@ rd@@ like@@ wo@@ rd like@@ Manuel@@ wo@@ rd"

    end
end


text:
using JetPack, Jets, Test

@testset "JopDiagonal, non constant diagonal, T=$T" for T in (Float64, Complex{Float64})
    diagonal = rand(T,512)
    A = JopDiagonal(diagonal)
    m = rand(domain(A))
    @test A*m ≈ diagonal .* m
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopDiagonal, non constant real diagonal, complex spaces" begin
    diagonal = rand(512)
    A = JopDiagonal(JetSpace(Complex{Float64}, 512), diagonal)
    m = rand(domain(A))
    @test eltype(m) <: Complex
    @test A*m ≈ diagonal .* m
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopDiagonal, constant diagonal" begin
    diagonal = rand(Float64)
    A = JopDiagonal(JetSpace(Float64,512), diagonal)
    m = rand(domain(A))
    @test A*m ≈ diagonal .* m
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end


text:
mutable struct ExtendedDMD{D,O,C} <: abstractKoopmanOperator
    koopman::D
    output::O
    basis::C
end

# Make the struct callable for transformations
(m::ExtendedDMD)(u; p = []) = m.basis(u, p = p)

# Some nice functions
LinearAlgebra.eigen(m::ExtendedDMD) = eigen(m.koopman)
LinearAlgebra.eigvals(m::ExtendedDMD) = eigvals(m.koopman)
LinearAlgebra.eigvecs(m::ExtendedDMD) = eigvecs(m.koopman)

function ExtendedDMD(X::AbstractArray, Ψ::abstractBasis; p::AbstractArray = [],  B::AbstractArray = reshape([], 0,0), dt::T = 1.0)  where T <: Real
    return ExtendedDMD(X[:, 1:end-1], X[:, 2:end], Ψ, p = p, B = B, dt = dt)
end

function ExtendedDMD(X::AbstractArray, Y::AbstractArray, Ψ::abstractBasis; p::AbstractArray = [], B::AbstractArray = reshape([], 0,0), dt::T = 1.0) where T <: Real
    @assert dt >= zero(typeof(dt)) "Provide positive dt"
    @assert size(X)[2] .== size(Y)[2] "Provide consistent dimensions for data"
    @assert size(Y)[1] .<= size(Y)[2] "Provide consistent dimensions for data"

    # Based upon William et.al. , A Data-Driven Approximation of the Koopman operator

    # Number of states and measurements
    N,M = size(X)

    # Compute the transformed data
    Ψ₀ = hcat([Ψ(xi, p = p) for xi in eachcol(X)]...)
    Ψ₁ = hcat([Ψ(xi, p = p) for xi in eachcol(Y)]...)
    Op = ExactDMD(Ψ₀, Ψ₁) # Initial guess based upon the basis

    # Transform back to states
    if isempty(B)
        B = X*pinv(Ψ₀)
    end

    # TODO Maybe reduce the observable space here
    return ExtendedDMD(Op, B, Ψ)
end

function update!(m::ExtendedDMD, x::AbstractArray, y::AbstractArray; p::AbstractArray = [], dt::T = 0.0, threshold::Float64 = 1e-3) where T <: Real
    Ψ₀ = m.basis(x, p = p)
    Ψ₁ = m.basis(y, p = p)
    update!(m.koopman, Ψ₀, Ψ₁, dt = dt, threshold = threshold)
    return
end

# TODO How to implement continouos time dynamics?
# We would need ∂Ψ/∂x or ∂Ψ/∂t
function dynamics(m::ExtendedDMD)
    # Create a set of nonlinear eqs
    p_ = m.output*m.koopman.Ã
    function dudt_(du, u, p, t)
        mul!(du,p_,m.basis(u, p = p))
    end
end

function linear_dynamics(m::ExtendedDMD; discrete::Bool = true)
    return dynamics(m.koopman, discrete = discrete)
end

# Reduction for basis
# This is a fairly naive approach of doing this
function reduce_basis(m::ExtendedDMD; threshold = 1e-5)
    b = m.output*m.koopman.Ã
    inds = sum(abs, b, dims = 1) .> threshold
    return Basis(m.basis.basis[vec(inds)], variables(m.basis), parameters = parameters(m.basis))
end


text:
##
using Plots
gr()
x = range(-2, stop=2, length=100)
y = range(sqrt(2), stop=2, length=100)
f(x,y) = x * y - x - y + 1

plot(x,y,f,st=:surface)

##
plot(rand(3))
plot!(rand(5),bordercolor=:white,
gridcolor=:white
)

##
gr()
x = -1.5:0.1:1.5
y = -1:0.1:1
z(x,y) = 3.434

plt = surface(x, y, z, fillcolor=:blue, α=0.8, grid=false, )

##
using Plots
gr()
f(x,y) = (3x + y^2) * abs(sin(x) + cos(y))

x = 0:0.1:20
y = 0:0.1:10
z = [f(i, j) for i in x, j in y]' # この転置を忘れるとデータが矛盾し、グラフが変になるので要注意．

plot(x,y,z, st=:surface, camera=(30, 50), ratio=:equal, proj=:"3d")

##
using PyPlot

thetas = range(0, stop=pi,   length=50)
phis   = range(0, stop=pi / 2, length=50)

X(r,theta,phi) = r * sin(theta) * sin(phi)
Y(r,theta,phi) = r * sin(theta) * cos(phi)
Z(r,theta,phi) = r * cos(theta)

xs = [X(1, theta, phi) for theta in thetas, phi in phis]
ys = [Y(1, theta, phi) for theta in thetas, phi in phis]
zs = [Z(1, theta, phi) for theta in thetas, phi in phis]

plot_surface(xs, ys, zs, alpha=1)

gcf()

##
using PyCall
plot = pyimport("matplotlib.pyplot")
mpl3 = pyimport("mpl_toolkits.mplot3d")

using PyPlot

x = -1.5:0.1:1.5
y = -1:0.1:1
z(x,y) = 3.434

plot_surface(x,y,z)

##

using PyPlot

x = -1.5:0.1:1.5
y = -1:0.1:1
f(x,y) = 3.434

z = f.(x', y)
fig = surf(x, y, z)

##
using WhiskerTypes
a = HydroWhisker(1, 2, 3, 4, 5, 6, 7)


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

include("model_NetworkSecurityGroupListResult.jl")
include("model_NetworkSecurityGroupPropertiesFormat.jl")
include("model_Resource.jl")
include("model_SecurityRuleListResult.jl")
include("model_SecurityRulePropertiesFormat.jl")
include("model_SubResource.jl")
include("model_NetworkInterface.jl")
include("model_NetworkSecurityGroup.jl")
include("model_SecurityRule.jl")
include("model_Subnet.jl")


text:
module adapters_postgresql_placeholder_test

using Test # @test
using Octo.Adapters.PostgreSQL # Repo from to_sql Schema.model Octo.PlaceHolder SELECT FROM WHERE

Repo.debug_sql()

Repo.connect(
    adapter = Octo.Adapters.PostgreSQL,
    dbname = "postgresqltest",
    user = "postgres",
)

struct User
end
Schema.model(User, table_name="users")

Repo.execute([DROP TABLE IF EXISTS User])
Repo.execute(Raw("""CREATE TABLE IF NOT EXISTS users (
                     ID SERIAL,
                     name VARCHAR(255),
                     salary FLOAT(8),
                     PRIMARY KEY (ID) )
                 """))

❔ = Octo.PlaceHolder

# https://en.wikipedia.org/wiki/SQL_injection

u = from(User)
userName = """ ' OR '1'='1 """ # Incorrectly filtered escape characters
df = Repo.query([SELECT * FROM u WHERE u.name == ❔ AND u.salary > ❔], [userName, 2000])
@test size(df) == (0,)

@test to_sql([Octo.PlaceHolder]) == "\$1"

Repo.disconnect()

end # module adapters_postgresql_placeholder_test


text:
using FFTW, Jets, JetPackTransforms, Test

@testset "fft, 1D, complex" begin
    n = 512
    m = rand(n) + im * rand(n)
    d = rand(n) + im * rand(n)
    A = JopFft(ComplexF64, n)

    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
    expected = fft(m) / sqrt(n)
    observed = A * m
    @test expected ≈ observed
    expected = bfft(d) / sqrt(n)
    observed = A' * d
    @test expected ≈ observed
end

@testset "fft, alternative constructor" begin
    R = JetSpace(Float64,512)
    A = JopFft(R)
    B = JopFft(Float64,512)
    m = rand(R)
    @test A*m ≈ B*m
end

@testset "fft, 1D, real" begin
    n = 512
    spDom = JetSpace(Float64, n)
    spRng = symspace(JopFft, Float64, 1, n)
    m = rand(spDom)
    d = rand(spRng)
    d[1] = real(d[1])
    d[length(spRng)] = d[1]
    A = JopFft(Float64, n)
    @test range(A) == spRng
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "fft, 2D" begin
    n1, n2 = 128, 256
    A = JopFft(ComplexF64, n1, n2)
    m = rand(domain(A))
    d = rand(range(A))
    #=
    On windows causes a segfault... presumably a bug in BLAS?
    =#
    if !Sys.iswindows()
        lhs, rhs = dot_product_test(A, m, d)
        @test lhs ≈ rhs
    end
end

@testset "fft, 2D, real->complex" begin
    n1, n2 = 128, 256
    spDom = JetSpace(Float64, n1, n2)
    spRng = symspace(JopFft, Float64, 1, n1, n2)
    m = rand(spDom)
    d = rand(spRng)
    A = JopFft(Float64,n1,n2)
    @test spRng == range(A)
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs

    m = rand(spDom)
    m_roundtrip = A' * (A * m)
    @test m ≈ m_roundtrip

    d = A * m
    d_check = rfft(m) / sqrt(length(m))
    @test parent(d) ≈ d_check

    d_copy = copy(d)
    mm = A' * d
    @test d_copy ≈ d
    mm_check = brfft(parent(d), n1) / sqrt(length(m))
    @test mm_check ≈ mm
    @test mm_check ≈ m
    @test mm ≈ m
end

@testset "fft, 1D transform of 2D array" begin
    n1, n2 = 128, 256
    A = JopFft(ComplexF64, n1, n2; dims=(1,))
    m = rand(domain(A))
    d = A*m
    d_expected = similar(d)
    for i2 = 1:n2
        d_expected[:,i2] = fft(vec(m[:,i2])) / sqrt(n1)
    end
    @test d ≈ d_expected
    a = A'*d
    a_expected = similar(a)
    for i2 = 1:n2
        a_expected[:,i2] = bfft(vec(d[:,i2])) / sqrt(n1)
    end
    @test a ≈ a_expected
    if !Sys.iswindows()
        lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
        @test lhs ≈ rhs
    end
end

@testset "fft, 1D transform of 2D array, real->complex" begin
    n1, n2 = 128, 256
    A = JopFft(Float64, n1, n2; dims=(1,))
    m = rand(domain(A))
    d = A*m
    d_expected = similar(parent(d))
    for i2 = 1:n2
        d_expected[:,i2] = rfft(vec(m[:,i2])) / sqrt(n1)
    end
    @test parent(d) ≈ d_expected
    a = A'*d
    a_expected = similar(a)
    for i2 = 1:n2
        a_expected[:,i2] = brfft(vec(parent(d)[:,i2]),n1) / sqrt(n1)
    end
    @test a ≈ a_expected
    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
    @test lhs ≈ rhs
end

@testset "fft, 1D transform of 2D array, 2nd dim" begin
    n1, n2 = 128, 256
    A = JopFft(ComplexF64, n1, n2; dims=(2,))
    m = rand(domain(A))
    d = A*m
    d_expected = similar(d)
    for i1 = 1:n1
        d_expected[i1,:] = fft(vec(m[i1,:])) / sqrt(n2)
    end
    @test d ≈ d_expected
    a = A'*d
    a_expected = similar(a)
    for i1 = 1:n1
        a_expected[i1,:] = bfft(vec(d_expected[i1,:])) / sqrt(n2)
    end
    @test a ≈ a_expected
    if !Sys.iswindows()
        lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
        @test lhs ≈ rhs
    end
end

@testset "fft, 1D transform of 2D array, 2nd dim, real->complex" begin
    n1, n2 = 128, 256
    A = JopFft(Float64, n1, n2; dims=(2,))
    m = rand(domain(A))
    d = A*m

    d_expected = similar(parent(d))
    for i1 = 1:n1
        d_expected[i1,:] = rfft(vec(m[i1,:])) / sqrt(n2)
    end
    @test parent(d) ≈ d_expected
    a = A'*d
    a_expected = similar(a)
    for i1 = 1:n1
        a_expected[i1,:] = brfft(vec(d_expected[i1,:]), n2) / sqrt(n2)
    end
    @test a ≈ a_expected
    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
    @test lhs ≈ rhs
end

@testset "fft, 2D of 3D array, real to complex along first two dims" begin
    n1,n2,n3 = 128,256,10
    A = JopFft(Float64,n1,n2,n3;dims=(1,2))
    m = rand(domain(A))
    d = A*m
    d_expected = similar(parent(d))
    for i3=1:n3
        d_expected[:,:,i3] = rfft(m[:,:,i3]) / sqrt(n1*n2)
    end
    @test parent(d) ≈ d_expected
    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
    @test lhs ≈ rhs
end


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

mutable struct IpTag <: SwaggerModel
    ipTagType::Any # spec type: Union{ Nothing, String } # spec name: ipTagType
    tag::Any # spec type: Union{ Nothing, String } # spec name: tag

    function IpTag(;ipTagType=nothing, tag=nothing)
        o = new()
        validate_property(IpTag, Symbol("ipTagType"), ipTagType)
        setfield!(o, Symbol("ipTagType"), ipTagType)
        validate_property(IpTag, Symbol("tag"), tag)
        setfield!(o, Symbol("tag"), tag)
        o
    end
end # type IpTag

const _property_map_IpTag = Dict{Symbol,Symbol}(Symbol("ipTagType")=>Symbol("ipTagType"), Symbol("tag")=>Symbol("tag"))
const _property_types_IpTag = Dict{Symbol,String}(Symbol("ipTagType")=>"String", Symbol("tag")=>"String")
Base.propertynames(::Type{ IpTag }) = collect(keys(_property_map_IpTag))
Swagger.property_type(::Type{ IpTag }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IpTag[name]))}
Swagger.field_name(::Type{ IpTag }, property_name::Symbol) =  _property_map_IpTag[property_name]

function check_required(o::IpTag)
    true
end

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


text:
using JuMP
using CutPruners
using StructJuMP
using StructDualDynProg
using Compat, Compat.LinearAlgebra, Test
using StochOptInterface
const SOI = StochOptInterface

# solver independent tests
include("comp.jl")
include("stopcrit.jl")
include("paths.jl")

# load a solver
include("solvers.jl")

# solver dependent tests
include("optimize_stock.jl")
include("hydro_thermal_scheduling.jl")
include("prob5.2.jl")


text:
using Revise
using CUDA
using StaticArrays
using StatsBase: sample
using BenchmarkTools

nbins = 10
items = Int32(1e6)
hist = zeros(Float32, nbins)
x = ones(Float32, items)
idx = Int64.(rand(1:nbins, items))

hist_gpu = CuArray(hist)
x_gpu = CuArray(x)
idx_gpu = CuArray(idx)

hist .- Array(hist_gpu)
sum(hist) - sum(Array(hist_gpu))

# base kernel
function kernel!(h::CuDeviceVector{T}, x::CuDeviceVector{T}, xid::CuDeviceVector{S}) where {T,S}

    nbins = size(h, 1)
    it = threadIdx().x
    ib = blockIdx().x
    id = blockDim().x

    shared = @cuDynamicSharedMem(T, nbins)
    fill!(shared, 0)
    fill!(h, 0)
    sync_threads()

    i_tot = size(x, 1)
    iter = 0
    while iter * id < i_tot
        i = it + id * iter
        if i <= size(xid, 1)
            @inbounds k = Base._to_linear_index(h, xid[i])
            @inbounds CUDA.atomic_add!(pointer(shared, k), x[i])
        end
        iter += 1
    end
    sync_threads()
    # loop to cover cases where nbins > nthreads
    for i in 1:(nbins - 1) ÷ id + 1
        bin_id = it + id * (i - 1)
        if bin_id <= nbins
            @inbounds CUDA.atomic_add!(pointer(h, bin_id), shared[bin_id])
        end
    end
    return nothing
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist!(h::AbstractVector{T}, x::AbstractVector{T}, xid::AbstractVector{S}; MAX_THREADS=256) where {T,S}
    threads = min(MAX_THREADS, size(xid, 1))
    @cuda blocks = 1 threads = threads shmem = sizeof(T) * size(h,1) kernel!(h, x, xid)
    return
end

@btime CUDA.@sync hist!($hist_gpu, $x_gpu, $idx_gpu, MAX_THREADS=1024)


text:
const IIF = IncrementalInference
const KDE = KernelDensityEstimate

const NothingUnion{T} = Union{Nothing, T}

const CTs = CoordinateTransformations
const TU = TransformUtils

vectoarr2(v) = reshape(v, length(v),1)


text:
using CartesianGeneticProgramming
using PyCall
using Cambrian
using ArgParse
import Cambrian.mutate
import Random

"""
Demonstrates CGP on Discrete gym problems like MountainCar, CartPole, Acrobot

Requires the Python package gym, which can be installed with Conda.jl:
> Conda.add("gym")

Pybullet environments can be used if the "--pybullet" flag is provided. These
require separate installation:
> Conda.add("pybullet")
"""

s = ArgParseSettings()
@add_arg_table s begin
    "--cfg"
    help = "configuration script"
    default = "cfg/gym.yaml"
    "--env"
    help = "environment"
    default = "CartPole-v1"
    "--seed"
    help = "random seed"
    arg_type = Int
    default = 0
    "--pybullet"
    help = "use a pybullet env"
    action = :store_true
    "--ind"
    help = "individual for evaluation"
    arg_type = String
    default = ""

end
args = parse_args(ARGS, s)

function play_env(ind::CGPInd, env_name::String; seed::Int=0, render::Bool=false)
    env = gym.make(env_name)
    env.seed(seed)
    if render
        env.render(mode="human")
    end
    obs = env.reset()
    total_reward = 0.0
    done = false
    max_obs = Float64(max(-minimum(env.observation_space.low),
                          maximum(env.observation_space.high)))
    while ~done
        action = process(ind, obs ./ max_obs)
        if hasproperty(env.action_space, :n)
            # discrete env, use argmax (python 0-based indexing)
            action = argmax(action) - 1
        else
            # continuous env, normalize outputs
            h = env.action_space.high
            l = env.action_space.low
            action = ((action .+ 1) ./ 2) .* (h - l) .+ l
        end
        obs, reward, done, _ = env.step(action)
        if render
            env.render(mode="human")
        end
        total_reward += reward
    end
    env.close()
    [total_reward]
end

if args["pybullet"]
    pybullet_envs = pyimport("pybullet_envs")
end
gym = pyimport("gym")
env = gym.make(args["env"])
obs = env.reset()
n_in = length(env.observation_space.sample())
n_out = length(env.action_space.sample())
if hasproperty(env.action_space, :n)
    n_out = env.action_space.n
end
env.close()

cfg = get_config(args["cfg"]; n_in=n_in, n_out=n_out)
Random.seed!(args["seed"])

if length(args["ind"]) > 0
    ind = CGPInd(cfg, read(args["ind"], String))
    reward = play_env(ind, args["env"]; seed=args["seed"], render=true)
    println(reward)
else
    mutate(i::CGPInd) = goldman_mutate(cfg, i)
    fit(i::CGPInd) = play_env(i, args["env"])
    e = CGPEvolution(cfg, fit)
    run!(e)
end


text:
# lhs = left hand side; rhs = right hand side
# src = source, dst = destination
# order is: lhs, rhs, src, dst
function Arrow(lhs, rhs, src, dst; edge_fun = DirectedEdge, cross = false)
    if cross
        edges = edge_fun(src, dst)
    else
        edges = edge_fun.(src, dst)
    end
    edges = isa(edges, VecOrMat) ? edges : [edges]
    edges = isa(edges, Matrix) ? vec(edges) : edges
    edges = isa(lhs, Arrow) ? vcat(lhs.edges..., edges...) : edges
    edges = isa(rhs, Arrow) ? vcat(edges..., rhs.edges...) : edges
    Arrow(edges, lhs, rhs)
end

→(lhs, rhs) = Arrow(lhs, rhs, keep(lhs, Right, Node), keep(rhs, Left, AbstractNode))
⇒(lhs, rhs) = Arrow(lhs, rhs, keep(lhs, Right, Node), keep(rhs, Left, AbstractNode); cross = true)
←(lhs, rhs) = Arrow(lhs, rhs, keep(rhs, Left, Node), keep(lhs, Right, AbstractNode))
⇐(lhs, rhs) = Arrow(lhs, rhs, keep(rhs, Left, Node), keep(lhs, Right, AbstractNode); cross = true)
↔(lhs, rhs) = Arrow(lhs, rhs, keep(lhs, Right, Node), keep(rhs, Left, AbstractNode); edge_fun = UndirectedEdge)
⇔(lhs, rhs) = Arrow(lhs, rhs, keep(lhs, Right, Node), keep(rhs, Left, AbstractNode); edge_fun = UndirectedEdge, cross = true)

# if either lhs or rhs are a single node, broadcast arrow (barrow) behave like cross arrow (carrow)
barrows = [:→, :←, :↔]
carrows = [:⇒, :⇐, :⇔]
for (barrow, carrow) in zip(barrows,  carrows)
    @eval $barrow(lhs::Union{AbstractNode, Symbol}, rhs::Union{AbstractNode, Symbol}) = $carrow(lhs, rhs)
    @eval $barrow(lhs::Union{AbstractNode, Symbol}, rhs) = $carrow(lhs, rhs)
    @eval $barrow(lhs, rhs::Union{AbstractNode, Symbol}) = $carrow(lhs, rhs)
end


text:
abstract type FloatInterval{T<:AbstractFloat} <: Uniform{T} end

abstract type CloseOpen{ T<:AbstractFloat} <: FloatInterval{T} end
abstract type OpenClose{ T<:AbstractFloat} <: FloatInterval{T} end
abstract type CloseClose{T<:AbstractFloat} <: FloatInterval{T} end
abstract type OpenOpen{  T<:AbstractFloat} <: FloatInterval{T} end

struct CloseOpen12{T<:AbstractFloat} <: CloseOpen{T} end # interval [1,2)

struct CloseOpen01{ T<:AbstractFloat} <: CloseOpen{T}  end # interval [0,1)
struct OpenClose01{ T<:AbstractFloat} <: OpenClose{T}  end # interval (0,1]
struct CloseClose01{T<:AbstractFloat} <: CloseClose{T} end # interval [0,1]
struct OpenOpen01{  T<:AbstractFloat} <: OpenOpen{T}   end # interval (0,1)

struct CloseOpenAB{T<:AbstractFloat} <: CloseOpen{T} # interval [a,b)
    a::T
    b::T

    CloseOpenAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct OpenCloseAB{T<:AbstractFloat} <: OpenClose{T} # interval (a,b]
    a::T
    b::T

    OpenCloseAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct CloseCloseAB{T<:AbstractFloat} <: CloseClose{T} # interval [a,b]
    a::T
    b::T

    CloseCloseAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct OpenOpenAB{T<:AbstractFloat} <: OpenOpen{T} # interval (a,b)
    a::T
    b::T

    OpenOpenAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

check_interval(a, b) = a >= b && throw(ArgumentError("invalid interval specification"))

const FloatInterval_64 = FloatInterval{Float64}
const CloseOpen01_64   = CloseOpen01{Float64}
const CloseOpen12_64   = CloseOpen12{Float64}

CloseOpen01(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen01{T}()
CloseOpen12(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen12{T}()

CloseOpen(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen01{T}()
CloseOpen(a::T, b::T) where {T<:AbstractFloat} = CloseOpenAB{T}(a, b)

OpenClose(::Type{T}=Float64) where {T<:AbstractFloat} = OpenClose01{T}()
OpenClose(a::T, b::T) where {T<:AbstractFloat} = OpenCloseAB{T}(a, b)

CloseClose(::Type{T}=Float64) where {T<:AbstractFloat} = CloseClose01{T}()
CloseClose(a::T, b::T) where {T<:AbstractFloat} = CloseCloseAB{T}(a, b)

OpenOpen(::Type{T}=Float64) where {T<:AbstractFloat} = OpenOpen01{T}()
OpenOpen(a::T, b::T) where {T<:AbstractFloat} = OpenOpenAB{T}(a, b)

# convenience functions

CloseOpen(a, b) = CloseOpen(promote(a, b)...)
CloseOpen(a::T, b::T) where {T} = CloseOpen(AbstractFloat(a), AbstractFloat(b))

OpenClose(a, b) = OpenClose(promote(a, b)...)
OpenClose(a::T, b::T) where {T} = OpenClose(AbstractFloat(a), AbstractFloat(b))

CloseClose(a, b) = CloseClose(promote(a, b)...)
CloseClose(a::T, b::T) where {T} = CloseClose(AbstractFloat(a), AbstractFloat(b))

OpenOpen(a, b) = OpenOpen(promote(a, b)...)
OpenOpen(a::T, b::T) where {T} = OpenOpen(AbstractFloat(a), AbstractFloat(b))

## sampling

### fall-back on Random definitions

for CO in (:CloseOpen01, :CloseOpen12)
    @eval begin
        Sampler(RNG::Type{<:AbstractRNG}, ::$CO{T}, n::Repetition) where {T} =
            Sampler(RNG, Random.$CO{T}(), n)

        Sampler(::Type{<:AbstractRNG}, ::$CO{BigFloat}, ::Repetition) =
            Random.SamplerBigFloat{Random.$CO{BigFloat}}(precision(BigFloat))
    end
end

### new intervals 01

# TODO: optimize for BigFloat

for CO = (:OpenClose01, :OpenOpen01, :CloseClose01)
    @eval Sampler(RNG::Type{<:AbstractRNG}, I::$CO{T}, n::Repetition) where {T} =
              SamplerSimple(I, CloseOpen01(T))
end

rand(r::AbstractRNG, sp::SamplerSimple{OpenClose01{T}}) where {T} =
    one(T) - rand(r, sp.data)

rand(r::AbstractRNG, sp::SamplerSimple{OpenOpen01{T}}) where {T} =
    while true
        x = rand(r, sp.data)
        x != zero(T) && return x
    end

# optimizations (TODO: optimize for BigFloat too)

rand(r::AbstractRNG, sp::SamplerSimple{OpenOpen01{Float64}}) =
    reinterpret(Float64, reinterpret(UInt64, rand(r, sp.data)) | 0x0000000000000001)

rand(r::AbstractRNG, sp::SamplerSimple{OpenOpen01{Float32}}) =
    reinterpret(Float32, reinterpret(UInt32, rand(r, sp.data)) | 0x00000001)

rand(r::AbstractRNG, sp::SamplerSimple{OpenOpen01{Float16}}) =
    reinterpret(Float16, reinterpret(UInt16, rand(r, sp.data)) | 0x0001)

# prevfloat(T(2)) - 1 for IEEEFloat
upper01(::Type{Float64}) = 0.9999999999999998
upper01(::Type{Float32}) = 0.9999999f0
upper01(::Type{Float16}) = Float16(0.999)
upper01(::Type{BigFloat}) = prevfloat(one(BigFloat))

rand(r::AbstractRNG, sp::SamplerSimple{CloseClose01{T}}) where {T} =
    rand(r, sp.data) / upper01(T)

### CloseOpenAB

for (CO, CO01) = (CloseOpenAB => CloseOpen01,
                  OpenCloseAB => OpenClose01,
                  CloseCloseAB => CloseClose01,
                  OpenOpenAB => OpenOpen01)

    @eval Sampler(RNG::Type{<:AbstractRNG}, d::$CO{T}, n::Repetition) where {T} =
        SamplerTag{$CO{T}}((a=d.a, d=d.b - d.a, sp=Sampler(RNG, $CO01{T}(), n)))

    @eval rand(rng::AbstractRNG, sp::SamplerTag{$CO{T}}) where {T} =
        sp.data.a + sp.data.d  * rand(rng, sp.data.sp)
end


text:
using MAT
using MNIST

dataLocation = "data/mnist67.scale.1k.mat"

function update_labels(y)
    for i in 1:size(y,1)
        if y[i] % 2 == 0
             y[i] = -1
        else
             y[i] = 1
        end
    end
end

function getTestData()
    println("Loading testdata.")
    X, y = testdata()
    X = X'
    # normalize data
    X = X / 255
    update_labels(y)
    println("Loaded testdata.")
    return (X,y)
end

function getOldData()
    println("Loading dataset.")
    file = matopen(dataLocation)
    tmp = read(file, "X")

    X = full(tmp)
    y = read(file,"y")
    println("Loaded dataset.")
    return (X,y)
end

function getFullData()
    println("Loading dataset.")
    X, y = traindata()
    X = X'
    # normalize data
    X = X / 255
    update_labels(y)
    println("Loaded dataset.")
    return (X,y)
end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
# ------------------------------------------------------------
# UsesRegularization: Algorithm uses some regularization method
# ------------------------------------------------------------
@define_trait UsesRegularization = begin
    IsRegularized  # Algorithm uses the regularized decomposition method of Ruszczyński
    HasTrustRegion # Algorithm uses the trust-region method of Linderoth/Wright
    HasLevels      # Algorithm uses the level set method of Lemarcheral
end

@define_traitfn UsesRegularization init_solver!(lshaped::AbstractLShapedSolver) = begin
    function init_solver!(lshaped::AbstractLShapedSolver,!UsesRegularization)
        nothing
    end
end

@define_traitfn UsesRegularization log_regularization!(lshaped::AbstractLShapedSolver) = begin
    function log_regularization!(lshaped::AbstractLShapedSolver,!UsesRegularization)
        nothing
    end
end

@define_traitfn UsesRegularization log_regularization!(lshaped::AbstractLShapedSolver,t::Integer) = begin
    function log_regularization!(lshaped::AbstractLShapedSolver,t::Integer,!UsesRegularization)
        nothing
    end
end

@define_traitfn UsesRegularization take_step!(lshaped::AbstractLShapedSolver) = begin
    function take_step!(lshaped::AbstractLShapedSolver,!UsesRegularization)
        nothing
    end
end

@define_traitfn UsesRegularization process_cut!(lshaped::AbstractLShapedSolver,cut::HyperPlane) = begin
    function process_cut!(lshaped::AbstractLShapedSolver,cut::HyperPlane,!UsesRegularization)
        nothing
    end

    function process_cut!(lshaped::AbstractLShapedSolver,cut::HyperPlane,UsesRegularization)
        nothing
    end
end

@implement_traitfn function process_cut!(lshaped::AbstractLShapedSolver,cut::HyperPlane,IsRegularized)
    push!(lshaped.committee,cut)
    nothing
end

@define_traitfn UsesRegularization remove_inactive!(lshaped::AbstractLShapedSolver) = begin
    function remove_inactive!(lshaped::AbstractLShapedSolver,UsesRegularization)
        inactive = find(c->!active(lshaped,c),lshaped.committee)
        diff = length(lshaped.committee) - length(lshaped.structuredmodel.linconstr) - lshaped.nscenarios
        if isempty(inactive) || diff <= 0
            return false
        end
        if diff <= length(inactive)
            inactive = inactive[1:diff]
        end
        append!(lshaped.inactive,lshaped.committee[inactive])
        deleteat!(lshaped.committee,inactive)
        delconstrs!(lshaped.mastersolver.lqmodel,inactive)
        return true
    end
end

@define_traitfn UsesRegularization queueViolated!(lshaped::AbstractLShapedSolver) = begin
    function queueViolated!(lshaped::AbstractLShapedSolver,UsesRegularization)
        violating = find(c->violated(lshaped,c),lshaped.inactive)
        if isempty(violating)
            return false
        end
        gaps = map(c->gap(lshaped,c),lshaped.inactive[violating])
        for (c,g) in zip(lshaped.inactive[violating],gaps)
            enqueue!(lshaped.violating,c,g)
        end
        deleteat!(lshaped.inactive,violating)
        return true
    end
end

@define_traitfn UsesRegularization project!(lshaped::AbstractLShapedSolver) = begin
    function project!(lshaped::AbstractLShapedSolver,!UsesRegularization)
        nothing
    end

    function project!(lshaped::AbstractLShapedSolver,UsesRegularization)
        nothing
    end
end

# Is Regularized
# ------------------------------------------------------------
@define_traitfn IsRegularized update_objective!(lshaped::AbstractLShapedSolver)

@implement_traitfn function init_solver!(lshaped::AbstractLShapedSolver,IsRegularized)
    if lshaped.parameters.autotune
        σ̅ = norm(lshaped.x)
        σ̲ = 1.0
        σ = mean([σ̅,σ̲])
        @pack lshaped.parameters = σ,σ̅,σ̲
    end
    lshaped.solverdata.σ = lshaped.parameters.σ
    push!(lshaped.σ_history,lshaped.solverdata.σ)

    update_objective!(lshaped)
end

@implement_traitfn function log_regularization!(lshaped::AbstractLShapedSolver,IsRegularized)
    @unpack Q̃,σ = lshaped.solverdata
    push!(lshaped.Q̃_history,Q̃)
    push!(lshaped.σ_history,σ)
end

@implement_traitfn function log_regularization!(lshaped::AbstractLShapedSolver,t::Integer,IsRegularized)
    @unpack Q̃,σ = lshaped.solverdata
    lshaped.Q̃_history[t] = Q̃
    lshaped.σ_history[t] = σ
end

@implement_traitfn function take_step!(lshaped::AbstractLShapedSolver,IsRegularized)
    @unpack Q,Q̃,θ = lshaped.solverdata
    @unpack τ,σ̅,σ̲ = lshaped.parameters
    if Q + τ*(1+abs(Q)) <= Q̃
        lshaped.ξ[:] = lshaped.x[:]
        lshaped.solverdata.Q̃ = Q
        if abs(Q - Q̃) <= 0.5*(Q̃-θ)
            # Enlarge the trust-region radius
            lshaped.solverdata.σ = min(σ̅,lshaped.solverdata.σ*2)
            update_objective!(lshaped)
        end
        lshaped.solverdata.major_iterations += 1
    else
        lshaped.solverdata.σ *= 0.5
        lshaped.solverdata.σ = max(σ̲,lshaped.solverdata.σ)
        update_objective!(lshaped)
        lshaped.solverdata.minor_iterations += 1
    end
    nothing
end

@implement_traitfn function update_objective!(lshaped::AbstractLShapedSolver,IsRegularized)
    # Linear regularizer penalty
    c = copy(lshaped.c)
    c -= (1/lshaped.solverdata.σ)*lshaped.ξ
    append!(c,fill(1.0,lshaped.nscenarios))
    setobj!(lshaped.mastersolver.lqmodel,c)

    # Quadratic regularizer penalty
    qidx = collect(1:length(lshaped.ξ)+lshaped.nscenarios)
    qval = fill(1/lshaped.solverdata.σ,length(lshaped.ξ))
    append!(qval,zeros(lshaped.nscenarios))
    if applicable(setquadobj!,lshaped.mastersolver.lqmodel,qidx,qidx,qval)
        setquadobj!(lshaped.mastersolver.lqmodel,qidx,qidx,qval)
    else
        error("The regularized decomposition algorithm requires a solver that handles quadratic objectives")
    end
end

# HasTrustRegion
# ------------------------------------------------------------
@define_traitfn HasTrustRegion set_trustregion!(lshaped::AbstractLShapedSolver)
@define_traitfn HasTrustRegion enlarge_trustregion!(lshaped::AbstractLShapedSolver)
@define_traitfn HasTrustRegion reduce_trustregion!(lshaped::AbstractLShapedSolver)

@implement_traitfn function init_solver!(lshaped::AbstractLShapedSolver,HasTrustRegion)
    if lshaped.parameters.autotune
        Δ = max(1.0,0.01*norm(lshaped.x,Inf))
        Δ̅ = 0.05*norm(lshaped.x,Inf)
        @pack lshaped.parameters = Δ,Δ̅
    end
    lshaped.solverdata.Δ = lshaped.parameters.Δ
    push!(lshaped.Δ_history,lshaped.solverdata.Δ)

    set_trustregion!(lshaped)
end

@implement_traitfn function log_regularization!(lshaped::AbstractLShapedSolver,HasTrustRegion)
    @unpack Q̃,Δ = lshaped.solverdata
    push!(lshaped.Q̃_history,Q̃)
    push!(lshaped.Δ_history,Δ)
end

@implement_traitfn function log_regularization!(lshaped::AbstractLShapedSolver,t::Integer,HasTrustRegion)
    @unpack Q̃,Δ = lshaped.solverdata
    lshaped.Q̃_history[t] = Q̃
    lshaped.Δ_history[t] = Δ
end

@implement_traitfn function take_step!(lshaped::AbstractLShapedSolver,HasTrustRegion)
    @unpack Q,Q̃,θ = lshaped.solverdata
    @unpack γ = lshaped.parameters
    if Q <= Q̃ - γ*abs(Q̃-θ)
        lshaped.solverdata.cΔ = 0
        lshaped.ξ[:] = lshaped.x[:]
        lshaped.solverdata.Q̃ = Q
        enlarge_trustregion!(lshaped)
        lshaped.solverdata.major_iterations += 1
    else
        reduce_trustregion!(lshaped)
        lshaped.solverdata.minor_iterations += 1
    end
    nothing
end

@implement_traitfn function set_trustregion!(lshaped::AbstractLShapedSolver,HasTrustRegion)
    l = max.(lshaped.structuredmodel.colLower, lshaped.ξ-lshaped.solverdata.Δ)
    append!(l,fill(-Inf,lshaped.nscenarios))
    u = min.(lshaped.structuredmodel.colUpper, lshaped.ξ+lshaped.solverdata.Δ)
    append!(u,fill(Inf,lshaped.nscenarios))
    setvarLB!(lshaped.mastersolver.lqmodel,l)
    setvarUB!(lshaped.mastersolver.lqmodel,u)
end

@implement_traitfn function enlarge_trustregion!(lshaped::AbstractLShapedSolver,HasTrustRegion)
    @unpack Q,Q̃,θ = lshaped.solverdata
    @unpack τ,Δ̅ = lshaped.parameters
    if abs(Q - Q̃) <= 0.5*(Q̃-θ) && norm(lshaped.ξ-lshaped.x,Inf) - lshaped.solverdata.Δ <= τ
        # Enlarge the trust-region radius
        lshaped.solverdata.Δ = min(Δ̅,2*lshaped.solverdata.Δ)
        set_trustregion!(lshaped)
        return true
    else
        return false
    end
end

@implement_traitfn function reduce_trustregion!(lshaped::AbstractLShapedSolver,HasTrustRegion)
    @unpack Q,Q̃,θ = lshaped.solverdata
    ρ = min(1,lshaped.solverdata.Δ)*(Q-Q̃)/(Q̃-θ)
    if ρ > 0
        lshaped.solverdata.cΔ += 1
    end
    if ρ > 3 || (lshaped.solverdata.cΔ >= 3 && 1 < ρ <= 3)
        # Reduce the trust-region radius
        lshaped.solverdata.cΔ = 0
        lshaped.solverdata.Δ = (1/min(ρ,4))*lshaped.solverdata.Δ
        set_trustregion!(lshaped)
        return true
    else
        return false
    end
end

# Has Levels
# ------------------------------------------------------------

@implement_traitfn function init_solver!(lshaped::AbstractLShapedSolver,HasLevels)
    # θs
    for i = 1:lshaped.nscenarios
        addvar!(lshaped.projectionsolver.lqmodel,-Inf,Inf,1.0)
    end
end

@implement_traitfn function log_regularization!(lshaped::AbstractLShapedSolver,HasLevels)
    @unpack Q̃ = lshaped.solverdata
    push!(lshaped.Q̃_history,Q̃)
end

@implement_traitfn function log_regularization!(lshaped::AbstractLShapedSolver,t::Integer,HasLevels)
    @unpack Q̃ = lshaped.solverdata
    lshaped.Q̃_history[t] = Q̃
end

@implement_traitfn function take_step!(lshaped::AbstractLShapedSolver,HasLevels)
    @unpack Q,Q̃ = lshaped.solverdata
    @unpack τ = lshaped.parameters
    if Q + τ*(1+abs(Q)) <= Q̃
        lshaped.solverdata.Q̃ = Q
    end
    nothing
end

@implement_traitfn function process_cut!(lshaped::AbstractLShapedSolver,cut::HyperPlane{OptimalityCut},HasLevels)
    #addconstr!(lshaped.projectionsolver.lqmodel,lowlevel(cut)...) TODO: Rewrite with MathOptInterface
end

@implement_traitfn function project!(lshaped::AbstractLShapedSolver,HasLevels)
    @unpack θ,Q̃ = lshaped.solverdata
    @unpack λ = lshaped.parameters
    # Copy current master problem (TODO: Rewrite with MathOptInterface)
    lshaped.projectionsolver.lqmodel = copy(lshaped.mastersolver.lqmodel)
    # Update level
    c = sparse(getobj(lshaped.projectionsolver.lqmodel))
    L = (1-λ)*Q̃ + λ*θ
    addconstr!(lshaped.projectionsolver.lqmodel,c.nzind,c.nzval,-Inf,L)
e
    # Update regularizer
    q = -copy(lshaped.ξ)
    append!(q,zeros(lshaped.nscenarios))
    setobj!(lshaped.projectionsolver.lqmodel,q)
    # Quadratic regularizer penalty
    qidx = collect(1:length(lshaped.ξ)+lshaped.nscenarios)
    qval = ones(length(lshaped.ξ))
    append!(qval,zeros(lshaped.nscenarios))
    if applicable(setquadobj!,lshaped.projectionsolver.lqmodel,qidx,qidx,qval)
        setquadobj!(lshaped.projectionsolver.lqmodel,qidx,qidx,qval)
    else
        error("The level set algorithm requires a solver that handles quadratic objectives")
    end
    # Solve projection problem
    lshaped.projectionsolver(lshaped.mastervector)
    if status(lshaped.projectionsolver) == :Infeasible
        error("Projection problem is infeasible, aborting procedure.")
    end
    # Update master solution
    ncols = lshaped.structuredmodel.numCols
    x = getsolution(lshaped.projectionsolver)
    lshaped.x[1:ncols] = x[1:ncols]
    lshaped.ξ[:] = lshaped.x[:]
end


text:
using Documenter
using PeaceVote
makedocs(sitename="PeaceVote.jl",pages = ["index.md"])

deploydocs(
     repo = "github.com/PeaceFounder/PeaceVote.jl.git",
 )


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Test SparseArray" begin
    n_row = 100
    n_col = 110

    # Initial a SparseArray with size(10, 11)
    T_e = Int64
    T = DS.SparseArray{T_e}
    sa = T(n_row, n_col)
    @test sa isa T
    @test size(sa) == (n_row, n_col)
    @test length(sa) == 0
    @test eltype(sa) == T_e

    # SetIndex test: SetIndex data into SparseArray
    sa[5, 4] = 54
    @test sa[5, 4] == 54
    @test length(sa) == 1

    # SetIndex test: SetIndex data into SparseArray at same place
    sa[5, 4] = 5454
    @test sa[5, 4] == 5454
    @test length(sa) == 1

    # SetIndex test: SetIndex 0 into SparseArray and length shouldn't increase
    sa[6, 5] = 0
    @test sa[6, 5] == 0
    @test length(sa) == 1

    # SetIndex test: SetIndex many data into SparseArray
    sa = T(n_row, n_col)
    len = 0
    passed_check_value = []
    passed_check_length = []
    for i = 1:n_row
        for j = 1:n_col
            len += 1
            sa[i, j] = len
            push!(passed_check_value, sa[i, j] == len)
            push!(passed_check_length, length(sa) == len)
        end
    end
    @test all(passed_check_value)
    @test all(passed_check_length)

    # SetIndex test: SetIndex many 0 into SparseArray
    len = length(sa)
    passed_check_length = []
    for i = 2:n_row
        for j = 2:n_col
            sa[i, j] = 0
            len -= 1
            push!(passed_check_length, length(sa) == len)
        end
    end
    @test all(passed_check_length)

    new_sa = T(10, 5)
    new_sa[3, 3] = 53
    @test repr(new_sa) == "SparseArray{Int64}(10, 5)[\n\t[3, 3]: 53\n]"
end


text:
Base.@kwdef mutable struct NodeOptions
    rate::Float64 = 100
end

Base.@kwdef mutable struct NodeFlags
    did_error::Threads.Atomic{Bool} = Threads.Atomic{Bool}(false)
    is_running::Threads.Atomic{Bool} = Threads.Atomic{Bool}(false)
    should_finish::Threads.Atomic{Bool} = Threads.Atomic{Bool}(false)
end

"""
    NodeIO

Describes the input/output mechanisms for the node. Each node should store this type
internally and add the necessary I/O mechanisms inside of the `setupIO!(::NodeIO, ::Node)`
method.

I/O mechanisms are added to a `NodeIO` object via [`add_publisher!`](@ref) and
[`add_subscriber!`](@ref).
"""
struct NodeIO
    ctx::Union{Nothing,ZMQ.Context}
    pubs::Vector{PublishedMessage}
    subs::Vector{SubscribedMessage}
    opts::NodeOptions
    flags::NodeFlags

    function NodeIO(ctx::ZMQ.Context = ZMQ.context(); opts...)
        new(
            ctx,
            PublishedMessage[],
            SubscribedMessage[],
            NodeOptions(; opts...),
            NodeFlags(),
        )
    end
end

"""
    add_publisher!(nodeIO, msg, args...)

Adds / registers a publisher to `nodeIO`. This method should only be called once
per unique message, across all nodes in the network, since each message should only
ever have one publisher. The `msg` can be any `ProtoBuf.ProtoType` message (usually
generated using `ProtoBuf.protoc`). Since this is stored as an abstract `ProtoBuf.ProtoType`
type internally, the user should store the original type inside the node.  The remaining
arguments are passed directly to the constructor for [`Publisher`](@ref).

This function adds a new [`PublishedMessage`](@ref) to `nodeIO.pubs`. During the `compute`
    method, the user should modify the original concrete `msg` stored in the node. The
    data can then be published by calling `publish` on the corresponding `PublishedMessage`.

# Example
Inside of the node constructor:

    ...
    test_msg = TestMsg(x = 1, y = 2, z= 3)
    ...

Inside of `setupIO!`:

    ...
    ctx = ZMQ.Context()
    ipaddr = ip"127.0.0.1"
    port = 5001
    add_publisher!(nodeIO, node.test_msg, ctx, ipaddr, port, name="TestMsg_publisher")
    ...

Inside of `compute`:

    ...
    node.test_msg.x = 1  # modify the message as desired
    publish(getIO(node).pubs[1])  # or whichever is the correct index
    ...

"""
function add_publisher!(nodeio::NodeIO, msg::MercuryMessage, pub::Publisher)
    push!(nodeio.pubs, PublishedMessage(msg, pub))
end

"""
    add_subscriber!(nodeIO, msg, args...)

Adds / registers a subscriber to `nodeIO`. The `msg` can be any
`ProtoBuf.ProtoType` message (usually generated using `ProtoBuf.protoc`). Since
this is stored as an abstract `ProtoBuf.ProtoType` type internally, the user
should store the original type inside the node.  The remaining arguments are
passed directly to the constructor for [`Subscriber`](@ref).

This function adds a new [`SubscribedMessage`](@ref) to `nodeIO.subs`. A
separate asynchronous task is created for each subscriber when the node is
launched.  During the `compute` method, the user can access the latest data by
reading from the message stored in their node. To avoid data races and minimize
synchronization, it's usually best practice to obtain the lock on the message
(stored in `SubscribedMessage`) and copy the data to a local variable (likely
also stored in the node) that can be used by the rest of the `compute` method
without worrying about the data being overwritted by the ongoing subscriber
task.

# Example
In the node constructor:

    ...
    test_msg = TestMessage(x = 0, y = 0, z = 0)
    ...

In `setupIO!`:

    ...
    ctx = ZMQ.Context()
    ipaddr = ip"127.0.0.1"
    port = 5001
    add_subscriber!(nodeIO, node.test_msg, ctx, ipaddr, port, name = "TestMsg_subscriber")
    ...

In `compute`:

    ...
    testmsg = getIO(node).subs[1]  # or whichever is the correct index
    lock(testmsg.lock) do
        node.local_test_msg = node.test_msg  # or convert to a different type
    end
    # use node.local_test_msg in the rest of the code
    ...
"""
function add_subscriber!(nodeio::NodeIO, msg::MercuryMessage, sub::Subscriber)
    push!(nodeio.subs, SubscribedMessage(msg, sub))
end

"""
    Node

A independent process that communicates with other processes via pub/sub ZMQ channels.
The process is assumed to run indefinately.

# Defining a new Node
Each node should contain a `NodeData` element, which stores a list of the publishers
and subscribers and some other associated data.

The publisher and subscribers for the node should be "registered" with the `NodeData`
using the `add_publisher!` and `add_subscriber!` methods. This allows the subscribers to
be automatically launched as separate tasks when launching the nodes.

The constructor for the node should initialize any variables and register the needed
publishers and subscribers with `NodeData`.

Each loop of the process will call the `compute` method once, which needs to be
implemented by the user. A lock for each subscriber task is created in `NodeData.sub_locks`.
It's recommended that the user obtains the lock and copies the data into a local variable
for internal use by the `compute` function.

# Launching the node
The blocking process that runs the node indefinately is called via `run(node)`. It's
recommended that this is launched asynchronously via

    node_task = @task run(node)
    schedule(node_task)
"""
abstract type Node end

##############################
# REQUIRED INTERFACE
##############################
compute(::Node)::Nothing =
    error("The `compute` method hasn't been implemented for your node yet!")

# NOTE: This method may be automatically defined using codegen and the TOML file
#       in the future
setupIO!(::Node, ::NodeIO) =
    error("The `setupIO` method hasn't been implemented for your node yet!")

##############################
# OPTIONAL INTERFACE
##############################
startup(::Node)::Nothing = nothing
finishup(::Node)::Nothing = nothing
getcontext(node::Node)::Union{Nothing,ZMQ.Context} = getIO(node).ctx
getIO(node::Node)::NodeIO = node.nodeio

function getname(::T) where {T<:Node}
    # Note that this only works well when each node is only instantiated once
    return string(T)
end

##############################
# PROVIDED INTERFACE (don't change)
##############################
getoptions(node::Node) = getIO(node).opts
getflags(node::Node) = getIO(node).flags
getrate(node::Node)::Float64 = getoptions(node).rate
function shouldnodefinish(node::Node)::Bool
    return getflags(node).should_finish[]
end

function stopnode(node::Node; timeout = 1.0)
    getflags(node).should_finish[] = true
    t_start = time()
    while (time() - t_start < timeout)
        yield()
        if !(getflags(node).is_running[])
            return true
        end
    end
    # If it doesn't stop in time, forcefully close the node
    @warn "Node timed out. Forcefully closing the node."
    closeall(node)
    return false
end

publishers(node::Node) = getIO(node).pubs
subscribers(node::Node) = getIO(node).subs

for pubsub in ((:publisher, :pubs), (:subscriber, :subs))
    @eval $(Symbol("get", pubsub[1]))(node::Node, index::Integer) =
        getIO(node).$(pubsub[2])[index]
    @eval function $(Symbol("get", pubsub[1]))(node::Node, name::String)
        index = findfirst(getIO(node).$(pubsub[2])) do msg
            getname(msg) == name
        end
        if !isnothing(index)
            return $(Symbol("get", pubsub[1]))(node, index)
        end
        return nothing
    end
    @eval $(Symbol("num", pubsub[1], "s"))(node::Node) = length(getIO(node).$(pubsub[2]))
end

"""
    numsubscribers(node)

Get the number of ZMQ subscribers attached to the node
"""
numsubscribers

"""
    numpublishers(node)

Get the number of ZMQ publishers attached to the node
"""
numpublishers

"""
    getsubscriber(node, index)
    getsubscriber(node, name)

Get a  [`SubscribedMessage`](@ref) attached to `node`, either by it's integer index or it's name.
"""
getsubscriber

"""
    getpublisher(node, index)
    getpublisher(node, name)

Get a  [`PublishedMessage`](@ref) attached to `node`, either by it's integer index or it's name.
"""
getpublisher

"""
    launch(node)

Run the main loop of the node indefinately. This method automatically sets up any necessary
subscriber tasks and then calls the `compute` method at a fixed rate.

This method should typically be wrapped in an `@async` or `@spawn` call.
"""
function launch(node::Node)
    try
        rate = getrate(node)
        lrl = LoopRateLimiter(rate)

        # Run any necessary startup
        startup(node)

        getflags(node).is_running[] = true

        @rate while !shouldnodefinish(node)
            # Check the subscribers for new messages
            poll_subscribers(node)

            # Check the subscribers for new messages
            compute(node)

            GC.gc(false)
            yield()
        end lrl
        @info "Closing node $(getname(node))"
        finishup(node)
        closeall(node)
    catch err
        if err isa InterruptException
            @info "Closing node $(getname(node)). Got Keyboard Interrupt."
        else
            @warn "Closing node $(getname(node)). Closed with error."
            @error "Node failed" exception = (err, catch_backtrace())

            Base.display_error(err)
            Base.show_exception_stack(err, stacktrace(catch_backtrace()))
            getflags(node).did_error[] = true
            rethrow(err)
        end
        finishup(node)
        closeall(node)
    end
    getflags(node).is_running[] = false
end

function poll_subscribers(node::Node)
    nodeio = getIO(node)

    for submsg in nodeio.subs
        receive(submsg)
    end
end

function closeall(node::Node)
    @info "Closing down all of $(getname(node))'s NodeIO connections"

    nodeio = getIO(node)
    # Close publishers and subscribers
    for submsg in nodeio.subs
        close(submsg.sub)
    end
    for pubmsg in nodeio.pubs
        close(pubmsg.pub)
    end

    return nothing
end

function node_sockets_are_open(node::Node)
    nodeio = getIO(node)
    return (
        all([isopen(submsg.sub) for submsg in nodeio.subs]) && all([isopen(pubmsg.pub) for pubmsg in nodeio.pubs])
    )
end

#! format: off
function printstatus(node::Node)
    is_running = getflags(node).is_running[]
    printstyled("Node name: ", bold=true); println(getname(node))
    printstyled("  Is running? ", bold=true); println(is_running)
    if !is_running
        printstyled("  Did error? ", bold=true); println(getflags(node).did_error[])
    end
    if (numpublishers(node) > 0)
        printstyled("  Publishers:\n", bold=true)
        for pub in publishers(node)
            printstatus(pub, indent=4)
        end
    end
    if (numsubscribers(node) > 0)
        printstyled("  Subscribers:\n", bold=true)
        for sub in subscribers(node)
            printstatus(sub, indent=4)
        end
    end
end
#! format: on


text:
using Latexify
using LaTeXStrings

test_array = ["x/y * d" :x ; :( (t_sub_sub - x)^(2*p) ) 3//4 ]
@test latexinline.(test_array) == [
    L"$\frac{x}{y} \cdot d$"                          L"$x$"         ;
    L"$\left( t_{sub\_sub} - x \right)^{2 \cdot p}$"  L"$\frac{3}{4}$"
    ]

# @test_throws ErrorException latexify("x/y"; bad_kwarg="should error")


text:
using ImportMacros, Test, Pkg

# Add the test-package
Pkg.develop(PackageSpec(path = joinpath(@__DIR__, "A")))

@testset "@import" begin
    # A.B.C
    @test (@import A.B.C as abc) == nothing
    @import A.B.C as abc
    @test_throws UndefVarError baz()
    @test_throws UndefVarError C.baz()
    @test_throws UndefVarError Baz()
    @test_throws UndefVarError C.Baz()
    @test abc.baz()
    @test abc.Baz()

    # A.B
    @test (@import A.B as ab) == nothing
    @import A.B as ab
    @test_throws UndefVarError bar()
    @test_throws UndefVarError B.bar()
    @test_throws UndefVarError Bar()
    @test_throws UndefVarError B.Bar()
    @test ab.bar()
    @test ab.Bar()

    # A
    @test (@import A as a) == nothing
    @import A as a
    @test_throws UndefVarError foo()
    @test_throws UndefVarError A.foo()
    @test_throws UndefVarError Foo()
    @test_throws UndefVarError A.Foo()
    @test a.foo()
    @test a.Foo()
end # testset

@testset "@import with specific objects" begin
    # A.B.C
    @test (@import A.B.C.baz as bz) == nothing
    @import A.B.C.baz as bz
    @import A.B.C.Baz as Bz
    @import A.B.C.@baz as @bz
    @test_throws UndefVarError baz()
    @test_throws UndefVarError C.baz()
    @test_throws UndefVarError Baz()
    @test_throws UndefVarError C.Baz()
    @test_throws UndefVarError @eval $(Symbol("@baz"))
    @test_throws UndefVarError @eval C.$(Symbol("@baz"))
    @test bz()
    @test Bz()
    @test @eval (@bz 1 2 3) == (1, 2, 3)

    # A.B
    @test (@import A.B.bar as b) == nothing
    @import A.B.bar as b
    @import A.B.Bar as Br
    @import A.B.@bar as @b
    @test_throws UndefVarError bar()
    @test_throws UndefVarError B.bar()
    @test_throws UndefVarError Bar()
    @test_throws UndefVarError B.Bar()
    @test_throws UndefVarError @eval $(Symbol("@bar"))
    @test_throws UndefVarError @eval B.$(Symbol("@bar"))
    @test b()
    @test Br()
    @test @eval (@b 1 2) == (1, 2)

    # A
    @test (@import A.foo as f) == nothing
    @import A.foo as f
    @import A.Foo as F
    @import A.@foo as @f
    @import A.@Foo as @F
    @test_throws UndefVarError foo()
    @test_throws UndefVarError A.foo()
    @test_throws UndefVarError Foo()
    @test_throws UndefVarError A.Foo()
    @test_throws UndefVarError @eval $(Symbol("@foo"))
    @test_throws UndefVarError @eval A.$(Symbol("@foo"))
    @test_throws UndefVarError @eval $(Symbol("@Foo"))
    @test_throws UndefVarError @eval A.$(Symbol("@Foo"))
    @test f()
    @test F()
    @test @eval (@f 1) == 1
    @test @eval @F
end

@testset "@using" begin
    # A.B.C
    @test (@using A.B.C as abc) == nothing
    @using A.B.C as abc
    @test baz()
    @test C.baz()
    @test_throws UndefVarError Baz()
    @test C.Baz()
    @test abc.baz()
    @test abc.Baz()

    # A.B
    @test (@using A.B as ab) == nothing
    @using A.B as ab
    @test bar()
    @test B.bar()
    @test_throws UndefVarError Bar()
    @test B.Bar()
    @test ab.bar()
    @test ab.Bar()

    # A
    @test (@using A as a) == nothing
    @using A as a
    @test foo()
    @test A.foo()
    @test_throws UndefVarError Foo()
    @test A.Foo()
    @test a.foo()
    @test a.Foo()
end # testset


text:
using Druid
using Test

@testset "Aggregations" begin
    include("aggregations.jl")
end

@testset "Datasources" begin
    include("datasources.jl")
end

@testset "Dimensionspecs" begin
    include("dimensionspecs.jl")
end

@testset "Filters" begin
    include("filters.jl")
end

@testset "Granularities" begin
    include("granularities.jl")
end

@testset "Havingspecs" begin
    include("having.jl")
end

@testset "Intervals" begin
    include("intervals.jl")
end

@testset "PostAggregations" begin
    include("post_aggregations.jl")
end

@testset "VirtualColumns" begin
    include("virtualcolumns.jl")
end

@testset "Queries" begin
    include("queries.jl")
end


text:
# -------------- CH 11: POLICY GRADIENT ESTIMATION --------------

export  simulate,
        FiniteDifferenceGradient,
        gradient,
        RegressionGradient,
        RewardToGoGradient,
        BaselineSubtractionGradient

"""
    function simulate(𝒫::MDP, s, π, d)

Generate depth-'d' sequence (trajectory `τ`) of state-action-reward tuples following policy `π` starting from state `s`.
"""
function simulate(𝒫::MDP, s, π, d)
    τ = []
    # iterate to depth d
    for i = 1:d
        # follow policy from current state
        a = π(s)
        # simulate dynamics for current state-action
        s′, r = 𝒫.TR(s,a)
        # save trajectory step
        push!(τ, (s,a,r))
        # move to new state (and repeat)
        s = s′
    end
    # return record of state-action-reward tuples
    return τ
end

"""
    struct FiniteDifferenceGradient
        𝒫 # problem
        b # initial state distribution
        d # depth
        m # number of samples
        δ # step size
    end
"""
struct FiniteDifferenceGradient
    𝒫 # problem
    b # initial state distribution
    d # depth
    m # number of samples
    δ # step size
end

@doc raw"""
    function gradient(M::FiniteDifferenceGradient, π, θ)

Calculate policy gradient for policy `π` at parameterization `θ`, with respect to `θ`.

Returns: $\nabla_{\theta} U(\pi_{\theta})$
"""
function gradient(M::FiniteDifferenceGradient, π, θ)

    𝒫, b, d, m, δ, γ, n = M.𝒫, M.b, M.d, M.m, M.δ, M.𝒫.γ, length(θ)

    # helper function: generate δ⋅(one-hot vector length n, index i-hot)
    Δθ(i) = [i == k ? δ : 0.0 for k in 1:n]

    # get trajectory τ's discounted (depth-d) reward.
    R(τ) = sum(r*γ^(k-1) for (k, (s,a,r)) in enumerate(τ))

    # calculate mean of discounted reward over m rollouts using policy π and parameterization θ, each starting from random state s ∼ b.
    U(θ) = mean(R(simulate(𝒫, rand(b), s->π(θ, s), d)) for i in 1:m)

    # Calculate finite difference of U(θ) for each dimension, using mean of rollouts to estimate policy value function.
    ΔU = [U(θ + Δθ(i)) - U(θ) for i in 1:n]

    return ΔU ./ δ
end

# ----- Regression Gradient ------
#
# Estimate gradient using m > 2n samples of finite differences in random directions.
"""
    struct RegressionGradient
        𝒫 # problem
        b # initial state distribution
        d # depth
        m # number of samples
        δ # step size
    end
"""
struct RegressionGradient
    𝒫 # problem
    b # initial state distribution
    d # depth
    m # number of samples
    δ # step size
end

"""
    function gradient(M::RegressionGradient, π, θ)

Estimate policy gradient centered at `θ` using `m` random perturbations and then using least-squares over estimated finite-differences of each.

# Example

```julia
julia> using Random

julia> using LinearAlgebra

julia> f(x) = x^2 + 1e-2*randn()

julia> m = 20

julia> δ = 1e-2

julia> ΔX = [δ.*randn() for i = 1:m]

julia> x0 = 2.0

julia> ΔF = [f(x0 + Δx) - f(x0) for Δx in ΔX]

julia> pinv(ΔX) * ΔF       # should be around f′(x==2) == 2×(x==2) == 4
```
"""
function gradient(M::RegressionGradient, π, θ)
    𝒫, b, d, m, δ, γ = M.𝒫, M.b, M.d, M.m, M.δ, M.𝒫.γ
    # m random perturbations of θ from uniform distriubtion over radius δ hypersphere.
    ΔΘ = [δ.*normalize(randn(length(θ)), 2) for i = 1:m]
    # calculate discounted reward from trajectory τ
    R(τ) = sum(r*γ^(k-1) for (k, (s,a,r)) in enumerate(τ))
    # calculate discounted reward from rollout using policy π and parameterization θ, starting from random state s ∼ b.
    U(θ) = R(simulate(𝒫, rand(b), s->π(θ,s), d))
    # finite differenes
    ΔU = [U(θ + Δθ) - U(θ) for Δθ in ΔΘ]
    # direct solve least-squares solution of gradient using pseudoinverse
    return pinv(reduce(hcat, ΔΘ)') * ΔU
end

"""
    struct LikelihoodRatioGradient
        𝒫 # problem
        b # initial state distribution
        d # depth
        m # number of samples
        ∇logπ # gradient of log likelihood
    end
"""
struct LikelihoodRatioGradient
    𝒫 # problem
    b # initial state distribution
    d # depth
    m # number of samples
    ∇logπ   # gradient of log likelihood
            #  ie: gradient w.r.t. θ of π(θ,a,s),
            #   where π(θ,a,s) = probability(a | s; θ)
end

"""
    function gradient(M::LikelihoodRatioGradient, π, θ)

"""
function gradient(M::LikelihoodRatioGradient, π, θ)
    𝒫, b, d, m, ∇logπ, γ = M.𝒫, M.b, M.d, M.m, M.∇logπ, M.𝒫.γ

    # πθ: stochastic/non-deterministic function of s
    πθ(s) = π(θ, s)

    # Calculate total discounted reward from trajectory τ
    R(τ) = sum(r*γ^(k-1) for (k, (s,a,r)) in enumerate(τ))

    # ∇₀log[p₀(τ)] × R(τ) = ∑ₖ₌₁(∇₀log[π(θ,a⁽ᵏ⁾,s⁽ᵏ⁾)]) × R(τ)
    ∇U(τ) = sum(∇logπ(θ, a, s) for (s,a) in τ)*R(τ)

    # monte carlo estimate of ∇₀U(θ) = 𝔼ₜ[∇₀log(p₀(τ)R(τ))] ∼ τ
    return mean(∇U(simulate(𝒫, rand(b), πθ, d)) for i in 1:m)
end

"""
    struct RewardToGoGradient
        𝒫 # problem
        b # initial state distribution
        d # depth
        m # number of samples
        ∇logπ # gradient of log likelihood
    end
"""
struct RewardToGoGradient
    𝒫 # problem
    b # initial state distribution
    d # depth
    m # number of samples
    ∇logπ # gradient of log likelihood
end

"""
    function gradient(M::RewardToGoGradient, π, θ)

"""
function gradient(M::RewardToGoGradient, π, θ)
    𝒫, b, d, m, ∇logπ, γ = M.𝒫, M.b, M.d, M.m, M.∇logπ, M.𝒫.γ
    πθ(s) = π(θ, s)

    # see eq. 11.24: == γᵏ⁻¹ × r_to-go⁽ᵏ⁾
    R(τ, k) = sum(r*γ^(l-1) for (l,(s,a,r)) in zip(k:d, τ[k:end]))
    # see eq. 11.24
    ∇U(τ) = sum(∇logπ(θ, a, s)*R(τ,k) for (k, (s,a,r)) in enumerate(τ))

    # monte carlo approximation to 𝔼xpectation over τ
    return mean(∇U(simulate(𝒫, rand(b), πθ, d)) for i in 1:m)
end

"""
    struct BaselineSubtractionGradient
        𝒫 # problem
        b # initial state distribution
        d # depth
        m # number of samples
        ∇logπ # gradient of log likelihood
    end
"""
struct BaselineSubtractionGradient
    𝒫 # problem
    b # initial state distribution
    d # depth
    m # number of samples
    ∇logπ # gradient of log likelihood
end

"""
    function gradient(M::BaselineSubtractionGradient, π, θ)

"""
function gradient(M::BaselineSubtractionGradient, π, θ)
    𝒫, b, d, m, ∇logπ, γ = M.𝒫, M.b, M.d, M.m, M.∇logπ, M.𝒫.γ
    πθ(s) = π(θ, s)

    l(a, s, k) = ∇logπ(θ, a, s)*γ^(k-1)
    # r_to-go
    R(τ, k) = sum(r*γ^(j-1) for (j,(s,a,r)) in enumerate(τ[k:end]))

    numer(τ) = sum(l(a,s,k).^2*R(τ,k) for (k,(s,a,r)) in enumerate(τ))

    denom(τ) = sum(l(a,s,k).^2 for (k,(s,a)) in enumerate(τ))
    # 11.43
    base(τ) = numer(τ) ./ denom(τ)

    trajs = [simulate(𝒫, rand(b), πθ, d) for i in 1:m]
    # 11.43
    rbase = mean(base(τ) for τ in trajs)
    # 11.28
    ∇U(τ) = sum(l(a,s,k).*(R(τ,k).-rbase) for (k,(s,a,r)) in enumerate(τ))
    # 11.28
    return mean(∇U(τ) for τ in trajs)
end


text:
using NMF
using Test
using Random
using LinearAlgebra
using StatsBase

tests = ["utils",
         "initialization",
         "spa",
         "multupd",
         "alspgrad",
         "coorddesc",
         "greedycd",
         "interf"]

println("Running tests:")
for t in tests
    tp = "$t.jl"
    println("* $tp ...")
    include(tp)
end


text:
const beale = let res_init=zeros(3), jac_init=zeros(3,2), x_init=[0.1, 0.1]

    function res(x, r)
        r[1] = 1.5 - x[1] * (1.0 - x[2])
        r[2] = 2.25 - x[1] * (1.0 - x[2]^2)
        r[3] = 2.625 - x[1] * (1.0 -  x[2]^3)
        return r
    end

    function jac(x, jac)
        jac[1, 1] = x[2] - 1.0
        jac[1, 2] = x[1]
        jac[2, 1] = x[2] ^ 2 - 1.0
        jac[2, 2] = 2.0 * x[1] * x[2]
        jac[3, 1] = x[2] ^ 3 - 1.0
        jac[3, 2] = 3.0 * x[1] * (x[2] ^ 2)
        return jac
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("beale", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end


text:
using FRGRealTime
using HCubature
using Plots

δk=0.02
ppfunps_delta(p0, k, kprim, m, T)

plot(p0 -> 2*ppfunps_delta(p0, 10.0, 10.0, -2.0, 145.0), 19.7, 19.9)
plot(
    p0 -> FRGRealTime.delta2_intcosthqs(p0, 0.001, 10.0, 10.0, -2.0, 145.0),
    19.7, 19.9
)

plot(
    p0 -> FRGRealTime.F4tildekSmooth_intcosthqs(p0,0.001,10.0,-2.0,145.0,400.0,0.1,0.00001),
    19.0, 20.9
)

function testf(p0, k)
    -hquadrature(
        x -> FRGRealTime.delta2Smooth_intcosthqs(p0, 0.001, k, x, -2.0, 145.0,0.00,0.01),
        k,
        400.0,initdiv=500,rtol=1e-8,atol=1e-8,
    )[1]
end

plot!(p0->testf(p0,10.0),19.0, 20.9)

plot(
    p0 -> FRGRealTime.delta2_intcosthqs(p0, 0.001, 10.0, 10.0, -2.0, 145.0),
    1.0, 30.9
)

plot(
    p0 -> FRGRealTime.PvdkF1Tildeps(p0, 0.001, 10.0, 10.0, -2.0, 145.0),
    30.0,
    40.0,
)

plot(
    p0 -> FRGRealTime.PvdkF1Tildeps(p0, 10.0, 10.0, -2.0, 145.0),
    30.0,
    40.0,
)

plot!(
    p0 -> FRGRealTime.PvdkF1Tildeps2(p0, 10.0, 10.0, -2.0, 145.0),
    30.0,
    40.0,
)

plot(
    p0 ->
        FRGRealTime.flowpp_intcostheqs(p0, 0.0000000, 10.0, 10.0, -50.0, 145.0),
    2 * Epi(10.0, -50),
    20.0,
)

plot!(
    p0 -> 2 * FRGRealTime.ppfunps(p0, 10.0, 10.0, -50.0, 145.0),
    2 * Epi(10.0, -50),
    20.0,
)

testf()=2 * Epi(k, m)

function flowppintp0qs(p0, qsmax, k, m, T)
    quadgk(
        ps -> flowppintp0(p0, ps, k, m, T),
        0.0,
        qsmax,
        atol = 1e-4,
        rtol = 1e-4,
    )
end

plot(p0->FRGRealTime.flowpm(p0, p0, 10.0, -50.0, 100.0),0.0000001,0.1)

plot(k->FRGRealTime.loopfunpm(0.1, 0.1, k, -10.0, 100.0),4.0,20.0)

plot(p0->FRGRealTime.loopfunpm(p0, p0+0.1, 4.0, -10.0, 100.0),0.0,0.00000001)

using Plots,FiniteDifferences

plot(p0->central_fdm(5,1)(x->FRGRealTime.deltasumkfix(
    x,
    130.001,
    10.0,
    145.0,
    400.0,
    1.0,
    800.0,
    m2fun,
    lamfun,
),p0),1.0,800.0)

plot(p0->FRGRealTime.deltasumkfix(
    p0,
    0.001,
    80.0,
    145.0,
    400.0,
    1.0,
    800.0,
    m2fun,
    lamfun,
),250.0,260.0)

using HCubature

hquadrature(
    x -> FRGRealTime.deltafun(x - 10.0,0.1),
    3.0,
    20.0,
    rtol = 1e-12,
    atol = 1e-12,
    initdiv = 20,
)[1]

hquadrature(k->hquadrature(
    x -> FRGRealTime.deltadxfun(x - 10.0,0.01),
    k,
    20.0)[1],1.0,20.0)[1]

plot(
    p0 -> FRGRealTime.deltasumkfix(
        p0,
        0.00001,
        1.0,
        145.0,
        4.0,
        1.0,
        150.0,
        m2fun,
        lamfun,
    ),38.0,42.0
)

plot(
    p0 -> FRGRealTime.VImintqs(
        p0,
        0.0001,
        50.0,
        100.0,
        4.0,
        1.0,
        200.0,
        m2funnew,
        lamfunnew,
    )+FRGRealTime.VImintqs_delta1(
        p0,
        0.0001,
        50.0,
        100.0,
        4.0,
        1.0,
        200.0,
        m2fun,
        lamfun,
    ),52.0,60.0
)

plot(p0->central_fdm(5,1)(x->FRGRealTime.VImintqs(
    p0,
    0.0001,
    x,
    50.0,
    100.0,
    4.0,
    1.0,
    200.0,
    m2fun,
    lamfun,
),Epi(50.0,100.0)),2.0,20.0)

plot!(
    p0 -> central_fdm(3,1)(x->FRGRealTime.VImintqsSimple(
        p0,
        0.0001,
        x,
        50.0,
        100.0,
        4.0,
        10.0^2,
        -10.0,
        200.0,
    ),Epi(50.0,100.0)),
    2.0,
    20.0,
)

using Plots
plot(
    p0 -> FRGRealTime.VImintqs(
        p0,
        0.0001,
        100.0,
        40.0,
        4.0,
        1.0,
        1600.0,
        m2fun,
        lamfun,
    ),
    1.0,
    800.0,label="I2, p0-q0",legend=:topleft,
)

sqrt(50^2+m2fun(50))
xaxis!("p0")

using QuadGK

function testf2(f;initdiv,kwarg...)
    quadgk(f,0.0,10.0;order=initdiv,kwarg...)
end