Fix findmax and findmin with iterables

[nalimilan]

The state is not guaranteed to be equivalent to the index of the element.
Fixes #14085

I wonder whether iterstate could be removed, as it doesn't appear to be used anywhere after this change, and it's potentially misleading.

[tkelman]

ref #13203, @timholy is this okay?

also you'll need to rebase any pr's you opened in the last day to fix appveyor nevermind looks okay now

[timholy]

This could just be i and you could get rid of iold.

[nalimilan]

Of course. I fixed the code a bit too directly.

[timholy]

Yes, since you're now using an integer counter, it's fine to get rid of iterstate.

[nalimilan]

@timholy I wonder why you didn't write the code that way in the first place? Did you have a more complex case in mind? Or for performance? My version sounds so terribly simple. :-)

I've updated the PR.

[nalimilan]

Actually, this PR loses the support for array indexing, it will only return the linear index. It looks like a simple solution would be to have the current method as a fallback for generic iterables, and add a specialized method for AbstractArray using for i in eachindex(A) and returning the index for the maximum. Does that sound OK?

We also need to make a decision as regards backporting, as this is technically an API break -- even if leaving it broken for ranges isn't acceptable either.

[timholy]

Nah, I think your version is just plain better 😄. I suspect I was wrestling with the idea of whether we should be returning the state (e.g., a CartesianIndex) rather than an integer index, since in some cases a[i] will be an error if i is an integer. But this would be breaking change, so I should have done the rewrite the way you did.

Thanks for doing this!

[timholy]

I think findmax has always returned an integer index, so for now backward compatibility wins (esp. for backports).

But you're right that a better solution would be nice.

[nalimilan]

@timholy OK, cool. This commit is good for backporting. But I wonder whether we shouldn't improve the behavior in master and return a cartesian index for arrays. That sounds more natural than a linear index (and indeed ?ind2sub gives precisely this example). Adding these lines could be enough:

    if isa(a, AbstractArray)
        return (m, ind2sub(size(a), mi))
    else
        return (m, mi)
    end

[timholy]

That wouldn't be type-stable, but it wouldn't be hard to make it so. In general I like the idea of returning a native iterator. One issue, though, is that this simple solution will break if/when we switch to supporting arrays with arbitrary indexes (i.e., not just 1:size(A, d) but also allow indexes which range, e.g., over -5:5).

[nalimilan]

Why wouldn't it be type-stable? The behavior would depend only on the type of the array. Regarding arrays with arbitrary indexes, we would have to take into account not only the size of the array, but the array itself (or its type) so that a specialized method could translate the index correctly. I can add a, say, readable_index(a, i) helper method which would also differ from ind2sub in that it would be a no-op for LinearSlow arrays, for which eachindex already returns cartesian indices.

I've also discovered one possible failure case: the one in which the iterable is also indexable, but the indexes differ from the rank of the elements. The only case I know of at the moment is strings, for which findmax doesn't sound very useful. But in an ideal world we would dispatch on an Indexable trait, in which case we would iterate over indexes rather than using the iteration protocol.

[timholy]

While this may change someday, isa does not currently get evaluated at compile time:

julia> function foo(x)
           if isa(x, AbstractArray)
               return 1
           else
               return 1.0
           end
       end
foo (generic function with 1 method)

julia> @code_warntype foo(rand(2,2))
Variables:
  x::Array{Float64,2}

Body:
  begin  # none, line 2:
      unless (Main.isa)(x::Array{Float64,2},Main.AbstractArray)::Bool goto 0 # none, line 3:
      return 1
      goto 1
      0:  # none, line 5:
      return 1.0
      1: 
  end::Union{Float64,Int64}

But if you define it as

foo(x::AbstractArray) = 1
foo(x) = 1.0

then it is type-stable.

[timholy]

Regarding arrays with arbitrary indexes, we would have to take into account not only the size of the array, but the array itself (or its type) so that a specialized method could translate the index correctly. I can add a, say, readable_index(a, i) helper method which would also differ from ind2sub in that it would be a no-op for LinearSlow arrays, for which eachindex already returns cartesian indices.

This is getting pretty close to the implementation we have now---the only difference is that I held onto the state and called sub2ind at the end, and you're proposing holding onto an integer and calling ind2sub at the end.

Maybe a better fix would be to leave the current one in place and write a special exception for Ranges? Really, for findmax there's no point iterating over each element of the range, because you can compute it directly (the maximum value is either first(r) or last(r)).

[nalimilan]

I don't think we can get away with this by making an exception for ranges: the fallback implementation should work in all cases and not rely on assumptions that are not in the API contract of iterables to support all possible custom types. (Of course, a more efficient implementation for ranges would be useful for other purposes.)

My proposal differs from the current implementation in two ways:

1. It would only apply the AbstractArray-specific behavior to those types.
2. It would return cartesian indices even for LinearFast arrays.

[timholy]

I guess there's still a part of me that wonders whether the state is the more fundamental object, as you yourself pointed out above for strings. But your plan makes a lot of sense. I have no objections.

[nalimilan]

But can I always call sub2ind, passing it the index obtained via eachindex, or do I need a special function for LinearSlow arrays?

(As regards the iteration state, we might consider requiring that it is equal to the index for indexable types.)

[timholy]

It would be a no-op for LinearSlow arrays, or perhaps it should extract the tuple from inside the CartesianIndex. However, if one is passing back the state, there is no need to call sub2ind---that's only if you want to pass back an integer.

(As regards the iteration state, we might consider requiring that it is equal to the index for indexable types.)

I wonder about that too. Anyone else want to chime in about whether this would be overly restrictive? Obviously Ranges would need some restructuring.

[toivoh]

I'm not sure how BitArray iteration is implemented right now, but I've got a feeling that nice iteration state could involve a (perhaps partially shifted) piece of the packed storage.

I don't think we should be in the business of prescribing the format of iteration state. Better in that case to introduce a function that will translate it into an index or similar, based on the type that is iterated over.

[StefanKarpinski]

+1 to what @toivoh said.

[timholy]

Sounds like we shouldn't assume that the state can be used as an index. So I'd say go ahead with your plan, @nalimilan, and we'll deal with difficulties from generalization as (or if) they arise.

[nalimilan]

OK. Could you tell me how I can extract the index tuple from the CartesianIndex?

[timholy]

julia/base/multidimensional.jl

Line 18 in 9c90d2c

I::NTuple{N,Int}

[nalimilan]

Thanks, I confused this simple structure with the more complex iterator state.

Though there are more issues to tackle:

• If we return a (cartesian index) tuple of indices for matrices and above, should keep returning a scalar for vectors? The inconsistency doesn't sound great to me (cf. e.g. zip() returns scalar instead of tuple when passed a single argument in 0.4 #13979), yet returning a tuple with one element isn't exactly what I'd call an improved usability.
• The same question applies to findmin(A, region). Should we keep returning an array of scalars when a single dimension is passed as region? Should an array of tuples be returned when several dimensions are passed?

At least always returning a linear index was simple. :-)

[timholy]

I don't know whether this makes it more or less complicated, but there's findn for the give-me-cartesian-indexes case. So we could just return linear indexes and call it good enough. Or we could be ambitious and clean this mess up, if only we figure out how it should be done 😄.

[timholy]

Maybe one could specify the desired output type, e.g.,

findmax(Int, a)
findmax(CartesianIndex, a)

[nalimilan]

Maybe one could specify the desired output type, e.g.,

That would be a possibility (and could also lead to replace findn(x) with find(CartesianIndex, x)), but that still doesn't answer the question about what to return by default. The current solution of returning a linear index is consistent with find, and it's predictable. Its drawback is that it's not necessarily the most efficient for LinearSlow arrays, though I'm not sure it matters in practice (maybe it would be more important for find?).

So I can see two solutions, to apply identically to find, findmax and findmin:

1. Keep returning a linear index as now and with this PR.
2. Return whatever eachindex returns by default.

In both cases, we can optionally support passing Int or CartesianIndex as the first argument to request a specific index type.

We can also consider in this discussion findmax(x, region), which currently returns a linear index, where it could also usefully return the index of the maximum on the dimension over which the reduction is performed. find(CartesianIndex, x, region) could work here too.

[timholy]

Yes, I also favor returning a linear index by default.

[nalimilan]

So, let's say we go with (and backport) this PR, and will tackle the question of findmax(CartesianIndex, x) later. I'll merge it if nobody objects.

[timholy]

👍