text:
#using Reduce, SpecialFunctions, Base.Cartesian, Sobol
#Reduce.Rational(false)

@generated function marginalize(x, ::Tuple{Val{N}}) where N
    quote
        out = zero(eltype(x))
        @nexprs $N i -> @inbounds out += x[i] *  marginal(Val{i-1}())
        out
    end
end

@generated function marginalize(x, ::Tuple{Val{N1},Val{N2}}) where {N1, N2}
    quote
        out = zero(eltype(x))
        @nexprs $N1 i -> begin
            @nexprs $N2 j -> begin
                out += x[j,i] * marginalize(Val{i-1}()) * marginalize(Val{j-1}())
            end
        end
        out
    end
end
@generated function marginalize(x, ::Tuple{Val{N1},Val{N2},Val{N3}}) where {N1, N2, N3}
    quote
        out = zero(eltype(x))
        @nexprs $N1 i -> begin
            @nexprs $N2 j -> begin
                @nexprs $N3 k -> begin
                    out += x[k,j,i] * marginalize(Val{i-1}()) * marginalize(Val{j-1}()) * marginalize(Val{k-1}())
                end
            end
        end
        out
    end
end
@generated function marginalize(x, ::Tuple{Val{N1},Val{N2},Val{N3},Val{N4}}) where {N1, N2, N3, N4}
    quote
        out = zero(eltype(x))
        @nexprs $N1 i -> begin
            @nexprs $N2 j -> begin
                @nexprs $N3 k -> begin
                    @nexprs $N4 l -> begin
                        out += x[l,k,j,i] * marginalize(Val{i-1}()) * marginalize(Val{j-1}()) * marginalize(Val{k-1}()) * marginalize(Val{l-1}())
                    end
                end
            end
        end
        out
    end
end

@generated function marginal_gen(::Val{N}, ::Type{T} = Float32) where {N,T}
    f = rcall(:(int( x^$N * exp(-x^2/2), x, -4.8, 4.8)))
    :(convert(T,$f))
end
@generated function marginal(::Val{N}, ::Type{T} = Float32) where {N,T}
    marginal_gen(Val{N}(), T)
end
function marginal(i::I, ::Type{T} = Float32) where {I <: Integer, T}
    out = @nif 16 d -> (i == 2(d-1)) d -> ( marginal(Val{2(d-1)}(), T) ) d -> ( zero(T) )
    out
end

@generated function uni_poly(x, coef_v, ::Val{N}) where N
    Nm1 = N - 1
    quote
        @inbounds out = coef_v[1]
        @nexprs $Nm1 i -> @inbounds out += x^i * coef_v[i+1]
        out * exp(-x^2/2)
    end
end

function coef(i::I, ::Type{T} = Float32) where {I <: Integer, T}
    out = @nif 160 d -> (i <= 2(d-1)) d -> ( 1/marginal(Val{2(d-1)}(), T) ) d -> ( zero(T) )
    out
end
function coef(p, T = Float32)
    c = similar(p, T)
    for i ∈ eachindex(c)
        c[i] = prod(x -> coef(x, T), p[i])
    end
    c
end

"""
a is the Array (nested VectorOfArray) form of the polynomial.
v is a vector of coefficients, typically created by solving a system of equations.
"""
poly(::Val{0}, n::NTuple{ns,I} = (32,), ::Type{S} = Float32) where {ns,S,I<:Integer} = fill(zero(S), n...)
@generated function poly(::Val{N}, n::NTuple{ns,I} = (32,), ::Type{S} = Float32) where {N,ns,S,I<:Integer}
    Nm1 = N - ns
    v_last = Symbol("v_", Nm1)
    quote
        v_0 = fill(zero($S), n...)#In preparation for the deprecation of zeros.
        @nextract $Nm1 v d -> [ v_{d-1} ]
        $v_last
    end
end

@generated function poly(v_coefs::Vector{T}, poly_exponents::Vector{NTuple{p,Cuint}}, n::NTuple{ns,I} = (32,), ::Type{S} = Float32) where {ns,p,T,S,I<:Integer}
    N = p - ns
    v_last = Symbol("v_", N)
    quote
        v_0 = poly(Val{$p}(), n, S)
        for i ∈ eachindex(v_coefs)
            ex = poly_exponents[i]
            @nextract $N v d -> begin
                for j ∈ length(v_{d-1}):ex[d]
                    push!(v_{d-1}, poly(Val{$N - d}(), n, S))
                end
                v_{d-1}[ ex[d] + 1 ]
            end
            (@nref $ns $v_last d -> ex[d + $N] + 1) = v_coefs[i]
        end
        v_0
    end
end
#const PolyDist{1,T} = Vector{T}
#const PolyDist{2,T} = Vector{Vector{T}}
struct PolyDist{N,T,V} <: AbstractArray{N,T}
    v::V
end
function PolyDist(β::Vector{T}, poly::Vector{NTuple{N,I}}) where {T,N,I<:Integer}

end

@generated function mss4(x, y, ::Val{S}, ::Val{N}) where {S,N}
    quote
        T = eltype(x)
        out = zero(T)
        @nexprs $N i -> @inbounds out += sum( vload( Vec{$S,T}, x, $S*(i-1)+1 ) * vload( Vec{$S,T}, y, $S*(i-1)+1 ) )
        out
    end
end
@generated function mssuneven(x,y,::Val{S},::Val{N},::Val{R}) where {S,N,R}
    rn = 1 + S*N
    quote
        T = eltype(x)
        out = zero(T)
        @nexprs $N i -> @inbounds out += sum( vload( Vec{$S,T}, x, $S*(i-1)+1 ) * vload( Vec{$S,T}, y, $S*(i-1)+1 ) )
        @inbounds out += sum( vload( Vec{$R,T}, x, $rn ) * vload( Vec{$R,T}, y, $rn ) )
        out
    end
end
@generated function mss(x, y, ::Val{N}, ::Val{S}) where {S,N}
    n, r = divrem(N,S)
    r == 0 ? :(mss4(x, y, Val{$S}(), Val{$n}())) : :(mssuneven(x,y, Val{$S}(), Val{$n}(), Val{$r}()))
end

Base.length(x::PolyDist) = length(x.v)
@generated function Base.size(x_0::PolyDist{N}) where N
    quote
        @nextract $N x i -> x_{i-1}[1]
        @ntuple $N i -> length(x_{i-1})
    end
end
Base.getindex(x::PolyDist, i) = x.v[i]
@generated function Base.getindex(x_0::PolyDist{N}, I::Vararg{Int, N}) where N
    quote
        @nextract $N x i -> x_{i-1}[I[i]]
    end
end
@generated function Base.setindex!(x_0::PolyDist{N}, v, I::Vararg{Int, N}) where N
    Nm1 = N - 1
    x_last = Symbol("x_", Nm1)
    quote
        @nextract $Nm1 x i -> x_{i-1}[I[i]]
        $x_last[I[$N]] = v
    end

end
function Base.setindex!(x_0::PolyDist{N,T,V}, v::Vs, i::I) where {N,T,Vs,V <: AbstractArray{Vs},I <: Integer}
    x_0.v[i] = v
end

@generated function marginalize_last(v_0::PolyDist{N,T}, s::NTuple{ns,Val{I}} = (Val{32}(),), nextshape = (32,)) where {N,T,ns,I}
    Nm1 = N - ns
    Nm2 = Nm1 - 1
    v_last = Symbol("v_", Nm1)
    m_out = Symbol("m_", Nm1)
    quote
        $m_out = poly(Val{$Nm1}(), nextshape )
        @nloops $Nm1 i d -> begin
            1:length(v_{$Nm1-d})
        end d -> begin
            d > 1 && begin
                m_{d-1} = poly(Val{d-1}(), nextshape )
                push!(m_d, m_{d-1})
            end
            v_{$N-d} = v_{$Nm1-d}[i_{d}]
        end begin
            m_1[i_1] = marginalize($v_last, s )
        end
        $m_out
    end
end
@generated function marginalize_last2(v_0::PolyDist{N,T}, s::NTuple{ns,Val{I}} = (Val{32}(),), nextshape = (32,)) where {N,T,ns,I}
    Nm1 = N - ns
    v_last = Symbol("v_", Nm1)
    m_out = Symbol("m_", Nm1)
    quote
        local $m_out
        @nloops $Nm1 i d -> begin
            1:length(v_{$Nm1-d})
        end d -> begin
            m_d = d == $Nm1 ? poly(Val{$Nm1}(), nextshape ) : m_{1+d}[i_{d}]
            v_{$N-d} = v_{$Nm1-d}[i_{d}]
        end begin
            m_1[i_1] = marginalize($v_last, s )
        end
        $m_out
    end
end

@generated function receive_and_marginalize_last(v_0::PolyDist{N,T}, message, s = (Val{32}(),) ) where {N,T}
    Nm1 = N - length(s)
    v_last = Symbol("v_", Nm1)
    m_out = Symbol("m_", Nm1)
    plan = plan_rfft(message)
    inv_plan = inv(plan)
    mess_fft = plan * message
    v_fft = similar(mess_fft)
    quote
        @nloops $Nm1 i d -> begin
            1:length(v_{$Nm1-d})
        end d -> begin
            m_{d} = d == $Nm1 ? poly(Val{$Nm1}()) : m_{1+d}[i_{d}]
            v_{$N-d} = v_{$Nm1-d}[i_{d}]
        end begin
            A_mul_B!($v_fft, $plan, $v_last)
            $v_fft .*= $mess_fft
            A_mul_B!($v_last, $inv_plan, $v_fft)
            $m_1[i_1] = marginalize($v_last, $(s...) )
        end
        $m_out
    end
end

@generated function max_degree_poly(::Val{N}, l::Int, T = Int) where N
    quote
        j_0 = l
        s_0 = 0
        out = Vector{NTuple{N,T}}(0)
        @nloops $N i p -> begin
            0:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p
            j_{$N-p+1} = l - s_{$N-p+1}
        end begin
            push!(out, (@ntuple $N j -> T(i_{j})) )
        end
        out
    end
end
function tupGridd(x::NTuple{N,T}) where {N,T}
    ntuple( i -> begin
        x[i] > 0 ? :( $(Symbol('a' + i - 1, "_", x[i])) - $(Symbol('a' + i - 1, "_", x[i]-1)) ) : Symbol('a' + i - 1, "_", x[i])
     end, Val{N}())
end

@generated function lnorm_poly(::Val{N}, d::Int, l=0.5, T = Int) where N
    quote
        j_0 = d
        s_0 = 0
        reciprical = inv(l)
        dl = d^l
        out = Vector{NTuple{N,T}}(0)
        @nloops $N i p -> begin
            0:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p^l
            j_{$N-p+1} = floor(Int, (dl - s_{$N-p+1})^reciprical)
        end begin
            push!(out, (@ntuple $N j -> T(i_{j})) )
        end
        out
    end
end

@generated function polycount(::Val{N}, l::Int) where N
    seq = ( 1:1, 1:3, 1:9, 1:14 )
    step = ( 1:1, 2:3, 4:9, 10:14 )
    quote
        l == 1 && return 1
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = 0
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                p == $N ? $seq[i_p] : $step[i_p]
            end begin
                out += 1
            end
        end
        out
    end
end
@generated function polycount2(::Val{N}, seq = [ -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ], l = length(seq)-1) where N
    quote
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = 0
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                if p == $N
                    0:seq[i_p+1]
                else
                    1+seq[i_p]:seq[i_p+1]
                end
            end begin
                out += 1
            end
        end
        out
    end
end
@generated function polycount2(seq::NTuple{N,Vector{T}} = ([ -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ],[ -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ]), l = length(seq[1])-1) where {N,T <: Integer}
    quote
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = 0
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                if p == $N
                    0:seq[p][i_p+1]
                else
                    1+seq[p][i_p]:seq[p][i_p+1]
                end
            end begin
                out += 1
            end
        end
        out
    end
end

@generated function pcvec(::Val{N}, seq = [ -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ], l = length(seq)-2) where N
    quote
        l += 1
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = Vector{NTuple{$N,Cuint}}(0)
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                if p == $N
                    0:seq[i_p+1]
                else
                    1+seq[i_p]:seq[i_p+1]
                end
            end begin
                push!(out, (@ntuple $N k -> Cuint(k)) )
            end
        end
        out
    end
end

([-1,0,1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18],
[-1,0,1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18])
function poly_basis_gen(::Val{2})
    p = pcvec2(Val{2}(), (-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17), 192)
    p, coef(p)
end
function poly_basis_gen(::Val{3})
    p = pcvec2(Val{2}(), ( -1, 0, 1, 2, 3, 4, 5, 6,  8, 10,  14 ), 256)
    p, coef(p)
end
function poly_basis_gen(::Val{4})
    p = pcvec2(Val{3}(), ( -1, 0, 1, 3,4, 10, 15 ), 320)
    p, coef(p)
end
function poly_basis_gen(::Val{5})
    p = pcvec2((
        [ -1, 0, 1, 3,11, 14 ],
        [ -1, 0, 1, 3,10, 15 ],
        [ -1, 0, 1, 3,10, 15 ],
        [ -1, 0, 1, 3,10, 15 ],
        [ -1, 0, 1, 3,10, 15]), 384)
    p, coef(p)
end
function poly_basis_gen(::Val{6})
    p = pcvec2((
        [ -1, 0, 1, 3,5, 15 ],
        [ -1, 0, 1, 3,5, 15 ],
        [ -1, 0, 1, 3,5, 15 ],
        [ -1, 0, 1, 3,5, 14 ],
        [ -1, 0, 1, 3,6, 14 ],
        [ -1, 0, 1, 3,6, 14 ]), 448)
    p, coef(p)
end
function poly_basis_gen(::Val{7})
    p = pcvec2((
        [ -1, 0, 1, 3,14, 14 ],
        [ -1, 0, 1, 3,14, 15 ],
        [ -1, 0, 1, 3,14, 15 ],
        [ -1, 0, 1, 3,14, 15 ],
        [ -1, 0, 1, 3,14, 15 ],
        [ -1, 0, 1, 3,13, 15 ],
        [ -1, 0, 1, 3,13, 15 ]), 1024)
    p, coef(p)
end
function poly_basis_gen(::Val{8})
    p = pcvec2((
        [ -1, 0, 1, 3,9, 15 ],
        [ -1, 0, 1, 3,9, 15 ],
        [ -1, 0, 1, 3,8, 15 ],
        [ -1, 0, 1, 3,8, 14 ],
        [ -1, 0, 1, 3,9, 14 ],
        [ -1, 0, 1, 3,9, 14 ],
        [ -1, 0, 1, 3,9, 14 ],
        [ -1, 0, 1, 3,9, 14 ]), 1052)
    p, coef(p)
end
function poly_basis_gen(::Val{9})
    p = pcvec2((
        [ -1, 0, 1, 3,4, 15 ],
        [ -1, 0, 1, 3,4, 15 ],
        [ -1, 0, 1, 3,4, 15 ],
        [ -1, 0, 1, 3,4, 14 ],
        [ -1, 0, 1, 3,5, 14 ],
        [ -1, 0, 1, 3,5, 14 ],
        [ -1, 0, 1, 3,5, 14 ],
        [ -1, 0, 1, 3,5, 14 ],
        [ -1, 0, 1, 3,5, 14 ]), 1280)
    p, coef(p)
end

function poly_basis_gen(::Val{10})
    p = pcvec2((
        [ -1, 0, 1, 3, 12, 15 ],
        [ -1, 0, 1, 3, 12, 15 ],
        [ -1, 0, 1, 3, 12, 15 ],
        [ -1, 0, 1, 2, 12, 15 ],
        [ -1, 0, 1, 2, 13, 15 ],
        [ -1, 0, 1, 2, 13, 15 ],
        [ -1, 0, 1, 2, 13, 15 ],
        [ -1, 0, 1, 2, 13, 15 ],
        [ -1, 0, 1, 2, 13, 15 ],
        [ -1, 0, 1, 2, 13, 14 ]), 2112)
    p, coef(p)
end
function poly_basis_gen(::Val{11})
    p = pcvec2((
        [ -1, 0, 1, 2, 11, 14 ],
        [ -1, 0, 1, 2, 11, 14 ],
        [ -1, 0, 1, 2, 11, 15 ],
        [ -1, 0, 1, 2, 11, 15 ],
        [ -1, 0, 1, 2, 11, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 15 ]), 2304)
    p, coef(p)
end
function poly_basis_gen(::Val{12})
    p = pcvec2((
        [ -1, 0, 1, 2, 8, 14 ],
        [ -1, 0, 1, 2, 8, 14 ],
        [ -1, 0, 1, 2, 7, 15 ],
        [ -1, 0, 1, 2, 7, 15 ],
        [ -1, 0, 1, 2, 7, 15 ],
        [ -1, 0, 1, 2, 7, 15 ],
        [ -1, 0, 1, 2, 7, 15 ],
        [ -1, 0, 1, 2, 7, 15 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ]), 2496)
    p, coef(p)
end
function poly_basis_gen(::Val{13})
    p = pcvec2((
        [ -1, 0, 1, 2, 9, 15 ],
        [ -1, 0, 1, 2, 9, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 15 ],
        [ -1, 0, 1, 2, 10, 14 ],
        [ -1, 0, 1, 2, 10, 14 ],
        [ -1, 0, 1, 2, 10, 14 ],
        [ -1, 0, 1, 2, 10, 14 ],
        [ -1, 0, 1, 2, 10, 14 ],
        [ -1, 0, 1, 2, 10, 14 ],
        [ -1, 0, 1, 2, 10, 14 ]), 3584)
    p, coef(p)
end
function poly_basis_gen(::Val{14})
    p = pcvec2((
        [ -1, 0, 1, 2, 6, 15 ],
        [ -1, 0, 1, 2, 6, 15 ],
        [ -1, 0, 1, 2, 6, 15 ],
        [ -1, 0, 1, 2, 6, 14 ],
        [ -1, 0, 1, 2, 6, 14 ],
        [ -1, 0, 1, 2, 6, 14 ],
        [ -1, 0, 1, 2, 6, 14 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ],
        [ -1, 0, 1, 2, 7, 14 ]), 3840)
    p, coef(p)
end
function poly_basis_gen(::Val{15})
    p = pcvec2((
        [ -1, 0, 1, 2, 4, 14 ],
        [ -1, 0, 1, 2, 4, 14 ],
        [ -1, 0, 1, 2, 4, 14 ],
        [ -1, 0, 1, 2, 4, 14 ],
        [ -1, 0, 1, 2, 4, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ],
        [ -1, 0, 1, 2, 3, 14 ]), 4096)
    p, coef(p)
end
poly_basis_gen(::Val{N}) where N = throw("Size $N not yet implemented.")
@generated function pcvec2(::Val{N}, seq, total_size = (N+1)*64, ::Type{I} = Int) where {N,I<:Integer}
    quote
        l = length(seq) - 1
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = Vector{NTuple{$N,I}}(undef, 0)
        sizehint!(out, total_size)
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                if p == $N
                    0:seq[i_p+1]
                else
                    1+seq[i_p]:seq[i_p+1]
                end
            end begin
                push!(out, (@ntuple $N p -> I(k_p)) )
            end
        end
        out
    end
end
@generated function pcvec2(seq::NTuple{N,Vector{T}}, total_size = (N+1)*64) where {N, T<: Integer}
    quote
        l = length(seq[1]) - 1
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = Vector{NTuple{$N,T}}(undef, 0)
        sizehint!(out, total_size)
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                if p == $N
                    0:seq[p][i_p+1]
                else
                    1+seq[p][i_p]:seq[p][i_p+1]
                end
            end begin
                push!(out, (@ntuple $N k -> T(k)) )
            end
        end
        out
    end
end

@generated function pcvec3(::Val{N}) where N
    seq = ( (-1,0,1,2,3,4,5,6,7,8,9,10,11,12,12,13,13,13,13,13,13,14,14,15), )[N-1]
    l = length(seq) - 2
    quote
        l += 1
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = Vector{NTuple{$N,Int}}(undef, 0)
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                if p == $N
                    0:seq[i_p+1]
                else
                    1+seq[i_p]:seq[i_p+1]
                end
            end begin
                push!(out, (@ntuple $N k) )
            end
        end
        out
    end
end
@generated function smolcount(::Val{N}, l::Int, λ = 1) where N
    quote
        j_0 = l
        lλ = l^λ
        λi = 1/λ
        l_0 = 1
        s_0 = 0
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p^λ - 1
            j_{$N-p+1} = round(Int, (lλ - s_{$N-p+1})^λi)
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            println( ( @ntuple $N i ) )
        end
    end
end


text:
# This file is based on a part of Julia LinearAlgebra/src/special.jl. License is MIT: https://julialang.org/license

# Methods operating on different special matrix types

# Usually, reducedim_initarray calls similar, which yields a sparse matrix for a
# Diagonal/Bidiagonal/Tridiagonal/SymTridiagonal matrix. However, reducedim should
# yield a dense vector to increase performance.
Base.reducedim_initarray(A::Union{Bidiagonal,Tridiagonal,SymTridiagonal}, region, init, ::Type{R}) where {R} = fill(convert(R, init), Base.reduced_indices(A,region))

# Interconversion between special matrix types

# conversions from Diagonal to other special matrix types
Bidiagonal(A::Diagonal) = Bidiagonal(A.diag, fill!(similar(A.diag, length(A.diag)-1), 0), :U)
SymTridiagonal(A::Diagonal) = SymTridiagonal(A.diag, fill!(similar(A.diag, length(A.diag)-1), 0))
Tridiagonal(A::Diagonal) = Tridiagonal(fill!(similar(A.diag, length(A.diag)-1), 0), A.diag,
                                       fill!(similar(A.diag, length(A.diag)-1), 0))

# conversions from Bidiagonal to other special matrix types
Diagonal(A::Bidiagonal) = Diagonal(A.dv)
SymTridiagonal(A::Bidiagonal) =
    iszero(A.ev) ? SymTridiagonal(A.dv, A.ev) :
        throw(ArgumentError("matrix cannot be represented as SymTridiagonal"))
Tridiagonal(A::Bidiagonal) =
    Tridiagonal(A.uplo == 'U' ? fill!(similar(A.ev), 0) : A.ev, A.dv,
                A.uplo == 'U' ? A.ev : fill!(similar(A.ev), 0))

# conversions from SymTridiagonal to other special matrix types
Diagonal(A::SymTridiagonal) = Diagonal(A.dv)
Bidiagonal(A::SymTridiagonal) =
    iszero(A.ev) ? Bidiagonal(A.dv, A.ev, :U) :
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))
Tridiagonal(A::SymTridiagonal) =
    Tridiagonal(copy(A.ev), A.dv, A.ev)

# conversions from Tridiagonal to other special matrix types
Diagonal(A::Tridiagonal) = Diagonal(A.d)
Bidiagonal(A::Tridiagonal) =
    iszero(A.dl) ? Bidiagonal(A.d, A.du, :U) :
    iszero(A.du) ? Bidiagonal(A.d, A.dl, :L) :
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))

# conversions from AbstractTriangular to special matrix types
Bidiagonal(A::AbstractTriangular) =
    LinearAlgebra.isbanded(A, 0, 1) ? Bidiagonal(diag(A, 0), diag(A,  1), :U) : # is upper bidiagonal
    LinearAlgebra.isbanded(A, -1, 0) ? Bidiagonal(diag(A, 0), diag(A, -1), :L) : # is lower bidiagonal
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))

const ConvertibleSpecialMatrix = Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal,AbstractTriangular}

convert(T::Type{<:Diagonal},       m::Union{Bidiagonal,SymTridiagonal,Tridiagonal}) = m isa T ? m :
    isdiag(m) ? T(m) : throw(ArgumentError("matrix cannot be represented as Diagonal"))
convert(T::Type{<:SymTridiagonal}, m::ConvertibleSpecialMatrix) = m isa T ? m :
    issymmetric(m) && LinearAlgebra.isbanded(m, -1, 1) ? T(m) : throw(ArgumentError("matrix cannot be represented as SymTridiagonal"))
convert(T::Type{<:Tridiagonal},    m::ConvertibleSpecialMatrix) = m isa T ? m :
    LinearAlgebra.isbanded(m, -1, 1) ? T(m) : throw(ArgumentError("matrix cannot be represented as Tridiagonal"))

convert(T::Type{<:LowerTriangular}, m::Bidiagonal) = m isa T ? m :
    istril(m) ? T(m) : throw(ArgumentError("matrix cannot be represented as LowerTriangular"))
convert(T::Type{<:UpperTriangular}, m::Bidiagonal) = m isa T ? m :
    istriu(m) ? T(m) : throw(ArgumentError("matrix cannot be represented as UpperTriangular"))

# specialized +/- for structured matrices. If these are removed, it falls
# back to broadcasting which has ~2-10x speed regressions.
# For the other structure matrix pairs, broadcasting works well.

# For structured matrix types with different non-zero diagonals the underlying
# representations must be promoted to the same type.
# For example, in Diagonal + Bidiagonal only the main diagonal is touched so
# the off diagonal could be a different type after the operation resulting in
# an error. See issue #28994

function (+)(A::Bidiagonal, B::Diagonal)
    newdv = A.dv + B.diag
    Bidiagonal(newdv, A.ev, A.uplo)
end

function (-)(A::Bidiagonal, B::Diagonal)
    newdv = A.dv - B.diag
    Bidiagonal(newdv, A.ev, A.uplo)
end

function (+)(A::Diagonal, B::Bidiagonal)
    newdv = A.diag + B.dv
    Bidiagonal(newdv, B.ev, B.uplo)
end

function (-)(A::Diagonal, B::Bidiagonal)
    newdv = A.diag-B.dv
    Bidiagonal(newdv, -B.ev, B.uplo)
end

function (+)(A::Diagonal, B::SymTridiagonal)
    newdv = A.diag+B.dv
    SymTridiagonal(A.diag+B.dv, B.ev)
end

function (-)(A::Diagonal, B::SymTridiagonal)
    newdv = A.diag-B.dv
    SymTridiagonal(newdv, -B.ev)
end

function (+)(A::SymTridiagonal, B::Diagonal)
    newdv = A.dv+B.diag
    SymTridiagonal(newdv, A.ev)
end

function (-)(A::SymTridiagonal, B::Diagonal)
    newdv = A.dv-B.diag
    SymTridiagonal(newdv, A.ev)
end

# this set doesn't have the aforementioned problem

+(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl+B.ev, A.d+B.dv, A.du+B.ev)
-(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl-B.ev, A.d-B.dv, A.du-B.ev)
+(A::SymTridiagonal, B::Tridiagonal) = Tridiagonal(A.ev+B.dl, A.dv+B.d, A.ev+B.du)
-(A::SymTridiagonal, B::Tridiagonal) = Tridiagonal(A.ev-B.dl, A.dv-B.d, A.ev-B.du)

function (+)(A::Diagonal, B::Tridiagonal)
    newdv = A.diag+B.d
    Tridiagonal(B.dl, newdv, B.du)
end

function (-)(A::Diagonal, B::Tridiagonal)
    newdv = A.diag-B.d
    Tridiagonal(-B.dl, newdv, -B.du)
end

function (+)(A::Tridiagonal, B::Diagonal)
    newdv = A.d+B.diag
    Tridiagonal(A.dl, newdv, A.du)
end

function (-)(A::Tridiagonal, B::Diagonal)
    newdv = A.d-B.diag
    Tridiagonal(A.dl, newdv, A.du)
end

function (+)(A::Bidiagonal, B::Tridiagonal)
    newdv = A.dv+B.d
    Tridiagonal((A.uplo == 'U' ? (B.dl, newdv, A.ev+B.du) : (A.ev+B.dl, newdv, B.du))...)
end

function (-)(A::Bidiagonal, B::Tridiagonal)
    newdv = A.dv-B.d
    Tridiagonal((A.uplo == 'U' ? (-B.dl, newdv, A.ev-B.du) : (A.ev-B.dl, newdv, -B.du))...)
end

function (+)(A::Tridiagonal, B::Bidiagonal)
    newdv = A.d+B.dv
    Tridiagonal((B.uplo == 'U' ? (A.dl, newdv, A.du+B.ev) : (A.dl+B.ev, newdv, A.du))...)
end

function (-)(A::Tridiagonal, B::Bidiagonal)
    newdv = A.d-B.dv
    Tridiagonal((B.uplo == 'U' ? (A.dl, newdv, A.du-B.ev) : (A.dl-B.ev, newdv, A.du))...)
end

function (+)(A::Bidiagonal, B::SymTridiagonal)
    newdv = A.dv+B.dv
    Tridiagonal((A.uplo == 'U' ? (B.ev, A.dv+B.dv, A.ev+B.ev) : (A.ev+B.ev, A.dv+B.dv, B.ev))...)
end

function (-)(A::Bidiagonal, B::SymTridiagonal)
    newdv = A.dv-B.dv
    Tridiagonal((A.uplo == 'U' ? (-B.ev, newdv, A.ev-B.ev) : (A.ev-B.ev, newdv, -B.ev))...)
end

function (+)(A::SymTridiagonal, B::Bidiagonal)
    newdv = A.dv+B.dv
    Tridiagonal((B.uplo == 'U' ? (A.ev, newdv, A.ev+B.ev) : (A.ev+B.ev, newdv, A.ev))...)
end

function (-)(A::SymTridiagonal, B::Bidiagonal)
    newdv = A.dv-B.dv
    Tridiagonal((B.uplo == 'U' ? (A.ev, newdv, A.ev-B.ev) : (A.ev-B.ev, newdv, A.ev))...)
end

# fixing uniform scaling problems from #28994
# {<:Number} is required due to the test case from PR #27289 where eltype is a matrix.

function (+)(A::Tridiagonal{<:Number}, B::UniformScaling)
    newd = A.d .+ B.λ
    Tridiagonal(A.dl, newd, A.du)
end

function (+)(A::SymTridiagonal{<:Number}, B::UniformScaling)
    newdv = A.dv .+ B.λ
    SymTridiagonal(newdv, A.ev)
end

function (+)(A::Bidiagonal{<:Number}, B::UniformScaling)
    newdv = A.dv .+ B.λ
    Bidiagonal(newdv, A.ev, A.uplo)
end

(+)(A::UniformScaling, B::Tridiagonal{<:Number}) = Tridiagonal(B.dl, A.λ .+ B.d, B.du)
(+)(A::UniformScaling, B::SymTridiagonal{<:Number}) = SymTridiagonal(A.λ .+ B.dv, B.ev)

function (+)(A::UniformScaling, B::Bidiagonal{<:Number})
    newdv = A.λ .+ B.dv
    Bidiagonal(newdv, B.ev, B.uplo)
end

function (-)(A::UniformScaling, B::Tridiagonal{<:Number})
    newd = A.λ .- B.d
    Tridiagonal(-B.dl, newd, -B.du)
end

function (-)(A::UniformScaling, B::SymTridiagonal{<:Number})
    newdv = A.λ .- B.dv
    SymTridiagonal(newdv, -B.ev)
end

function (-)(A::UniformScaling, B::Bidiagonal{<:Number})
    newdv = A.λ .- B.dv
    Bidiagonal(newdv, -B.ev, B.uplo)
end

# fill[stored]! methods
fillstored!(A::Bidiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)
fillstored!(A::Tridiagonal, x) = (fill!(A.dl, x); fill!(A.d, x); fill!(A.du, x); A)
fillstored!(A::SymTridiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)

_small_enough(A::Bidiagonal) = size(A, 1) <= 1
_small_enough(A::Tridiagonal) = size(A, 1) <= 2
_small_enough(A::SymTridiagonal) = size(A, 1) <= 2

function fill!(A::Union{Bidiagonal,Tridiagonal,SymTridiagonal}, x)
    xT = convert(eltype(A), x)
    (iszero(xT) || _small_enough(A)) && return fillstored!(A, xT)
    throw(ArgumentError("array of type $(typeof(A)) and size $(size(A)) can
    not be filled with $x, since some of its entries are constrained."))
end

one(A::Bidiagonal{T}) where T = Bidiagonal(fill!(similar(A.dv, typeof(one(T))), one(T)), fill!(similar(A.ev, typeof(one(T))), zero(one(T))), A.uplo)
one(A::Tridiagonal{T}) where T = Tridiagonal(fill!(similar(A.du, typeof(one(T))), zero(one(T))), fill!(similar(A.d, typeof(one(T))), one(T)), fill!(similar(A.dl, typeof(one(T))), zero(one(T))))
one(A::SymTridiagonal{T}) where T = SymTridiagonal(fill!(similar(A.dv, typeof(one(T))), one(T)), fill!(similar(A.ev, typeof(one(T))), zero(one(T))))
# equals and approx equals methods for structured matrices
# SymTridiagonal == Tridiagonal is already defined in tridiag.jl

# SymTridiagonal and Bidiagonal have the same field names
==(A::Diagonal, B::Union{SymTridiagonal, Bidiagonal}) = iszero(B.ev) && A.diag == B.dv
==(B::Bidiagonal, A::Diagonal) = A == B

==(A::Diagonal, B::Tridiagonal) = iszero(B.dl) && iszero(B.du) && A.diag == B.d
==(B::Tridiagonal, A::Diagonal) = A == B

function ==(A::Bidiagonal, B::Tridiagonal)
    if A.uplo == 'U'
        return iszero(B.dl) && A.dv == B.d && A.ev == B.du
    else
        return iszero(B.du) && A.dv == B.d && A.ev == B.dl
    end
end
==(B::Tridiagonal, A::Bidiagonal) = A == B

==(A::Bidiagonal, B::SymTridiagonal) = iszero(B.ev) && iszero(A.ev) && A.dv == B.dv
==(B::SymTridiagonal, A::Bidiagonal) = A == B

import LinearAlgebra: TypeFuncs
LinearAlgebra.isstructurepreserving(::Union{typeof(abs),typeof(big)}, ::Union{Tridiagonal,SymTridiagonal,Bidiagonal}) = true
LinearAlgebra.isstructurepreserving(::TypeFuncs, ::Union{Tridiagonal,SymTridiagonal,Bidiagonal}) = true
LinearAlgebra.isstructurepreserving(::TypeFuncs, ::Ref{<:Type}, ::Union{Tridiagonal,SymTridiagonal,Bidiagonal}) = true

ext: jl


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   WDCFILES.txt
#
#   The WDC files contains the Kp and Ap indices.
#
#   For more information, see:
#
#       https://www.gfz-potsdam.de/en/kp-index/
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

################################################################################
#                       Private Structures and Variables
################################################################################

"""
Structure to store the interpolations of the data in WDC files.

# Fields

* `Kp`: Kp index.
* `Ap`: Ap index.

"""
struct _WDC_Structure
    Kp::_space_indices_itp_constant{SVector{8,Float64}}
    Ap::_space_indices_itp_constant{SVector{8,Float64}}
end

# Remote files: *.wdc
#
# This set will be configured in the function `init_space_indices()`.
_wdcfiles = RemoteFileSet(".wdc files", Dict{Symbol,RemoteFile}())

# Optional variable that will store the WDC data.
@OptionalData(_wdc_data, _WDC_Structure,
              "Run `init_space_indices()` with `:wdcfiles` in `enabled_files` array to initialize required data.")

################################################################################
#                               Public Functions
################################################################################

export get_Kp, get_Ap

#                                   Getters
# ==============================================================================

"""
    function get_Kp(JD::Number)

Return the Kp index at Julian Day `JD`.

"""
function get_Kp(JD::Number)
    Kp_day = get(SatelliteToolbox._wdc_data).Kp(JD)

    # Get the hour of the day and return the appropriate Kp.
    y, m, d, h, min, sec = JDtoDate(JD)

    return Kp_day[ floor(Int, h/3) + 1 ]
end

"""
    function get_Ap(JD::Number; mean::Tuple{Int} = (), daily = false)

Return the Ap index.

If `mean` is a tuple of two integers `(hi, hf)`, then the average between `hi`
and `hf` previous hours will be computed.

If `mean` is empty and `daily` is `true`, then the day average will be computed.

If `mean` keyword is empty, and `daily` keyword is `false`, then the Ap at
Julian day `JD` will be computed.

By default, `mean` is empty and `daily` is `false`.

"""
function get_Ap(JD::Number; mean::Tuple = (), daily = false)
    # Check if we must compute the mean of previous hours.
    if isempty(mean)
        Ap_day = get(SatelliteToolbox._wdc_data).Ap(JD)

        # Check if we must compute the daily mean.
        if daily
            return sum(Ap_day)/8
        else
            # Get the hour of the day and return the appropriate Ap.
            y, m, d, h, min, sec = JDtoDate(JD)

            return Ap_day[ floor(Int, h/3) + 1 ]
        end
    else
        # Check the inputs.
        (length(mean) != 2) && @error "The keyword `mean` must be empty or a tuple with exactly 2 integers."
        hi = mean[1]
        hf = mean[2]
        (hi > hf) && @error "The first argument of the keyword `mean` must be lower than the second."

        # Assemble the vector with the previous hours that will be averaged.
        hv = hi:3:hf

        # Compute the mean.
        Ap_sum = 0
        for h in hv
            Ap_sum += get_Ap(JD - h/24; mean = (), daily = false)
        end

        return Ap_sum/length(hv)
    end
end

################################################################################
#                              Private Functions
################################################################################

"""
    function _init_wdcfiles(;force_download = false, local_dir = nothing, wdcfiles_oldest_year = year(now())-3)

Initialize the data in the WDC files by creating `_wdcfiles_data`. The
initialization process is composed of:

1. Download the files, if it is necessary;
2. Parse the files;
3. Create the interpolations and the structures.

If the keyword `force_download` is `true`, then the files will always be
downloaded.

The user can also specify a location for the directory with the WDC files using
the keyword `local_dir`. If it is `nothing`, which is the default, then the
file will be downloaded.

The user can select what is the oldest year in which the data will be downloaded
by the keyword `wdcfiles_oldest_year`. By default, it will download the data
from 3 previous years.

The user can select what is the newest year in which the data will be downloaded
by the keyword `wdcfiles_newest_year`. It it is `nothing`, which is the default,
then it is set to the current year.

"""
function _init_wdcfiles(;force_download = false, local_dir = nothing,
                        wdcfiles_oldest_year = year(now())-3,
                        wdcfiles_newest_year = nothing)

    years     = Int[]
    filepaths = String[]

    if local_dir == nothing
        (wdcfiles_newest_year == nothing) && (wdcfiles_newest_year = year(now()))
        _prepare_wdc_remote_files(wdcfiles_oldest_year, wdcfiles_newest_year)
        download(_wdcfiles; force = force_download)

        # Get the files available and sort them by the year.
        for (sym,wdcfile) in _wdcfiles.files
            #
            # The year must not be obtained by the data inside the file,
            # because it contains only 2 digits and will break in 2032.
            # We will obtain the year by the symbol of the remote file. The
            # symbol name is:
            #
            #       kpYYYY
            #
            # where `YYYY` is the year.
            push!(years, parse(Int, String(sym)[3:6]))
            push!(filepaths, path(wdcfile))
        end
    else
        # If the user provided a directory, check what files are available.
        # Notice that the name must be the same as the ones online.
        for (root, dirs, files) in walkdir(local_dir)
            for file in files
                if occursin(r"^kp[1-2][0-9][0-9][0-9].wdc$", file)
                    year = parse(Int, file[3:6])

                    # Check if the year is not older than the oldest year.
                    if year >= wdcfiles_oldest_year
                        @info "Found WDC file `$file` related to the year `$year`."
                        push!(filepaths, joinpath(root, file))
                        push!(years,     year)
                    end
                end
            end
        end
    end

    p = sortperm(years)

    push!(_wdc_data,       _parse_wdcfiles(filepaths[p], years[p]))

    nothing
end

"""
    function _parse_wdcfiles(filepaths::Vector{String}, years::Vector{Int})

Parse the WDC files with paths in `filepaths` related to the years in `years`.

**Notice that the files must be sorted by the year!**

"""
function _parse_wdcfiles(filepaths::Vector{String}, years::Vector{Int})
    # Allocate the raw data.
    JD = Float64[]
    Kp = SVector{8,Float64}[]
    Ap = SVector{8,Float64}[]

    for (filepath, year) in zip(filepaths, years)

        open(filepath) do file
            # Read each line.
            for ln in eachline(file)
                # Get the Julian Day.
                month = parse(Int, ln[3:4])
                day   = parse(Int, ln[5:6])

                # The JD of the data will be computed at noon. Hence, we will be
                # able to use the nearest-neighbor algorithm in the
                # interpolations.
                JD_k  = DatetoJD(year, month, day, 12, 0, 0)

                # Get the vector of Kps and Aps.
                Ap_k = zeros(Float64,8)
                Kp_k = zeros(Float64,8)

                for i = 1:8
                    Kp_k[i] = parse(Int, ln[2(i-1) + 13:2(i-1) + 14])/10
                    Ap_k[i] = parse(Int, ln[3(i-1) + 32:3(i-1) + 34])
                end

                # Add data to the vector.
                push!(JD, JD_k)
                push!(Kp, SVector{8,Float64}(Kp_k))
                push!(Ap, SVector{8,Float64}(Ap_k))
            end
        end
    end

    # Create the interpolations for each parameter.
    knots    = (JD,)

    # Create the interpolations.
    itp_Kp = interpolate(knots, Kp, Gridded(Constant()))
    itp_Ap = interpolate(knots, Ap, Gridded(Constant()))

    _WDC_Structure(itp_Kp, itp_Ap)
end

"""
    function _prepare_wdc_remote_files(oldest_year::Number, newest_year::Number)

Configure all the WDC remote files between `newest_year` and `oldest_year`.
Notice that previous years will never be updated whereas the current year will
be updated daily.

If `oldest_year` is greater than current year, then only the files from the
current year will be downloaded.

If `newest_year` is smaller than `oldest_year`, then only the files from the
`oldest_year` will be downloaded.

This function modifies the global variable `_wdcfiles`.

"""
function _prepare_wdc_remote_files(oldest_year::Number, newest_year::Number)
    # Get the current year.
    current_year = year(now())

    # If `oldest_year` is greater than current year, then consider only the
    # current year.
    (oldest_year > current_year) && (oldest_year = current_year)
    (newest_year < oldest_year)  && (newest_year = oldest_year)
    (newest_year > current_year) && (newest_year = current_year)

    # For the current year, we must update the remote file every day. Otherwise,
    # we do not need to update at all.
    for y = oldest_year:newest_year
		filename = "kp$y"
        sym = Symbol(filename)
        file_y = @RemoteFile("ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/wdc/$filename.wdc",
                             file="$filename.wdc",
                             updates= (y == current_year) ? :daily : :never)

        merge!(_wdcfiles.files, Dict(sym => file_y))
    end

    nothing
end


text:
for f in (:seed_zeros, :seed_position, :seed_strength)
  @eval $f(v::Vector{<:Point},cfg::ComplexGradientConfig) =
      Point.($f(circulation,v,cfg)...)
end

for f in (:seed_zeros, :seed_position, :seed_strength)
  @eval $f(v::Vector{<:Blob},cfg::ComplexGradientConfig) =
      Blob.($f(circulation,v,cfg)...,Elements.blobradius(v))
end


text:
#
## Update active bots

function check_for_new_bots!(config_paths, users)
    config_files = glob(joinpath(Pinguu.base_dir[2:end], "users/*/*/config.json"))
    api = JSON.parsefile("usr/pinguu/api.json")
    telegram = JSON.parsefile(joinpath(Pinguu.base_dir[2:end], "telegram.json"))
    for file in config_files
        tmp = JSON.parsefile(file)
        if "false" == tmp["dry_mode"] && tmp["ws"] == 1 && tmp["status"] != "dead"
            !haskey(config_paths, tmp["symbol"]) ? config_paths[tmp["symbol"]] = [] : ""
            bot_path = joinpath("/", dirname(file))
            user_id = split(split(file, "pinguu_data/users/")[end], "/")[1]
            if haskey(api, user_id) && haskey(tmp, "strategy") && tmp["pureEMA"] == 0
                if !haskey(tmp, "symbols")
                    if !(bot_path in config_paths[tmp["symbol"]])
                        push!(config_paths[tmp["symbol"]], bot_path)
                    end
                else
                    for (k,v) in tmp["symbols"]
                        if !(bot_path in config_paths[k])
                            push!(config_paths[k], bot_path)
                        end
                    end
                end
            end
        end
    end
end
checkForNewBots!(config_paths, users) = @task check_for_new_bots!(config_paths, users)
export check_for_new_bots!

#
## Update active bots (function)

function check_streams(config_paths, binance_websocket_api_manager, data, client, users)
    #
    ## Put currently needed streams in dict
    tmp_markets = Dict()

    #
    ## Get current needed streams
    for (k,v) in config_paths
        for bot_path in v
            tmp = JSON.parsefile(joinpath(bot_path, "config.json"))
            if !haskey(tmp_markets, tmp["symbol"]) && tmp["status"] != "dead"
                tmp_markets[tmp["symbol"]] = []
            end
            if tmp["status"] != "dead"
                for condition in tmp["strategy"]["BUY"]["conditions"]
                    kline = tmp["strategy"]["BUY"]["indicators"][condition]["kline"]
                    if !haskey(tmp, "symbols")
                        if !(kline in tmp_markets[tmp["symbol"]])
                            push!(tmp_markets[tmp["symbol"]], kline)
                        end
                    else
                        for (k1,v1) in tmp["symbols"]
                            if !(kline in tmp_markets[k1])
                                push!(tmp_markets[k1], kline)
                            end
                        end
                    end
                end
            end
        end
    end
#     @info("Need to update streams?")

    invests = glob(joinpath(Pinguu.base_dir[2:end], "invest/users/*/filled/*/*/*"))
    for invest in invests
        symbol = split(basename(invest), "_")[1]
        if !haskey(tmp_markets, symbol)
            tmp_markets[symbol] = ["kline_1m"]
        else
            if !("kline_1m" in tmp_markets[symbol])
                push!(tmp_markets[symbol], "kline_1m")
            end
        end
    end
    #
    ## Check if a new stream has to start
    for (k,v) in tmp_markets
        for kline in v
            ## Check if stream is running
            create = true
            for (k1,v1) in binance_websocket_api_manager.stream_list
                if uppercase(v1["markets"][1]) == k && v1["channels"][1] == kline && v1["status"] in ["running", "starting", "restarting"]
                    create = false
                elseif uppercase(v1["markets"][1]) == k && v1["channels"][1] == kline
                    create = false
                    binance_websocket_api_manager.set_restart_request()
                end
            end
            if create
                binance_websocket_api_manager.create_stream([kline], [k], stream_label=k * " | " * split(kline, "_")[end], output="UnicornFy")
                !haskey(data, k) ? data[k] = Dict() : ""
                datapoints = 200
                data[k][kline] = Pinguu.get_klines(client, k, kline, datapoints)
            end
        end
    end

    #
    ## Check if a stream has to stop
    for (k,v) in binance_websocket_api_manager.stream_list
        if uppercase(v["markets"][1]) != "!USERDATA"
            kline = v["channels"][1]
            market = uppercase(v["markets"][1])
            stop_this_stream = false
            if haskey(tmp_markets, market)
                if !(kline in tmp_markets[market])
                    stop_this_stream = true
                end
            else
                stop_this_stream = true
            end
            if stop_this_stream
                binance_websocket_api_manager.stop_stream(k)
                sleep(10)
                binance_websocket_api_manager.delete_stream_from_stream_list(k)
            end
        #
        ## Check for user streams
        elseif uppercase(v["markets"][1]) == "!USERDATA"
            for (k,v) in users
                if v["api_secret"] == v["api_secret"] && v["api_key"] == v["api_key"]
                end
            end
        end
    end
#     @info("Check user streams")
    #
    ## Check for user streams
    for (k,v) in users
        create = true
        for (k1,v1) in binance_websocket_api_manager.stream_list
            if uppercase(v1["markets"][1]) == "!USERDATA"
                if v["api_secret"] == v1["api_secret"] && v["api_key"] == v1["api_key"]
                    create = false
                end
            end
        end
        if create
            binance_websocket_api_manager.create_stream("arr", "!userData", stream_label="User id: " * k,  api_key=v["api_key"], api_secret=v["api_secret"], output="UnicornFy")
        end
    end

end
checkStreams(config_paths, binance_websocket_api_manager, data, client, users) = @task check_streams(config_paths, binance_websocket_api_manager, data, client, users)

#
## Store data task

function store_data(data, symbol, kline_type)
    data_dict = Dict()
    data_dict["time"] = data[symbol][kline_type].close_time
    data_dict["price"] = data[symbol][kline_type].close
    open(joinpath(Pinguu.base_dir, "ticker/" * symbol * "_" * kline_type * "_graph.json"), "w") do f
        JSON.print(f, data_dict, 4)
    end
    if kline_type == "kline_1m"
        tmp = [data[symbol][kline_type].close_time, data[symbol][kline_type].open, data[symbol][kline_type].high, data[symbol][kline_type].low, data[symbol][kline_type].close]
        tmp_data = []
        for i in 1:length(tmp[1])
            push!(tmp_data, [1000 .* datetime2unix.(tmp[1][i]) .+ 1, tmp[2][i], tmp[3][i], tmp[4][i], tmp[5][i]])
        end
        open(joinpath(Pinguu.base_dir, "ticker/" * symbol * "_" * kline_type * "_candle.json"), "w") do f
            JSON.print(f, tmp_data, 4)
        end
    end
end

storedata(data, symbol, kline_type) = @task store_data(data, symbol, kline_type)

#
## Print Task (for notebook)

function print_task(data, kline_type)
    IJulia.clear_output(true)
    plot_symbol = []
    for (k,v) in data
        push!(plot_symbol, plot(v[kline_type].close_time[end-60:end], v[kline_type].close[end-60:end], label="", title=k))
    end
    p = plot(plot_symbol..., size=(1500,800))
    display(p)
end

printtask(data, kline_type) = @task print_task(data, kline_type)

#
## Analyse Task

function analyze_data(data, symbol)
    results = Dict()
    for (k,v) in data[symbol]
        results[k] = Dict()
        #
        ## Bollinger Bands
        results[k]["BB"] = Dict()
        for n in 20:10:50
            results[k]["BB"][string(n)] = Dict()
            for sigma in 1.0:1.0:5.0
                results[k]["BB"][string(n)][string(sigma)] = Indicators.bbands(data[symbol][k].close, n=n, sigma=sigma);
            end
        end

        #
        ## EMA
        results[k]["EMA"] = Dict()
        for n in [8, 13, 21, 55]
            results[k]["EMA"][string(n)] = Indicators.ema(data[symbol][k].close, n=n)
        end
    end
    return results
end
# export analyze_data

function check_for_signal(data, results, symbol, strategy)
    if strategy["BUY"]["type"] == "active"
        signals = []
        for con in strategy["BUY"]["conditions"]
            v = strategy["BUY"]["indicators"][con]
            prices = data[symbol][v["kline"]].close[end-1:end]
            if v["type"] == "BB"

                ## Above
                if v["above"]["active"] == 1
                    n = string(v["above"]["n"])
                    sigma = string(float(v["above"]["sigma"]))
                    bb_high = results[v["kline"]]["BB"][n][sigma][:,3][end-1:end]
                    if prices[2] >= bb_high[2] && prices[1] < bb_high[1]
                        push!(signals, con)
                    end
                end

                ## Below
                if v["below"]["active"] == 1 && !(con in signals)
                    n = string(v["below"]["n"])
                    sigma = string(float(v["below"]["sigma"]))
                    bb_low = results[v["kline"]]["BB"][n][sigma][:,1][end-1:end]
                    if prices[2] < bb_low[2]
                        push!(signals, con)
                    end
                end
            elseif v["type"] == "EMA"
                levels = sort(v["levels"])
                level = []
                for lvl in levels
                    if string(lvl) in keys(results[v["kline"]]["EMA"])
                        push!(level, results[v["kline"]]["EMA"][string(lvl)][end-1:end])
                    end
                end
                if length(levels) == 4
                    if v["trigger"] == "transition"
                        if level[1][2] > level[2][2] > level[3][2] > level[4][2] && level[1][2] > level[2][2] > level[3][2] <= level[4][2]
                            push!(signals, con)
                        end
                    elseif v["trigger"] == "trend"
                        if level[1][2] > level[2][2] > level[3][2] > level[4][2]
                            push!(signals, con)
                        end
                    end
                end
            end
        end

        ## Decision
        if strategy["BUY"]["union"] == 0
            buy_signal = length(signals) > 0 ? true : false
        else
            buy_signal = length(signals) == length(strategy["BUY"]["conditions"]) ? true : false
        end
        buy_note = ""
        for s in signals
            buy_note *= s
        end
    end
    sell_signal = false
    sell_note = ""
    return buy_signal, sell_signal, buy_note, sell_note
end
# export check_for_signal

function analyse_task(data, symbol, config_paths)
    check_paths = Dict()
    for (symbol, paths) in config_paths
        for bot_path in paths
            config = JSON.parsefile(joinpath(bot_path, "config.json"))
            open_orders = JSON.parsefile(joinpath(bot_path, "open_orders.json"))
            if length(open_orders["NEW"]) == 0 && length(open_orders["BUY"]) == 0 &&
                length(open_orders["SELL"]) == 0
                check_paths[symbol] = paths
            end
            if isfile(joinpath(bot_path, "manual_order.json"))
                check_paths[symbol] = paths
            end
        end
    end
    for (symbol, paths) in check_paths
        results = analyze_data(data, symbol)
        for bot_path in paths
            config = JSON.parsefile(joinpath(bot_path, "config.json"))
            if config["status"] == "active"
                open_orders = JSON.parsefile(joinpath(bot_path, "open_orders.json"))

                if length(open_orders["NEW"]) == 0 && length(open_orders["BUY"]) == 0 && length(open_orders["SELL"]) == 0 && !isfile(joinpath(bot_path, "manual_order.json"))
                    buy, sell, buy_note, sell_note = check_for_signal(data, results, symbol, config["strategy"])
                    if buy
                        open_orders["NEW"]["time"] = now()
                        open_orders["NEW"]["trigger"] = buy_note
                        open(joinpath(bot_path, "open_orders.json"), "w") do f
                            JSON.print(f, open_orders, 4)
                        end
                    end
                    if sell
                        #
                        ## EMA Bot

                    end
                elseif length(open_orders["NEW"]) == 0 && length(open_orders["BUY"]) == 0 && length(open_orders["SELL"]) == 0 && isfile(joinpath(bot_path, "manual_order.json"))
                    tmp = JSON.parsefile(joinpath(bot_path, "manual_order.json"))
                    open_orders["NEW"]["time"] = tmp["time"]
                    open_orders["NEW"]["trigger"] = tmp["trigger"]
                    open(joinpath(bot_path, "open_orders.json"), "w") do f
                        JSON.print(f, open_orders, 4)
                    end
                    rm(joinpath(bot_path, "manual_order.json"))
                elseif isfile(joinpath(bot_path, "manual_order.json"))
                    rm(joinpath(bot_path, "manual_order.json"))
                end
            end
        end
    end
end
# export analyse_task

analysetask(data, symbol, config_paths) = @task analyse_task(data, symbol, config_paths)


text:
function load_motor_data()
    df = CSV.read("data/motor.csv");
    x = df.times;
    y = df.accel;
    # Scaling the x-variables to lie in [0,1]
    L = maximum(x) - minimum(x);
    t = (x .- minimum(x))/L;
    return x, y, t
end

function alpha(p::Int)
    a = zeros(p, 1);
    for i = 0:p-1
        for j = 0:i
            for k = i:p-1
                a[i+1] = a[i+1] + (-1.0)^(k-j)/((k+j+1)*
                        factorial(j)*
                        factorial(p-1)*
                        factorial(p-1-k)*
                        factorial(i-j)*
                        factorial(k-i));
            end
        end
    end
    return a;
end

function spline_kernel(t::Array{Float64}, p::Int)

    #TO-DO
    # Check inputs
    # Maybe the Vector(range) can be done more eleganly.

    if all(diff(t) .> 0)
        monotonic = 1;
    elseif all(diff(t) .< 0)
        monotonic = 0;
    else
        throw(DomainError("t is not strictly monotonic"));
    end

    if size(t,2) != 1 && size(t,1) == 1;
        t = t';
    end

    fp = factorial.(p-1:-1:0)
    a = alpha(p).*fp;
    if monotonic == 1;
        U = (repeat(t,1,p).^Vector(p-1:-1:0)')./fp';
        V = (repeat(t,1,p).^Vector(p:2*p-1)').*a';
    elseif monotonic == 0
        U = (repeat(t,1,p).^Vector(p:2*p-1)')./fp';
        V = (repeat(t,1,p).^Vector(p-1:-1:0)').*a';
    end

    return U,V;
end

function spline_kernel_matrix(U::Array{Float64}, V::Array{Float64})
    return tril(U*V') + triu(V*U',1);
end

function logp(Σ, T, y, σf, σn)
    m = rank(T);
    n = length(y);
    Ki = cholesky(σf^2*Σ  + σn^2*I);
    if LinearAlgebra.issuccess(Ki) == false
        throw(DomainError("Cholesky Factorization failed"));
    end
    A = T'*(Ki\T);
    C = (Ki\T)*(A\(T'*inv(Ki)));
    lp = -0.5*y'*(Ki\y)  +
          0.5*y'*(C*y)   +
         -0.5*logdet(Ki) +
         -0.5*logdet(A)  +
         -(n-m)*0.5*log(2*π);
    return lp;
end

function dlogp(Σ, T, y, σf, σn)
    m = rank(T);
    n = length(y);
    Ki = cholesky(σf^2*Σ  + σn^2*I);
    if LinearAlgebra.issuccess(Ki) == false
        throw(DomainError("Cholesky Factorization failed"));
    end
    A = T'*(Ki\T);
    # TO_DO: Solving for the identity?
    dKf = 2*σf*Σ;
    dAf = -T'*(Ki\dKf)*(Ki\T);
    dKn = 2*σn*I;
    dAn = -T'*(inv(Ki)*dKn)*(Ki\T);
    Pf = Ki\dKf*inv(Ki);
    Pn = Ki\(dKn*inv(Ki));
    G = T*(A\(T'*(Ki\y)));
    dlpf =  0.5*y'*(Pf*y)                +
            0.5*(-2*(y'*Pf)*G + G'*Pf*G) +
           -0.5*tr(Ki\dKf)               +
           -0.5*tr(A\dAf);
   dlpn =   0.5*y'*(Pn*y)                   +
            0.5*(-2*(y'*Pn)*G + G'*Pn*G)    +
           -0.5*tr(inv(Ki)*dKn)             +
           -0.5*tr(A\dAn);
    return [dlpf; dlpn];
end

## Auxillary functions
function create_T(t, m)
    T = ones(length(t), m);
    for ν = 2:m
        T[:, ν] = t.^(ν-1)/factorial(ν - 1);
    end
    return T;
end

## Leave-one-out (LOO), what to do with explicit basis functions?


text:
using JuLIP, Test, IPFitting, DataFrames
using JuLIP.Potentials: evaluate_d
using JuLIP.MLIPs: IPSuperBasis
using ACE1
using IPFitting: Dat, LsqDB
Lsq = IPFitting.Lsq
Err = IPFitting.Errors
using LinearAlgebra: norm
using Random

# generate random data
function generate_data(species, L, rmax, N, calc)
   data = Dat[]
   for n = 1:N
      at = bulk(species; cubic=true, pbc=true) * L
      rattle!(at, rand() * rmax)
      E = energy(calc, at)
      F = forces(calc, at)
      push!(data, Dat(at, "rand"; E = E, F = F))
   end
   return data
end

##

Random.seed!(1)
r0 = rnn(:Si)
calc = StillingerWeber()
data = generate_data(:Si, 2, 0.33*r0, 100, calc)

# 2 stands for 2 neighbours i.e. body-order 3
b3basis(deg) = rpi_basis(species = :Si, N = 2, maxdeg = deg,
                         r0 = r0, rcut = cutoff(calc), rin = 0.5*r0,
                         pin = 2)

##
err_erms = Float64[]
err_frms = Float64[]
degrees = [4, 8, 12, 16, 20]
for deg in degrees
   B = b3basis(deg)
   @show length(B)
   db = LsqDB("", B, data)
   IP, fitinfo = Lsq.lsqfit( db,
                         E0 = 0.0,
             weights = Dict("default" => Dict("E" => 100.0, "F" => 1.0)),
                         error_table = true )
   push!(err_erms, fitinfo["errors"]["relrmse"]["set"]["E"])
   push!(err_frms, fitinfo["errors"]["relrmse"]["set"]["F"])
end

##
df = DataFrame( :degrees => degrees,
                :relrms_E => err_erms,
                :relrms_F => err_frms )
display(df)

(@test minimum(err_erms) < 2e-5) |> println
(@test minimum(err_frms) < 5e-3) |> println


text:
function subgradient!(dd_model::DDModel)
    #get an upper bound
    alpha = dd_model.solver_data.alpha
    lower_bound = dd_model.current_lower_bound

    if dd_model.upper_bound_function != nothing
        dd_model.upper_bound = dd_model.upper_bound_function(dd_model)
        step = alpha*abs(lower_bound - dd_model.upper_bound)/(norm([dd_model.residual_equality;dd_model.residual_inequality])^2)
    else
        upper_bound = nothing
        step = alpha
    end

    lambda_equality_delta = step*dd_model.residuals_equality  #update multipliers
    lambda_inequality_delta = step*dd_model.residuals_inequality
    lambda_linkvar_delta = step*dd_model.residuals_linkvars

    return lambda_equality_delta,lambda_inequality_delta,lambda_linkvar_delta
end

#TODO
function cutting_plane!(dd_model::DDModel)
end


text:
using Plots
using Flux
using Distributions
using LinearAlgebra
using Random
# using ForwardDiff
# using StaticArrays
# using CSV
# using DataFrames
# using Base.Iterators: repeated
# using Statistics

function leapfrog(q, p, step_size, num_steps; masses = [])
    d = length(q)
    if d != length(p)
        error("Position vector and momentum vector must have equal length.")
    end
    # if d == 1
    #     p = [p]
    #     q = [q]
    # end

    # define masses of kinetic energy
    if isempty(masses)
        masses = d > 1 ? ones(d) : 1
    else
        if length(masses)!= d
            error("Length of masses does not equal number of parameters")
        end
    end

    # new state according to leapfrom for hamilton dynamics
    p -= step_size .* ∇U(q)[1] / 2
    for jump in 1:1:num_steps-1
        q += step_size .* p ./ masses
        p -= step_size .* ∇U(q)[1]
    end
    q += step_size .* p ./ masses
    p -= step_size .* ∇U(q)[1] / 2 # the second half for the kinetic variables
    # println("p=",-p)
    # println("q=",q)
    return q, -p

end

## example of leapfrog on h(q,p) = q^2/2 + p^2/2
H(q, p) = q^2 / 2  + p^2 / 2
U(q) = q^2 / 2
gradU(q) = q                # analytical gradient
∇U(q) = gradient(U, q)      # numerical gradient

q = range(-1, 1, length=100)
p = range(-1, 1, length=100)

H.(q,p)
g = ∇U.(q)

contour(p,q, (p,q) -> H.(q,p), levels=1, color=:blue); xlims!((-2,2)); ylims!((-2,2))

function test_leapfrog(ϵ, L)
    gradU(q) = q
    ∇U(q) = q

    q = 0.0
    p = 1.0
    P = zeros(L)
    Q = zeros(L)
    for i in 1:1:L
        q, p = leapfrog.(q,p,ϵ,i)
        P[i]=p
        Q[i]=q
    end
    return P,Q
end

P,Q = test_leapfrog(0.5,20)
scatter!(P,Q,legend=nothing)

function hmc(current_q, step_size, num_steps; masses=[])
    q = current_q
    d = length(q)
    # if d != length(p)
    #     error("Position vector and momentum vector must have equal length.")
    # end
    # if d == 1
    #     p = [p]
    #     q = [q]
    # end

    # define masses of kinetic energy
    if isempty(masses)
        masses = d > 1 ? ones(d) : 1
    else
        if length(masses) != d
            error("Length of masses does not equal number of parameters")
        end
    end
    p = rand(MvNormal(zeros(d), Diagonal(masses)))
    current_p = p

    # new state according to leapfrom for hamilton dynamics
    q, p = leapfrog(q, p, step_size, num_steps; masses = [])

    # Evaluate potential and kinetic energies at start and end of trajectory
    current_U = U(current_q)
    current_K = sum(current_p.^2 ./ (2*masses))
    proposed_U = U(q)
    proposed_K = sum(p.^2 / (2*masses))

    if  rand(Uniform(0,1)) < exp(current_U - proposed_U + current_K - proposed_K)
        return q
    else
        return current_q
    end
end

μ = zeros(2)
Σ = [1 0.8;0.8 1]
# μ, Σ = zeros(d), Diagonal(ones(d))
U(q) = -logpdf(MvNormal(μ,Σ), q)
∇U(q) = gradient(q->logpdf(MvNormal(μ,Σ), q), q)

function main(n=1000)
    sampled_values = zeros(n,2)
    sampled_values[1,:] = [0, -0.0]
    for i in 2:n
        sampled_values[i,:] = hmc(sampled_values[i-1,:], 0.1, 10)
    end
    return sampled_values
end
Random.seed!(123)
samples = main()

histogram(samples[:,2])

x_range=range(-3, stop=3, length=100)
y_range=range(-3, stop=3, length=100)
Z = [pdf(MvNormal(μ, Σ), [x,y]) for x in x_range, y in y_range]
contour(x_range, y_range, Z)
plot!(samples[:,1], samples[:,2], color="red")

using GeometryBasics
heatmap(x_range, y_range, Z, color=:deep)
plot!(Point2.(eachcol(samples')), color="red", label="samples")

anim = @gif for i=1:10:1000
    scatter!(samples[i,1], samples legend=nothing, color="blue")
end every 5

max_stars_repo_name: cmerkatas/HMC

max_issues_repo_name: cmerkatas/HMC

max_forks_repo_name: cmerkatas/HMC


text:
module Kombinator
  # TODO: introduce a graph problem in the abstract hierarchy? Graph-related functions could be defined on it.
  # TODO: document the generic functions.
  # TODO: refactor the package further to remove the old st_prim()-like names? Only solve(::CombinatorialInstance, ::CombinatorialAlgorithm)::CombinatorialSolution should be used. At least, don't export all these names.
  # TODO: rename matching to "perfect matching" and introduce "imperfect matching".
  # TODO: add a note to rather use LightGraph's algorithms instead of these in this package for anything series.
  # TODO: add a link to other packages like LightGraph for the corresponding algorithms, so that they seamlessly integrate with the others?
  # TODO: how to add a link to functions like budgeted_msets_lp_select, budgeted_msets_lp_all? They can be quite useful (but do not bring much in terms of performance)

  using LinearAlgebra
  using LightGraphs
  using JuMP

  import Base: values
  import LightGraphs: src, dst
  import JuMP: value

  import Hungarian
  import Munkres

  abstract type CombinatorialInstance end
  abstract type CombinatorialSolution end
  abstract type CombinatorialAlgorithm end

  # Define the algorithm types here, at the global level: some names might be shared by several algorithms.
  # This is also the reason why there is no supplementary abstraction layer (e.g. a base type for all
  # bipartite-matching algorithms).
  # First, generic-named algorithms (dynamic, greedy, etc.); then, the ones with a recognised name (think Bellman-Ford).
  struct DynamicProgramming <: CombinatorialAlgorithm; end # When the DP algorithm has no more specific name exists, like Bellman-Ford
  struct GreedyAlgorithm <: CombinatorialAlgorithm; end
  struct LagrangianAlgorithm <: CombinatorialAlgorithm; end # Usually, for a crude approximation with no guarantee
  struct LagrangianRefinementAlgorithm <: CombinatorialAlgorithm; end # Usually, for an approximation term
  struct IteratedLagrangianRefinementAlgorithm <: CombinatorialAlgorithm; end # Usually, for a constant approximation ratio

  struct BellmanFordAlgorithm <: CombinatorialAlgorithm; end
  struct HungarianAlgorithm <: CombinatorialAlgorithm; end
  struct PrimAlgorithm <: CombinatorialAlgorithm; end

  # The case of linear programming is usually a bit more complex: in some cases, there may be several formulations.
  # In that case, don't use the generic object here, but rather a problem-specific object that can specify the
  # formulation to use.
  abstract type AbstractLinearProgramming <: CombinatorialAlgorithm; end
  struct LinearProgramming <: AbstractLinearProgramming; end

  """
      approximation_ratio(::CombinatorialInstance, ::CombinatorialAlgorithm)

  Returns the approximation ratio for this algorithm when run on the input instance. When the problem is minimising
  a function (e.g., finding the minimum-cost path), the ratio is defined as one constant ``r \\le 1.0``
  (ideally, the lowest) such that

  ``f(x^\\star) \\leq f(x) \\leq r \\cdot f(x^\\star),``

  where ``f(x^\\star)`` is the cost of the optimum solution and ``f(x)`` the one of the returned solution. On the
  contrary, for a maximisation problem, the definition is reversed (with ``r \\le 1.0``):

  ``r \\cdot f(x^\\star) \\geq f(x) \\geq f(x^\\star).``

  For an exact algorithm, the ratio is always ``1.0``, for both minimisation and maximisation problems.

  The returned ratio might be constant, if the algorithm provides a constant ratio; if the ratio is not constant (i.e.
  instance-dependent), it may either be a worst-case value or a truly instance-dependent ratio. Depending on the
  algorithm, this behaviour might be tuneable.

  If the algorithm has no guarantee, it should return `NaN`.
  """
  function approximation_ratio(::CombinatorialInstance, ::CombinatorialAlgorithm)
    return 1.0
  end

  """
      approximation_term(::CombinatorialInstance, ::CombinatorialAlgorithm)

  Returns the approximation term for this algorithm when run on the input instance. When the problem is minimising
  a function (e.g., finding the minimum-cost path), the term is defined as one constant ``t \\ge 0.0``
  (ideally, the lowest) such that

  ``f(x^\\star) \\leq f(x) \\leq f(x^\\star) + t,``

  where ``f(x^\\star)`` is the cost of the optimum solution and ``f(x)`` the one of the returned solution. On the
  contrary, for a maximisation problem, the definition is reversed (with ``t \\ge 0.0``):

  ``f(x^\\star) \\leq f(x) \\leq f(x^\\star) - t.``

  For an exact algorithm, the term is always ``0.0``, for both minimisation and maximisation problems.

  The returned term might be constant, if the algorithm provides a constant term; if the term is not constant (i.e.
  instance-dependent), it may either be a worst-case value or a truly instance-dependent term. Depending on the
  algorithm, this behaviour might be tuneable.

  If the algorithm has no guarantee, it should return `NaN`.
  """
  function approximation_term(::CombinatorialInstance, ::CombinatorialAlgorithm)
    return 0.0
  end

  # Include the actual contents of the package.
  include("helpers.jl")

  include("BipartiteMatching/matching.jl")
  include("BipartiteMatching/matching_hungarian.jl")
  include("BipartiteMatching/matching_dp.jl")
  include("BipartiteMatching/matching_budgeted.jl")
  include("BipartiteMatching/matching_budgeted_dp.jl")
  include("BipartiteMatching/matching_budgeted_lagrangian.jl")

  include("ElementaryPath/ep.jl")
  include("ElementaryPath/ep_dp.jl")
  include("ElementaryPath/ep_budgeted.jl")
  include("ElementaryPath/ep_budgeted_dp.jl")

  include("MSet/mset.jl")
  include("MSet/mset_greedy.jl")
  include("MSet/mset_dp.jl")
  include("MSet/mset_lp.jl")
  include("MSet/mset_budgeted.jl")
  include("MSet/mset_budgeted_dp.jl")
  include("MSet/mset_budgeted_lp.jl")

  include("SpanningTree/st.jl")
  include("SpanningTree/st_prim.jl")
  include("SpanningTree/st_budgeted.jl")
  include("SpanningTree/st_budgeted_lagrangian.jl")

  # Export all symbols. Code copied from JuMP.
  const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

  for sym in names(@__MODULE__, all=true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_")
      continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") && Base.isidentifier(sym_string[2:end])))
      continue
    end
    @eval export $sym
  end
end


text:
using PyPlot
using HDF5, JLD
using Statistics

using SubmodularMaximization

pygui(false)
name = "compare_adaptive_solvers"
fig_path = "./fig/$name"
data_path = "./data/$name"
mkpath(fig_path)
mkpath(data_path)

########
# Params
########
num_trials = 100

num_events = 50
num_agents = 50
num_sensors = 10
nominal_area = 0.6

max_success_probability = 1.0
sensor_radius = sqrt(nominal_area / (num_agents * pi))
station_radius = 4*sensor_radius

acceptable_suboptimality = 0.4
suboptimality_per_agent = acceptable_suboptimality / num_agents
communication_range = 2 * station_radius

agent_specification = ProbabilisticAgentSpecification(max_success_probability,
                                                      sensor_radius,
                                                      station_radius,
                                                      num_sensors)

###################
# Generate problems
###################

# each problem has a different set of agents and events
problems = map(1:num_trials) do unused
  agents = generate_agents(agent_specification, num_agents)

  events = generate_events(num_events)
  f(x) = mean_detection_probability(x, events)

  problem = ExplicitPartitionProblem(f, agents)
end
@save "$data_path/partition_matroids" map(x->x.partition_matroid, problems)

############################################
# Do intermediate processing of each problem
############################################

# now analyze weights
println("Generating weight matrices")
weight_matrices = map(problems, 1:length(problems)) do problem, ii
  println("Weight matrix $ii")
  compute_weight_matrix(problem)
end
@save "$data_path/weights" weight_matrices

println("Generating total weights")
total_weights = map(weight_matrices, 1:length(problems)) do weight_matrix, ii
  println("Total weights $ii")
  total_weight(weight_matrix)
end

println("Generating edge sets")
edge_sets = map(weight_matrices, 1:length(problems)) do weight_matrix, ii
  println("Triangle $ii")
  extract_triangle(weight_matrix)
end

println("Generating local partition sizes")
local_partition_sizes = map(problems) do problem
  compute_local_num_partitions(suboptimality_per_agent, problem)
end
@save "$data_path/local_partition_sizes" local_partition_sizes

println("Generating global partition sizes")
global_partition_sizes = map(problems) do problem
  compute_global_num_partitions(suboptimality_per_agent, problem)
end
@save "$data_path/global_partition_sizes" global_partition_sizes

###########################
# Generate adaptive solvers
###########################

dag_solver_array = []

println("Generating global solvers")
global_adaptive_solvers = map(1:length(problems)) do ii
  generate_by_global_partition_size(num_agents, global_partition_sizes[ii])
end
push!(dag_solver_array, global_adaptive_solvers)

println("Generating global range solvers")
global_range_solvers = map(1:length(problems)) do ii
  RangeSolver(problems[ii], global_adaptive_solvers[ii], communication_range)
end
push!(dag_solver_array, global_range_solvers)

println("Generating local solvers")
local_adaptive_solvers = map(1:length(problems)) do ii
  generate_by_local_partition_size(local_partition_sizes[ii])
end
push!(dag_solver_array, local_adaptive_solvers)

println("Generating local range solvers")
local_range_solvers = map(1:length(problems)) do ii
  RangeSolver(problems[ii], local_adaptive_solvers[ii], communication_range)
end
push!(dag_solver_array, local_range_solvers)

###################################
# Evaluate weights of deleted edges
###################################
function deleted_weights(solvers)
  println("...evaluating solver")
  mean(map(1:length(problems)) do ii
    deleted_edge_weight(solvers[ii], weight_matrices[ii])
  end)
end

println("Evaluating deleted weights")
mean_deleted = map(deleted_weights, dag_solver_array)

################
# Obtain results
################

function evaluate_adaptive_solver(solver_instances, title)
  values = map(1:length(problems)) do ii
    solve_dag(solver_instances[ii], problems[ii]).value
  end

  println("Evaluating $title")
  push!(titles, title)
  push!(results, values)
end

titles = String[]
results = Array[]

println("Evaluating Myopic")
push!(titles, "Myopic")
push!(results, map(x->solve_myopic(x).value, problems))

println("Evaluating Random")
push!(titles, "Random")
push!(results, map(x->solve_random(x).value, problems))

evaluate_adaptive_solver(global_adaptive_solvers, "Global Adaptive")
evaluate_adaptive_solver(global_range_solvers, "Global Range")
evaluate_adaptive_solver(local_adaptive_solvers, "Local Adaptive")
evaluate_adaptive_solver(local_range_solvers, "Local Range")

println("Evaluating Sequential")
push!(titles, "Sequential")
push!(results, map(x->solve_sequential(x).value, problems))

results = hcat(results...)
@save "$data_path/results" results

################
# Generate plots
################

function compute_range_frequency(x::Array{Int64, 1})
  range = collect(minimum(x):maximum(x))

  frequencies = zeros(Int64, length(range))

  for value in x
    frequencies[value - range[1] + 1] += 1
  end

  (range, frequencies)
end

function frequency_bar(x::Array{Int64, 1})
  r, f = compute_range_frequency(x)
  bar(r, f)
end

figure()
boxplot(results, notch=false, vert=false)
yticks(1:length(titles), titles)
xlabel("Expect. num. ident. events")
tight_layout()

save_fig(fig_path, "results")

# plot histograms of edge weights
figure()
PyPlot.plt.hist(total_weights, 20)
tt = "Total Graph Weight Frequency"

save_fig(fig_path, tt)
title(tt)

figure()
PyPlot.plt.hist(vcat(edge_sets...), 20)
tt = "Edge Weight Frequency"
save_fig(fig_path, tt)
title(tt)

# plot histograms of partition sizes
figure()

frequency_bar(global_partition_sizes)
tt = "Global Partition Size Frequency"
ylabel("Freq.")
xlabel("Partition size, \$n_d\$")
tight_layout()

save_fig(fig_path, tt)
title(tt)

figure()
frequency_bar(vcat(local_partition_sizes...))
tt = "Local Partition Size Frequency"
ylabel("Freq.")
xlabel("Local partition size, \$k_i\$")
tight_layout()

save_fig(fig_path, tt)
title(tt)

# plot deleted weights
bar_order = [1, 2, 3, 4]
figure()
barh(1:4, mean_deleted[bar_order])
yticks(1:4, titles[2:5][bar_order])
tt = "Deleted Edge Weights"
xlabel("Cumulative deleted edge weight")
tight_layout()

save_fig(fig_path, tt)
title(tt)

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
# Stricker Klaus 2019 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - https://github.com/strickek/ReadWriteDlm2.jl

"""
## ReadWriteDlm2
`ReadWriteDlm2` functions `readdlm2()`, `writedlm2()`, `readcsv2()` and
`writecsv2()` are similar to those of stdlib.DelimitedFiles, but with additional
support for `Date`, `DateTime`, `Time`, `Complex`, `Rational`, `Missing` types
and special decimal marks.

### `readcsv2(), writecsv2()`:
For "decimal dot" users the functions `readcsv2()` and `writecsv2()` have the
respective defaults: Delimiter is `','` (fixed) and `decimal='.'`.

### `readdlm2(), writedlm2()`:
The basic idea of these functions is to support the "decimal comma countries".
They use `';'` as default delimiter and `','` as default decimal mark.
"Decimal dot" users of these functions need to define `decimal='.'`

### Detailed Documentation:
For more information about functionality and (keyword) arguments see `?help` for
`readdlm2()`, `writedlm2()`, `readcsv2()` and `writecsv2()`.
"""
module ReadWriteDlm2

using Dates
using DelimitedFiles
using DelimitedFiles: readdlm_string, val_opts

export readdlm2, writedlm2, readcsv2, writecsv2

# readdlm2() and readcsv2()
# =========================

"""

    dfregex(df::AbstractString, locale::AbstractString=\"english\")

Create a regex string `r\"^...\$\"` for the given `Date` or `DateTime`
`format`string `df`.

The regex groups are named according to the `format`string codes. `locale` is
used to calculate min and max length of month and day names (for codes: UuEe).
"""
function dfregex(df::AbstractString, locale::AbstractString="english")
    # calculate min and max string lengths of months and day_of_weeks names
    Ule = try extrema([length(Dates.monthname(i;locale=locale)) for i in 1:12])catch; (3, 9) end
    ule = try extrema([length(Dates.monthabbr(i;locale=locale)) for i in 1:12])catch; (3, 3) end
    Ele = try extrema([length(Dates.dayname(i;locale=locale)) for i in 1:7])catch; (6, 9) end
    ele = try extrema([length(Dates.dayabbr(i;locale=locale)) for i in 1:7])catch; (3, 3) end

    codechars = 'y', 'Y', 'm', 'u', 'e', 'U', 'E', 'd', 'H', 'M', 'S', 's', 'Z', 'z', '\\'
    r = "^ *"; repeat_count = 1; ldf = length(df); dotsec = false
    for i = 1:ldf
        repeat_next = ((i < ldf) && (df[(i + 1)] == df[i])) ? true : false
        ((df[i] == '.') && (i < ldf) && (df[(i + 1)] == 's')) && (dotsec = true)
        repeat_count = (((i > 2) && (df[(i - 2)] != '\\') && (df[(i - 1)] == df[i])) ||
                        ((i == 2) && (df[1] == df[2]))) ? (repeat_count + 1) : 1
        r = r * (
        ((i > 1) && (df[(i - 1)] == '\\')) ? string(df[i]) :
        ((df[i] == 'y') && (repeat_count < 5) && !repeat_next) ? "(?<y>\\d{1,4})" :
        ((df[i] == 'y') && (repeat_count > 4) && !repeat_next) ? "(?<y>\\d{1,$repeat_count})" :
        ((df[i] == 'Y') && (repeat_count < 5) && !repeat_next) ? "(?<y>\\d{1,4})" :   # new
        ((df[i] == 'Y') && (repeat_count > 4) && !repeat_next) ? "(?<y>\\d{1,$repeat_count})" :
        ((df[i] == 'm') && (repeat_count == 1) && !repeat_next) ? "(?<m>0?[1-9]|1[012])" :
        ((df[i] == 'm') && (repeat_count == 2) && !repeat_next) ? "(?<m>0[1-9]|1[012])" :
        ((df[i] == 'm') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<m>0[1-9]|1[012])" :
        ((df[i] == 'u') && (repeat_count == 1)) ? "(?<u>[A-Za-z\u00C0-\u017F]{$(ule[1]),$(ule[2])})" :
        ((df[i] == 'U') && (repeat_count == 1)) ? "(?<U>[A-Za-z\u00C0-\u017F]{$(Ule[1]),$(Ule[2])})" :
        ((df[i] == 'e') && (repeat_count == 1)) ? "(?<e>[A-Za-z\u00C0-\u017F]{$(ele[1]),$(ele[2])})" :
        ((df[i] == 'E') && (repeat_count == 1)) ? "(?<E>[A-Za-z\u00C0-\u017F]{$(Ele[1]),$(Ele[2])})" :
        ((df[i] == 'd') && (repeat_count == 1) && !repeat_next) ? "(?<d>0?[1-9]|[12]\\d|3[01])" :
        ((df[i] == 'd') && (repeat_count == 2) && !repeat_next) ? "(?<d>0[1-9]|[12]\\d|3[01])" :
        ((df[i] == 'd') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<d>0[1-9]|[12]\\d|3[01])" :
        ((df[i] == 'H') && (repeat_count == 1) && !repeat_next) ? "(?<H>0?\\d|1\\d|2[0-3])" :
        ((df[i] == 'H') && (repeat_count == 2) && !repeat_next) ? "(?<H>0\\d|1\\d|2[0-3])" :
        ((df[i] == 'H') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<H>0\\d|1\\d|2[0-3])" :
        ((df[i] == 'M') && (repeat_count == 1) && !repeat_next) ? "(?<M>\\d|[0-5]\\d)" :
        ((df[i] == 'M') && (repeat_count == 2) && !repeat_next) ? "(?<M>[0-5]\\d)" :
        ((df[i] == 'M') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<M>[0-5]\\d)" :
        ((df[i] == 'S') && (repeat_count == 1) && !repeat_next) ? "(?<S>\\d|[0-5]\\d)" :
        ((df[i] == 'S') && (repeat_count == 2) && !repeat_next) ? "(?<S>[0-5]\\d)" :
        ((df[i] == 'S') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<S>[0-5]\\d)" :
        ((df[i] == '.') && dotsec) ? "" :
        ((df[i] == '.')) ? "\\." :
        ((df[i] == 's') && (dotsec == true) && (repeat_count < 4) && !repeat_next) ? "(\\.(?<s>\\d{0,3}0{0,6}))?" :
        ((df[i] == 's') && (dotsec == true) && (repeat_count > 3) && !repeat_next) ? "(\\.(?<s>\\d{$(repeat_count)}))?" :
        ((df[i] == 's') && (dotsec == false) && (repeat_count < 4) && !repeat_next) ? "(?<s>\\d{3})?" :
        ((df[i] == 's') && (dotsec == false) && (repeat_count > 3) && !repeat_next) ? "(?<s>\\d{$(repeat_count)})?" :
        ((df[i] == 'z') && !repeat_next) ? "(?<z>[\\+|\\-]?(0\\d|1\\d|2[0-3]):?[0-5]\\d)" :
        ((df[i] == 'Z') && !repeat_next) ? "(?<Z>[A-Z]{3,14})" :
        in(df[i], codechars) ? "" : string(df[i])
        )
    end
    return Regex(r * " *" * string('$'))
end

"""

    parseothers(y::AbstractString, doparsetime::Bool, doparsecomplex::Bool, doparserational::Bool)

Parse string `y` for `Time`, `Complex` and `Rational` format and if match return the value.
Otherwise return the input string `y`.
"""
function parseothers(y, doparsetime, doparsecomplex, doparserational)

    if doparsetime # parse Time
        mt = match(r"^ *(0?\d|1\d|2[0-3])[:Hh]([0-5]?\d)(:([0-5]?\d)([\.,](\d{1,3})(\d{1,3})?(\d{1,3})?)?)? *$", y)
        if mt != nothing
            h = parse(Int, mt[1])
            mi = parse(Int, mt[2])
            (mt[4] == nothing) ? s = ms = us = ns = 0 :
            s  = parse(Int, lpad(string(mt[4]), 2, string(0))); (mt[6] == nothing) ? ms = us = ns = 0 :
            ms = parse(Int, rpad(string(mt[6]), 3, string(0))); (mt[7] == nothing) ? us = ns = 0 :
            us = parse(Int, rpad(string(mt[7]), 3, string(0))); (mt[8] == nothing) ? ns = 0 :
            ns = parse(Int, rpad(string(mt[8]), 3, string(0)))
            return Dates.Time(h, mi, s, ms, us, ns)
        end
    end

    if doparsecomplex # parse Complex
        mc = match(r"^ *(-?\d+(\.\d+)?([eE]-?\d+)?|(-?\d+)//(\d+)) ?([\+-]) ?(\d+(\.\d+)?([eE]-?\d+)?|(\d+)//(\d+))(\*im|\*i|\*j|im|i|j) *$", y)
        if mc != nothing
            real =
                ((mc[4] != nothing) && (mc[5] != nothing)) ? //(parse(Int, mc[4]), parse(Int, mc[5])) :
                ((mc[2] == nothing) && (mc[3] == nothing)) ? parse(Int, mc[1]) : parse(Float64, mc[1])
            imag =
                ((mc[10] != nothing) && (mc[11] != nothing)) ? //(parse(Int, mc[6]*mc[10]), parse(Int, mc[11])) :
                ((mc[8] == nothing) && (mc[9] == nothing)) ? parse(Int, mc[6]*mc[7]) : parse(Float64, mc[6]*mc[7])
            return complex(real, imag)
        end
    end

    if doparserational # parse Rational
        mr = match(r"^ *(-?\d+)//(-?\d+) *$", y)
        if mr != nothing
            nu = parse(Int, mr[1])
            de = parse(Int, mr[2])
            return //(nu, de)
        end
    end

    return y
end

"""

    readcsv2(source, T::Type=Any; opts...)

Equivalent to `readdlm2()` with delimiter `','` and `decimal='.'`. Default Type
`Any` activates parsing of `Bool`, `Int`, `Float64`, `Complex`, `Rational`,
`Missing`, `DateTime`, `Date` and `Time`.

# Code Example
```jldoctest
julia> using ReadWriteDlm2

julia> B = Any[1 complex(1.5,2.7);1.0 1//3];

julia> writecsv2("test.csv", B)

julia> readcsv2("test.csv")
2×2 Array{Any,2}:
 1    1.5+2.7im
 1.0    1//3
```
"""
readcsv2(input; opts...) =
    readdlm2auto(input, ',', Any, '\n', false; decimal='.', opts...)

readcsv2(input, T::Type; opts...) =
    readdlm2auto(input, ',', T, '\n', false; decimal='.', opts...)

"""

    readdlm2(source; options...)
    readdlm2(source, T::Type; options...)
    readdlm2(source, delim::AbstractChar; options...)
    readdlm2(source, delim::AbstractChar, T::Type; options...)
    readdlm2(source, delim::AbstractChar, eol::AbstractChar; options...)
    readdlm2(source, delim::AbstractChar, T::Type, eol::AbstractChar; options...)

Read a matrix from `source`. The `source` can be a text file, stream or byte
array. Each line (separated by `eol`, this is `'\\n'` by default) gives one row.
The columns are separated by `';'`, another `delim` can be defined.

Pre-processing of `source` with regex substitution changes the decimal marks
from `d,d` to `d.d`. For default `rs` the keyword argument `decimal=','` sets
the decimal Char in the `r`-string of `rs`. When a special regex substitution
tuple `rs=(r.., s..)` is defined, the argument `decimal` is not used.
Pre-processing can be switched off with: `rs=()`.

In addition to stdlib readdlm(), data is also parsed for `Dates` formats (ISO),
the `Time` format `\"HH:MM[:SS[.s{1,9}]]\"` and for complex and rational numbers.
To deactivate parsing dates/time set: `dfs=\"\", dtfs=\"\"`.
`locale` defines the language of day (`E`, `e`) and month (`U`, `u`) names.

The result will be a (heterogeneous) array of default type `Any`.
Other (abstract) types for the array elements could be defined.
If data is empty, a `0×0 Array{T,2}` is returned.
If `dfheader=true` instead of `header=true` is used, the returned tuple
(data::Array{T,2}, header::Array{Symbol,1}) fits for `DataFrames`.

# Additional Keyword Arguments

* `decimal=','`: Decimal mark Char used by default `rs`, irrelevant if `rs`-tuple is not the default one
* `rs=(r\"(\\d),(\\d)\", s\"\\1.\\2\")`: Regex (r,s)-tuple, the default change d,d to d.d if `decimal=','`
* `dtfs=\"yyyy-mm-ddTHH:MM:SS.s\"`: Format string for DateTime parsing
* `dfs=\"yyyy-mm-dd\"`: Format string for Date parsing
* `locale=\"english\"`: Language for parsing dates names
* `dfheader=false`: Return header in DataFrames format if `true`
* `missingstring=\"na\"`: How missing values are represented

Find more information about `readdlm()` functionality and (keyword) arguments -
 which are also supported by `readdlm2()` - in `help` for `readdlm()`.

# Code Example
```jldoctest
julia> using ReadWriteDlm2

julia> A = Any[1 1.2; "text" missing];

julia> writedlm2("test.csv", A)

julia> readdlm2("test.csv")
2×2 Array{Any,2}:
 1        1.2
  "text"   missing
```
"""
readdlm2(input; opts...) =
    readdlm2auto(input, ';', Any, '\n', true; opts...)

readdlm2(input, T::Type; opts...) =
    readdlm2auto(input, ';', T, '\n', false; opts...)

readdlm2(input, dlm::AbstractChar; opts...) =
    readdlm2auto(input, dlm, Any, '\n', true; opts...)

readdlm2(input, dlm::AbstractChar, T::Type; opts...) =
    readdlm2auto(input, dlm, T, '\n', false; opts...)

readdlm2(input, dlm::AbstractChar, eol::AbstractChar; opts...) =
    readdlm2auto(input, dlm, Any, eol, true; opts...)

readdlm2(input, dlm::AbstractChar, T::Type, eol::AbstractChar; opts...) =
    readdlm2auto(input, dlm, T, eol, false; opts...)

function readdlm2auto(input, dlm, T, eol, auto;
        decimal::AbstractChar=',',
        rs::Tuple=(r"(\d),(\d)", s"\1.\2"),
        dtfs::AbstractString="yyyy-mm-ddTHH:MM:SS.s",
        dfs::AbstractString="yyyy-mm-dd",
        locale::AbstractString="english",
        dfheader::Bool=false,
        missingstring::AbstractString="na",
        opts...)

    if dfheader == true
        opts = [opts...]
        push!(opts, :header => true)
    end

    ((!isempty(dtfs) && !occursin(Regex("[^YymdHMSs]"), dtfs)) ||
    (!isempty(dfs) && !occursin(Regex("[^YymdHMSs]"), dfs))) && info(
        """
        Format string for DateTime(`$dtfs`) or Date(`$dfs`)
        contains numeric code elements only. At least one non-numeric
        code element or character is needed for parsing dates.
        """)

    # "parsing-logic" - defaults
    doparsedatetime = !isempty(dtfs)
    doparsedate = !isempty(dfs)
    doparsetime = false
    doparsecomplex = true
    doparserational = true
    doparsemissing = false
    doparsenothing = false
    convertarray = true
    T2 = Any

    # "parsing-logic" for different T::Types
    ((typeintersect(DateTime, T)) == Union{}) && (doparsedatetime = false)
    ((typeintersect(Date, T)) == Union{}) && (doparsedate = false)
    doparsetime = ((doparsedatetime && doparsedate) || ((Time <: T) && !(Any <:T)))
    T <: Union{AbstractFloat, AbstractString, Char} && (T2 = T; convertarray = false)
    ((typeintersect(Complex, T)) == Union{}) && (doparsecomplex = false)
    ((typeintersect(Rational, T)) == Union{}) && (doparserational = false)
    Missing <: T && (doparsemissing = true)
    Nothing <: T && (doparsenothing = true)
    (Any <: T) && (convertarray = false)

    s = read(input, String)

    # empty input data - return empty array
    if (isempty(s) || (s == string(eol)))
        return Array{T2}(undef, 0, 0)
    end

    if (!isempty(rs) && (decimal != '.')) # do pre-processing of decimal mark

        # adopt dfs if decimal between two digits of formatstring
        if (rs == (r"(\d),(\d)", s"\1.\2")) && doparsedate
            drs = (Regex("([YymdHM])$decimal([YymdHM])"), s"\1.\2")
            dfs = replace(dfs, drs[1] => drs[2])
        end

        # adopt dtfs if decimal between two digits of formatstring
        if (rs == (r"(\d),(\d)", s"\1.\2")) && doparsedatetime
            drs = (Regex("([YymdHMSs])$decimal([YymdHMSs])"), s"\1.\2")
            dtfs = replace(dtfs, drs[1] => drs[2])
        end

        # Change default regex substitution Tuple if decimal != ','
        if ((rs == (r"(\d),(\d)", s"\1.\2")) && (decimal != ','))
            rs = (Regex("(\\d)$decimal(\\d)"), s"\1.\2")
        end

        # Error if decimal mark to replace is also the delim Char
        "1"*string(dlm)*"1" != replace("1"*string(dlm)*"1", rs[1] => rs[2]) &&
        error(
            """
            Error: Decimal mark to replace is also the delim Char.
            Pre-processing with decimal mark Regex substitution for
            `$(dlm)` (= delim!!) is not allowed - change rs/decimal or delim!
            """)

        # Regex substitution decimal
        s = replace(s, rs[1] => rs[2])

    end

    # Using stdlib DelimitedFiles internal function to read dlm-string
    z = readdlm_string(s, dlm, T2, eol, auto, val_opts(opts))

    if isa(z, Tuple)
         y, h = z
         if dfheader == true # format header to DataFrames format
             h = Symbol.(reshape(h, :))
             z = (y, h)
         end
    else
         y = z
    end

    # Formats, Regex for Date/DateTime parsing
    doparsedatetime && (dtdf = DateFormat(dtfs, locale); rdt = dfregex(dtfs, locale))
    doparsedate && (ddf = DateFormat(dfs, locale); rd = dfregex(dfs, locale))

    for i in eachindex(y)
        if isa(y[i], AbstractString)
            if doparsedatetime && occursin(rdt, y[i]) # parse DateTime
                try y[i] = DateTime(y[i], dtdf) catch; end
            elseif doparsedate && occursin(rd, y[i]) # parse Date
                try y[i] = Date(y[i], ddf) catch; end
            elseif doparsemissing && y[i] == missingstring # parse Missing
                try y[i] = missing catch; end
            elseif doparsenothing && y[i] == "nothing" # parse Nothing
                try y[i] = nothing catch; end
            else # parse Time, Complex and Rational
                try y[i] = parseothers(y[i], doparsetime, doparsecomplex, doparserational) catch; end
            end
        end
    end

    if convertarray
        isa(z, Tuple) ? z = (convert(Array{T}, y), h) : z = convert(Array{T}, z)
    end

    return z

end # End function readdlm2auto()

# writedlm2() and writecsv2()
# ===========================

"""

    floatformat(a, decimal::AbstractChar)

Convert Int or Float64 numbers to String and change decimal mark.
"""
function floatformat(a, decimal)
    a = string(a)
    (decimal != '.') && (a = replace(a, '.' => decimal))
    return a
end

"""

    timeformat(a, decimal::AbstractChar)

Convert Time to String, optional with change of decimal mark for secounds.
"""
function timeformat(a, decimal)
    a = string(a)
    (decimal != '.') && (a = replace(a, '.' => decimal))
    return a
end

"""

    Complexformat(a, decimal::AbstractChar, imsuffix::AbstractString)

Convert Complex number to String, optional change of decimal and/or imsuffix.
"""
function complexformat(a, decimal, imsuffix)
    a = string(a)
    (imsuffix != "im") && (a = string(split(a, "im")[1], imsuffix))
    (decimal != '.') && (a = replace(a, '.' => decimal))
    return a
end

"""

    writecsv2(f, A; opts...)

Equivalent to `writedlm2()` with fixed delimiter `','` and `decimal='.'`.

# Code Example
```jldoctest
julia> using ReadWriteDlm2

julia> B = Any[1 complex(1.5,2.7);1.0 1//3];

julia> writecsv2("test.csv", B)

julia> read("test.csv", String)
"1,1.5 + 2.7im\\n1.0,1//3\\n"
```
"""
writecsv2(f, a; opts...) =
    writedlm2auto(f, a, ','; decimal='.', opts...)

"""

    writedlm2(f, A; opts...)
    writedlm2(f, A, delim; opts...)

Write `A` (a vector, matrix, or an iterable collection of iterable rows) as
text to `f`(either a filename string or an IO stream). The columns will be
separated by `';'`, another `delim` (Char or String) can be defined.

By default, a pre-processing of values takes place. Before writing as strings,
decimal marks are changed from `'.'` to `','`. With the keyword argument
`decimal=` another decimal mark can be defined.
To switch off this pre-processing set: `decimal='.'`.

In `writedlm2()` the output format for `Date` and `DateTime` data can be
defined with format strings. Defaults are the ISO formats. Day (`E`, `e`)
and month (`U`, `u`) names are written in the `locale` language. For writing
`Complex` numbers the imaginary component suffix can be selected with the
`imsuffix=` keyword argument.

# Additional Keyword Arguments

* `decimal=','`: Character for writing decimal marks
* `dtfs=\"yyyy-mm-ddTHH:MM:SS.s\"`: DateTime write format
* `dfs=\"yyyy-mm-dd\"`: Date write format
* `locale=\"english\"`: Language for DateTime writing
* `imsuffix=\"im\"`: Complex Imag suffix `\"im\"`, `\"i\"` or `\"j\"`
* `missingstring=\"na\"`: How missing values are written

# Code Example
```jldoctest
julia> using ReadWriteDlm2, Dates

julia> A = Any[1 1.2; "text" Date(2017)];

julia> writedlm2("test.csv", A)

julia> read("test.csv", String)
"1;1,2\\ntext;2017-01-01\\n"
```
"""
writedlm2(io::IO, a; opts...) =
    writedlm2auto(io, a, ';'; opts...)

writedlm2(io::IO, a, dlm; opts...) =
    writedlm2auto(io, a, dlm; opts...)

writedlm2(f::AbstractString, a; opts...) =
    writedlm2auto(f, a, ';'; opts...)

writedlm2(f::AbstractString, a, dlm; opts...) =
    writedlm2auto(f, a, dlm; opts...)

function writedlm2auto(f, a, dlm;
        decimal::AbstractChar=',',
        dtfs::AbstractString="yyyy-mm-ddTHH:MM:SS.s",
        dfs::AbstractString="yyyy-mm-dd",
        locale::AbstractString="english",
        imsuffix::AbstractString="im",
        missingstring::AbstractString="na",
        opts...)

    ((!isempty(dtfs) && !occursin(Regex("[^YymdHMSs]"), dtfs)) ||
    (!isempty(dfs) && !occursin(Regex("[^YymdHMSs]"), dfs))) && info(
        """
        Format string for DateTime(`$dtfs`) or Date(`$dfs`)
        contains numeric code elements only. At least one non-numeric
        code element or character is needed for parsing dates.
        """)

    (string(dlm) == string(decimal)) && error(
        "Error: decimal = delim = ´$(dlm)´ - change decimal or delim!")

    ((imsuffix != "im") && (imsuffix != "i") && (imsuffix != "j")) && error(
    "Only `\"im\"`, `\"i\"` or `\"j\"` are valid arguments for keyword `imsuffix=`.")

    if isa(a, Union{Nothing, Missing, Number, TimeType})
         a = [a]  # create 1 element Array
    end

    if isa(a, AbstractArray)
        fdt = !isempty(dtfs)  # Bool: format DateTime
        dtdf = DateFormat(dtfs, locale)
        fd = !isempty(dfs)    # Bool: format Date
        ddf = DateFormat(dfs, locale)
        ft = (decimal != '.')   # Bool: format Time (change decimal)
        fc = ((imsuffix != "im") || (decimal != '.')) # Bool: format Complex

        # create b for manipulation/write - keep a unchanged
        b = similar(a, Any)
        for i in eachindex(a)
            b[i] =
            isa(a[i], AbstractFloat) ? floatformat(a[i], decimal) :
            isa(a[i], Missing) ? missingstring :
            isa(a[i], Nothing) ? "nothing" :
            isa(a[i], DateTime) && fdt ? Dates.format(a[i], dtdf) :
            isa(a[i], Date) && fd ? Dates.format(a[i], ddf) :
            isa(a[i], Time) && ft ? timeformat(a[i], decimal) :
            isa(a[i], Complex) && fc ? complexformat(a[i], decimal, imsuffix) :
            string(a[i])
        end
        else  # a is not a Number, TimeType or Array -> no preprocessing
        b = a
    end

    writedlm(f, b, dlm; opts...)

    end # End function writedlm2auto()

end # End module ReadWriteDlm2