Mailman 3 July 2011 - NumPy-Discussion

Are .M and .H removed in NumPy 1.6?
by Sturla Molden 12 Jul '11

12 Jul '11

After upgrading EPD, I just discovered that my ndarrays no longer have .M and .H attributes. Were they deprectated, or is my NumPy not working correctly? Sturla

2 2

Missing Values Discussion
by Travis Oliphant 11 Jul '11

11 Jul '11

Hi all, I want to first apologize for stepping into this discussion a bit late and for not being able to participate adequately. However, I want to offer a couple of perspectives, and my opinion about what we should do as well as clarify what I have instructed Mark to do as part of his summer work. First, the discussion has reminded me how valuable it is to get feedback from all points of view. While it does lengthen the process, it significantly enhances the result. I strongly hope we can continue the tradition of respectful discussion on this mailing list where people's views are treated with respect --- even if we don't always have the time to understand them in depth. I also really appreciate people taking the time to visit on the phone call with me as it gave me a chance to understand many opinions quickly and at least start to form a possibly useful opinion. Basically, because there is not consensus and in fact a strong and reasonable opposition to specific points, Mark's NEP as proposed cannot be accepted in its entirety right now. However, I believe an implementation of his NEP is useful and will be instructive in resolving the issues and so I have instructed him to spend Enthought time on the implementation. Any changes that need to be made to the API before it is accepted into a released form of NumPy can still be made even after most of the implementation is completed as far as I understand it. This is because most of the disagreement is about the specific ability to manipulate the masks independently of assigning missing data and the creation of an additional np.HIDE (np.IGNORE) concept at the Python level. Despite some powerful arguments on both sides of the discussion, I am confident that we can figure out an elegant solution that will work long term. My current opinion is that I am very favorable to making easy the use-case that has been repeatedly described of having "missing data" that is *always* missing and then having "hidden data" that you don't want to think about for a particular set of calculations (but you also don't want to through away by over-writing). I think it is important to make it easy to keep that data around without over-writing but also have the "idea" of that kind of missing data different than the idea of data you can't care about because it just isn't there. I also think it is important for the calculation infrastructure to have just one notion of "missing data" which Mark's NEP handles beautifully. It seems to me that some of the disagreement is one of perspective in that Mark articulates very well the position of "generic programming, make-opaque-the-implementation" perspective with a focus on the implications of missing data for calculations. Nathaniel and Matthew articulate well the perspective of "focusing" on the data object itself and the desire to keep separate the different ideas behind missing data that have been described --- as well as a powerfully described description of the NumPy tradition of exposing the raw data to the Python side without hiding too much of the implementation from the user. I think it's a healthy discussion. But, I would like to see Mark's code get completed so that we can start talking about code examples. Please don't interpret my instructing Mark to finish the code as "it's been decided". I simply think it's the best path forward to ultimately resolving the concerns. I would like to see an API worked out before summer's end --- and I'm hopeful everyone will be excited about what the resulting design is. I do think there is room for agreement in the present debate if we all remember to keep listening to each other. It takes a lot of effort to understand somebody else's point of view. I have been grateful to see evidence I see of that behavior multiple times (in Mark's revamping of the NEP, in Matthew Brett's re-statement of his interpretation of Mark's views, in Nathaniel's working hard to engage the dialogue even in the throes of finishing his PhD, and many other examples). It makes me very happy to be a part of this community. I look forward to times when I can send more thoughtful and technical emails than this one. All the best, -Travis

3 11

Problem with ufunc of a numpy.ndarray derived class
by Martin DRUON 11 Jul '11

11 Jul '11

Hi, I have a problem with the ufunc return type of a numpy.ndarray derived class. In fact, I subclass a numpy.ndarray using the tutorial : http://docs.scipy.org/doc/numpy/user/basics.subclassing.html But, for example, if I execute the "max" ufunc from my subclass, the return type differs from the return type of the numpy ufunc. This is my code (testSubclassNumpy.py), "copy/paste" from the tutorial : # -*- coding: utf-8 -*- import numpy class MySubClass(numpy.ndarray): def __new__(cls, input_array, info=None): obj = numpy.asarray(input_array).view(cls) obj.info = info return obj def __array_finalize__(self, obj): #print 'In __array_finalize__:' #print ' self is %s' % repr(self) #print ' obj is %s' % repr(obj) if obj is None: return self.info = getattr(obj, 'info', None) def __array_wrap__(self, out_arr, context=None): #print 'In __array_wrap__:' #print ' self is %s' % repr(self) #print ' arr is %s' % repr(out_arr) # then just call the parent return numpy.ndarray.__array_wrap__(self, out_arr, context) >>> import numpy >>> numpy.__version__ '1.6.0' >>> import testSubclassNumpy >>> a = numpy.random.random(size=(10,10)) >>> t = testSubclassNumpy.MySubClass(a) >>> type(a) >>> type(t) >>> a.max() 0.99207693069079683 >>> t.max() MySubClass(0.9920769306907968) >>> type(numpy.max(a)) >>> type(numpy.max(t)) This problem seems to be appeared with the latest version of numpy. Today, I use Python 2.7.2 + numpy 1.6.0 but I didn't have this problem with python 2.6.6 and numpy 1.5.1. Is it a bug ? or perhaps I have made mistake somewhere... thanks, Martin

1 0

FloatingPointError: overflow encountered in multiply
by Roy Lowrance 11 Jul '11

11 Jul '11

I have a 1D float64 array ts. I want to square each element, so I compute x = ts * ts I get a floating point overflow error. However, when I access each element separately and multiple, I get no error: for i in ts.shape[0]:

2 2

Post-doctoral position in image & data analysis for image-based high-content screening
by Thouis (Ray) Jones 11 Jul '11

11 Jul '11

Short summary: A funded post-doctoral position in data analysis, machine-learning, and statistics applied to biological image-based high-content screening is available at the BioPhenics platform of the Curie Institute (Paris, France). The position involves development and maintenance of tools in Python, Matlab, or R that will be used both to analyze current screens and as a contribution to a larger community. Expertise in machine learning or statistical analysis with large biological datasets is desirable, good knowledge of image processing tools is an asset. Full details at: http://dl.dropbox.com/u/16028921/PostDocApp_MachineLearning_BFX.pdf Thouis Jones Institut Curie

1 0

Re: [Numpy-discussion] gist gist: 1068264
by Nathaniel Smith 11 Jul '11

11 Jul '11

Hi Bruce, I'm replying on the list instead of on github, to make it easier for others to join in the discussion if they want. [For those joining in: this was a comment posted at https://gist.github.com/1068264 ] On Fri, Jul 8, 2011 at 10:36 AM, bsouthey wrote: > I presume missing float values could be addressed with one of the 'special' ranges such as 'Indeterminate' in IEEE 754 (http://babbage.cs.qc.edu/IEEE-754/References.xhtml). The outcome should be determined by the IEEE special operations. Right. An IEEE 754 double has IIRC about 2^53 distinct bit-patterns that all mean "not a number". A few of these are used to signal different invalid operations: In [20]: hex(np.asarray([np.nan]).view(dtype=np.uint64)[0]) Out[20]: '0x7ff8000000000000L' In [21]: hex(np.log([0]).view(dtype=np.uint64)[0]) Out[21]: '0xfff0000000000000L' In [22]: hex(np.divide([0.], [0,]).view(dtype=np.uint64)[0]) Out[22]: '0xfff8000000000000L' ...but that only accounts for, like, 10 of the 2^53 or something. The rest are simply unused. So what R does, and what we would do for dtype-style NAs, is just pick one of those (ideally the same one R uses), and declare that that is *not* not a number; it's NA. > So my real concern is handling integer arrays: > 1) How will you find where the missing values are in an array? If there is a variable that denotes missing values are present (NA_flags?) then do you have to duplicate code to avoid this searching when an array has no missing values? Each dtype has a bunch of C functions associated with it that say how to do comparisons, assignment, etc. In the miniNEP design, we add a new function to this list called 'isna', which every dtype that wants to support NAs has to define. Yes, this does mean that code which wants to treat NAs separately has to check for and call this function if it's present, but that seems to be inevitable... *all* of the dtype C functions are supposedly optional, so we have to check for them before calling them and do something sensible if they aren't defined. We could define a wrapper that calls the function if its defined, or else just fills the provided buffer with zeros (to mean "there are no NAs), and then code which wanted to avoid a special case could use that. But in general we probably do want to handle arrays that might have NAs differently from arrays which don't have NAs, because if there are no NAs present then it's quicker to skip the handling altogether. That's true for any NA implementation. > 2) What happens if a normal operation equates to that value: If you use max(np.int8), such as when adding 1 to an array with an element of 126 or when overflow occurs: >>>> np.arange(120,127, dtype=np.int8)+2 > array([ 122, 123, 124, 125, 126, 127, -128], dtype=int8) > The -128 corresponds to the missing element but is the second to last element now missing? This is worse if the overflow is larger. Yeah, in the design as written, overflow (among other things) can create accidental NAs. Which kind of sucks. There are a few options: -- Just live with it. -- We could add a flag like NPY_NA_AUTO_CHECK, and when this flag is set, the ufunc loop runs 'isna' on its output buffer before returning. If there are any NAs there that did not arise from NAs in the input, then it raises an error. (The reason we would want to make it a flag is that this checking is pointless for dtypes like NA-string, and mostly pointless for dtypes like NA-float.) Also, we'd only want to enable this if we were using the NPY_NA_AUTO_UFUNC ufunc-delegation logic, because if you registered a special ufunc loop *specifically for your NA-dtype*, then presumably it knows what it's doing. This would also allow such an NA-dtype-specific ufunc loop to return NAs on purpose if it wanted to. -- Use a dtype that adds a separate flag next to the actual integer to indicate NA-ness, instead of stealing one of the integer's values. So your NA-int8 would actually be 2 bytes, where the first byte was 1 to indicate NA, or 0 to indicate that the second byte contains an actual int8. If you do this with larger integers, say an int32, then you have a choice: you could store your int32 in 8 bytes, in which case arithmetic etc. is fast, but you waste a bit of memory. Or you could store your int32 in 5 bytes, in which case arithmetic etc. become somewhat slower, but you don't waste any memory. (This latter case would basically be like using an unaligned or byteswapped array in current numpy, in terms of mechanisms and speed.) -- Nothing in this design rules out a second implementation of NAs based on masking. Personally, as you know, I'm not a big fan, but if it were added anyway, then you could use that for your integers as well. A related issue is, of the many ways we *can* do integer NA-dtype, which one *should* we do by default. I don't have a strong opinion, really; I haven't heard anyone say that they have huge quantities of integer-plus-NA data that they want to manipulate and memory/speed/allowing the full range of values are all really important. (Maybe that's you?) In the design as written, they're all pretty trivial to implement (you just tweak a few magic numbers in the dtype structure), and probably we should support all of them via more-or-less exotic invocations of np.withNA. (E.g., 'np.withNA(np.int32, useflag=True, flagsize=1)' to get a 5-byte int32.) ...I kind of like that NPY_NA_AUTO_CHECK idea, it's pretty clean and would definitely make things safer. I think I'll add it. -- Nathaniel

2 2

code review request: masked iteration, fix 'where=' in ufuncs, documentation
by Mark Wiebe 10 Jul '11

10 Jul '11

>From the pull request: https://github.com/numpy/numpy/pull/108 The two flags NPY_ITER_WRITEMASKED and NPY_ITER_ARRAYMASK now fully work. This made it easy to fix the ufunc 'where=' bug. Also added documentation of the new iterator flags and inline functions for dealing with masks. Special thanks to Chuck, Eric, and Nathaniel for reviewing and providing feedback on my last pull request. -Mark

1 0

code review request: masked dtype transfers
by Mark Wiebe 10 Jul '11

10 Jul '11

I've just made pull request 105: https://github.com/numpy/numpy/pull/105 This adds public API PyArray_MaskedCopyInto and PyArray_MaskedMoveInto, which behave analogously to the corresponding unmasked functions. To expose this with a reasonable interface, I added a function np.copyto, which takes a 'where=' parameter just like the element-wise ufuncs. One thing which needs discussion is that I've flagged 'putmask' and PyArray_PutMask as deprecated, because 'copyto' PyArray_MaskedMoveInto handle what those functions do but in a more flexible fashion. If there are any objections to deprecating 'putmask' and PyArray_PutMask, please speak up! Thanks, Mark

4 14

non-integer powers in scimath.power
by Nathaniel Dosé 09 Jul '11

09 Jul '11

Is it true that numpy.lib.scimath.power() only accepts integer powers? >>> np.info(scimath.power) ... >>> Parameters >>> ---------- >>> x : array_like >>> The input value(s). >>> p : array_like of ints In any case, scimath.power() doesn't complain when you send in a non-int, and in fact gives a pretty correct answer: >>> scimath.power(-4., .5) (1.2246063538223773e-16+2j) Note I'm wondering about any non-integers, not just .5 Also note that the "normal" numpy power() function accepts non-ints: >>> np.info(np.power) ... >>> Parameters >>> ---------- >>> x1 : array_like >>> The bases. >>> x2 : array_like >>> The exponents. ... >>> np.power(4, .5) 2.0 (Sorry, I tried to follow the source for scimath.py to see what really happens at the C level, but lose the trail as soon as I hit numpy.core.umath...) -Nathaniel

1 0

potential bug in PyArray_MoveInto and PyArray_CopyInto?
by James Bergstra 08 Jul '11

08 Jul '11

In numpy 1.5.1, the functions PyArray_MoveInto and PyArray_CopyInto don't appear to treat strides correctly. Evidence: PyNumber_InPlaceAdd(dst, src), and modifies the correct subarray to which dst points. In the same context, PyArray_MoveInto(dst, src) modifies the first two rows of the underlying matrix instead of the first two columns. PyArray_CopyInto does the same. Is there something subtle going on here? James -- http://www-etud.iro.umontreal.ca/~bergstrj

2 3