feature request findmax(A,dims)

[floswald]

getting the maximum and its location in one pass is very useful. at current it's not supported along dimension d of an array.

A = rand(3,3,3)
now = mapslices(findmax,A,3)
#or
for i in 1:3
    for j in 1:3
         now[i,j] = findmax(A[i,j,:])
    end
end
then = findmax(A,3)

it's odd to have

max(A,dims)

and not findmax(A,dims)

[timholy]

If no one else beats me to it, I'll tackle it eventually, so I'm assigning it to myself.

[ivarne]

Up for grabs? This seems like a great start for anyone who wants to contribute to Julia, but needs a easy start.

[floswald]

Sure ill give it a shot. I just submit a pull request when im done?

sent from mobile
On May 2, 2014 11:26 AM, "Ivar Nesje" notifications@github.com wrote:

Up for grabs? This seems like a great start for anyone who wants to
contribute to Julia, but needs a easy start.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/6716#issuecomment-42014347
.

[timholy]

That would be awesome. @floswald, of course you're free to implement this any way you please, but if you'd like a recommendation: you're basically wanting to implement reduction with two output arrays, one containing the value of the maximum/minimum and the other the index at which it occurs. Look at gen_reduction_body in reducedim.jl. Depending on your experience, you may have some trouble parsing all those macros, but http://docs.julialang.org/en/latest/devdocs/cartesian/ should help you get started. I highly recommend doing cut/paste and using macroexpand to help understand what's happening in that code---it's actually really simple, as the brevity of the function presumably suggests, but the syntax is not pretty.

To wrap your "core" in a nice function, see the 3-line definition of maximum! (for example) at the bottom of the same file.

[floswald]

ok sounds great - i'm on it. i'm away for the weekend though, so if anyone needs this pronto, feel free to preempt me. thanks a lot for the pointer to Base.Cartesian. I had no idea about the existence of this - it's massive.

[floswald]

hi guys,
i'm a bit stuck here. I cannot run gen_reduction_body with macroexpand for some reason, even after allocating A and R, it tells me wrong number of arguments? If I understand correctly gen_reduction_body should be very close to what I want. I think the crucial lines for me to change are

    @nloops $N i A d->(j_d = sizeR_d==1 ? 1 : i_d) begin
    @inbounds (@nref $N R j) = ($F)((@nref $N R j), (@nref $N A i))

but what are they doing? I guess I would have (@nref $N R j) = findmax((@nref $N A i)) or so, but how do I know which dimension I'm reducing over? what does this do:

@nextract $N sizeR d->size(R,d)

sorry i couldnt find any docs on @nextract.
asking questions here may be highly inefficient (certainly compared to you doing it yourself), do you want me to post somewhere else?

[timholy]

For me this works:

using Base.Cartesian
macroexpand(Base.gen_reduction_body(2, max))

I guess I would have (@nref $N R j) = findmax((@nref $N A i)) or so

Since findmax allocates memory, you don't want to call it in the inner loop. A better plan would be to define your function to simply be a comparison operator, and maintain two reduced arrays, R will contain the value and I will contain the index. Do the equivalent of the following:

maxval = X[1]
maxind = 1
for i = 2:length(X)
    Xval = X[i]
    if Xval > maxval
        maxval = Xval
        maxind = i
    end
end

In "Cartesian" this looks like

Rval = @nref $N R j
Aval = @nref $N A i
if ($F)(Rval, Aval)
     (@nref $N R j), (@nref $N I j) = Aval, k
end

where k (the linear index) is an integer that gets incremented on every iteration.

sorry i couldnt find any docs on @nextract.

I didn't document @nextract because it's redundant with @nexprs:

@nextract $N sizeR d->size(R,d)

is the same as

@nexprs $N d->(sizeR_d = size(R,d))

(I wonder if we should deprecate @nextract.)

[floswald]

Hi @timholy
i just submitted a pull request. i just noted i should have submitted that not from my master branch of my julia fork but from branch "findmax(A,dim)" ...can't change that now. anyway, I implemented something quite close to your suggestion. The biggest problem was to figure out the linear indexing rule. I had some substantial help from my buddy @tlamadon with this. So what it does now is to hold one dimension d constant, looping over all dims "before" and "after" d, ordered according to size(A). I tested the results against mapslices(Base.findmax, A, d) and they all pass. The timing tests are quite encouraging, as you might have guessed.

[floswald]

still trying to understand your cartesian snippet:

Rval = @nref $N R j
Aval = @nref $N A i
if ($F)(Rval, Aval)
     (@nref $N R j), (@nref $N I j) = Aval, k
end

can I ask what ($F)(Rval, Aval) stands for?
is that a function which takes R[j_1,j_2,...],A[i_1,i_2,...] as arguments?
why is there an if clause?
how does k get increased? i was looking to get this to work

macroexpand(:(@nloops 2 i x begin P, k = @nlinear 2 x i end))

but seems outdated (or wrongly used). sorry if this becoming too much work for you.

[timholy]

F = < for example, if you're implementing findmax. (Or you could make F = > and reverse the order of the arguments, if that's more intuitive.) It's the same as the snippet under "Do the equivalent of the following" above.

For increasing k, just have a k += 1 somewhere in the body of the loop. Don't use @nlinear, that's not even in Base (and just incrementing k will obviate the need for @nlinear, and be much more efficient as well).