Cartesian Product iteration

[GlenHertz]

This is a follow-up to the discussion on the julia-users mailing list here.

There seemed to be good progress but I didn't see any code committed.

[JeffBezanson]

In fact, I've been thinking of adding this to the language as follows:

for I in X...
    A[I...] = B[I...]
end

where X is an n-tuple of iterables. This expands to a nest of n for loops and delivers a tuple of elements in I. Other syntax/behavior ideas welcome.

[StefanKarpinski]

I like both the idea and the syntax. People are inevitably going to ask for support for X being other kinds of collections. Is that feasible or can we not infer their size well enough?

[JeffBezanson]

It is feasible, since after all we can't infer the size of every tuple anyway. You just want to avoid spurious use of non-tuples, like {1:n,1:n} instead of (1:n,1:n). This can be implemented by calling the equivalent of gen_cartesian_map.

[timholy]

I like it a lot. Maybe an interesting use case would be to consider what it would look like to compute the discrete Laplacian of a multidimensional array this way. Presumably you'd have to do iteration over the components of I?

[johnmyleswhite]

+1

[timholy]

Would you be able to do math on the index tuples? For example, suppose I wanted to implement imfilter for arrays of arbitrary dimensions: I might want something like this: (edit: rats, stray keypress had me submit early, continuing to edit)

function imfilter(X,K)
    center = ceil([size(K)...]/2)
    Y = similar(X)
    for I in size(X)...    # Note: have to handle boundaries somehow!
        s = zero(eltype(X))
        for J in size(K)...
            s += K[J...]*X[I...+J...-center...]
        end
        Y[I...] = s
    end
end

[StefanKarpinski]

It would be something like ranges(X)... not size(I)... since that would iterate a single time with I equal to size(I). Also: man, that's a lot of dots.

[JeffBezanson]

We could define things like interior(X, 2) to get various useful index ranges for use with ....
A syntax like X[(I.+J.-center)...] might work.

[timholy]

Stefan, you're right of course.

Another thing to think about is that there may be powerful performance arguments to be able to reduce most array accesses to linear indexing. I seem to remember that 4d indexing was something like 20x slower than linear, when the accesses were in a pattern that allowed you to precalculate the outer dimensions. The imfilter example would probably be best implemented in terms of a list of coefficients and stride offsets.

But for boundary conditions you also have to keep track of the actual coordinates, and feed that into your linear indexing, so you do need cartesian iteration in parallel with linear indexing. If you imagine implementing reflecting, periodic, and NA boundary conditions (the latter meaning "just pretend any undefined elements don't exist, normalizing by the ones you have"), it seems like the sub2ind-style transformations could get fairly hairy. Certainly, this was a pain to handle in Grid.jl (which implements even more boundary conditions), although it eventually worked out.

[GunnarFarneback]

For the record there's an imfilter implementation for arbitrary dimensions at
GunnarFarneback@2c3d599

You definitely want to do the inner loop over the kernel as a linear loop and in many cases handle boundary conditions with pre-padding.

[timholy]

A bit off-topic, but @GunnarFarneback, your imfilter is certainly better than the one in Images.jl. I'd welcome a pull request, if you're willing to share.

While linear indexing when you can do it will always be faster, I suspect my 20x number needs revision in light of e550406. Nice work, Jeff.

[GunnarFarneback]

Sharing is no problem, deciding how to fit it in is more of one. But maybe it's time to get that implementation replaced.

Going back to topic, the bilateral filter may be a good concrete example of how to stress a Cartesian product language feature. In 2D a somewhat naive implementation looks like

function bilateral(in::Matrix{Float64}, kernel_size::Int,
                   sigma_s::Float64, sigma_i::Float64)
    k_s = -0.5 / sigma_s^2
    k_i = -0.5 / sigma_i^2
    kernel_radius = fld(kernel_size - 1, 2)
    out = similar(in)
    height = size(in, 1)
    width = size(in, 2)
    for x = 1:width
        for y = 1:height
            center_value = in[y, x]
            sum_w = 0.0
            sum = 0.0

            for dx = -kernel_radius:kernel_radius
                xi = x + dx
                if xi < 1 || xi > width
                    continue
                end
                for dy = -kernel_radius:kernel_radius
                    yi = y + dy
                    if yi < 1 || yi > height
                        continue
                    end

                    value = in[yi, xi]
                    diff = value - center_value
                    wi = exp(k_s * (dx^2 + dy^2) + k_i * diff^2)
                    sum_w += wi
                    sum += wi * value
                end
            end

            out[y, x] = sum / sum_w
        end
    end
    return out
end

It would be lovely if this could be written compactly for arbitrary dimensionality.

[timholy]

Just wanted to give another example here, showing how useful/necessary this will be for switching over to returning views.

It would be nice to write a sum function something like this:

function (+)(A::AbstractSubView, B::AbstractSubView)
    if size(A) != size(B)
        error("Dimension mismatch")
    end
    C = Array(promote(eltype(A), eltype(B)), size(A))
    for I in C...
        C[I...] = A[indexes(A,I.)...] + B[indexes(B,I.)...]
    end
    C
end

I suspect we won't be happy with the performance until this type of syntax generates code something like this:

# Implementation for 3d
function (+)(A::SubView, B::SubView)
    if size(A) != size(B)
        error("Dimension mismatch")
    end
    C = Array(promote(eltype(A), eltype(B)), size(A))
    s1 = size(C, 1)
    s2 = size(C, 2)
    strdC2 = s1
    strdC3 = s1*s2
    strdA2 = size(parent(A), 1)
    strdA3 = size(parent(A), 2)*strdA2
    strdB2 = size(parent(B), 1)
    strdB3 = size(parent(B), 2)*strdB2
    for i3 = 1:size(C, 3)
        o3 = (i3-1)*strdC2
        oA3 = (A.indexes[3][i3]-1)*strdA3
        oB3 = (B.indexes[3][i3]-1)*strdB3
        for i2 = 1:s2
            o2 = o3 + (i2-1)*strdC2
            oA2 = oA3 + (A.indexes[2][i2]-1)*strdA2
            oB2 = oB3 + (B.indexes[2][i2]-1)*strdB2
            for i1 = 1:s1
                C[i1+o2] = parent(A)[A.indexes[1][i1]+oA2] + parent(B)[B.indexes[1][i1]+oB2]
            end
        end
    end
end

Here I've not required that SubViews be strided, i.e., the indexes could be AbstractVectors.