`accumarray` makes anomalous calls to its function argument

Question

Short version: The function passed as the fourth argument to `accumarray` sometimes gets called with arguments that are not consistent with specifications encoded the first argument to `accumarray`. As a result, functions used as arguments to `accumarray` must test for what are, in effect, anomalous conditions. The question is: how can an a 1-expression anonymous function test for such anomalous conditions? And more generally: how can write anonymous functions that are robust to `accumarray`'s undocumented behavior? Full version: The code below is a drastically distilled version of a problem that ate up most of my workday today. First some definitions: ``` idxs = [1:3 1:3 1:3]'; vals0 = [1 4 6 3 5 7 6 Inf 2]'; vals1 = [1 Inf 6 3 5 7 6 4 2]'; anon = @(x) max(x(~isinf(x))); ``` Note `vals1` is obtained from `vals0` by swapping elements 2 and 8. The "anonymous" function `anon` computes the maximum among the non-infinite elements of its input. Given these definitions, the two calls below…

Accepted Answer

Short answer

The fourth input argument of `accumarray`, `anon` in this case, must return a scalar for any input.

Long answer (and discussion about index sorting)

Consider the output when the indexes are sorted:

```
>> [idxsSorted,sortInds] = sort(idxs)
>> accumarray(idxsSorted, vals0(sortInds), [], anon)
ans =
     6
     5
     7
>> accumarray(idxsSorted, vals1(sortInds), [], anon)
ans =
     6
     5
     7
```

Now, all the documentation has to say about this is the following:

If the subscripts in subs are not sorted, fun should not depend on the order of the values in its input data.

How does this relate the trouble with `anon`? It is a clue, as this forces `anon` to be called for the complete set of values for a given `idx` rather than a subset/subarray, as Luis Mendo suggested.

Consider how `accumarray` would work for a non-sorted list of indexes and values:

```
>> [idxs vals0 vals1]
ans =
     1     1     1
     2     4   Inf
     3     6     6
     1     3     3
     2     5     5
     3     7     7
     1     6     6
     2   Inf     4
     3     2     2
```

For both `vals0` and `vals1`, the `Inf` belongs to the set where `idxs` equals 2. Since `idxs` is not sorted, it does not process all values for `idxs=2` in one shot, at first. The actual algorithm (implementation) is opaque, but it seems to start by assuming that `idxs` is sorted, processing each single-valued block of the first argument. This is verifiable by putting a breakpoint in `fun`, the function reference by fourth input argument. When it encounters a 1 in `idxs` for the second time, it seems to start over, but with subsequent calls to `fun` containing all the values for a given index. Presumably `accumarray` calls some implementation of `unique` to fully-segment `idxs` (incidentally, order is not preserved). As kjo suggests, this is the point where `accumarray` actually processes the inputs as described in the documentation, following steps 1-5 here ("Find out how many unique indices there are..."). As a result, it crashes for `vals1`, when `anon(Inf)` is called, but not for `vals0`, which instead calls `anon(4)` on the first try.

However, even if it followed those steps exactly on the first go, it would not necessarily be robust if a complete subarray of values contained just `Inf`s (consider that `anon([Inf Inf Inf]`) returns an empty matrix too). It is a requirement, although an understated one, that `fun` must return a scalar. What is not clear from the documentation is that it must return a scalar, for any inputs, not just what is expected based on the high-level description of the algorithm.

Workaround:

```
anon = @(x) max([x(~isinf(x));-Inf]);
```

`accumarray` makes anomalous calls to its function argument

Solution

Problem