text:
# @generated function static_tuple(t::SDTuple{N,X,P}) where {N,X,P}
#     i = 0
#     sdt = Expr(:tuple)
#     Xv = tuple_type_to_value_tuple(X)
#     for n ∈ 1:N
#         push!(sdt.args, Xv[n] == -1 ? Expr(:ref, :x, (i += 1)) : Expr(:call, Expr(:curly, :StaticInt, Xv[n])))
#     end
#     q = Expr(:block, Expr(:meta, :inline))
#     i > 0 && push!(q.args, :(x = t.x))
#     push!(q.args, sdt)
#     q
# end

@inline mulsizeof(::Type{T}, x::Number) where {T} = vmul_fast(sizeof(T), x)
@generated mulsizeof(::Type{T}, ::StaticInt{N}) where {T,N} = Expr(:call, Expr(:curly, :StaticInt, N*sizeof(T)))
@inline mulsizeof(::Type{T}, ::Tuple{}) where {T} = ()
@inline mulsizeof(::Type{T}, x::Tuple{X}) where {T,X} = (mulsizeof(T, first(x)), )
@inline mulsizeof(::Type{T}, x::Tuple) where {T} = (mulsizeof(T, first(x)), mulsizeof(T, Base.tail(x))...)

@inline bytestrides(A::AbstractArray{T}) where {T} = mulsizeof(T, ArrayInterface.strides(A))

# @generated function mulsizeof(::Type{T}, ::Type{S}) where {T, N, S <: Tuple{Vararg{Any,N}}}
#     Smul = Expr(:curly, :Tuple)
#     for n in 1:N
#         Sₙ = S.parameters[n]
#         push!(Smul.args, Sₙ == -1 ? -1 : sizeof(T) * Sₙ)
#     end
#     Smul
# end
# @inline function mulsizeof(::Type{T}, t::SDTuple{N,X,P}) where {X,N,P,T}
#     SDTuple{N,mulsizeof(T, X),P}(mulsizeof(T, t.x))
# end

@inline memory_reference(A::AbstractArray) = memory_reference(device(A), A)
@inline memory_reference(::CPUPointer, A) = pointer(A)
@inline function memory_reference(::ArrayInterface.CheckParent, A)
    P = parent(A)
    if P === A
        memory_reference(ArrayInterface.CPUIndex(), A)
    else
        memory_reference(device(P), P)
    end
end
@inline memory_reference(::ArrayInterface.CPUIndex, A) = throw("Memory access for $(typeof(A)) not implemented yet.")

"""
  abstract type AbstractStridedPointer{T,N,C,B,R,X,O} end

T: element type
N: dimensionality
C: contiguous dim
B: batch size
R: rank of strides
X: strides
O: offsets
"""
abstract type AbstractStridedPointer{T<:NativeTypes,N,C,B,R,X<:Tuple{Vararg{Any,N}},O<:Tuple{Vararg{Any,N}}} end

@inline Base.eltype(::AbstractStridedPointer{T}) where {T} = T

struct StridedPointer{T,N,C,B,R,X,O} <: AbstractStridedPointer{T,N,C,B,R,X,O}
    p::Ptr{T}
    strd::X
    offsets::O
end
@inline StridedPointer{T,N,C,B,R}(ptr::Ptr{T}, strd::X, o::O) where {T,N,C,B,R,X,O} = StridedPointer{T,N,C,B,R,X,O}(ptr, strd, o)
@inline StridedPointer{T,N,C,B,R,X}(ptr::Ptr{T}, strd::X, o::O) where {T,N,C,B,R,X,O} = StridedPointer{T,N,C,B,R,X,O}(ptr, strd, o)

@inline function stridedpointer(A::AbstractArray{T}) where {T <: NativeTypes}
    stridedpointer(memory_reference(A), contiguous_axis(A), contiguous_batch_size(A), stride_rank(A), bytestrides(A), offsets(A))
end
@inline function stridedpointer(
    ptr::Ptr{T}, ::Contiguous{C}, ::ContiguousBatch{B}, ::StrideRank{R}, strd::X, offsets::O
) where {T<:NativeTypes,C,B,R,N,X<:Tuple{Vararg{Any,N}},O<:Tuple{Vararg{Any,N}}}
    StridedPointer{T,N,C,B,R,X,O}(ptr, strd, offsets)
end
@inline Base.strides(ptr::AbstractStridedPointer) = ptr.strd
@inline ArrayInterface.offsets(ptr::AbstractStridedPointer) = ptr.offsets
@inline ArrayInterface.contiguous_axis_indicator(ptr::AbstractStridedPointer{T,N,C}) where {T,N,C} = contiguous_axis_indicator(Contiguous{C}(), Val{N}())

@generated function zerotuple(::Val{N}) where {N}
    t = Expr(:tuple);
    for n in 1:N
        push!(t.args, :(Zero()))
    end
    Expr(:block, Expr(:meta,:inline), t)
end

@inline function zero_offsets(sptr::StridedPointer{T,N,C,B,R}) where {T,N,C,B,R}
    StridedPointer{T,N,C,B,R}(sptr.p, sptr.strd, zerotuple(Val{N}()))
end
@inline zstridedpointer(A) = zero_offsets(stridedpointer(A))

@inline function Base.similar(sptr::StridedPointer{T,N,C,B,R,X,O}, ptr::Ptr{T}) where {T,N,C,B,R,X,O}
    StridedPointer{T,N,C,B,R,X,O}(ptr, sptr.strd, sptr.offsets)
end
# @inline noalias!(p::StridedPointer) = similar(p, noalias!(pointer(p)))
@inline function similar_no_offset(sptr::StridedPointer{T,N,C,B,R,X,O}, ptr::Ptr{T}) where {T,N,C,B,R,X,O}
    StridedPointer{T,N,C,B,R,X}(ptr, sptr.strd, zerotuple(Val{N}()))
end

@inline Base.pointer(ptr::StridedPointer) = ptr.p
Base.unsafe_convert(::Type{Ptr{T}}, ptr::AbstractStridedPointer{T}) where {T} = pointer(ptr)
# Shouldn't need to special case Array
# function stridedpointer(A::Array{T,N}) where {T,N}
#     StridedPointer{T,1,0,ntuple(identity,Val{N}()),ntuple(n -> isone(n) ? 1 : -1, Val{N}()), N, N-1}(pointer(A), Base.tail(strides(A)))
# end

@inline vload(ptr::AbstractStridedPointer) = vload(pointer(ptr))
@inline vstore!(ptr::AbstractStridedPointer{T}, v::T) where {T} = vstore!(pointer(ptr), v)

@generated function nopromote_axis_indicator(::AbstractStridedPointer{<:Any,N}) where {N}
    t = Expr(:tuple); foreach(n -> push!(t.args, Expr(:call, Expr(:curly, :Val, true))), 1:N)
    Expr(:block, Expr(:meta, :inline), t)
end

@inline vload(ptr::AbstractStridedPointer{T,0}, i::Tuple{}, ::Val{A}) where {A,T} = vload(pointer(ptr), Val{A}())
@inline gep(ptr::AbstractStridedPointer{T,0}, i::Tuple{}) where {T} = pointer(ptr)

# Fast compile path?
@inline function vload(ptr::AbstractStridedPointer{T,N,C,B,R,X,NTuple{N,StaticInt{0}}}, i::Tuple{Vararg{Any,N}}, ::Val{A}) where {A,T,N,C,B,R,X}
    vload(pointer(ptr), tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), Val{A}())
end
@inline function vload(ptr::AbstractStridedPointer{T,N,C,B,R,X,NTuple{N,StaticInt{0}}}, i::Tuple{Vararg{Any,N}}, m, ::Val{A}) where {A,T,N,C,B,R,X}
    vload(pointer(ptr), tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}())
end
@inline function vload(ptr::AbstractStridedPointer{T,N,C,B,R,X,O}, i::Tuple{Vararg{Any,N}}, ::Val{A}) where {A,T,N,C,B,R,X,O}
    vload(pointer(ptr), tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), contiguous_axis_indicator(ptr)), Val{A}())
end
@inline function vload(ptr::AbstractStridedPointer{T,N,C,B,R,X,O}, i::Tuple{Vararg{Any,N}}, m, ::Val{A}) where {A,T,N,C,B,R,X,O}
    vload(pointer(ptr), tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}())
end
@inline function vload(ptr::AbstractStridedPointer{T}, i::Tuple{I}, ::Val{A}) where {A,T,I}
    vload(pointer(ptr), tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), Val{A}())
end
@inline function vload(ptr::AbstractStridedPointer{T}, i::Tuple{I}, m, ::Val{A}) where {A,T,I}
    vload(pointer(ptr), tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}())
end
# Ambiguity: 1-dimensional + 1-dim index -> Cartesian (offset) indexing
@inline function vload(ptr::AbstractStridedPointer{T,1,C,B,R,X,O}, i::Tuple{I}, ::Val{A}) where {A,T,I,C,B,R,X,O}
    vload(pointer(ptr), tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), contiguous_axis_indicator(ptr)), Val{A}())
end
@inline function vload(ptr::AbstractStridedPointer{T,1,C,B,R,X,O}, i::Tuple{I}, m, ::Val{A}) where {A,T,I,C,B,R,X,O}
    vload(pointer(ptr), tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}())
end
# Ambiguity: 1-dimensional + 1-dim index -> Cartesian (offset) indexing
@inline function vload(ptr::AbstractStridedPointer{T,1,C,B,R,X,Tuple{StaticInt{0}}}, i::Tuple{I}, ::Val{A}) where {A,T,I,C,B,R,X}
    vload(pointer(ptr), tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), Val{A}())
end
@inline function vload(ptr::AbstractStridedPointer{T,1,C,B,R,X,Tuple{StaticInt{0}}}, i::Tuple{I}, m, ::Val{A}) where {A,T,I,C,B,R,X}
    vload(pointer(ptr), tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}())
end

# align, noalias, nontemporal
@inline function vstore!(
    ptr::AbstractStridedPointer{T,N,C,B,R,X,NTuple{N,StaticInt{0}}}, v, i::Tuple{Vararg{Any,N}}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,N,C,B,R,X,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T,N,C,B,R,X,NTuple{N,StaticInt{0}}}, v, i::Tuple{Vararg{Any,N}}, m, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,N,C,B,R,X,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T,N,C,B,R,X,O}, v, i::Tuple{Vararg{Any,N}}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,N,C,B,R,X,O,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), contiguous_axis_indicator(ptr)), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T,N,C,B,R,X,O}, v, i::Tuple{Vararg{Any,N}}, m, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,N,C,B,R,X,O,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T}, v, i::Tuple{I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,I,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T}, v, i::Tuple{I}, m, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,I,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T,1,C,B,R,X,O}, v, i::Tuple{I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,I,C,B,R,X,O,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), contiguous_axis_indicator(ptr)), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T,1,C,B,R,X,O}, v, i::Tuple{I}, m, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,I,C,B,R,X,O,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T,1,C,B,R,X,Tuple{StaticInt{0}}}, v, i::Tuple{I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,I,C,B,R,X,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::AbstractStridedPointer{T,1,C,B,R,X,Tuple{StaticInt{0}}}, v, i::Tuple{I}, m, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,I,C,B,R,X,A,S,NT}
    vstore!(pointer(ptr), v, tdot(ptr, i, strides(ptr), contiguous_axis_indicator(ptr)), m, Val{A}(), Val{S}(), Val{NT}())
end
@inline function gep(ptr::AbstractStridedPointer{T,N,C,B,R,X,NTuple{N,StaticInt{0}}}, i::Tuple{Vararg{Any,N}}) where {T,N,C,B,R,X}
    gep(pointer(ptr), tdot(ptr, i, strides(ptr), nopromote_axis_indicator(ptr)))
end
@inline function gep(ptr::AbstractStridedPointer{T,N,C,B,R,X,O}, i::Tuple{Vararg{Any,N}}) where {T,N,C,B,R,X,O}
    gep(pointer(ptr), tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), nopromote_axis_indicator(ptr)))
end
@inline function gep(ptr::AbstractStridedPointer{T}, i::Tuple{I}) where {T, I}
    gep(pointer(ptr), tdot(ptr, i, strides(ptr), nopromote_axis_indicator(ptr)))
end
@inline function gep(ptr::AbstractStridedPointer{T,1,C,B,R,X,O}, i::Tuple{I}) where {T, I,C,B,R,X,O}
    gep(pointer(ptr), tdot(ptr, map(vsub_fast, i, offsets(ptr)), strides(ptr), nopromote_axis_indicator(ptr)))
end
@inline function gep(ptr::AbstractStridedPointer{T,1,C,B,R,X,Tuple{StaticInt{0}}}, i::Tuple{I}) where {T, I,C,B,R,X}
    gep(pointer(ptr), tdot(ptr, i, strides(ptr), nopromote_axis_indicator(ptr)))
end

struct StridedBitPointer{N,C,B,R,X,O} <: AbstractStridedPointer{Bit,N,C,B,R,X,O}
    p::Ptr{Bit}
    strd::X
    offsets::NTuple{N,Int}
end
function StridedBitPointer{N,C,B,R}(p::Ptr{Bit}, strd::X, offsets::O) where {N,C,B,R,X,O}
    StridedBitPointer{N,C,B,R,X,O}(p, strd, offsets)
end
@inline Base.pointer(p::StridedBitPointer) = p.p
# @inline stridedpointer(A::BitVector) = StridedBitPointer{1,1,0,(1,)}(Base.unsafe_convert(Ptr{Bit}, pointer(A.chunks)), (StaticInt{1}(),), (StaticInt{1}(),))
@inline stridedpointer(A::BitVector) = StridedBitPointer{1,1,0,(1,)}(Base.unsafe_convert(Ptr{Bit}, pointer(A.chunks)), (StaticInt{1}(),), (1,))
@generated function stridedpointer(A::BitArray{N}) where {N}
    q = quote;
        s = size(A)
        @assert iszero(s[1] & 7) "For performance reasons, `BitArray{N}` where `N > 1` are required to have a multiple of 8 rows.";
    end
    sone = :(StaticInt{1}());
    strd = Expr(:tuple, sone, :s_2); offsets = Expr(:tuple, sone, sone);
    last_stride = next_stride = :s_2
    push!(q.args, :(s_2 = size(A,1))); # >>> 3
    R = Expr(:tuple, 1, 2)
    for n ∈ 3:N
        next_stride = Symbol(:s_, n)
        push!(q.args, Expr(:(=), next_stride, Expr(:call, :(*), Expr(:ref, :s, n-1), last_stride)))
        push!(strd.args, next_stride)
        # push!(offsets.args, :(StaticInt{1}()))
        push!(offsets.args, 1)
        last_stride = next_stride
        push!(R.args, n)
    end
    push!(q.args, :(StridedBitPointer{$N,1,0,$R}(Base.unsafe_convert(Ptr{Bit}, pointer(A.chunks)), $strd, $offsets)))
    q
end
@inline function similar_no_offset(sptr::StridedBitPointer{N,C,B,R,X,O}, ptr::Ptr{Bit}) where {N,C,B,R,X,O}
    StridedBitPointer{N,C,B,R}(ptr, sptr.strd, zerotuple(Val{N}()))
end

@generated function gesp(ptr::StridedBitPointer{N,C,B,R}, i::Tuple{Vararg{Any,N}}) where {N,C,B,R}
    quote
        $(Expr(:meta,:inline))
        offs = ptr.offsets
        StridedBitPointer{$N,$C,$B,$R}(ptr.p, ptr.strd, Base.Cartesian.@ntuple $N n -> vsub_fast(offs[n], i[n]))
    end
end
@generated function pointerforcomparison(p::StridedBitPointer{N}) where {N}
    inds = Expr(:tuple); foreach(_ -> push!(inds.args, :(Zero())), 1:N)
    Expr(:block, Expr(:meta,:inline), Expr(:call, :gep, :p, inds))
end
# @inline tdot(ptr::StridedBitPointer, a, b, c) = tdot(Bool, a, b, c) >>> StaticInt(3)

# There is probably a smarter way to do indexing adjustment here.
# The reasoning for the current approach of geping for Zero() on extracted inds
# and for offsets otherwise is best demonstrated witht his motivational example:
#
# A = OffsetArray(rand(10,11), 6:15, 5:15);
# for i in 6:15
    # s += A[i,i]
# end
# first access is at zero-based index
# (first(6:16) - ArrayInterface.offsets(a)[1]) * ArrayInterface.strides(A)[1] + (first(6:16) - ArrayInterface.offsets(a)[2]) * ArrayInterface.strides(A)[2]
# equal to
#  (6 - 6)*1 + (6 - 5)*10 = 10
# i.e., the 1-based index 11.
# So now we want to adjust the offsets and pointer's value
# so that indexing it with a single `i` (after summing strides) is correct.
# Moving to 0-based indexing is easiest for combining strides. So we gep by 0 on these inds.
# E.g., gep(stridedpointer(A), (0,0))
# ptr += (0 - 6)*1 + (0 - 5)*10 = -56
# new stride = 1 + 10 = 11
# new_offset = 0
# now if we access the 6th element
# (6 - new_offse) * new_stride
# ptr += (6 - 0) * 11 = 66
# cumulative:
# ptr = pointer(A) + 66 - 56  = pointer(A) + 10
# so initial load is of pointer(A) + 10 -> the 11th element w/ 1-based indexing

function double_index_quote(C,B,R::NTuple{N,Int},I1,I2,typ) where {N}
    # place into position of second arg
    J1 = I1 + 1; J2 = I2 + 1;
    @assert (J1 != B) & (J2 != B)
    Cnew = ((C == J1) | (C == J2)) ? -1 : (C - (J1 < C))
    push!(typ.args, Cnew); push!(typ.args, B);
    strd = Expr(:tuple); offs = Expr(:tuple);
    inds = Expr(:tuple); Rtup = Expr(:tuple)
    for n in 1:N
        if n == J1
            push!(inds.args, :(Zero()))
        elseif n == J2
            push!(strd.args, Expr(:call, :vadd_fast, Expr(:ref, :strd, J1), Expr(:ref, :strd, J2)))
            push!(offs.args, :(Zero()))
            push!(inds.args, :(Zero()))
            push!(Rtup.args, max(R[J1], R[J2]))
        else
            push!(strd.args, Expr(:ref, :strd, n))
            push!(offs.args, Expr(:ref, :offs, n))
            push!(inds.args, Expr(:ref, :offs, n))
            push!(Rtup.args, R[n])
        end
    end
    push!(typ.args, Rtup)
    gepedptr = Expr(:call, :gep, :ptr, inds)
    newptr = Expr(:call, typ, gepedptr, strd, offs)
    Expr(:block, Expr(:meta,:inline), :(strd = ptr.strd), :(offs = ptr.offsets), newptr)
end
@generated function double_index(ptr::StridedPointer{T,N,C,B,R}, ::Val{I1}, ::Val{I2}) where {T,N,C,B,R,I1,I2}
    double_index_quote(C,B,R,I1,I2, Expr(:curly, :StridedPointer, :T, N - 1))
end
@generated function double_index(ptr::StridedBitPointer{N,C,B,R}, ::Val{I1}, ::Val{I2}) where {N,C,B,R,I1,I2}
    double_index_quote(C,B,R,I1,I2, Expr(:curly, :StridedBitPointer, N - 1))
end

@inline stridedpointer(ptr::AbstractStridedPointer) = ptr

struct FastRange{T,F,S,O}# <: AbstractRange{T}
    f::F
    s::S
    offset::O
end
FastRange{T}(f::F,s::S,o::O) where {T,F,S,O} = FastRange{T,F,S,O}(f,s,o)
@inline function stridedpointer(r::AbstractRange{T}) where {T}
    FastRange{T}(ArrayInterface.static_first(r), ArrayInterface.static_step(r), One())
end
@inline function gesp(r::FastRange{T}, i::Tuple{I}) where {I,T}
    ii = first(i) - r.offset
    f = r.f
    s = r.s
    FastRange{T}(f + ii * s, r.s, Zero())
end
@inline vload(r::FastRange{T}, i::Tuple{I}) where {T,I} = convert(T, r.f) + convert(T, r.s) * (first(i) - convert(T, r.offset))
@inline vload(r::FastRange, i::Tuple, m::Mask) = vload(r, i)
@inline vload(r::FastRange, i::Tuple, m::Bool) = vload(r, i)
# @inline Base.getindex(r::FastRange, i::Integer) = vload(r, (i,))
@inline Base.eltype(::FastRange{T}) where {T} = T


text:
using Polynomials: polyint
using GeometricIntegrators.BasisFunctions: legendre_polynomial

@testset "$(rpad("Legendre Basis Tests",80))" begin

    b = LegendreBasis(Float64, 5)

    @test evaluate(b, 1, 0.0) == 1.0
    @test evaluate(b, 1, 0.5) == 1.0
    @test evaluate(b, 1, 1.0) == 1.0

    @test evaluate(b, 2, 0.5) == 0.0
    @test evaluate(b, 3, 0.5) != 0.0
    @test evaluate(b, 4, 0.5) == 0.0
    @test evaluate(b, 5, 0.5) != 0.0

    @test derivative(b, 1, 0.0) == 0.0
    @test derivative(b, 1, 0.5) == 0.0
    @test derivative(b, 1, 1.0) == 0.0

    @test derivative(b, 2, 0.0) == derivative(b, 2, 0.5) == derivative(b, 2, 1.0)

    @test integral(b, 1, 0.0) == 0.0
    @test integral(b, 1, 0.5) == 0.5
    @test integral(b, 1, 1.0) == 1.0

    @test integral(b, 2, 0.0) == 0.0
    @test integral(b, 2, 1.0) == 0.0

    for q in 1:10, p in 1:10
        f1, p1 = legendre_polynomial(q)
        f2, p2 = legendre_polynomial(p)
        int = polyint(p1 * p2)
        @test f1 * f2 * (int(1.0) - int(0.0)) ≈ (q == p ? 1.0 : 0.0) atol=1e-3
    end

end

ext: jl


text:
using SeaPearl

"""
Proposition for using SeaPearl CP solver for the Eternity II puzzle
Reference https://en.wikipedia.org/wiki/Eternity_II_puzzle
"""

include("IOmanager.jl")

struct InputData
    n    ::Int #number of lines
    m    ::Int #number of columns
    pieces ::Matrix{Int} #the n*m pieces
end

struct OutputDataEternityII
    nb_sols ::Int
    orientation::Array{Int, 4}# dims = (nb_sols, n, m, 5) where five corresponds to (id, u,r, d, l) (u=upper edge, ...)
end

"""
    model_eternity2(input_file; order=[1,2,3,4], limit=nothing)

return the SeaPearl model for to the solve_eternity2 problem using SeaPearl.AllDifferent  and SeaPearl.TableConstraint (without solving it)

# Arguments
- `input_file`: file containing the pieces of the game as well as the dimensions
- 'order' : Vector, giving the order of edges for the IO manager. example : [1,4,2,3] means it is given as [up,down,left,right]
- 'limit' : Int, giving the number of solutions after which it will stop searching. if nothing given, it will lookk for all the solutions
"""
function model_eternity2(input_file; order=[1,2,3,4], limit=nothing)
    trailer = SeaPearl.Trailer()
    model = SeaPearl.CPModel(trailer)
    model.limit.numberOfSolutions = limit
    inputData = getInputData(input_file;order=order)
    n = inputData.n
    m = inputData.m
    pieces = inputData.pieces

    table = Matrix{Int}(undef, 5, 4*n*m)

    maxi = maximum(pieces)

    for i = 1:n*m
        for k = 1:4
            table[1,4*(i-1) + k] = i
            for j = 2:5
                table[j,4*(i-1) + k] = pieces[i, (j+k+1)%4+1]
            end
        end
    end

    id = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)
    u = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#up
    r = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#right
    d = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#down
    l = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#left

    for i = 1:n, j=1:m
        id[i,j] = SeaPearl.IntVar(1, n*m, "id_"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, id[i,j]; branchable=true)
        u[i,j] = SeaPearl.IntVar(0, maxi, "u_"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, u[i,j]; branchable=false)
        r[i,j] = SeaPearl.IntVar(0, maxi, "r_"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, r[i,j]; branchable=false)
        d[i,j] = SeaPearl.IntVar(0, maxi, "d_"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, d[i,j]; branchable=false)
        l[i,j] = SeaPearl.IntVar(0, maxi, "l_"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, l[i,j]; branchable=false)

        vars = SeaPearl.AbstractIntVar[id[i,j], u[i,j], r[i,j],d[i,j], l[i,j]]
        push!(model.constraints, SeaPearl.TableConstraint(vars, table, trailer))

        if (j==m) push!(model.constraints, SeaPearl.EqualConstant(r[i,j], 0, trailer)) end
        if (j==1) push!(model.constraints, SeaPearl.EqualConstant(l[i,j], 0, trailer)) end
        if (i==1) push!(model.constraints, SeaPearl.EqualConstant(u[i,j], 0, trailer)) end
        if (i==n) push!(model.constraints, SeaPearl.EqualConstant(d[i,j], 0, trailer)) end
    end
    for i = 1:n, j=1:m
        if (j<m) push!(model.constraints, SeaPearl.Equal(r[i,j], l[i,j+1], trailer)) end
        if (i<n) push!(model.constraints, SeaPearl.Equal(d[i,j], u[i+1,j], trailer)) end
    end

    #breaking some symmetries

    #if count(==(2),count(==(0),pieces,dims=2))==4
        #index = findfirst(==(2),vec(count(==(0),pieces,dims=2)))
        #push!(model.constraints,SeaPearl.EqualConstant(id[1,1], index, trailer))
        #SeaPearl.assign!(id[1,1].domain, index)
    #end

    push!(model.constraints, SeaPearl.AllDifferent(id, trailer))
    model.adhocInfo = Dict([("n", n), ("m", m)])
    return model
end

"""
    model_eternity2_fast(input_file; order=[1,2,3,4], limit=nothing)

return the SeaPearl model for to the solve_eternity2 problem, using SeaPearl.AllDifferent  and SeaPearl.TableConstraint (without solving it). It uses less variables than v1. To be used with print_fast

# Arguments
- `input_file`: file containing the pieces of the game as well as the dimensions
- 'order' : Vector, giving the order of edges for the IO manager. example : [1,4,2,3] means it is given as [up,down,left,right]
- 'limit' : Int, giving the number of solutions after which it will stop searching. if nothing given, it will lookk for all the solutions
"""
function model_eternity2_fast(input_file; order=[1,2,3,4], limit=nothing)
    trailer = SeaPearl.Trailer()
    model = SeaPearl.CPModel(trailer)
    model.limit.numberOfSolutions = limit
    inputData = getInputData(input_file;order=order)
    n = inputData.n
    m = inputData.m
    pieces = inputData.pieces

    table = Matrix{Int}(undef, 5, 4*n*m)

    maxi = maximum(pieces)

    for i = 1:n*m
        for k = 1:4
            table[1,4*(i-1) + k] = i
            for j = 2:5
                table[j,4*(i-1) + k] = pieces[i, (j+k+1)%4+1]
            end
        end
    end

    id = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)
    u = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#up
    r = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#right
    d = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#down
    l = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#left
    src_v = Matrix{SeaPearl.AbstractIntVar}(undef,n,m+1) #horizontal edges
    src_h = Matrix{SeaPearl.AbstractIntVar}(undef,n+1,m) #vertical edges

    for i = 1:n
        src_v[i,m+1] = SeaPearl.IntVar(0, maxi, "src_v"*string(i)*string(m+1), trailer)
        SeaPearl.addVariable!(model, src_v[i,m+1]; branchable=false)
    end
    for j =1:m
        src_h[n+1,j] = SeaPearl.IntVar(0, maxi, "src_h"*string(n+1)*string(j), trailer)
        SeaPearl.addVariable!(model, src_h[n+1,j]; branchable=false)
    end

    for i = 1:n, j=1:m
        id[i,j] = SeaPearl.IntVar(1, n*m, "id_"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, id[i,j]; branchable=true)

        src_v[i,j] = SeaPearl.IntVar(0, maxi, "src_v"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, src_v[i,j]; branchable=false)
        src_h[i,j] = SeaPearl.IntVar(0, maxi, "src_h"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, src_h[i,j]; branchable=false)
    end

    for i = 1:n, j=1:m
        u[i,j] = src_h[i,j]
        d[i,j] = src_h[i+1,j]
        l[i,j] = src_v[i,j]
        r[i,j] = src_v[i,j+1]

        vars = SeaPearl.AbstractIntVar[id[i,j], u[i,j], r[i,j],d[i,j], l[i,j]]
        push!(model.constraints, SeaPearl.TableConstraint(vars, table, trailer))

        if (j==m) push!(model.constraints, SeaPearl.EqualConstant(r[i,j], 0, trailer)) end
        if (j==1) push!(model.constraints, SeaPearl.EqualConstant(l[i,j], 0, trailer)) end
        if (i==1) push!(model.constraints, SeaPearl.EqualConstant(u[i,j], 0, trailer)) end
        if (i==n) push!(model.constraints, SeaPearl.EqualConstant(d[i,j], 0, trailer)) end
    end

    #breaking some symmetries

    #if count(==(2),count(==(0),pieces,dims=2))==4
        #index = findfirst(==(2),vec(count(==(0),pieces,dims=2)))
        #push!(model.constraints,SeaPearl.EqualConstant(id[1,1], index, trailer))
        #SeaPearl.assign!(id[1,1].domain, index)
    #end

    push!(model.constraints, SeaPearl.AllDifferent(id, trailer))
    model.adhocInfo = Dict([("n", n), ("m", m)])
    return model
end

"""
    model_eternity2_rotation(input_file; order=[1,2,3,4], limit=nothing)

return the SeaPearl model for to the solve_eternity2 problem, using SeaPearl.AllDifferent  and SeaPearl.TableConstraint (without solving it). This time, it branches on id+orientation, which could be more conveninent for learning, but it adds variables.

# Arguments
- `input_file`: file containing the pieces of the game as well as the dimensions
- 'order' : Vector, giving the order of edges for the IO manager. example : [1,4,2,3] means it is given as [up,down,left,right]
- 'limit' : Int, giving the number of solutions after which it will stop searching. if nothing given, it will lookk for all the solutions
"""
function model_eternity2_rotation(input_file; order=[1,2,3,4], limit=nothing)
    trailer = SeaPearl.Trailer()
    model = SeaPearl.CPModel(trailer)
    model.limit.numberOfSolutions = limit
    inputData = getInputData(input_file;order=order)
    n = inputData.n
    m = inputData.m
    pieces = inputData.pieces

    table = Matrix{Int}(undef, 6, 4*n*m) #n*m pieces with for different orientation, 6 for orientation + id + u + r + d +

    maxi = maximum(pieces)
    orientation = Matrix{SeaPearl.AbstractIntVar}(undef, n,m)

    for i = 1:n*m
        for k = 1:4
            table[1,4*(i-1) + k] = i
            table[6,4*(i-1) + k] = 4*(i-1) + k
            for j = 2:5
                table[j,4*(i-1) + k] = pieces[i, (j+k+1)%4+1]
            end
        end
    end

    id = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)
    u = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#up
    r = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#right
    d = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#down
    l = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#left

    src_v = Matrix{SeaPearl.AbstractIntVar}(undef,n,m+1) #vertical edges
    src_h = Matrix{SeaPearl.AbstractIntVar}(undef,n+1,m) #horizontal edges

    for i = 1:n
        src_v[i,m+1] = SeaPearl.IntVar(0, maxi, "src_v"*string(i)*string(m+1), trailer)
        SeaPearl.addVariable!(model, src_v[i,m+1]; branchable=false)
    end
    for j =1:m
        src_h[n+1,j] = SeaPearl.IntVar(0, maxi, "src_h"*string(n+1)*string(j), trailer)
        SeaPearl.addVariable!(model, src_h[n+1,j]; branchable=false)
    end

    for i = 1:n, j=1:m
        orientation[i,j] = SeaPearl.IntVar(1, 4*n*m, "orientation_"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, orientation[i,j]; branchable=true)
        id[i,j] = SeaPearl.IntVar(1, n*m, "id_"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, id[i,j]; branchable=false)
        src_v[i,j] = SeaPearl.IntVar(0, maxi, "src_v"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, src_v[i,j]; branchable=false)
        src_h[i,j] = SeaPearl.IntVar(0, maxi, "src_h"*string(i)*string(j), trailer)
        SeaPearl.addVariable!(model, src_h[i,j]; branchable=false)
    end

    for i = 1:n, j=1:m
        u[i,j] = src_h[i,j]
        d[i,j] = src_h[i+1,j]
        l[i,j] = src_v[i,j]
        r[i,j] = src_v[i,j+1]

        vars = SeaPearl.AbstractIntVar[id[i,j], u[i,j], r[i,j],d[i,j], l[i,j], orientation[i,j]]
        push!(model.constraints, SeaPearl.TableConstraint(vars, table, trailer))

        if (j==m) push!(model.constraints, SeaPearl.EqualConstant(r[i,j], 0, trailer)) end
        if (j==1) push!(model.constraints, SeaPearl.EqualConstant(l[i,j], 0, trailer)) end
        if (i==1) push!(model.constraints, SeaPearl.EqualConstant(u[i,j], 0, trailer)) end
        if (i==n) push!(model.constraints, SeaPearl.EqualConstant(d[i,j], 0, trailer)) end
    end

    #breaking some symmetries

    #if count(==(2),count(==(0),pieces,dims=2))==4
        #index = findfirst(==(2),vec(count(==(0),pieces,dims=2)))
        #push!(model.constraints,SeaPearl.EqualConstant(id[1,1], index, trailer))
        #SeaPearl.assign!(id[1,1].domain, index)
    #end

    push!(model.constraints, SeaPearl.AllDifferent(id, trailer))
    model.adhocInfo = Dict([("n", n), ("m", m)])
    return model
end

"""
    solve_eternity2(input_file; order=[1,2,3,4], variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic(), limit=nothing)

Solve the SeaPearl model for to the eternity2 problem, using SeaPearl.MinDomainVariableSelection heuristique
and  SeaPearl.AllDifferent and SeaPearl.TableConstraint, and the function model_AIS

# Arguments
- `n::Int`: dimension
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
- 'order' : Vector, giving the order of edges for the IO manager. example : [1,4,2,3] means it is given as [up,down,left,right]
- 'limit' : Int, giving the number of solutions after which it will stop searching. if nothing given, it will lookk for all the solutions
- 'modeling' : modeling funciton, fast by default
"""
function solve_eternity2(input_file; order=[1,2,3,4], variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic(), limit=1,modeling=model_eternity2_fast)
    model = modeling(input_file; order=order,variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection, valueSelection=valueSelection)
    return model
end

"""
    solve_eternity2(model::SeaPearl.CPModel;benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())

Solve the SeaPearl model for to the eternity2 problem, using an existing model (from model_eternity2)

# Arguments
- `model::SeaPearl.CPModel`: model (from model_AIS)
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
"""
function solve_eternity2(model::SeaPearl.CPModel; variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection, valueSelection=valueSelection)
    return model
end

"""
    outputFromSeaPearl(model::SeaPearl.CPModel; optimality=false)

Return the results of the eternity2 problem as type OutputDataEternity, taking the solved model as argument

# Arguments
- `model::SeaPearl.CPModel`: needs the model to be already solved (by solve_eternity2)
"""
function outputFromSeaPearl(model::SeaPearl.CPModel)
    solutions = model.statistics.solutions
    nb_sols = length(solutions)
    n = model.adhocInfo["n"]
    m = model.adhocInfo["m"]
    orientation = Array{Int, 4}(undef, nb_sols, n,m,5)
    for (ind,sol) in enumerate(solutions)
        for i in 1:n, j in 1:m
            orientation[ind, i, j,1]= sol["id_"*string(i)*string(j)]
            orientation[ind, i, j, 2] = sol["u_"*string(i)*string(j)]
            orientation[ind, i, j, 3] = sol["r_"*string(i)*string(j)]
            orientation[ind, i, j, 4] = sol["d_"*string(i)*string(j)]
            orientation[ind, i, j, 5] = sol["l_"*string(i)*string(j)]
        end
    end
    return OutputDataEternityII(nb_sols, orientation)
end

function outputFromSeaPearl_fast_orientation(model::SeaPearl.CPModel; optimality=false)
    solutions = model.statistics.solutions
    nb_sols = length(solutions)
    n = model.adhocInfo["n"]
    m = model.adhocInfo["m"]
    orientation = Array{Int, 4}(undef, nb_sols, n,m,5)
    for (ind,sol) in enumerate(solutions)
        for i in 1:n, j in 1:m
            orientation[ind, i, j,1]= sol["id_"*string(i)*string(j)]
            orientation[ind, i, j, 2] = sol["src_h"*string(i)*string(j)]
            orientation[ind, i, j, 3] = sol["src_v"*string(i)*string(j+1)]
            orientation[ind, i, j, 4] = sol["src_h"*string(i+1)*string(j)]
            orientation[ind, i, j, 5] = sol["src_v"*string(i)*string(j)]
        end
    end
    return OutputDataEternityII(nb_sols, orientation)
end

"""
    print_eternity2(sol::Array{Int,3})

print a solution

# Arguments
- 'sol::Array{Int,3}' : one solution from OutputDataEternityII
"""
function print_eternity2(sol::Array{Int,3})
    id = sol[:,:,1]
    u = sol[:,:,2]
    r = sol[:,:,3]
    d = sol[:,:,4]
    l = sol[:,:,5]
    n = size(id,1)
    m = size(id,2)
    for k in 1:9*m
        print("-")
    end
    println()
    for i in 1:n
        print("|")
        for j in 1:m
            printstyled("   "*string(u[i,j],pad=2)*"   ", color=u[i,j])
            print("|")
        end
        println()
        print("|")
        for j in 1:m
            printstyled(string(l[i,j],pad=2),color=l[i,j])
            printstyled(" "*string(id[i,j],pad=2)*" ")
            printstyled(string(r[i,j],pad=2),color=r[i,j])
            print("|")
        end
        println()
        print("|")
        for j in 1:m
            printstyled("   "*string(d[i,j],pad=2)*"   ", color=d[i,j])
            print("|")
        end
        println()
        for k in 1:9*m
            print("-")
        end
        println()
    end
end

"""
    print_eternity2(output::OutputDataEternityII;limit=1)

print a certain number of solution

# Arguments
- 'output::OutputDataEternityII' : output from outputFromSeaPearl2
- 'limit::Int' : number of solutions to print
"""
function print_eternity2(output::OutputDataEternityII;limit=1)
    println("Let's show "*string(limit)*" solutions :")
    for i in 1:limit
        sol = output.orientation[i,:,:,:]
        println("Solution "*string(i))
        print_eternity2(sol)
    end
end

"""
    print_eternity2(model::SeaPearl.CPModel;limit=1)

print a certain number of solution

# Arguments
- 'model::SeaPearl.CPModel' : model solved
- 'limit::Int' : number of solutions to print
"""
function print_eternity2(model::SeaPearl.CPModel;limit=1)
    output = outputFromSeaPearl(model)
    print_eternity2(output;limit=limit)
end

function print_eternity2_fast_orientation(model::SeaPearl.CPModel;limit=1)
    output = outputFromSeaPearl_fast_orientation(model)
    print_eternity2(output;limit=limit)
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
abstract type AbstractHESN <: AbstractEchoStateNetwork end

struct HESN{T, S<:AbstractArray{T}, I, B, F, N, M, U} <: AbstractHESN
    res_size::I
    in_size::I
    out_size::I
    train_data::S
    physics_model_data::S
    prior_model::M
    u0::U
    tspan::Tuple
    datasize::I
    dt::T
    alpha::T
    nla_type::N
    activation::F
    W::S
    W_in::S
    states::S
    extended_states::B
end

function HESN(W::AbstractArray{T},
        train_data::AbstractArray{T},
        prior_model::Any,
        u0::AbstractArray{T},
        tspan::Tuple,
        datasize::Int,
        W_in::AbstractArray{T};
        activation::Any = tanh,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    #Create physics data with one extra step ahead extra
    trange = collect(range(tspan[1], tspan[2], length = datasize))
    dt = trange[2]-trange[1]
    tsteps = push!(trange, dt + trange[end])
    tspan_new = (tspan[1], dt+tspan[2])
    physics_model_data = prior_model(u0, tspan_new, tsteps)
    physics_informed_data = vcat(train_data, physics_model_data[:, 1:end-1])

    in_size = size(physics_informed_data, 1)
    out_size = size(train_data, 1)
    res_size = size(W, 1)

    if size(W_in, 1) != res_size
        throw(DimensionMismatch("size(W_in, 1) must be equal to size(W, 1)"))
    elseif size(W_in, 2) != in_size
        throw(DimensionMismatch("size(W_in, 2) must be equal to in_size"))
    end

    states = states_matrix(W, W_in, physics_informed_data, alpha, activation, extended_states)

    return HESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(activation),
        typeof(nla_type),
        typeof(prior_model),
        typeof(u0)}(res_size, in_size, out_size, train_data, physics_model_data,
        prior_model, u0, tspan, datasize, dt, alpha, nla_type, activation, W, W_in, states, extended_states)
end

"""
    HESNpredict(esn::AbstractLeakyESN, predict_len::Int, prior_data::AbstractArray{Float64}, model_size::Int, W_out::AbstractArray{Float64})

Return prediction the a starting after the training time using HESN model.
"""

function HESNpredict(hesn::AbstractHESN,
    predict_len::Int,
    W_out::AbstractArray{Float64})

    output = zeros(Float64, hesn.out_size, predict_len)
    x = hesn.states[:, end]
    predict_tsteps = [hesn.tspan[2]+hesn.dt]
    [append!(predict_tsteps, predict_tsteps[end]+hesn.dt) for i in 1:predict_len]
    tspan_new = (hesn.tspan[2]+hesn.dt, predict_tsteps[end])
    u0 = hesn.physics_model_data[:, end]
    physics_prediction_data = hesn.prior_model(u0, tspan_new, predict_tsteps)[:, 2:end]

    for i=1:predict_len
        x_new = nla(hesn.nla_type, x)
        x_new = vcat(x_new, physics_prediction_data[:, i]) #<-- append states of prior model at current tstep
        out = (W_out*x_new) #<-- prediction w/ reservoir state & prior model values given solved W_out
        output[:, i] = out
        out = vcat(out, physics_prediction_data[:, i]) #<-- append states of prior model for input of next prediction
        x = leaky_fixed_rnn(hesn.activation, hesn.alpha, hesn.W, hesn.W_in, x, out)
    end
    return output
end


text:
module TestRunner

export  get_tests_structure, run_all_tests,
        get_tests_structure_as_json, get_tests_results_as_json,
        children, line, name, result, details, stacktrace,
        RESULT, test_success, test_failure, test_error, test_pending, test_not_run

using FactCheck
using Compat
abstract TestStructureNode

@enum RESULT test_success test_failure test_error test_pending test_not_run

type RootNode <: TestStructureNode
  children::Vector{TestStructureNode}
end

type FactsCollectionNode <: TestStructureNode
    line::Int
    name::AbstractString
    children::Vector{TestStructureNode}
end

type ContextNode <: TestStructureNode
    line::Int
    name::AbstractString
    children::Vector{TestStructureNode}
end

type FactNode <: TestStructureNode
    line::Int
    name::AbstractString
    result::RESULT
    details::AbstractString
    stacktrace::AbstractString
end

FactNode(line::Int, name::AbstractString, result::RESULT) = FactNode(line, name, result, "")
FactNode(line::Int, name::AbstractString, result::RESULT, details::AbstractString) = FactNode(line, name, result, details, "")
FactNode(line::Int, name::AbstractString, result::FactCheck.Success) = FactNode(line, name, test_success, _get_details(line, result))
FactNode(line::Int, name::AbstractString, result::FactCheck.Pending) = FactNode(line, name, test_pending, _get_details(line, result))
FactNode(line::Int, name::AbstractString, result::FactCheck.Failure) = FactNode(line, name, test_failure, _get_details(line, result))
FactNode(line::Int, name::AbstractString, result::FactCheck.Error)   = FactNode(line, name, test_error,   _get_details(line, result), _get_stacktrace(result))
FactNode(line::Int, name::AbstractString) = FactNode(line, name, test_not_run)

_get_details(line::Int, result::FactCheck.Result) = _replace_line_number(line, sprint(show, result)) |> _strip_ascii_escape_codes
_replace_line_number(line::Int, s::AbstractString) = replace(s, r"line:(-?\d+)", "line:$line")
_strip_ascii_escape_codes(s::AbstractString) = replace(s, r"\x1b[^m]*m", "")
_get_stacktrace(result::FactCheck.Error) = sprint(showerror, result.err, result.backtrace)

get_tests_structure(test_file_path::AbstractString) = test_file_path |> _get_file_content |> _get_tests_structure |> RootNode

_get_file_content(test_file_path::AbstractString) = test_file_path |> readstring |> (content -> "begin\n" * content * "\nend") |> parse |> _fixLineNumbers!

function _fixLineNumbers!(expression_tree_node::Expr)
  map!(_fixLineNumbers!, expression_tree_node.args)
  expression_tree_node
end
_fixLineNumbers!(line_number_node::LineNumberNode) = LineNumberNode(line_number_node.file, line_number_node.line - 1)
_fixLineNumbers!(other) = other
_fixLineNumbers(expression_tree_node::Expr) = _fixLineNumbers!(deepcopy(expression_tree_node))

function _get_tests_structure(expression_tree_node::Expr, tests_results::Vector{FactCheck.Result} = FactCheck.Result[] , line::Int = 0)
    results = Vector{TestStructureNode}()
    for child_node in expression_tree_node.args
        line = get(_try_get_line(child_node), line)
        _procces_expression_tree_node(child_node, results, tests_results, line)
    end
    results
end

_try_get_line(node::LineNumberNode) = Nullable{Int}(node.line)
_try_get_line(_) = Nullable{Int}()

function _procces_expression_tree_node(tree_node::Expr, results::Vector{TestStructureNode}, tests_results::Vector{FactCheck.Result}, line::Int)
  length(tree_node.args)>0 || return

  get_name(expr::Expr) = length(expr.args)>2?expr.args[3]:""
  get_children(expr::Expr) = _get_tests_structure(expr.args[2], tests_results, line)
  get_fact_result!(rs::Vector{FactCheck.Result}) =  isempty(rs)?test_not_run:pop!(rs)
  dump_pending_test_result!(rs::Vector{FactCheck.Result}) = isempty(rs) || pop!(rs)
  create_custom_pending_details(expr::Expr) = "$(expr.args[2].args[1]) $(expr.args[2].head) $(expr.args[2].args[2])"

  if tree_node.args[1] == :facts
      node = FactsCollectionNode(line, get_name(tree_node), get_children(tree_node))
      push!(results, node)
  elseif tree_node.args[1] == :context
      node = ContextNode(line, get_name(tree_node), get_children(tree_node))
      push!(results, node)
  elseif  tree_node.args[1] == Symbol("@fact")
     node =  FactNode(line, get_name(tree_node), get_fact_result!(tests_results))
     push!(results, node)
  elseif tree_node.args[1] == Symbol("@pending")
      dump_pending_test_result!(tests_results)
      node = FactNode(line, get_name(tree_node), test_pending, create_custom_pending_details(tree_node) )
      push!(results, node)
  else
      append!(results, _get_tests_structure(tree_node, tests_results, line))
  end
end

_procces_expression_tree_node(_...) = nothing

function run_all_tests(test_file_path::AbstractString)
  results = _get_results(test_file_path)
  content = _get_file_content(test_file_path)
  _get_tests_structure(content, results) |> RootNode
end

function _get_results(test_file_path::AbstractString)
  original_allresults  = FactCheck.allresults
  FactCheck.clear_results()
  original_STDOUT = STDOUT
  (out_read, out_write) = redirect_stdout()

  evalfile(test_file_path)

  close(out_write)
  close(out_read)
  redirect_stdout(original_STDOUT)
  results = copy(FactCheck.allresults)
  empty!(FactCheck.allresults)
  append!(FactCheck.allresults, original_allresults)
  results |> reverse
end

children(_::FactNode) = Vector{TestStructureNode}()
children(node::TestStructureNode) =  node.children
line(_::RootNode) = 0
line(node::TestStructureNode) = node.line
name(_::RootNode) = "root"
name(node::TestStructureNode) = node.name
result(node::FactNode) = node.result
details(node::FactNode) = node.details
stacktrace(node::FactNode) = node.stacktrace

include("TestRunnerJSON.jl")

get_tests_structure_as_json(test_file_path::AbstractString) = test_file_path |> get_tests_structure |> json

get_tests_results_as_json(test_file_path::AbstractString) = test_file_path |> run_all_tests |> json

end # module


text:
__precompile__()

module PointClouds

import Base:
    show, keys, haskey,
    getindex, setindex!,
    vcat, length, endof

import NearestNeighbors:
    knn, inrange

using NearestNeighbors
using StaticArrays

export PointCloud,
    # Point cloud data access
    positions,
    normals,
    # Spatial indexing
    knn,
    inrange,
    # Data handling
    split_cloud,
    # Adding columns
    add_normals!

include("cloud.jl")

end # module


text:
using Sigma

# Create a uniformly distributed random variable
x = uniform(0,1)

# Sample 100 values distributed according to x
rand(x, 100)

# Compute the probability that x^2 is greater than 0.6
prob(x^2 > 0.6)

# Create another uniformly distributed random variable (independent from y)
y = uniform(0,1)

# Compute the probability that x^2 is greater than 0.6, under the condition that
# x + y < 1
prob(x^2 > 0.6, x + y < 1)

# Sample from x conditioned on x + y < 1
rand(x, x + y < 1)


text:
pwd() == joinpath(@__DIR__, "bin") && cd(@__DIR__) # allow starting app from bin/ dir

using Venti
Venti.main()


text:
@testset "generate" begin
    @testset "rand_graph" begin
        n, m = 10, 20
        m2 = m ÷ 2
        x = rand(3, n)
        e = rand(4, m2)

        g = rand_graph(n, m, ndata=x, edata=e, graph_type=GRAPH_T)
        @test g.num_nodes == n
        @test g.num_edges == m
        @test g.ndata.x === x
        if GRAPH_T == :coo
            s, t = edge_index(g)
            @test s[1:m2] == t[m2+1:end]
            @test t[1:m2] == s[m2+1:end]
            @test g.edata.e[:,1:m2] == e
            @test g.edata.e[:,m2+1:end] == e
        end

        g = rand_graph(n, m, bidirected=false, seed=17, graph_type=GRAPH_T)
        @test g.num_nodes == n
        @test g.num_edges == m

        g2 = rand_graph(n, m, bidirected=false, seed=17, graph_type=GRAPH_T)
        @test edge_index(g2) == edge_index(g)
    end

    @testset "knn_graph" begin
        n, k = 10, 3
        x = rand(3, n)
        g = knn_graph(x, k; graph_type=GRAPH_T)
        @test g.num_nodes == 10
        @test g.num_edges == n*k
        @test degree(g, dir=:in) == fill(k, n)
        @test has_self_loops(g) == false

        g = knn_graph(x, k; dir=:out, self_loops=true, graph_type=GRAPH_T)
        @test g.num_nodes == 10
        @test g.num_edges == n*k
        @test degree(g, dir=:out) == fill(k, n)
        @test has_self_loops(g) == true

        graph_indicator = [1,1,1,1,1,2,2,2,2,2]
        g = knn_graph(x, k; graph_indicator, graph_type=GRAPH_T)
        @test g.num_graphs == 2
        s, t = edge_index(g)
        ne = n*k÷2
        @test all(1 .<= s[1:ne] .<= 5)
        @test all(1 .<= t[1:ne] .<= 5)
        @test all(6 .<= s[ne+1:end] .<= 10)
        @test all(6 .<= t[ne+1:end] .<= 10)
    end
end


text:
using NumericalMethodsforEngineers
using Test

a = [16.; 4.; 5.; 8.; -4.; 22.]
kdiag = [1; 3; 6]

amat = fromskyline(kdiag, a)

b = [4., 2., 5.]

nfix = 1
no = [2;]
val = [5.0;]
penalty = 1.0e20

a[kdiag[no]] .+= penalty
b[no] = a[kdiag[no]] .* val

sparin!(a, kdiag)
spabac!(a, b, kdiag)
d = copy(b)
#------------------------------------

a = amat
b = [4., 2., 5.]

a[2, 2] += 1.0e20
b[2] = a[2, 2] * 5.0

a = sparse(a)

F = cholesky(a)
y = F.L \ b
c = F.U \ y

@test round.(c, digits=4) == round.(a \ b, digits=4)
@test round.(c, digits=4) == round.(d, digits=4)


text:
#=
    AbstractGUI
      AbstractContext, AbsractContent
         AbstractConstrain, AbstractConstruct, AbstractContemper

                                                   AbstractGUI
                            AbstractContext                                     AbstractContent
                   AbstractConstrainContext                             AbstractConstrainContent
                   AbstractConstructContext                             AbstractConstructContent
                   AbstractContemperContext                             AbstractContemperContent

=#

#=
    An Algebra of Content and Context

    Jeffrey Sarnoff 2019-Oct-09

    Context pervades, informs, reflects  Content
    Content suffuses, regards, mapsfrom  Context

    meld   ::  [A, A] => A
    whole   :: [B, B] => B
    pervade :: [A, B] => B
    inform  :: [A, B] => B
    suffuse :: [B, A] => A
    regard  :: [B, A] => A

=#

mutable struct HtmlElement
    htmlelem::Symbol
    hyperelem::Function
    cssclass::MaybeTuple
    cssid::MaybeString
    content::MaybeContent
    context::MaybeContext

    function HtmlElement(elem::Symbol, cssclass::MaybeNTuple{N,String}=nothing; id::MaybeString=nothing,
                         content::MaybeContent=nothing, context::MaybeContext=nothing) where {N}
       helem = HyperHtml[elem]
       return new(elem, helem, cssclass, id, content, context)
    end

    function HtmlElement(elem::Symbol, cssclass::String; id::MaybeString=nothing,
                         content::MaybeContent=nothing, context::MaybeContext=nothing)
       helem = HyperHtml[elem]
       return new(elem, helem, (cssclass,), id, content, context)
    end
end

Base.string(io::IO, x::HtmlElement) = string(io, hyperhtml(x))
Base.string(x::HtmlElement) = string(hyperhtml(x))

render(io::IO, x::HtmlElement) = Hyperscript.render(io, Hyperscript.Pretty(hyperhtml(x)))
render(x::HtmlElement) = Hyperscript.render(stdout, Hyperscript.Pretty(hyperhtml(x)))

hyperhtml(x::Nothing) = nothing
hyperhtml(x::Context{T}) where {T} = hyperhtml(x.value)
hyperhtml(x::Content{T}) where {T} = hyperhtml(x.value)
hyperhtml(x::String) = x

function hyperhtml(x::HtmlElement)
    context = hyperhtml(x.context)
    content = hyperhtml(x.content)
    result =
        if !isnothing(x.cssclass)
            classes = join(x.cssclass, " ")
            if !isnothing(x.cssid)
                x.hyperelem(classes, :id=>x.cssid)
            else
                x.hyperelem(classes)
            end
        elseif !isnothing(x.cssid)
            x.hyperelem(:id=>x.cssid)
        else
            x.hyperelem()
        end
    if !isnothing(context)
        if !isnothing(content)
            result(content,context)
        else
            result(context)
        end
    elseif !isnothing(content)
        result(content)
    else
        result
    end
end

#=
asection = HtmlElement(:section, ("asection",)) # HtmlElement(:section, GUI.section_, ("asection",), nothing)
adiv=HtmlElement(:div, "adiv") # HtmlElement(:div, GUI.div_, ("adiv",), nothing)
alink=HtmlElement(:a, ("alink")) # HtmlElement(:a, GUI.a_, ("alink",), nothing)

hyperhtml(asection) #

hyperhtml(adiv) #

hyperhtml(alink) #

hyperhtml(adiv)(hyperhtml(alink)) #

hyperhtml(adiv)(hyperhtml(alink)("url")) #
url

hyperhtml(adiv),hyperhtml(alink)("url") # (

, url)
hyperhtml(asection)(hyperhtml(adiv),hyperhtml(alink)("url")) #

url

bdiv = HtmlElement(:div, "bdiv", context=Context(adiv)) # HtmlElement(:div, GUI.div_, ("bdiv",), nothing, nothing, Context{HTML_Element}(HTML_Element(:div, GUI.div_, ("adiv",), nothing)))
hyperhtml(bdiv) #

hyperhtml(bdiv)("b") #

b

cdiv = HtmlElement(:div, "cdiv", content=GUI.Content("inside cdiv"), context=GUI.Context(alink)) # HtmlElement(:div, GUI.div_, ("ddiv",), nothing, Content{HtmlElement}(HtmlElement(:a, GUI.a_, ("alink",), nothing, nothing, nothing)), Context{String}("abcd"))
hyperhtml(cdiv) #
inside cdiv

ddiv = HtmlElement(:div, "ddiv", content=GUI.Content(alink), context=GUI.Context("abcd")) # HtmlElement(:div, GUI.div_, ("ddiv",), nothing, Content{HtmlElement}(HtmlElement(:a, GUI.a_, ("alink",), nothing, nothing, nothing)), Context{String}("abcd"))
hyperhtml(ddiv) #
abcd

=#


text:
## Delta rule
function delta_rule!(Q, s::Int, a::Int, r, α)
    Q[a,s] += α*(r - Q[a,s])
end

function delta_rule!(Q, a::Int, r::T, α) where T <: Real
    Q[a] += α*(r - Q[a])
end

function delta_rule!(Q, s::Int, r::T, α) where T <: AbstractArray
    q = @view(Q[:,s])
    @. q += α * (r - q)
end

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
# plot recipes, see https://docs.juliaplots.org/latest/recipes/ and
# https://github.com/JuliaPlots/RecipesBase.jl

mutable struct WithBeam
    system::Vector{<:Element}
    beam::Beam
end

# type recipe, e.g. for `plot(WithBeam(system, beam))`
@recipe f(::Type{WithBeam}, data::WithBeam) = begin
    seriestype := :shape
    linecolor --> color(data.beam.λ)
    fillcolor --> color(data.beam.λ)
    fillalpha --> 0.1
    label --> string(data.beam.λ)
    aspect_ratio --> :equal
    xlabel --> "Distance along Beam Axis"
    ylabel --> "1/e^2 Boundary"
    # choose a very large number of discretizations to have some
    # chance of approximating minimum waist radii decently
    ds = discretize(data.system, 200)
    N = length(ds) + 1
    Tw = typeof(float(spotsize(data.beam)))
    Tz = typeof(float(location(data.beam)))
    ws = Vector{Tw}(undef, N)
    zs = Vector{Tz}(undef, N)
    ws[1] = spotsize(data.beam)
    zs[1] = location(data.beam)
    beam = data.beam
    for i = 1:length(ds)
        beam = transform(ds[i], beam)
        ws[i+1] = spotsize(beam)
        zs[i+1] = location(beam)
    end
    xs, ys = vcat(zs, reverse(zs)), vcat(ws, (-1.0) .* reverse(ws))
    [(xs[i], ys[i]) for i in 1:length(xs)]
end

const _λ_hue_data = [
    380e-9 260 # "far" UV
    440e-9 250 # violet
    485e-9 180 # blue
    523e-9 120 # blue-green
    572e-9 60 # yellow
    620e-9 20
    700e-9 0 # (infra) red
]
const hue_of_λ = Interpolations.LinearInterpolation(
    _λ_hue_data[:,1], _λ_hue_data[:,2]
)

"""

Return a (very approximately) correct color from the package
[Colors](http://docs.juliaplots.org/latest/colors/) for a given
wavelength (specified in SI units, i.e. in the range from
approximately `380e-9` to `700e-9`).

"""
function color(λ)
    # very rough transformation of λ to a color via package Colors
    λmin = 380e-9
    λmax = 700e-9
    if λ < λmin
        # ultraviolet (far spectral violet, dark)
        return Colors.HSL(260.0, 1.0, 0.25)
    elseif λ > λmax
        # infrared (red, dark)
        return Colors.HSL(0, 1.0, 0.25)
    end
    color = Colors.HSL(hue_of_λ(λ), 1.0, 0.5)
end

# user recipe, e.g. for `plot(system, beam)`
@recipe f(system::Vector{<:Element}, beam::Beam) =
    WithBeam(system, beam)


text:
using DataFrames
using CSV
using Plots
#%%
data_jacobi = open("jacobi_p.dat")
df_jacobi=CSV.read(data_jacobi, header=false, delim=" ")
x_jacobi=df_jacobi[1]
y_jacobi=df_jacobi[2]
z_jacobi=reshape(df_jacobi[3],(101,101))

data_gs = open("gauss_seidel_p.dat")
df_gs=CSV.read(data_gs, header=false, delim=" ")
x_gs=df_gs[1]
y_gs=df_gs[2]
z_gs=transpose(reshape(df_gs[3],(101,101)))

data_s = open("sor_p.dat")
df_s=CSV.read(data_s, header=false, delim=" ")
x_s=df_s[1]
y_s=df_s[2]
z_s=transpose(reshape(df_s[3],(101,101)))

data_sol = open("sol_p.dat")
df_sol=CSV.read(data_sol, header=false, delim=" ")
x_sol=df_sol[1]
y_sol=df_sol[2]
z_sol=transpose(reshape(df_sol[3],(101,101)))
length(x_sol)/101
x=y_sol[1:101]
y=x
#%%
plot(x,y,z_jacobi, fill=true, label="jacobi")
xlabel!("x")
ylabel!("y")
title!("jacobi")
plot!(dpi=600)
savefig("jacobi_p.png")
# xaxis!(:log10)
# yaxis!(:log10)
#%%

plot(x,y,z_gs, fill=true, label="gauss_seidel")
xlabel!("x")
ylabel!("y")
title!("gauss seidel")
plot!(dpi=600)
savefig("gs_p.png")

#%%

plot(x,y,z_s, fill=true, label="sor")
xlabel!("x")
ylabel!("y")
title!("sor")
plot!(dpi=600)
savefig("sor_p.png")

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
function print_eigenvector(vec, precision)
	s = size(vec, 1)
	for i=1:s
		print(round(vec[i], precision), "\t")
	end
	println()
end

function eigenvalues_ham(hamiltonian, precision, print, Braak_form)
	size_ham = size(hamiltonian)
	λ = eigfact!(view(hamiltonian, 1:size_ham[1]-2, 1:size_ham[2]-2))
	states = Array{String}(length(λ[:values]))
	#print_matrix(hamiltonian, 2)
	if print
		println("\n Eigenvalues:")
		for i=1:size(λ[:values], 1)
			@printf("#%d |     ", i)
			@printf("%.2f \t current state |", round(λ[:values][i], precision))
			str = ""
			if Braak_form
				if i%2 == 0
					str = str * "+ "* join(convert(Int, floor((i-1)/2))) * ">"
				else
					str = str * "- " * join(convert(Int, floor(i/2))) * ">"
				end
			else
				if i%2 == 0
					str = str * "e, "* convert(Int, floor((i-1)/2))*">"
				else
					str = str * "g, " * convert(Int, floor(i/2)) *">"
				end
			end
			println(str)
			states[i] = str
		end
		#println("Extra λ |  ", round(hamiltonian[size_ham[1], size_ham[2]], precision))

		println("\n Eigenvectors:")
		for i=1:size(λ[:vectors], 2)
			@printf("#%d |     ", i)
			print_eigenvector(λ[:vectors][:, i], precision)
		end
	end
	return λ#(λ, states)
end

function sort_eigs(eigs, unsorted_eigs, ε)
	incr_ratio = eigs[:, 2] - eigs[:, 1]
	new_incr = unsorted_eigs - eigs[:, 2]
	for i=1:size(new_incr, 1)
		if abs(new_incr[i]-incr_ratio[i]) > ε ### out acceptable bounds...
			#info("qui")
			for j=i+1:size(new_incr, 1)
				if abs(unsorted_eigs[j] - eigs[i,2]-incr_ratio[i]) <= ε ##each other very close
					temp = unsorted_eigs[j]
					unsorted_eigs[j] = unsorted_eigs[i]
					unsorted_eigs[i] = temp
					new_incr = unsorted_eigs - eigs[:, 2]
				end
			end
		end
	end
	return unsorted_eigs
end

function get_corresponding_eigenvalue(hamiltonian, λ, str)
	size_ham = size(hamiltonian)
	level = eval(parse(str[search(str,',')+1:end-1]))
  # Grabbig qubit state
  index = -1
  if str[search(str,'|')+1] == 'g'
    index = level*2+1
  else
   index = level*2+2
  end
	test_hamiltonian = deepcopy(hamiltonian)
	test_hamiltonian[index, index] = 0.0
	λ_test = eigfact!(sub(test_hamiltonian, 1:size_ham[1]-1, 1:size_ham[2]-1))
	for x=1: size(λ, 1)
		found = false
		for y =1: size(λ_test[:values], 1)
			if round(λ[x], 5)==round(λ_test[:values][y, 1], 5)
				found = true
			end
		end
		found == false && return λ[x]
	end
	return NaN
end


text:
typeof(2 + 3) #> Int64
:(2 + 3) #> :(2 + 3)
typeof(:(2 + 3)) #> Expr
exp = :(2 + 3)
names(exp)
#3-element Array{Symbol,1}:
# :head
# :args
# :typ
exp.head #> :call
exp.args #> 3-element Array{Any,1}:  :+ 2 3
exp.typ #> Any
2 + 3 == + (2, 3) #> true

exp2 = :(2 + a * b - c)
exp2.args #> 3-element Array{Any,1}: :- :(2 + a * b) :c
:(2 + a * b).args #> 3-element Array{Any,1}:  :+ 2  :(a * b)
:(a * b).args #> 3-element Array{Any,1}: :* :a :b

x = 5 #> 5
:x #> :x

quote
   a = 42
   b = a^2
   a - b
end
# quote  # none, line 2:
#     a = 42 # line 3:
#     b = a^2 # line 4:
#     a - b
# end
:(a = 42; b = a^2; a - b)
# returns:
# quote
#     a = 42
#     begin
#         b = a^2
#         a - b
#     end
# end

type Expr
  head::Symbol
  args::Array{Any,1}
  typ
end

dump(:(2 + 2))
# Expr
#   head: Symbol call
#   args: Array(Any,(3,))
#     1: Symbol +
#     2: Int64 2
#     3: Int64 2
#   typ: Any

dump(:(2 + a * b - c))
# Expr
#   head: Symbol call
#   args: Array(Any,(3,))
#     1: Symbol -
#     2: Expr
#       head: Symbol call
#       args: Array(Any,(3,))
#         1: Symbol +
#         2: Int64 2
#         3: Expr
#           head: Symbol call
#           args: Array(Any,(3,))
#           typ: Any
#       typ: Any
#     3: Symbol c
#   typ: Any

symbol('\'') #> :'
symbol("'") #> :'

max_stars_repo_path: Module 1/Chapter07/Old/expressions.jl

max_issues_repo_path: Module 1/Chapter07/Old/expressions.jl

max_forks_repo_path: Module 1/Chapter07/Old/expressions.jl


text:
module Specs

import ..CryptoGroups: order, a, b, cofactor, modulus, bitlength

include("spec.jl")
include("field_specs.jl")

include("legacy.jl")
include("primitives.jl")

include("curve_constants.jl")
include("modp_constants.jl")

export ECP, EC2N, Koblitz, MODP

end


text:
const evecURL = "http://api.eve-central.com/api/"

evec_filename(url::AbstractString) = joinpath(cacheDir, replace(replace(url, evecURL, ""), urlReplaceKeys, c->urlReplace[c[1]]) * ".xml")

function get_evec(endpoint::AbstractString; query=Dict{AbstractString, Any}(), timeoutHours::Float64 = 6.0)
	url = evecURL * endpoint * "?" * Requests.format_query_str(query)
	fileName = evec_filename(url)
	data = mtime(fileName) + timeoutHours * 60 * 60 > time()? xml_read(fileName) : nothing
	if data == nothing
		headers = Dict("User-Agent"=>userAgent)
		info("GET from $url")
		resp = Requests.get(url; headers=headers)
		if resp.status != 200
			info("Response status: $(resp.status)")
			return nothing
		end
		xml_write(fileName, resp.data)
		xdoc = LightXML.parse_string(resp.data)
		data = LightXML.root(xdoc)
	end
	return data
end


text:
using LinearAlgebra
using Statistics
using LinearAlgebra.BLAS
import StatsBase.percentile

const Mat{T} = AbstractArray{T, 2}
const Tup{T} = Tuple{Mat{T}, Mat{T}}

abstract type AbstractGBM{T, I} end

"""
General Boltzmann Machine, including parameter matrix (element type T)
and a mask matrix (element type I) for decimation.
"""
mutable struct GBM{T, I} <: AbstractGBM{T, I}
    W::Matrix{T}
    b::Matrix{T}
    W_mask::Array{I, 2}
    PLrange::Dict
end

"""
    GBM(n_unit::Int)

Construct GBM with number of unit `n_unit`.

 * T - type of GBM parameters [weights and biases;] (default: Float64)
 * I - type of Mask-matrix for decimation (default: Int64)
"""
function GBM(n_unit::Int)
    GBM{Float64, Int64}(
        rand(n_unit, n_unit),
        rand(n_unit, 1),
        ones(Int64, (n_unit, n_unit)) - Matrix(I, n_unit, n_unit),
        Dict("min" => 0.0, "max" => 0.0)
    )
end

## gradient calculation

"""
    calc_grad(gmb::GBM, X::AbstractArray, ctx::Dict)

Calculate gradeint of log-pseudo-likelihood with respect to GBM parameter.

# References
- A. Hyvärinen, Neural Comput. **18**, 2283 (2006)
"""
function calc_grad(gbm::GBM{T}, X::Mat{T}, ctx::Dict)::Tup{T} where T

    alpha = T(ctx[:alpha])
    beta = T(ctx[:beta])

    dev = X' - tanh.(beta * (gbm.W * X' .+ gbm.b))
    dW = @getarray(ctx, :dW_buf, size(gbm.W), similar(gbm.W))

    # same as: dW = beta .* dev * X / size(X, 1)
    gemm!('N', 'N', beta / size(X, 1), dev, X, T(0.0), dW)
    symmetrize!(dW)
    axpy!(-alpha, sign.(gbm.W), dW)

    # gradient for bias params.
    db = beta * mean(dev, dims=2)

    return (dW, db)
end

"""
    update(gbm::GBM, grads::Tuple, ctx::Dict)

Update GBM parameters with SGD method.

Update weight parameters of GBM using provided tuple `dtheta = (dW, db)`
of parameter gradients. Optimal method is Stochastic Gradient Descent [SGD].
"""
function update(gbm::GBM{T}, grads::Tup{T}, ctx::Dict) where T
    dW, db = grads
    lr = T(ctx[:lr])

    # same as: gbm.W += lr * dW
    axpy!(lr, dW, gbm.W)
    axpy!(lr, db, gbm.b)

    symmetrize!(gbm.W)
    gbm.W .*= gbm.W_mask
end

## Pseudo likelihood calculation

"""
    pseudo_likelihood(gbm::GBM, X::AbstractArray, ctx::Dict)

Calculate pseudo-likelihood[PL] with spin-example data `X`.
"""
function pseudo_likelihood(gbm::GBM{T}, X::Mat{T}, ctx::Dict)::T where T
    alpha = T(ctx[:alpha])
    beta = T(ctx[:beta])
    H = tr(X * gbm.W * X') + sum(X * gbm.b)
    tmp = sum(log.(cosh.(beta * (gbm.W * X' .+ gbm.b))))
    PL = (beta * H - tmp) / size(X, 1)
    return PL - alpha * sum(abs.(gbm.W))
end

"""
    fit(gbm::GBM, X::AbstractArray; delta=1e-1 [, opts])

Fit the weight parameters of GBM to data `X`.

# Arguments
- `delta::Float`: Threshold for thinning out parameters.

Options that can be provided in the `opts` dictionary:

* `:alpha` - Decay rate of the L1 regularization term. (default: 0.01)
* `:beta` - Inverse temperature. (default: 1.0)
* `:lr` - Learning rate. (default: 1.0)
* `:epochs` - Maximum permissible number of step in SGD. (default: 20)
* `:n_iter` - Number of learning repetitions. (default: 100)

# References
- E. Aurell & M. Ekeberg, Phys. Rev. Lett. **108**, 090201 (2012).

# Examples
```jldoctest
julia> using InverseIsing

julia> samples = [1 -1 -1;] # Spin configuration.

julia> model = GBM(3) # Set the number of units.

julia> fit(model, samples)

julia> W = infer(model); output = decode(W)

julia> output
OrderedCollections.OrderedDict{Tuple{Int64,Int64},Int64} with 3 entries:
  (1, 2) => -1
  (1, 3) => -1
  (2, 3) => 1
```
    The above example means that the interaction between (1, 2) and (1, 3) is
    antiferromagnetic bond and only (2, 3) is ferromagnetic bond.

"""
function fit(gbm::GBM{T}, X::Mat; delta::T=1e-1, opts...) where {T, I}

    ctx = copy(Dict{Any, Any}(opts))
    option_checker!(ctx)
    epochs = ctx[:epochs]
    max_step = ctx[:n_iter]

    X_train = trans_type(T, X)
    variable_check(X_train)

    for n in 1:max_step

        old_likelihood = pseudo_likelihood(gbm, X_train, ctx)
        update(gbm, calc_grad(gbm, X_train, ctx), ctx)

        if n % epochs == 0
            gbm.W_mask[abs.(gbm.W) .< delta] .= 0
            gbm.W_mask .*= gbm.W_mask'
        end

        if isapprox(pseudo_likelihood(gbm, X_train, ctx), old_likelihood)
            break
        end
    end
    return ctx
end

### sub-utility functions ###
"""
    coef(gbm::GBM)

Get weight marix of a trained GBM.

# Examples
```jldoctest
julia> model = GBM(3)

julia> coef(model)
3×3 Array{Float64,2}:
 0.86902   0.4058     0.568231
 0.700701  0.0892165  0.399372
 0.224095  0.234413   0.380341
```
"""
function coef(gbm::GBM{T}) where T
    W = copy(gbm.W)
    # translate -0.0 to 0.0
    W[W .== -T(0.0)] .= T(0.0)
    return W
end

"""
    weights(gbm::GBM)

Same as coef. See also: [`coef`](@ref)
"""
weights = coef

"""
    infer(gbm::GBM)

Infer GBM parameter matrix using trained weight.

# Examples
```jldoctest
julia> using InverseIsing

julia> model = GBM(3)

julia> infer(model)
3×3 Array{Int64,2}:
 1  1  1
 1  1  1
 1  1  1

julia> fit(model, [1 -1 -1;])

julia> infer(model)
3×3 Array{Int64,2}:
  0  -1  -1
 -1   0   1
 -1   1   0
```
"""
function infer(gbm::GBM{T, I}) where {T, I}
    pred = heaviside(coef(gbm))
    return convert(Array{I}, pred)
end

max_stars_repo_name: oilneck/InverseIsing

max_issues_repo_name: oilneck/InverseIsing

max_forks_repo_name: oilneck/InverseIsing


text:
module AppendedTriangulationsTests

using Test
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using Gridap.Integration
using LinearAlgebra: ⋅
using Gridap.CellData
using Gridap.FESpaces

domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(domain,partition)

ncells = num_cells(model)
nin = ceil(Int,2*ncells/3)
cell_to_mask = fill(false,ncells)
cell_to_mask[1:nin] .= true

grid = get_grid(model)

trian_in = RestrictedTriangulation(grid,cell_to_mask)
trian_out = RestrictedTriangulation(grid,collect(Bool, .! cell_to_mask))

trian = lazy_append(trian_out,trian_in)
test_triangulation(trian)

order = 1
quad = CellQuadrature(trian,2*order)
quad_in = CellQuadrature(trian_in,2*order)
quad_out = CellQuadrature(trian_out,2*order)

q = get_coordinates(quad)
w = get_weights(quad)
@test isa(q,AppendedArray)
@test isa(w,AppendedArray)

V = TestFESpace(model=model,valuetype=Float64,order=order,reffe=:Lagrangian,conformity=:H1)

u(x) = x[1]+x[2]

_v = interpolate(u,V)
v = restrict(_v,trian)

e = u - v
el2 = sqrt(sum(integrate(e*e,trian,quad)))
@test el2 < 1.0e-8

_dv = get_cell_basis(V)
dv = restrict(_dv,trian)

cellmat =  integrate(∇(dv)⋅∇(dv),trian,quad)
@test isa(cellmat,AppendedArray)
@test isa(cellmat.a,CompressedArray)
@test isa(cellmat.b,CompressedArray)

end # module


text:
# Returns index of next zero (or error if none is found)
# pointer must point to first byte where the search begins
# This can be SIMD'd but it's way fast anyway.
function next_zero(p::Ptr{UInt8}, i::UInt32, lastindex::UInt32)::UInt32
    while i ≤ lastindex
        unsafe_load(p) === 0x00 && return i
        i += UInt32(1)
        p += 1
    end
    throw(LibDeflateError(UNTERMINATED_NULL_STRING))
end

struct SizedMemory
    ptr::Ptr{UInt8}
    len::UInt
end

SizedMemory(x) = SizedMemory(pointer(x), sizeof(x))
Base.length(x::SizedMemory) = x.len
Base.pointer(x::SizedMemory) = x.ptr

"Check if there are any 0x00 bytes in a block of memory"
function any_zeros(mem::SizedMemory)
    i = UInt(1)
    while i ≤ length(mem)
        unsafe_load(pointer(mem), i % Int) === 0x00 && return true
        i += UInt(1)
    end
    return false
end

# +---+---+---+---+==================================+
# |SI1|SI2|  LEN  |... LEN bytes of subfield data ...|
# +---+---+---+---+==================================+
"""
    GzipExtraField

Data structure for gzip extra data. Fields:

* `tag::NTuple{2, UInt8}` two-byte tag
* `data::Union{Nothing, UnitRange{UInt32}}` location of subfield data in original vector,
or `nothing` if empty.
"""
struct GzipExtraField
    tag::Tuple{UInt8, UInt8} # (SI1, SI2)
    data::Union{Nothing, UnitRange{UInt32}}
end

# The pointer points to the first byte of the first field
function parse_fields!(
	fields::Vector{GzipExtraField},
	ptr::Ptr{UInt8},
	index::UInt32,
remaining_bytes::UInt16
)
	empty!(fields)
    while !iszero(remaining_bytes)
        field = parse_extra_field(ptr, index, remaining_bytes)
        push!(fields, field)

        # We zero the range field on an empty subfield, so we take
        # that possibility into account
        field_len = field.data === nothing ? UInt16(0) : length(field.data) % UInt16
        total_len = field_len + UInt16(4)
        remaining_bytes -= total_len
        ptr += total_len
        index += total_len
    end
    return fields
end

# The pointer points to the first byte of the first field
function parse_fields(ptr::Ptr{UInt8}, index::UInt32, remaining_bytes::UInt16)
	parse_fields!(GzipExtraField[], ptr, index, remaining_bytes)
end

# The pointer points to the first byte of the extra fields
function parse_extra_field(ptr::Ptr{UInt8}, index::UInt32, remaining_bytes::UInt16)
    remaining_bytes < 4 && throw(LibDeflateError(EXTRA_DATA_TOO_LONG))
    s1 = unsafe_load(ptr)
    s2 = unsafe_load(ptr + 1)
    iszero(s2) && throw(LibDeflateError(EXTRA_DATA_INVALID)) # not allowed
    field_len = ltoh(unsafe_load(Ptr{UInt16}(ptr + 2)))
    field_len + 4 > remaining_bytes && throw(LibDeflateError(EXTRA_DATA_TOO_LONG))

    # If the field is empty, we use a Nothing to convey that
    range = if iszero(field_len)
        nothing
    else
        index+UInt32(4):index+UInt32(4)-UInt32(1)+field_len
    end
    return GzipExtraField((s1, s2), range)
end

"""
    is_valid_extra_data(ptr::Ptr{UInt8}, remaining_bytes::UInt16)

Check if the chunk of bytes pointed to by `ptr` and `remaining_bytes`
onward represent valid gzip metadata for the "extra" field.
"""
function is_valid_extra_data(ptr::Ptr{UInt8}, remaining_bytes::UInt16)
    while !iszero(remaining_bytes)
        # First four bytes: S1, S2, field_len
        remaining_bytes < 4 && return false
        # S2 must not be zero
        iszero(unsafe_load(ptr + 1)) && return false
        field_len = ltoh(unsafe_load(Ptr{UInt16}(ptr + 2)))
        remaining_bytes < field_len + 4 && return false
        ptr += 4 + field_len
        remaining_bytes -= UInt16(4) + field_len
    end
    return true
end

"""
    GzipHeader

Struct representing a gzip header. It has the following fields:
* `mtime::UInt32`: Modification time of file
* `filename::Union{Nothing, UnitRange{32}}` index of filename in header
* `comment::Union{Nothing, UnitRange{32}}` index of comment in header
* `extra::Union{Nothing, Vector{GzipExtraField}}` Extra gzip fields, if applicable.
"""
struct GzipHeader
    mtime::UInt32
    filename::Union{Nothing, UnitRange{UInt32}}
    comment::Union{Nothing, UnitRange{UInt32}}
    extra::Union{Nothing, Vector{GzipExtraField}}
end

"""
    parse_gzip_header

See `unsafe_parse_gzip_header`
"""
function parse_gzip_header(
    in::Vector{UInt8},
    max_len::UInt,
    extra_data::Union{Vector{GzipExtraField}, Nothing}=nothing
)
    return unsafe_parse_gzip_header(pointer(in), max_len, extra_data)
end

"""
    unsafe_parse_gzip_header(
        in_ptr::Ptr{UInt8}, max_len::UInt,
        extra_data::Union{Vector{GzipExtraField}, Nothing}=nothing
    )

Parse the input data, returning a `GzipHeader` object, or erroring.
The parser will not read more than `max_len` bytes. If a vector of gzip
extra data is passed, it will not allocate a new vector, but overwrite the given one.
"""
function unsafe_parse_gzip_header(
    in_ptr::Ptr{UInt8},
    max_len::UInt, # maximum length of header
    extra_data::Union{Vector{GzipExtraField}, Nothing}=nothing
)
    # header is at least 10 bytes
    max_len > 9 || throw(LibDeflateError(INPUT_DATA_TOO_SHORT))
    # Bytes 1 - 10. Check first four bytes, skip rest
    # +---+---+---+---+---+---+---+---+---+---+
    # |ID1|ID2|CM |FLG|     MTIME     |XFL|OS | (more-->)
    # +---+---+---+---+---+---+---+---+---+---+
    ptr = in_ptr - UInt(1) # zero-indexed pointer
    header = ltoh(unsafe_load(Ptr{UInt32}(ptr + 1)))
    header & 0x00ffffff == 0x00088b1f || throw(LibDeflateError(BAD_HEADER))
    FLAG_HCRC =    !iszero(header & 0x02000000)
    FLAG_EXTRA =   !iszero(header & 0x04000000)
    FLAG_NAME =    !iszero(header & 0x08000000)
    FLAG_COMMENT = !iszero(header & 0x10000000)
    mtime = ltoh(unsafe_load(Ptr{UInt32}(ptr + 5)))

    # 32-bit index because this library only works with 32-bit buffers anyway
    # (skip MTIME, XFL, OS), they're not useful anyway
    index = UInt32(11)

    extra = nothing
    if FLAG_EXTRA
        # +---+---+=================================+
        # | XLEN  |...XLEN bytes of "extra field"...| (more-->)
        # +---+---+=================================+
        extra_len = ltoh(unsafe_load(Ptr{UInt16}(ptr + index)))
        extra_vector = if extra_data === nothing
            GzipExtraField[]
        else
            extra_data
        end
        extra = parse_fields!(extra_vector, ptr + index + 2, index + UInt32(2), extra_len)
        index += extra_len + UInt32(2)
        index > max_len && throw(LibDeflateError(UNTERMINATED_NULL_STRING))
    end

    filename = nothing
    if FLAG_NAME
        # +=========================================+
        # |...original file name, zero-terminated...| (more-->)
        # +=========================================+
        zero_pos = next_zero(ptr + index, index, max_len % UInt32)
        zero_pos > max_len && throw(LibDeflateError(UNTERMINATED_NULL_STRING))
        filename = index:zero_pos - one(UInt32)
        index = zero_pos + one(UInt32)
    end

    # Skip comment
    comment = nothing
    if FLAG_COMMENT
        zero_pos = next_zero(ptr + index, index, max_len % UInt32)
        zero_pos > max_len && throw(LibDeflateError(UNTERMINATED_NULL_STRING))
        comment = index:zero_pos - one(UInt32)
        index = zero_pos + one(UInt32)
    end

    # Verify header CRC16, if present
    if FLAG_HCRC
        # Lower 16 bits of crc32 up to, not including, this index
        # +---+---+
        # | CRC16 |
        # +---+---+
        crc_obs_16 = unsafe_crc32(ptr + one(UInt), index - one(UInt)) % UInt16
        crc_exp_16 = ltoh(unsafe_load(Ptr{UInt16}(ptr + index)))
        crc_obs_16 == crc_exp_16 || throw(LibDeflateError(HEADER_CRC16_CHECKSUM_DOES_NOT_MATCH))
        index += UInt32(2)
        index > max_len && throw(LibDeflateError(UNTERMINATED_NULL_STRING))
    end

    return (index - UInt32(1), GzipHeader(mtime, filename, comment, extra))
end

"""
    GzipDecompressResult

Result of `LibDeflate`'s gzip decompression on byte vector. The fields `filename`
and `comment` specify the location of gzip feature data in the input vector.
When not applicable (e.g. the `comemnt` field is not applicable for gzip files
without the `FCOMMENT` flag), these fields are zeroed out.

It has the following fields:
* `len::UInt32` length of decompressed data
* `mtime::UInt32` timestamp of original data, or zero
* `filename::UnitRange{UInt32}` location of filename (or zero)
* `comment::UnitRange{UInt32}` location of gzip comment (or zero)
* `extra::Union{Nothing, Vector{GzipExtraField}}` gzip extra data (or `nothing`)
"""
struct GzipDecompressResult
    len::UInt32 # length of decompressed data
    header::GzipHeader
end

"""
    gzip_decompress!(::Decompressor, out::Vector{UInt8}, in::Vector{UInt8},
        max_len=typemax(Int))

Gzip decompress the input data into `out`, and resize `out` to fit. Throws an error
if `out` would be resized to larger than `max_len`. Returns a `GzipDecompressResult`.

See also: [`unsafe_gzip_decompress!`](@ref)
"""
function gzip_decompress!(
    decompressor::Decompressor,
    out_data::Vector{UInt8},
    in_data::Union{Vector{UInt8}, String, SubString{String}},
    extra_data::Union{Vector{GzipExtraField}, Nothing}=nothing;
    max_len::Integer=typemax(Int),
)
    GC.@preserve in_data out_data begin
        result = unsafe_gzip_decompress!(
            decompressor, out_data, UInt(max_len),
            pointer(in_data), sizeof(in_data) % UInt, extra_data
        )
    end

    length(out_data) == result.len || resize!(out_data, result.len)
    return result
end

"""
    unsafe_gzip_decompress!(decompressor::Decompressor, outdata::Vector{UInt8},
        max_outlen::UInt, in_ptr::Ptr{UInt8}, len::UInt,
        extra_data::Union{GzipExtraField, Nothing})

Use the `Decompressor` to decompress gzip data at `in_ptr` and `len` bytes forward
into `outdata`. If there is not enough room at `outdata`, resize `outdata`, except
if it would be bigger than `max_outlen`, in that case throw an error.
If `extra_data` is not `nothing`, reuse the passed vector

Return a `GzipDecompressResult`

See also: [`gzip_decompress!`](@ref)
"""
function unsafe_gzip_decompress!(
    decompressor::Decompressor,
    out_data::Vector{UInt8},
    max_outlen::UInt,
    in_ptr::Ptr{UInt8},
    len::UInt,
    extra_data::Union{Vector{GzipExtraField}, Nothing}=nothing
)
    # We need to have at least 2 + 4 + 4 bytes left after header
    nonheader_min_len = 2 + 4 + 4

    # First decompress header
    header_len, header = unsafe_parse_gzip_header(in_ptr, len - nonheader_min_len, extra_data)

    # Skip to end to check crc32 and data len
    # +---+---+---+---+---+---+---+---+
    # |     CRC32     |     ISIZE     | END OF FILE
    # +---+---+---+---+---+---+---+---+

    compressed_len = len - UInt(8) - header_len
    uncompressed_size = ltoh(unsafe_load(Ptr{UInt32}(in_ptr + len - UInt(4))))
    uncompressed_size > max_outlen && throw(LibDeflateError(OUTPUT_DATA_TOO_LONG))
    length(out_data) < uncompressed_size && resize!(out_data, uncompressed_size)

    # Now DEFLATE decompress
    unsafe_decompress!(Base.HasLength(), decompressor, pointer(out_data), uncompressed_size,
    in_ptr + header_len, compressed_len)

    # Check for CRC checksum and validate it
    crc_exp = ltoh(unsafe_load(Ptr{UInt32}(in_ptr + len - UInt(8))))
    crc_obs = unsafe_crc32(pointer(out_data), uncompressed_size % Int)
    crc_exp == crc_obs || throw(LibDeflateError(PAYLOAD_CRC132_CHECKSUM_DOES_NOT_MATCH))

    GzipDecompressResult(uncompressed_size, header)
end

"Computes maximal output length of a gzip compression"
function max_out_len(input_len::UInt, comment_len::UInt, filename_len::UInt, extra_len::UInt16, header_crc::Bool)
    # Taken from libdeflate source code
    # with slight modifications
    static = 10 + 8 + 9 # header + footer + padding

    n_chunks = max(cld(input_len, 10000), 1)
    len = static + input_len + 5 * n_chunks # 5 byte overhead per chunk
    len += comment_len + !iszero(comment_len) # incl. null byte
    len += filename_len + !iszero(filename_len) # incl. null byte
    len += extra_len + 2 * !iszero(extra_len) # incl. 2-byte leader
    len += 2*header_crc
    return len
end

"""
    gzip_compress!(
        compressor::Compressor,
        output::Vector{UInt8},
        input::Vector{UInt8},
        comment::Union{String, SubString{String}, Vector{UInt8}, Nothing}=nothing,
        filename::Union{String, SubString{String}, Vector{UInt8}, Nothing}=nothing,
        extra::Union{String, SubString{String}, Vector{UInt8}, Nothing}=nothing,
        header_crc::Bool=false
    )

Gzip compress `input` into `output` and resizing output to fit. Returns `output`.

Adds optional data `comment`, `filename`, `extra` and `header_crc`:
* `comment` and `filename` must not include the byte `0x00`.
* `extra` must be at most `typemax(UInt16)` bytes long.
* `header_crc` is true, add the header CRC checksum.

See also: [`unsafe_gzip_compress!`](@ref)
"""
function gzip_compress!(
    compressor::Compressor,
    output::Vector{UInt8},
    input::Union{String, SubString{String}, Vector{UInt8}};
    comment::Union{String, SubString{String}, Vector{UInt8}, Nothing}=nothing,
    filename::Union{String, SubString{String}, Vector{UInt8}, Nothing}=nothing,
    extra::Union{String, SubString{String}, Vector{UInt8}, Nothing}=nothing,
    header_crc::Bool=false
)
    # Resize output to maximal possible length
    maxlen = max_out_len(
        sizeof(input) % UInt,
        comment === nothing ? UInt(0) : sizeof(comment) % UInt,
        filename === nothing ? UInt(0) : sizeof(filename) % UInt,
        extra === nothing ? UInt16(0) : sizeof(extra) % UInt16,
        header_crc
    )
    # We add 8 extra bytes to make sure Libdeflate don't error due to off-by-one errors
    resize!(output, maxlen + 8)

    GC.@preserve output input comment filename extra begin
        mem_comment = comment === nothing ? nothing : SizedMemory(comment)
        mem_filename = filename === nothing ? nothing : SizedMemory(filename)
        mem_extra = extra === nothing ? nothing : SizedMemory(extra)

        n_bytes = unsafe_gzip_compress!(
            compressor,
            pointer(output),
            sizeof(output) % UInt,
            pointer(input),
            sizeof(input) % UInt,
            mem_comment,
            mem_filename,
            mem_extra,
            header_crc
        )
    end

    resize!(output, n_bytes % UInt)
    return output
end

"""
    unsafe_gzip_compress!(
        compressor::Compressor,
        out_ptr::Ptr{UInt8}, out_len::UInt,
        in_ptr::Ptr{UInt8}, in_len::UInt,
        comment::Union{SizedMemory, Nothing},
        filename::Union{SizedMemory, Nothing},
        extra::Union{SizedMemory, Nothing},
        header_crc::Bool
    )::Int

Use the `Compressor` to gzip compress input `in_len` bytes from `in_ptr`, into `out_ptr`.
If the resulting gzip data could be longer than `out_len`, throw an error.
Optionally, include gzip comment, filename or extra data. These should be represented
by a `SizedMemory` object, or `nothing` if it should be skipped.
* `comment` and `filename` must not include the byte `0x00`.
* `extra` must be at most `typemax(UInt16)` bytes long.
If `header_crc` is true, add the header CRC checksum.

Returns the number of bytes written to `out_ptr`.

See also: [`gzip_compress!`](@ref)
"""
function unsafe_gzip_compress!(
    compressor::Compressor,
    out_ptr::Ptr{UInt8},
    out_len::UInt,
    in_ptr::Ptr{UInt8},
    in_len::UInt,
    comment::Union{SizedMemory, Nothing},
    filename::Union{SizedMemory, Nothing},
    extra::Union{SizedMemory, Nothing},
    header_crc::Bool,
)
    # Check output len is long enough
    max_out_len(
        in_len,
        comment === nothing ? UInt(0) : length(comment),
        filename === nothing ? UInt(0) : length(filename),
        if extra === nothing
            UInt16(0)
        else
            # No more than typemax(UInt16) bytes for extra field
            length(extra) > typemax(UInt16) && throw(LibDeflateError(EXTRA_DATA_TOO_LONG))
            length(extra) % UInt16
        end,
        header_crc
    ) > out_len && throw(LibDeflateError(OUTPUT_DATA_TOO_LONG))

    # Write first four bytes - magix number, compression type, flags
    header = 0x00088b1f
    if comment !== nothing
        # Check for absence of zero byte
        any_zeros(comment) && throw(LibDeflateError(UNTERMINATED_NULL_STRING))
        header |= 0x10000000
    end
    if filename !== nothing
        # Check for absence of zero byte
        any_zeros(filename) && throw(LibDeflateError(UNTERMINATED_NULL_STRING))
        header |= 0x08000000
    end
    if extra !== nothing
        # Validate extra data
        is_valid_extra_data(pointer(extra), length(extra) % UInt16) || throw(LibDeflateError(EXTRA_DATA_INVALID))
        header |= 0x04000000
    end
    header = ifelse(header_crc, header | 0x02000000, header)
    ptr = out_ptr - 1
    unsafe_store!(Ptr{UInt32}(ptr + 1), htol(header))

    # Add system time (take lower 32 bits if it overflows)
    unsafe_store!(Ptr{UInt32}(ptr + 5), htol(unsafe_trunc(UInt32, time())))

    # Add system (unknown) and XFL (zero)
    unsafe_store!(Ptr{UInt16}(ptr + 9), htol(0x00ff))

    index = UInt(11)

    # Add in extra data
    if extra !== nothing
        unsafe_store!(Ptr{UInt16}(ptr + index), htol(length(extra) % UInt16))
        unsafe_copyto!(ptr + index + 2, pointer(extra), length(extra))
        index += UInt(2) + length(extra)
    end

    # Add in filename
    if filename !== nothing
        unsafe_copyto!(ptr + index, pointer(filename), length(filename))
        index += length(filename) + one(UInt) # null byte
        unsafe_store!(ptr + index - 1, 0x00)
    end

    # Add in comment
    if comment !== nothing
        unsafe_copyto!(ptr + index, pointer(comment), length(comment))
        index += length(comment) + one(UInt) # null byte
        unsafe_store!(ptr + index - 1, 0x00)
    end

    # Add in CRC16
    if header_crc
        header_crc = unsafe_crc32(ptr + one(UInt), index - one(UInt)) % UInt16
        unsafe_store!(Ptr{UInt16}(ptr + index), htol(header_crc))
        index += UInt(2)
    end

    # Add in compressed data
    remaining_outdata = out_len - index + 1 - 8 # tail
    n_compressed = unsafe_compress!(compressor, ptr + index, remaining_outdata, in_ptr, in_len)
    index += n_compressed

    # Add in crc32 of uncompressed data
    crc = unsafe_crc32(in_ptr, in_len)
    unsafe_store!(Ptr{UInt32}(ptr + index), htol(crc))
    index += 4

    # Add in isize (uncompressed size)
    unsafe_store!(Ptr{UInt32}(ptr + index), htol(in_len % UInt32))
    return (index + 3) % Int # 4 bytes isize - off-by-one
end


text:
if !@isdefined(interpret)
    include("boolfuck.jl")
end

bffibo = brainfucktoboolfuck(",>+>>>>++++++++++++++++++++++++++++++++++++++++++++>++++++++++++++++++++++++++++++++<<<<<<[>[>>>>>>+>+<<<<<<<-]>>>>>>>[<<<<<<<+>>>>>>>-]<[>++++++++++[-<-[>>+>+<<<-]>>>[<<<+>>>-]+<[>[-]<[-]]>[<<[>>>+<<<-]>>[-]]<<]>>>[>>+>+<<<-]>>>[<<<+>>>-]+<[>[-]<[-]]>[<<+>>[-]]<<<<<<<]>>>>>[++++++++++++++++++++++++++++++++++++++++++++++++.[-]]++++++++++<[->-<]>++++++++++++++++++++++++++++++++++++++++++++++++.[-]<<<<<<<<<<<<[>>>+>+<<<<-]>>>>[<<<<+>>>>-]<-[>>.>.<<<[-]]<<[>>+>+<<<-]>>>[<<<+>>>-]<<[<+>-]>[<+>-]<<<-]");
input = join(Char[10])

@time interpret(bffibo, input)
@time f = compile(bffibo)
@time f(input)

# Baseline - equivalent function in Julia
# 0.018314 seconds (53.82 k allocations: 3.214 MiB, 99.81% compilation time)
# 0.000014 seconds (31 allocations: 1.453 KiB)

# Time taken for fibonacci(10) in Boolfuck

# First interpreter
# 17.954638 seconds (6.12 M allocations: 200.450 MiB, 0.24% gc time, 4.89% compilation time)
# 16.866883 seconds (3.55 M allocations: 54.287 MiB, 0.05% gc time)

# Second interpreter avoiding excess memory use
# 15.507475 seconds (861.41 k allocations: 47.568 MiB, 0.06% gc time, 2.14% compilation time)
# 15.666484 seconds (90 allocations: 168.312 KiB)

# Third interpreter parser groups multiple +s and <>s
# 5.436422 seconds (147.41 k allocations: 8.672 MiB, 0.16% gc time, 1.89% compilation time)
# 5.266540 seconds (87 allocations: 168.172 KiB)

# First compiler - compile boolfuck to Julia
# Oneshot - 14.027388 seconds (3.92 M allocations: 183.405 MiB, 0.51% gc time, 97.50% compilation time)
# Compile only - 0.522731 seconds (509.95 k allocations: 30.057 MiB, 22.37% compilation time)
# Cold start - 13.808873 seconds (3.75 M allocations: 173.838 MiB, 1.44% gc time, 99.96% compilation time)
# Warm start - 0.004169 seconds (66 allocations: 7.188 KiB)

# Second compiler - change from gotos to while loops in compiled code

using Plots
function plotgraph()
    plot(
        [@elapsed f(String(Char[i])) for i = 1:50],
        label = "Compiled",
        xaxis = "n",
        yaxis = ("Time to calculate fib(n) in seconds, log10 scale", :log10)
    )
    plot!(
        [@elapsed interpret(bffibo, String(Char[i])) for i = 1:10],
        label = "Interpreted"
    )
    plot!(plot_title = "Time taken to calculate Fibonacci numbers")
    savefig("boolfuck-plot.svg")
end

max_stars_repo_name: stellartux/Esolangs

max_issues_repo_name: stellartux/Esolangs

max_forks_repo_name: stellartux/Esolangs


text:
using LinearAlgebra, SparseArrays,Statistics,DataFrames,DataFramesMeta

#To install Laplacians you need to put in Julia's Prompt "add Laplacians#master"
#For some reason Pkg.add() doesn't work with this.
#Pkg.add("MATLAB")
#using Laplacians
#using Pkg
#include(string(Pkg.dir("Laplacians") , "/src/matlabSolvers.jl"))
using Test
using Parameters
using Random

#Full Network
const y_full = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.464866 ; 0.730622 ; 0.619239; 0.753429; 0.0863208; 0.372084 ;  0.958089]
const id_full = [1; 1; 2;2 ; 3; 4;4 ;5;5 ;6;6]
const firmid_full = [1;2;1;1;1;1;2;2;2;2;3]
const obs_id_full = collect(1:11)

#Pruned Network
const y_p = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.464866 ; 0.730622 ; 0.619239; 0.753429; 0.0863208]
const id_p = [1; 1; 2;2 ; 3; 4;4 ;5;5]
const firmid_p = [1;2;1;1;1;1;2;2;2]
const obs_id_p = collect(1:9)

#Drop-Single Network
const y_lo = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.730622 ; 0.619239; 0.753429; 0.0863208]
const id_lo = [1; 1; 2;2 ; 3;3 ;4;4]
const firmid_lo = [1;2;1;1;1;2;2;2]
const obs_id_lo1 = [1;2;3;4;6;7;8;9]
const obs_id_lo = collect(1:8)

#Leave-Out Network Matchid
const match_id = [1; 2 ; 3;3 ; 4; 5; 6; 6]
@testset "PruneNetwork" begin
    @test find_connected_set(y_full,id_full,firmid_full; verbose=true) == (obs_id = obs_id_full , y = y_full , id = id_full, firmid = firmid_full )

    @test prunning_connected_set(y_full,id_full,firmid_full, obs_id_full; verbose=true) == (obs_id = obs_id_p , y = y_p , id = id_p, firmid = firmid_p )

    @test drop_single_obs(y_p,id_p,firmid_p,obs_id_p) == (obs_id = obs_id_lo1 , y = y_lo , id = id_lo, firmid = firmid_lo )

end

@testset "Movers & Matches" begin
    @test  compute_matchid(firmid_lo,id_lo) == [1; 2 ; 3;3 ; 4; 5; 6; 6]
    @test  compute_movers(id_lo, firmid_lo) == (movers = Bool[1,1,0,0,1,1,0,0], T= [2, 2, 2, 2, 2, 2, 2, 2])
end

Xtest = sparse([1;2;3;4;5;6;7;8;1;3;4;5],[1;1;2;2;3;3;4;4;5;5;5;5],[1.0;1.0;1.0;1.0;1.0;1.0;1.0;1.0;-1.0;-1.0;-1.0;-1.0])
xxtest = Xtest'*Xtest
compute_sol = approxcholSolver(xxtest;verbose=true)
settings_exact = Settings(leverage_algorithm = ExactAlgorithm(), person_effects=true, cov_effects=true)
settings_JLA = Settings(leverage_algorithm = JLAAlgorithm(num_simulations=3), person_effects=true, cov_effects=true)

@testset "Clustering" begin
    @test check_clustering(id_lo).nnz_2 == 12

    @test check_clustering(match_id).nnz_2 ==10

    @test index_constr(obs_id_lo,id_lo, obs_id_lo) == [obs_id_lo  obs_id_lo  obs_id_lo   obs_id_lo]

    @test index_constr(obs_id_lo,id_lo, match_id) == [obs_id_lo  obs_id_lo  [1;2;3;3;4;5;6;6]   obs_id_lo]

end

Random.seed!(1234)
rademacher1 = rand(1,8)  .> 0.5
rademacher1 = rademacher1 - (rademacher1 .== 0)

@testset "Solvers" begin
    @test compute_sol( [Xtest[1,:]...] ;verbose=true) ≈ [ 0.25, -0.5, -0.25, 0.0, -0.5]
    @test compute_sol([rademacher1*Xtest...];verbose=true) ≈ [ 1.0, 0.0, 1.0, -1.0, 0.0]
end

Random.seed!(1234)
const K=0
@testset "LambdaMatrices" begin
    @test    do_Pii(Xtest,obs_id_lo) == sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.75,0.75,0.5,0.5,0.75,0.75,0.5,0.5] )
    @test    eff_res(settings_exact.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_exact).Lambda_P == sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.75,0.75,0.5,0.5,0.75,0.75,0.5,0.5] )
    @test    eff_res(settings_exact.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_exact).Lambda_B_fe == sparse( [1,2,5,6], [1,2,5,6], [0.5,0.5,0.5,0.5],8,8)
    @test    eff_res(settings_exact.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_exact).Lambda_B_pe == sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.625,0.625,0.375,0.375,0.625,0.625,0.375,0.375] )
    @test    eff_res(settings_exact.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_exact).Lambda_B_cov == sparse( [1,2,5,6], [1,2,5,6], [-0.25,-0.25,-0.25,-0.25],8,8)

    # For now, just test that this executes
    @test eff_res(settings_JLA.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_JLA).Lambda_P != nothing

    # Original test
    @test_broken    eff_res(settings_JLA.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_JLA).Lambda_P  ==  sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.666667,0.666667,0.5,0.5,0.666667,0.666667,0.5,0.5] )
end

@testset "LeaveOut" begin
    #  @test  leave_out_estimation(y,id,firmid,nothing, settings_exact) ≈ [-0.0042 , 0.0019, 0.0037]
    #  @test  leave_out_estimation(y,id,firmid,nothing, settings_JLA) ≈ [-0.0042 , 0.0019, 0.0037]
end


text:
@testset "74.search-a-2d-matrix.jl" begin
    @test search_matrix([[1, 3, 5, 7], [10, 11, 16, 20], [23, 30, 34, 60]], 3) == true
    @test search_matrix([[1, 3, 5, 7], [10, 11, 16, 20], [23, 30, 34, 60]], 13) == false
end


text:
module Toolz
include("InternalList.jl")
using .List: head, tail, cons, nil, reverse, linkedlist

export ($), isCapitalized

($)(f, a) = f(a)
isCapitalized(s :: AbstractString) :: Bool = !isempty(s) && isuppercase(s[1])
end