Light Cone constructor for SpatioTemporalEmbedding about timeseriesprediction.jl HOT 5 CLOSED

function indices_within(radius, dimension)
    #Return indices β within hypersphere with radius
    #Floor radius to nearest integer. Does not lose points
    r = floor(Int,radius)
    #Hypercube of indices
    hypercube = CartesianIndices((repeat([-r:r], dimension)...,))
    #Select subset of hc which is in Hypersphere

    βs = [ β for β ∈ hypercube if norm(β.I) <= radius ]
end


function TimeseriesPrediction.SpatioTemporalEmbedding(
    s::AbstractArray{<:AbstractArray{T,Φ}}, ΔT, τ, speed, boundary::BC
    ) where {T,Φ, BC<:AbstractBoundaryCondition}
    if !(BC <: PeriodicBoundary) && typeof(boundary.c) != T
        throw(ArgumentError(
        "Boundary value must be same element type as the timeseries data."))
    end

    τs = Int[]
    βs = CartesianIndex{Φ}[]
    for delay = ΔT:-τ:1 # Backwards for forward embedding 1:τ:ΔT
        radius = max(speed*delay, 1.) # No radius lower than 1
        β = indices_within(radius, Φ)
        push!(βs, β...)
        push!(τs, repeat([ΔT - delay], length(β))...)
    end
    X = length(τs)
    return SpatioTemporalEmbedding{X}(τs, βs, boundary, size(s[1]))
end



function visualize_STE(em::STE{2})
    figure()
    grid()
    axis("equal")
    r = maximum([max(β.I...) for β in em.β])
    xticks(-r:r)
    yticks(-r:r)
    for τ in em.τ
        points = [ β.I for (delay,β) in zip(em.τ,em.β) if delay == τ]
        x = map(xy -> xy[1], points)
        y = map(xy -> xy[2], points)
        scatter(x,y, c="C$τ", s=100)

        #Plot Circle
        radius = maximum(norm.(points))
        x = radius*cos.(0:0.1:2π)
        y = radius*sin.(0:0.1:2π)
        plot(x,y, c="C$τ")
    end
    scatter([0],[0], color="black", s=200)
end


s = [rand(10,10) for i =1:5];
em = SpatioTemporalEmbedding(s, 5,1,1,PeriodicBoundary())

visualize_STE(em)

x[t] = (s[t], s[t+τ], ... , s[t+D*τ])

x[t] = (s[t], s[t+τ[1]], ..., s[t+τ[D]])