For various reasons, I want to create a static sparse vector (SSV) type. In contrast to what is done (for matrices) in LuxurySparse.jl, I want the sparsity structure to be part of the data type.
I do that by using a recursive datatype SSV parametrized by an index I that has a value for this index, and a next containing the rest of the vector.
struct EmptySSV{L} end
struct SSV{I,V,N} # Static Sparse Vector, recursive, in which the sparsity pattern is part of the type
value::V
next ::N
end
SSV{I}(value,next) where{I} = SSV{I,typeof(value),typeof(next)}(value,next)
function StaticVector(l,i,v) # general purpose, type unstable constructor
out = EmptySSV{l}()
for j = reverse(eachindex((i)))
out = SSV{i[j]}(v[j],out)
end
return out
end
for OP∈(:(+),:(-))
@eval function Base.$OP(v1::SSV{I1,V1,N1},v2::SSV{I2,V2,N2}) where{I1,V1,N1,I2,V2,N2}
if I1==I2 SSV{I1}($OP(v1.value,v2.value),$OP(v1.next,v2.next))
elseif I1< I2 SSV{I1}( v1.value ,$OP(v1.next,v2 ))
elseif I1> I2 SSV{I2}( v2.value ,$OP(v1 ,v2.next))
end
end
@eval Base.$OP(v::SSV{I,V,N} ,e::EmptySSV{L }) where{I,V,N,L} = SSV{I}(v.value,$OP(v.next,e ))
@eval Base.$OP(e::EmptySSV{L },v::SSV{I,V,N }) where{I,V,N,L} = SSV{I}(v.value,$OP(e ,v.next))
@eval Base.$OP( ::EmptySSV{L }, ::EmptySSV{L }) where{L } = EmptySSV{L}()
@eval Base.$OP( ::EmptySSV{L1}, ::EmptySSV{L2}) where{L1,L2 } = error("Attempt to add vectors of ≠ lengths")
end
# Tests
v = StaticVector(10,[3,5],[3.,5.])
w = StaticVector(10,[2,5],[2.,5.])
@show v[3]
@show v[5]
#@show v[1] # throws "Not stored"
#@show v[11] # throws "Out of bounds"
@show v+w===StaticVector(10,[2,3,5],[2.,3.,10.])
Everything goes fine until I try to make it possible to have vectorised operations on these vectors. I tried to overload broadcast:
function Base.Broadcast.broadcast(op,v1::SSV{I1,V1,N1},v2::SSV{I2,V2,N2}) where{I1,V1,N1,I2,V2,N2}
if I1==I2 SSV{I1}(broadcast(op,v1.value,v2.value),broadcast(op,v1.next,v2.next))
elseif I1< I2 SSV{I1}( v1.value ,broadcast(op,v1.next,v2 ))
elseif I1> I2 SSV{I2}( v2.value ,broadcast(op,v1 ,v2.next))
end
end
Base.Broadcast.broadcast(op,v::SSV{I,V,N} ,e::EmptySSV{L }) where{I,V,N,L} = SSV{I}(v.value,broadcast(op,v.next,e ))
Base.Broadcast.broadcast(op,e::EmptySSV{L },v::SSV{I,V,N }) where{I,V,N,L} = SSV{I}(v.value,broadcast(op,e ,v.next))
Base.Broadcast.broadcast(op, ::EmptySSV{L }, ::EmptySSV{L }) where{L } = EmptySSV{L}()
Base.Broadcast.broadcast(op, ::EmptySSV{L1}, ::EmptySSV{L2}) where{L1,L2 } = error("Attempt to add vectors of ≠ lengths")
but my methods do not seem to be called, I get an error that no method matching iterate(::SSV{3, Float64, SSV{5, Float64, EmptySSV{10}}}) which indeed I did not define, but did not call from my broadcast. So I was not supposed to overload broadcast maybe?
How are broadcast and iterate connected? I do not believe Julia would be able, with me implementing ìterate alone, to do the sort of recursive braiding I do in my + and - methods, and so I am concerned that I do not have all the bits in place…
Looking for inspiration, I find that SparseArray.jl has functions like _bcast_binary_map. … but called by what?
Anyone knows how things look under the hood?
) my intention to have the whole interface (though I need