Create a single API/test/benchmarking but with multiple independent implementations

[d1manson]

I've created a separate branch for working on this.

So far I've spent a few hours trying to write a fairly comprehensive readme plus I've done some very basic work to split up the implementations into separate files. Note though that I've made basically no attempt at actually enforcing the API yet, or even checking if the code runs ok...it doesn't!

Obviously there is still quite a bit of work to be done.

In terms of agreeing on a single API, I think the main questions are:

1. What to do with dtypes? It's always nice to give the user full control over this, but it's also great when the program is able to automatically choose the most sensible option. For example, summing int8's probably need to be upgraded to int32 or int64, maybe even float..the program could guess what is best based on some heuristic, although this would add extra time - the simplest heuristic would be to assume the worst case i.e. all values are already maximum at the current bitdepth and all will be added to the same group, thus if len(idx)=100 and vals.dtype=int8 then we assume the sum is going to need log2(100*255) bits. This calculations is specific to sum, mean var std have different requirements, and min max first last all any are easier to guess. The other question is what to do about the dtype of the fillvalue - I have chosen to enforce bools for the fillvalue of any all but, the user may legitimately want some form of nan or pseudo-nan (I remember reading a discussion about efficient nans for non-floats somewhere in pandas docs or maybe it was to do with numpy's MaskedArray ..anyway there's no simple answer to that question). If the fillvalue does require using float in the output, it may still be more efficient to do the main computation with some for of int and then only upgrade to float at the end. Perhaps the easiest thing to do is upgrade all the ambiguous stuff to at least float32 and then encourage the user to request a simpler dtype if they fell that would suffice...but then what about detecting overflows etc..I'm not sure what error handling of that kind is built into numpy?

2. Support for non-flat vals, matched to idx Although matlab allows this, I would prefer not to. Assuming both idx and vals have the same memory layout, it costs nothing for the user to flatten them, and it isn't too ugly syntactically. By enforcing that it then makes it simpler to identify when 2D idx is being used to request multidimensional output...this was always a bit confusing in matlab...see next point...

3. Support for multidimensional output, using 2D idx This is available in matlab and is a very helpful feature. I use it a lot for 2d binning, but it could be n-dimensional. I feel strongly that we should provide it, at least in the API (i.e. throw NotImplemented if need be..but it's actually very easy to do once numpy is involved!).

4. Support for scalar vals - I wrote a small section on this in the readme. I think it's not really that useful in general, but helpful for sum of vals=1...you might argue that the user could just go and use bincount, but Matlab accepts scalars and this version of accumarray can help with multidimensional bincounting, so it's probably a feature worth having. If need be, you could put a simple hack at the top of the function, which expands the scalar out to be a full vector..or throw NotImplemented.

5. Support broadcasting vals to idx This is complicated to implement and causes ambiguities if you also want to support 2D idx for multidimensional output (disucssed above). Matlab does not offer this, and I've never felt the need to try and use/implement it, so I strongly feel we should explicitly choose not to do it. It also means that the ufunc.at optimizations (if they ever happen) can be kept simpler....to be honest I'm not even sure whether broadcasting makes any sense in the accumarray context?!

6. Support for contiguous/sorted mode. I haven't mentioned this in the readme and I haven't implemented it yet, but I agree that something along these lines would be nice. There are a variety of things you might want, for example if the user says that idx is pre-sorted, various optimisations may be possible (i.e. the code still treats the inputs the same way, but is now able to assume that groups are already laid out contiguously in memory). Alternatively the user may want your contiguous concept, where idx is no longer an actual index, but now just a general label for contiguous sections. A variation on this would be to say that sections where idx=0, should be treated as belonging to a common "null" group, but all other contiguous blocks should be treated as being labeled 1, 2, ..., n. ..this can work nicely if idx is a bool array, as each true block is treated as a new label, but all the false blocks are grouped together...but sometimes the user will want a way to have two different groups adjacent to one another (i.e. no false element between them)..which you can't do with a simple bool array....anyway, the point is that there's a lot you could do here to make the user happy, some of which will offer extra optimizations, and some of which would just allow the user to be lazier.

7. User-specified size and bounds checking - I think it's definitely worth allowing the user to specify the size of the output, Matlab does this. I often need to do this when I'm doing analysis because I use accumarray twice, the second time with a subset of the original data, and I want the outputs to be the same size. I'm not sure exactly how much control you can give the user over bounds-checking...I guess with your weave stuff you could completely disable all checking if the user really wanted. If that's something you want to offer we should encapsulate it in a proper kwarg, however I think in general people will expect basic checking as I think all numpy code does. edit 0: by the way I don't like your choice of "incontiguous" as default...the english is _non_contiguous, but it's a bit of a confusing choice of default. I also haven't looked into the details what your downscaled mode offers...I'm not sure where it fits in to the picture.

---

edit 1:

Some of the above stuff probably warrants common scripts for numpy and weave, though pure python solutions would have to be kept separate (and in many cases may be simple/not applicable at all.). I'm not sure of the best way to deal with that given that we want to try and keep stuff as separate as possible...there may need to be some kind of compromise.

One thing that should be fairly easy to agree on is the order of args, I vote for following Matlab exactly and then adding on all the extra kwargs after that (order, dtype, mode).

Also, we should agree on the aliasing rules and exactly what the nan- versions of functions are supposed to do. I tried to describe that fairly explicitly in the readme...let me know if you disagree and/or want to add stuff.

Regarding the name of the main two arguments, I definitely prefer idx and vals, but if we permit a range of modes then idx is no longer necessarily such a great choice...but I'll stick with it for now!

Regarding the weave stuff, is there a good place to point people to in order to get it installed quickly...so far I haven't said anything about this in the readme...it's probably worth giving people a one-sentence explanation as to exactly why weave is so good.

I think you also mentioned a numba implementation..what did you say the status of that was? And what about a Cython one?..although I'm not sure there's going to be much to be gained from Cython. I have a small amount of experience with Cython and a tiny amount of experience with numba, but I don't feel like I really know enough to compare the pros and cons of each...if you are going to provide all these implementations, we should have a proper explanation as to why you might want each of them....I guess one of the advantages worth mentioning is that (some/all? of these things) work nicely on the GPU...though I think the simplicity of going to the GPU varies greatly?

In terms of making the C versions more useful to people, might it be worth providing some binaries..I guess that's a pretty big job, but it might increase uptake.

---

edit 2:

Ok, I didn't realise that the the non-numpy implementation did actually use numpy. ...I've now re-written it to be truly pure-python. It's relatively simple, though obviously not exactly performant.

edit 3:

I think it may be worth trying to wrap pandas. It looks like pandas uses cython under the hood for its groupby- see ..pandas/ib.pyx and ..core/groupby.py

---

edit 4:

Right, I've implemented a simple pandas version. I didn't spend any time trying to squeeze the best performance out of pandas, so it may be possible to do better. I've put benchmarking stats in the readme, but the basic story is that it doesn't beat the optimized numpy version, except where numpy has to rely on ufunc.at (ie. min max prod).

I've also implemented a ufunc-at version, for use with benchmarking.

[ml31415]

Wow, great job so far! Sorry for the delay, I was on a conference over the weekend. I'll go through step by step, hope I don't miss anything.

1. dtypes
We should do our best to provide some reasonable defaults, and it has to be the same default for all implementations, so this needs to go into a separate function. On the other hand, I guess it's fair to just take int64/float64 for sum/prod, and let the user choose himself to reduce this, if he knows the result to be smaller than that. Any heuristic will require inspection of the data beforehand. The better the assumption shall be, the more time costly it is, and the shorter, the more error prone to overflows it will be. I wouldn't want to do that.

For bool and int, I wouldn't bother with a separate nan type. If in doubt, the user can run an anynan over the same data, and use this as a mask afterwards. End of story. In general, the fillvalue must be representable with the selected output dtype. So the fillvalue should also be an input for the shared dtype selection function.

2. Matching vals to idx
Full acc. KISS principle.

3. ndim
Yes, should be in the API.

4. Scalar vals
The simple hack sounds reasonable. Not sure if I would bother spending much time in writing any optimized code for that rare case.

5. Broadcasting
KISS principle, let's forget this.

6. Contiguous mode
The downscaled thing is mainly doing a compression of the outputs, skipping all empty fields, by running unique over accmap first. It could be supported just like the scalar vals, by simply modifying accmap before acutally running the main action. It's a nicety, and I also won't bother to throw it out.

For the contiguous mode, this was my initial implementation, requireing an already sorted input. After hacking around a bit longer and improving the implementation for unsorted accmaps, it turned out, the speed advantages were actually not that great as I had thought, and I could get comparable speeds for arbitrary accmaps. It's basically one pointer indirection, that can get skipped, when accessing the data. Looks like the prefetching etc. is that great, that this is not a big cost anymore. As it is mainly a subcase of the unsorted map. It would make sense to throw that out completely.

Renaming incontiguous to something better, I'm totally fine with that. I'm not a native speaker, so I beg my pardon for the naming.

For your other suggestions on this, I'd say KISS principle again. If the user wants some 0-group, simply set some 0-values in the accmap, before running accumarray over it.

7. Bounds checking
As the functions in the inner loop, adding simple values are really simple, the whole process is mainly memory bound. So adding a bounds check is not a great overhead. E.g. the nansum and sum run in nearly the same speed, even if nansum has an additional val != val comparison inside. So for the sake of simplicity, and as segfaults are really ugly, I wouldn't consider a possibility to turn bounds checking off. It should be safe to use the library, even if there were some spooky values in your accmap.

So, I'm totally fine adding the sz argument.

8. Input variable naming
vals, idx, accmap, a ... I like speaking variable names, when ever possible. idx is quite generic, used million of times anywhere in any code, it doesn't describe well what it is meant for. Iterating over it idx[i] ... what do I get? accmap is more speaking. It's a map, how to accumulate vals, accumulation map. I'd clearly prefer that over idx. Matlab calls it subscripts, that may be another good option. vals vs a, both don't give any information, so they are equally bad I'd say. Though I wouldn't have any better suggestion on that right now, I'm fine with both names. Finding something better couldn't hurt though.

It may be a good idea to change fillvalue to fill_value, as used in numpy e.g. for masked arrays.

[ml31415]

For getting scipy weave to run, it's actually a single line with linux, not more than "apt-get install python-scipy" and everything works. I also know, what kind of pain in the ass it can be, to get toolchains or math libraries installed in windows. I haven't used windows for ages, and can't help too much with setting up scipy for it. The documentation there should cover that I hope. What weave generally does, is creating C-code on the fly, depending on the input data type, and storing the compilation results for reuse. A working C compiler is necessary for all this to work, so there is no way to provide binaries, that would run without a working toolchain.

I also wouldn't worry at all about the purepy implementation. Any sane person, trying to do some math with python will install numpy. If not - not my problem. No need to reinvent the wheel, just for a demonstration version of accumarray. If some things don't work there, so be it.

In terms of aliasing, you had the clearer approach. We should have a shared function, of converting any func-keyword-input into a unique function name string, with which every implementation has to deal itself then. Your implementation is a good start for that.

[d1manson]

Your English is pretty much perfect (I actually spent a year in Moscow as an EFL teacher so I should know..it's certainly a lot lot better than my Russian!).

Ok...

I think we need to flesh out the accumarray_utils.py with the common numpy stuff...but bear in mind that I think it is a nice feature to support a fully non-numpy implementation...there are some environments where it might be a pain to get numpy working, but where it's nice to have some simple accumarray for smallish datasets..anyway the code is pretty short and simple, so I don't think it's a big deal...it's also quite fun to benchmark it against numpy!

1. dtypes bool - agreed. sum - it's pretty simple to get a safe bound just using len(vals) and the number of bits in vals, for example 30 000 int16s can be safely summed into an int32, or 60 000 uint16s...but, yes, I think we should put all the logic in one utility function and work on it there..with regards dtype+fillvalue, my _get_minimum_dtype function is a basic first draft.

6. modes So, are you suggesting that the only thing we should support is something like this:

if <whatever_the_special_test_is>:
   idx_is_used = np.zeros(sz,dtype=bool)
   idx_is_used[idx] = True
   new_idx = np.cumsum(idx_is_used) - idx_is_used
   idx = new_idx[idx]
   sz = max(idx)

If that's what you mean, then could we treat it as a special case of fillvalue, e.g. fillvalue=None, means drop all the missing values...although for multidimensional indexing you could still end up with some empty spots..maybe we could just say that it's not supported for multidimensional output. The only other reason not to use fillvalue=None would be for the func='array' or func='sorted' cases, because there the user may actually want to use None as the fillvalue...but I don't think the user will be too annoyed if they cant have that.

Just to clarify, if the above is true, then there is no mode needed at all, right?

8. variable naming I really didn't like a, but I guess you are right that vals is equally meaningless!...not sure what else to suggest there. But for, the indices, how about dest_idx (for destination index - this is in the spirit of C/asm-like naming when you provide output pointers as input args)..I felt like "accumulation map" was a bit confusing - it sounded to me either like the name of a function/lambda or else the output of the accumulation (rather than the input to it).

I'm probably not going to do any more work on this in the next few days, so you can safely commit to the branch I was working on and it shouldn't be too confusing.

[d1manson]

9. benchmarking and testing I'm sure your stuff was good, but I was lazy and ended up basically starting from scratch in a separate file rather than trying to understand everything you were doing.
As discussed initially, we should aim to get testing to a pretty high standard.

[ml31415]

With start from scratch you mean the test_and_bench_simple module? I agree that testing and benchmarking should be separated. I'm also fine when the testing part with pytest will be up to me. I guess it could require some more comments. It's surely not too intuitive, what's going on there. For basic understanding I'd recommend pytest fixture and parametrization reads.

• https://pytest.org/latest/fixture.html
• https://pytest.org/latest/parametrize.html

accmap_all and accmap_compare are different fixtures. accmap_all simply returns different accum implementations for the tests. That are the first bunch of tests, testing basic functionality that has to be guaranteed by all implementations. This could also be achieved, by parametrizing all these tests separately with all implementations. Parametrizing the fixture just does this in a more convenient way.

accmap_compare provides a whole namespace with a bunch of testing ingredients, data with and without nan, depending on what shall be tested, a reference function etc. test_compare then runs all the listed functions in func_list against any combination of implementations provided by the accmap_compare fixture.

[ml31415]

fillvalue=None would be fine to me. This is what is done right now in that case:

accmap = np.unique(accmap, return_inverse=True)[1]

[ml31415]

Btw. please let's be more specific with the thrown errors in the future rework. It's really bad practice to only throw Exception. When I would like to catch any common problem with values, types, that might arise, and handle in a graceful way outside of accumarray, this is not possible right now. Catching any "harmless" error now requires having a catch all statement, which also catches any sort of other serious runtime errors, memory errors, io errors etc, which you're probably not able to handle in a graceful way.

[d1manson]

that all sounds good to me..

I stand by what I said above - I'm not going to touch anything for a couple of days at least so it should be easy for you to commit straight onto this branch

[ml31415]

Ok, I'll see where I get with this the next days.

I guess the last open thing is a final decision in terms of variable naming, I got the names initially from here: http://wiki.scipy.org/Cookbook/AccumarrayLike . So it's not really my invention. As it's probably the oldest public python accumarray implementation (dating back to 2010) and was promoted on scipy, it likely got some reach. I'd rather stick with that, then inventing something completely new, without good reason. dest_idx seems to me similarly descriptive, but also not really an improvement over accmap. Maybe one more suggestion how to see accmap in a more speaking way: Different values in accmap mark and group regions within a, just like different colors mark countries in an ordiary map.

a is besides being used in the cookbook quite commonly used throughout numpy and scipy for any array like input. cumprod, all, copy, zeros_like and most linalg, just to name a few, use a for that. Sticking with this naming scheme looks like the best idea to me, if there is no clearly better name. If we want to stress the fact, that we also support scalars as second argument, following this logic we might use x.

For some unclear reason I had taken the fillvalue without underscore, even if it's commonly used with it. This should be fixed into fill_value.

About the order= keyword, I guess this should not be added here. The conversion from C to F is cheap and easy, and can be done manually afterwards, if in doubt. This would save one more keyword. If I saw it correctly, there wasn't any logic implemented for it so far. For the input, C/F can be detected manually and acted upon accordingly. Maybe the output should just copy the ordering from the accmap?

impl_dict= shouldn't be part of the public API as you already wrote in the docstring. So I guess you don't mind prefixing an underscore there.

About the keyword order, I initially would have tended to a agree to keep matlab ordering. Thinking about it now, I guess this should be talked about again. I was once told, when I suggested some accumarray thing for numpy, that numpy is not matlab, and does not need to repeat it's quirks. func seems to me to be clearly the most used keyword argument, nearly being worth to be a mandatory, positional argument. size is something useful, but probably even less frequently used then fill_value. So while we could keep the order for fill_value and size in matlab style, we may have func as the first argument, just to allow it to be used without explicitly writing func= in front of it. And when we tend not to lean too close to the matlab implementation, it may be more speaking to use size= instead of sz=. All in all, the cookbook version syntax seems to be well thought through ...

If we could agree on all this, we would finally end up with:

def accumarray(accmap, a, func='sum', size=None, fill_value=0, dtype=None):
    ...

This are the conventions we talked about already, just to sum it up:

• fill_value=None triggers downscaling
• boundary checks always guaranteed
• a[accmap >= size] is treated as if it would be nan
• no special nan handling for bool / int
• no implicit matching of a to accmap

Are you fine with this?

[d1manson]

In no particular order...

I suspect that the cookbook accumarray hasn't been that widely adopted by people because it would have been improved upon already..and there's not that much activity on the SO question. When we get this repo sorted out we should endeavor to get it linked to from the cookbook page so people dont find that page and think that's the best thing out there.

C/F - I'm not sure if we're on the same page here... I have used it in the calls to ravel_multi_index and reshape at the start and end of the function respectively. I thought it would be nice for the user to be able to specify the layout of the ndarray that accumarray is creating from scratch. Admittedly the user can control this "manually" by passing in accmap[::-1,:] in order to flip the order of the output if need be, but that's a bit clumsy. Personally, my feeling is that order is a well understood/easily-ignored kwarg in numpy, so I don't feel like there's anything wrong with providing it here.

With the first two args, I guess you're right that a is used elsewhere in numpy, so maybe I should accept it. But I really don't like accmap! ufunc.at uses indices, which is basically what I am suggesting too. Here's some more options: group_idx grouping_idx target_idx out_idx at_idx or maybe a_label a_group a_idx a_index a_dest_idx..anything else?

I'm not too fussed about whether or not we match up with Matlab beyond the first two args, as long as we provide basically the same functionality. Python's superior argument passing syntax is basically doing a pretty good job already - my feeling is that different users are likely to use the other kwargs to different extents, so it's not that easy to provide a one size fits all optimal ordering...hence the thinking that maybe just sticking to Matlab's order is maybe not a terrible idea...but as I say, I'm not too fussed.

I'm not sure what you mean by "a[accmap >= size] is treated as nan"...where accmap>=size, the boundary check should throw an exception, surely? If you did treat those values as nans then when would you put them in the accumulated output?

To clarify the nans in bools and ints question...it's only relevant for the fill_value. And note that the user might want to choose +-Inf as the fill value as well as nan. My feeling on this may have changed a bit...perhaps we should upgrade to at least float32 if the user insists that they want non-finite values in their output array...both for bool and int. By default the fill_value is 0, so it would be the user's fault if they force the change in the dtype. We could always implement it as a second step at the end where you convert form int/bool to float and set the fill_value.

And what did you mean by the last point - "no implicit matching of a to accmap"?

[ml31415]

C/F is not a big deal in pure numpy operations to support, as numpy takes care of all iteration issues related to it. Having to take care about this in C is quite an annoyance, I wouldn't like to take care of without good reason. Everyone really taking care about the order can easily change the order afterwards imho. I had a bug in the C version once, not correctly recognizing F arrays, producing wrong results, so that's why I said autodetection and acting upon must be granted. Nothing to do for the numpy version here.

Let's pick group_idx then. Seems like a well descriptive choice.

What I meant with treat like nan was, ignore these values. If the size of the output is already given with size= and not determined by the biggest value of group_idx, there is no place in the output to accumulate anything there. It's like cropping the output. It's the current behaviour, nothing special. I just wanted to mention it here.

If dtype is specified, I guess we agree that specifying an invalid fill_value should simply raise a ValueError. If the user specifies dtype=int and fill_value=-inf, what else to do with it?

For more vague situations, fill_value=-inf for a bool operation, no dtype specified, there are points for both. I don't like too much implicit behaviour. Who knows, where the fill_value came from. Some calculation, that went wrong? A bool operation simply doesn't produce nan (or any other float) as an output. Neither should accumarray do, if not explicitly told to do so. Same with sums/prod on int and float values as fill_value. If the user insists on specifying +Inf as a fill_value, I'd expect at least the output dtype to be properly specified, so it's clear what should happen here. The analogy would be setting float values to an int. If it's a more or less sane operation, the value is rounded and set. If it's nan/inf, there is an error. Anyways, this is probably too much detail, we should probably skip this point, until there is some implementation for the ambiguation function.

No implicit matching, if the shapes don't fit together for a sane operation, there shouldn't be any fancy magic making it somehow work, but just a clear error.

[ml31415]

One more thing. I totally don't care about my name being presented anywhere in the code, readme whatever, or someone copying. This is such a small snippet, not worth to claim any authorship on anything imho. So I'd suggest scratching any name mentionings in the readme and the code. The git history is descriptive enough, if someone really cares. If you'd prefer to have it in the code itself, I'd rather suggest something explicit for that like contributors.txt.

[d1manson]

I don't like silently cropping for out-of-range indices. It's likely to hide bugs in users' code. If they really want to do that then they can do a simple masking before they pass the values in to accumarray - we could give an example of this in the docs. Alternatively you could provide a ignore_out_of_range kwarg, so that it's very clear they are potentially losing data. bincount uses the minlength approach, which lets the user specify the minimum size and then returns a larger array if necessary...that's also a bit weird but at least you don't silently drop data. ..I'm saying we should throw an exception for out of range.

I think I agree with everything you just said about dtypes, let's see what the code ends up looking like and discuss further if need be.

With the order I'm still not sure we're talking about the same thing- you can't autodetect the order of the output because it doesn't exist yet....you can autodetect the input, and you should, but I leave that up to you. Regarding the output, there are a few cases where the user would be grateful for being able to choose which dimensions are contiguous in memory:

• e.g. the next step after accumarray is to smooth along a subset of dimensions (e.g. 2D planes in a 3D stack, or 1D signals in a 2D stack.)..could also be more general image processing or peak detection.
• e.g. the next step after accumarray is another aggregation (e.g. firstly accum mean into 2D, then take mean of those means along a specific dimension - note the effect of doing it in two steps is to weight each bin in the intermediate result equally.)

I'm not sure how you were planning on doing the multidimensional indexing, but if you use ravel_multi_index and reshape then you actually don't need to worry about the order in C because you only have to deal with 1D indexing. If you're keen to write a single-pass loop over the data then I think it's fairly easy to inline the ravel_multi_index - it's really just a dot product of the index with the shape of the output.

Regarding authorship - I'm fine with removing all references in code, but it would be good to leave one line in the readme somewhere...not a big deal though.