The code that I used in the example is slightly different from the code I
showed, and also wrong. I added the condition that ir != jr when deciding
whether to record a value in the triplet form. However, that is the wrong
condition. I should check whether ioffset + ir is equal to joffset + jr.
function edist(a::AbstractMatrix, nblocks::Int, threshold)
N = size(a,2)
si = subinds(nblocks,N)
i = Int[]
j = Int[]
x = eltype(a)[]
for ii in si, jj in si
r = pairwise(Euclidean(), view(a,:,ii), view(a,:,jj))
ioffset = start(ii) - 1
joffset = start(jj) - 1
for ir in 1:size(r,1), jr in 1:size(r,2)
if (ie = ioffset + ir) != (je = joffset + jr) && (rij =
r[ir,jr]) < threshold
push!(i,ie)
push!(j,je)
push!(x,rij)
end
end
end
sparse(i,j,x)
end
It turns out that this does not affect the results but that was just a
matter of good fortune.