# Dot operation broadcast for a "StaticSparseVector"
**URL:** https://discourse.julialang.org/t/dot-operation-broadcast-for-a-staticsparsevector/66124
**Category:** Performance
**Tags:** broadcast, sparse, staticarrays
**Created:** [August 10, 2021, 11:27am UTC](https://discourse.julialang.org/t/dot-operation-broadcast-for-a-staticsparsevector/66124 "2021-08-10T11:27:29Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)
#### Post date: [August 10, 2021, 11:27am UTC](https://discourse.julialang.org/t/dot-operation-broadcast-for-a-staticsparsevector/66124/1 "2021-08-10T11:27:29Z")
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.
```julia
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`:
```julia
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](https://github.com/JuliaLang/julia/blob/master/stdlib/SparseArrays/src/sparsevector.jl) has functions like `_bcast_binary_map`. … but called by what?
Anyone knows how things look under the hood?
---
### Author:  [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)
#### Post date: [August 12, 2021, 4:54am UTC](https://discourse.julialang.org/t/dot-operation-broadcast-for-a-staticsparsevector/66124/2 "2021-08-12T04:54:04Z")
Maybe I can rephrase my question better: My SSV type has the peculiar property that `getindex` is not implemented for it. This is because this vector (by design) can combine elements of different types, and hence deciding what zero to return on indexing an element not explicitely stored is not possible. Nonsense? Well, this workhorse type is to be wrapped into more specific types…
**What do I need to implement to allow broadcasting on my own vector type?**
**Is it a requirement for broadcasting that indexing must be possible?**
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [August 12, 2021, 6:42am UTC](https://discourse.julialang.org/t/dot-operation-broadcast-for-a-staticsparsevector/66124/3 "2021-08-12T06:42:08Z")
There’s a list of functions required for broadcasting: [Interfaces · The Julia Language](https://docs.julialang.org/en/v1/manual/interfaces/#man-interfaces-broadcasting)
But if you can’t know the type of the zeros, how will broadcasting work? What happens with `isa.(Integer, x)`?
---
### Author:  [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)
#### Post date: [August 12, 2021, 6:50am UTC](https://discourse.julialang.org/t/dot-operation-broadcast-for-a-staticsparsevector/66124/4 "2021-08-12T06:50:29Z")
Hi Gustaphe,
I will study the link you gave me in fine details! Give me a few days…
Very good question indeed! Your question guides me: this should return a `StaticSparseVector` of booleans. A central issue: I do not want broadcasting to result in full vectors (unless… see below)
What if I add a `Vector` to a `StaticSparseVector`. I would return the same `Vector`, except with an additions _where relevant_, and so I never have to find out “what kind of zero”. Indeed, for efficiency, I do not want to spend time adding zeros.
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [August 12, 2021, 7:09am UTC](https://discourse.julialang.org/t/dot-operation-broadcast-for-a-staticsparsevector/66124/5 "2021-08-12T07:09:11Z")
The problem with that approach is, a sparse vector if booleans is just a representation of a vector where every element is `false` except the specified. But that doesn’t really make sense for you. How about
```julia
f(x) = x isa Integer || iszero(x)
all(f.(mySparseVector))
```
Should this be `true` or `false` for a sparse vector with only `Integer`s in it? Because intuitively I’d say `true`, but with your suggestion it’s either `false` or undefined.
---
### Author:  [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)
#### Post date: [August 12, 2021, 7:35am UTC](https://discourse.julialang.org/t/dot-operation-broadcast-for-a-staticsparsevector/66124/6 "2021-08-12T07:35:23Z")
Touché, excellent!
OK, my StaticSparseVector is a weird customer, and it is hereby renamed HeterogenousStaticSparseVector (HSSV). You point out that sparse and heterogenous put together make it impossible to implement the whole sparse vector interface, including, quite possibly, broadcasting. I’ll chew on this!
And indeed, this is not (or no longer!!! 😁) my intention to have the whole interface (though I need `a.*b`, and this operation has a well defined result…): HSSV will be wrapped into various other containers with a well defined `eltype`, to which we will be able to give a reasonable “AbstractVector” interface.
So why heterogenous? Well, you can do some cool stuff with it. One favorite example is implementing a matrix as sparse vector of sparse vector. How could I go static (stay off the heap) otherwise? A matrix is a row of columns, and each column has a different type (since sparsity pattern is part of the vector type). Of course for a matrix, you can (should!) implement this is other ways, but this flexibility will be valuable for forward automatic differentiation to higher orders with sparse partials, which are themselves forwarddiff objects with sparse partials…
Thank you for this discussion, this helps me a lot!