
http://www.mail-archive.com/numpy-discussion@scipy.org/msg17595.html
Prompted by the thread above, I decided to see what it would take to implement ufuncs with masking in C. I described the result here:
http://www.mail-archive.com/numpy-discussion@scipy.org/msg17698.html
Now I am starting a new thread. The present state of the work is now in github: http://github.com/efiring/numpy-work/tree/cfastma
I don't want to do any more until I have gotten some feedback from core developers. (And I would be delighted if someone wants to help with this, or take it over.)
1) The strategy I have started with is to make a full set of masked ufuncs alongside the existing ones, appending "_m" to their names. Only the binary ufuncs are implemented now, but the unary ufuncs can be handled similarly. Example:
multiply(x, y, out) # present ufunc: no change multiply_m(x, y, mask, out) # new
Where mask is True, the operation is skipped.
2) I have in mind the possibility of supporting two input masks and one output mask for binary operations. This would look like:
multiply_mm(x, y, maskx, masky, out, outmask)
outmask would be the logical_or of maskx and masky, and in the case of domained operations it would also be True where the arguments are outside the domain.
This form would provide the fastest support for masked arrays, but would also take quite a bit more work, and would expand the namespace even more. I'm not sure it's worth it.
3) I have not yet taken any steps to modify numpy.ma to take advantage of the new ufuncs, but I think that will be quite simple.
4) Likewise, to save time, I am now just borrowing the regular ufunc docstrings.
5) No tests yet, Stefan. They can be added as soon as there is agreement on API and general strategy.
6) The present implementation is based on conceptually small modifications of the existing numpy code generation system. It required a lot of cut and paste, and yields a lot of nearly duplicated code. There may be better ways to do it--especially if it turns out it needs to be redone in some modified form.
Eric

On Fri, May 15, 2009 at 7:48 PM, Eric Firing efiring@hawaii.edu wrote:
http://www.mail-archive.com/numpy-discussion@scipy.org/msg17595.html
Prompted by the thread above, I decided to see what it would take to implement ufuncs with masking in C. I described the result here:
http://www.mail-archive.com/numpy-discussion@scipy.org/msg17698.html
Now I am starting a new thread. The present state of the work is now in github: http://github.com/efiring/numpy-work/tree/cfastma
I don't want to do any more until I have gotten some feedback from core developers. (And I would be delighted if someone wants to help with this, or take it over.)
Here the if ... continue needs to follow the declaration:
if (*mp1) continue; float in1 = *(float *)ip1; float in2 = *(float *)ip2; *(float *)op1 = f(in1, in2);
I think this would be better as
if (!(*mp1)) { float in1 = *(float *)ip1; float in2 = *(float *)ip2; *(float *)op1 = f(in1, in2); }
But since this is actually a ternary function, you could define new functions, something like
double npy_add_m(double a, double b, double mask) { if (!mask) { return a + b; else { return a; } }
And use the currently existing loops. Well, you would have to add one for ternary functions.
Question, what about reduce? I don't think it is defined defined for ternary functions. Apart from reduce, why not just add, you already have the mask to tell you which results are invalid.
Chuck

Charles R Harris wrote:
On Fri, May 15, 2009 at 7:48 PM, Eric Firing <efiring@hawaii.edu mailto:efiring@hawaii.edu> wrote:
http://www.mail-archive.com/numpy-discussion@scipy.org/msg17595.html Prompted by the thread above, I decided to see what it would take to implement ufuncs with masking in C. I described the result here: http://www.mail-archive.com/numpy-discussion@scipy.org/msg17698.html Now I am starting a new thread. The present state of the work is now in github: http://github.com/efiring/numpy-work/tree/cfastma I don't want to do any more until I have gotten some feedback from core developers. (And I would be delighted if someone wants to help with this, or take it over.)
Chuck,
Thanks very much for the quick response.
Here the if ... continue needs to follow the declaration:
if (*mp1) continue; float in1 = *(float *)ip1; float in2 = *(float *)ip2; *(float *)op1 = f(in1, in2);
I was surprised to see the declarations inside the loop in the first place (this certainly is not ANSI-C), and I was also pleasantly surprised that letting them be after the conditional didn't seem to bother the compiler at all. Maybe that is a gcc extension.
I think this would be better as
if (!(*mp1)) { float in1 = *(float *)ip1; float in2 = *(float *)ip2; *(float *)op1 = f(in1, in2); }
I agree, and I thought of that originally--I think I did it with continue because it was easier to type it in, and it reduced the difference relative to the non-masked form.
But since this is actually a ternary function, you could define new functions, something like
double npy_add_m(double a, double b, double mask) { if (!mask) { return a + b; else { return a; } }
And use the currently existing loops. Well, you would have to add one for ternary functions.
That would incur the overhead of an extra function call for each element; I suspect it would slow it down a lot. My motivation is to make masked array overhead negligible, at least for medium to large arrays.
Also your suggestion above does not handle the case where an output argument is supplied; it would modify the output under the mask.
Question, what about reduce? I don't think it is defined defined for ternary functions. Apart from reduce, why not just add, you already have the mask to tell you which results are invalid.
You mean just do the operation and ignore the results under the mask? This is the way Pierre originally did it, if I remember correctly, but fairly recently people started objecting that they didn't want to disturb values in an output argument under a mask. So now ma jumps through hoops to satisfy this requirement, and it is consequently slow.
ufunc methods like reduce are supported only for the binary ops with one output, so they are automatically unavailable for the masked versions. To get around this would require subclassing the ufunc to make a masked version. This is probably the best way to go, but I suspect it is much more complicated than I can handle in the amount of time I can spend.
So maybe my proposed masked ufuncs are a slight abuse of the ufunc concept, or at least its present implementation. Unary functions with a mask, which I have not yet tried to implement, would actually be binary, so they would have reduce etc. methods that would not make any sense. Is there a way to disable (remove) the methods in this case?
Eric
Chuck
Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

On Sat, May 16, 2009 at 03:02, Eric Firing efiring@hawaii.edu wrote:
Charles R Harris wrote:
Here the if ... continue needs to follow the declaration:
if (*mp1) continue; float in1 = *(float *)ip1; float in2 = *(float *)ip2; *(float *)op1 = f(in1, in2);
I was surprised to see the declarations inside the loop in the first place (this certainly is not ANSI-C), and I was also pleasantly surprised that letting them be after the conditional didn't seem to bother the compiler at all. Maybe that is a gcc extension.
I believe they are a part of C99.

2009/5/16 Robert Kern robert.kern@gmail.com:
On Sat, May 16, 2009 at 03:02, Eric Firing efiring@hawaii.edu wrote:
Charles R Harris wrote:
Here the if ... continue needs to follow the declaration:
if (*mp1) continue; float in1 = *(float *)ip1; float in2 = *(float *)ip2; *(float *)op1 = f(in1, in2);
I was surprised to see the declarations inside the loop in the first place (this certainly is not ANSI-C), and I was also pleasantly surprised that letting them be after the conditional didn't seem to bother the compiler at all. Maybe that is a gcc extension.
I believe they are a part of C99.
Exactly (so not supported by Visual Studio).
Matthieu

Eric Firing wrote:
That would incur the overhead of an extra function call for each element; I suspect it would slow it down a lot. My motivation is to make masked array overhead negligible, at least for medium to large arrays.
You can use inline in that case - starting with numpy 1.3.0, inline can be used in C code in a portable way (it is a macro which points to compiler specific inline if the C99 inline is not available so it works even on broken compilers).
cheers,
David

On Sat, May 16, 2009 at 2:02 AM, Eric Firing efiring@hawaii.edu wrote:
Charles R Harris wrote:
On Fri, May 15, 2009 at 7:48 PM, Eric Firing <efiring@hawaii.edu mailto:efiring@hawaii.edu> wrote:
http://www.mail-archive.com/numpy-discussion@scipy.org/msg17595.html Prompted by the thread above, I decided to see what it would take to implement ufuncs with masking in C. I described the result here: http://www.mail-archive.com/numpy-discussion@scipy.org/msg17698.html Now I am starting a new thread. The present state of the work is now
in
github: http://github.com/efiring/numpy-work/tree/cfastma I don't want to do any more until I have gotten some feedback from
core
developers. (And I would be delighted if someone wants to help with this, or take it over.)
Chuck,
Thanks very much for the quick response.
Here the if ... continue needs to follow the declaration:
if (*mp1) continue; float in1 = *(float *)ip1; float in2 = *(float *)ip2; *(float *)op1 = f(in1, in2);
I was surprised to see the declarations inside the loop in the first place (this certainly is not ANSI-C), and I was also pleasantly surprised that letting them be after the conditional didn't seem to bother the compiler at all. Maybe that is a gcc extension.
Declarations at the top of a block have always been valid C.
I think this would be better as
if (!(*mp1)) { float in1 = *(float *)ip1; float in2 = *(float *)ip2; *(float *)op1 = f(in1, in2); }
I agree, and I thought of that originally--I think I did it with continue because it was easier to type it in, and it reduced the difference relative to the non-masked form.
But since this is actually a ternary function, you could define new functions, something like
double npy_add_m(double a, double b, double mask) { if (!mask) { return a + b; else { return a; } }
And use the currently existing loops. Well, you would have to add one for ternary functions.
That would incur the overhead of an extra function call for each element; I suspect it would slow it down a lot. My motivation is to make masked array overhead negligible, at least for medium to large arrays.
It overhead would be the same as it is now, the generic loops all use passed function pointers for functions like sin. Some functions like addition, which is intrinsic and not part of a library, are done in their own special loops that you will find further down in that file. The difficulty I see is that with the current machinery the mask will be converted to the same type as the added numbers and that could add some overhead.
Also your suggestion above does not handle the case where an output argument is supplied; it would modify the output under the mask.
Question, what about reduce? I don't think it is defined defined for ternary functions. Apart from reduce, why not just add, you already have the mask to tell you which results are invalid.
You mean just do the operation and ignore the results under the mask? This is the way Pierre originally did it, if I remember correctly, but fairly recently people started objecting that they didn't want to disturb values in an output argument under a mask. So now ma jumps through hoops to satisfy this requirement, and it is consequently slow.
OK. I'm not familiar with the uses of masked arrays.
ufunc methods like reduce are supported only for the binary ops with one output, so they are automatically unavailable for the masked versions. To get around this would require subclassing the ufunc to make a masked version. This is probably the best way to go, but I suspect it is much more complicated than I can handle in the amount of time I can spend.
I think reduce could be added for ternary functions, but it is a design decision how it should operate.
Chuck
participants (5)
-
Charles R Harris
-
David Cournapeau
-
Eric Firing
-
Matthieu Brucher
-
Robert Kern