Mailman 3 December 2011 - NumPy-Discussion

Fast Reading of ASCII files
by Chris.Barker Dec. 13, 2011

Dec. 13, 2011

Hi folks, This is a continuation of a conversation already started, but i gave it a new, more appropriate, thread and subject. On 12/6/11 2:13 PM, Wes McKinney wrote: > we should start talking > about building a *high performance* flat file loading solution with > good column type inference and sensible defaults, etc. ... > I personally don't > believe in sacrificing an order of magnitude of performance in the 90% > case for the 10% case-- so maybe it makes sense to have two functions > around: a superfast custom CSV reader for well-behaved data, and a > slower, but highly flexible, function like loadtable to fall back on. I've wanted this for ages, and have done some work towards it, but like others, only had the time for a my-use-case-specific solution. A few thoughts: * If we have a good, fast ascii (or unicode?) to array reader, hopefully it could be leveraged for use in the more complex cases. So that rather than genfromtxt() being written from scratch, it would be a wrapper around the lower-level reader. * key to performance is to have the text to number to numpy type happening in C -- if you read the text with python, then convert to numbers, then to numpy arrays, it's simple going to be slow. * I think we want a solution that can be adapted to arbitrary text files -- not just tabular, CSV-style data. I have a lot of those to read - and some thoughts about how. Efforts I have made so far, and what I've learned from them: 1) fromfile(): fromfile (for text) is nice and fast, but buggy, and a bit too limited. I've posted various notes about this in the past (and, I'm pretty sure a couple tickets). They key missing features are: a) no support form commented lines (this is a lessor need, I think) b) there can be only one delimiter, and newlines are treated as generic whitespace. What this means is that if you have whitespace-delimited file, you can read multiple lines, but if it is, for instance, comma-delimited, then you can only read one line at a time, killing performance. c) there are various bugs if the text is malformed, or doesn't quite match what you're asking for (ie.e reading integers, but the tet is float) -- mostly really limited error checking. I spent some time digging into the code, and found it to be really hard to track C code. And very hard to update. The core idea is pretty nice -- each dtype should know how to read itself form a text file, but the implementation is painful. The key issue is that for floats and ints, anyway, it relies on the C atoi and atof functions. However, there have been patches to these that handle NaN better, etc, for numpy, and I think a python patch as well. So the code calls the numpy atoi, which does some checks, then calls the python atoi, which then calls the C lib atoi (I think all that...) In any case, the core bugs are due to the fact that atoi and friends doesn't return an error code, so you have to check if the pointer has been incremented to see if the read was successful -- this error checking is not propagated through all those levels of calls. It got really ugly to try to fix! Also, the use of the C atoi() means that locales may only be handled in the default way -- i.e. no way to read european-style floats on a system with a US locale. My conclusion -- the current code is too much a mess to try to deal with and fix! I also think it's a mistake to have text file reading a special case of fromfile(), it really should be a separate issue, though that's a minor API question. 2) FileScanner: FileScanner is some code a wrote years ago as a C extension - it's limited, but does the job and is pretty fast. It essentially calls fscanf() as many times as it gets a successful scan, skipping all invalid text, then returning a numpy array. You can also specify how many numbers you want read from the file. It only supports floats. Travis O. asked it it could be included in Scipy way back when, but I suspect none of my code actually made it in. If I had to do it again, I might write something similar in Cython, though I am still using it. My Conclusions: I think what we need is something similar to MATLAB's fscanf(): what it does is take a C-style format string, and apply it to your file over an over again as many times as it can, and returns an array. What's nice about this is that it can be purposed to efficiently read a wide variety of text files fast. For numpy, I imagine something like: fromtextfile(f, dtype=np.float64, comment=None, shape=None): """ read data from a text file, returning a numpy array f: is a filename or file-like object comment: is a string of the comment signifier. Anything on a line after this string will be ignored. dytpe: is a numpy dtype that you want read from the file shape: is the shape of the resulting array. If shape==None, the file will be read until EOF or until there is read error. By default, if there are newlines in the file, a 2-d array will be returned, with the newline signifying a new row in the array. """ This is actually pretty straightforward. If it support compound dtypes, then you can read a pretty complex CSV file, once you've determined the dtype for your "record" (row). It is also really simple to use for the simple cases. But of course, the implementation could be a pain -- I've been thinking that you could get a lot of it by creating a mapping from numpy dtypes to fscanf() format strings, then simply use fscanf for the actual file reading. This would certainly be easy for the easy cases. (maybe you'd want to use sscanf, so you could have the same code scan strings as well as files) Ideally, each dtype would know how to read itself from a string, but as I said above, the code for that is currently pretty ugly, so it may be easier to keep it separate. Anyway, I'd be glad to help with this effort. -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker(a)noaa.gov

8 9

Re: [Numpy-discussion] NumPy-Discussion Digest, Vol 63, Issue 43
by Konrad Banachewicz Dec. 13, 2011

Dec. 13, 2011

U On 12/13/11, numpy-discussion-request(a)scipy.org <numpy-discussion-request(a)scipy.org> wrote: > Send NumPy-Discussion mailing list submissions to > numpy-discussion(a)scipy.org > > To subscribe or unsubscribe via the World Wide Web, visit > http://mail.scipy.org/mailman/listinfo/numpy-discussion > or, via email, send a message with subject or body 'help' to > numpy-discussion-request(a)scipy.org > > You can reach the person managing the list at > numpy-discussion-owner(a)scipy.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of NumPy-Discussion digest..." > > > Today's Topics: > > 1. Re: Fast Reading of ASCII files (Chris Barker) > 2. Re: Apparently non-deterministic behaviour of complex array > multiplication (kneil) > 3. Re: numpy.mean problems (Eraldo Pomponi) > 4. Re: Fast Reading of ASCII files (Bruce Southey) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Tue, 13 Dec 2011 10:08:44 -0800 > From: Chris Barker <chris.barker(a)noaa.gov> > Subject: Re: [Numpy-discussion] Fast Reading of ASCII files > To: denis <denis-bz-gg(a)t-online.de>, Discussion of Numerical Python > <numpy-discussion(a)scipy.org> > Message-ID: > <CALGmxEJt9Y0oaM1gkFSuFLwaBJNxfLk54x-N8+f8hT5VzjcVtQ(a)mail.gmail.com> > Content-Type: text/plain; charset="utf-8" > > NOTE: > > Let's keep this on the list. > > On Tue, Dec 13, 2011 at 9:19 AM, denis <denis-bz-gg(a)t-online.de> wrote: > >> Chris, >> unified, consistent save / load is a nice goal >> >> 1) header lines with date, pwd etc.: "where'd this come from ?" >> >> # (5, 5) svm.py bz/py/ml/svm 2011-12-13 Dec 11:56 -- automatic >> # 80.6 % correct -- user info >> 245 39 4 5 26 >> ... >> > I'm not sure I understand what you are expecting here: What would be > automatic? if itparses a datetime on the header, what would it do with it? > But anyway, this seems to me: > - very application specific -- this is for the users code to write > - not what we are talking about at this point anyway -- I think this > discussion is about a lower-level, does-the-simple-things-fast reader -- > that may or may not be able to form the basis of a higher-level fuller > featured reader. > > >> 2) read any CSVs: comma or blank-delimited, with/without column names, >> a la loadcsv() below >> > > yup -- though the column name reading would be part of a higher-level > reader as far as I'm concerned. > > >> 3) sparse or masked arrays ? >> >> sparse probably not, that seem pretty domain dependent to me -- though > hopefully one could build such a thing on top of the lower level reader. > Masked support would be good -- once we're convinced what the future of > masked arrays are in numpy. I was thinking that the masked array issue > would really be a higher-level feature -- it certainly could be if you need > to mask "special value" stype files (i.e. 9999), but we may have to build > it into the lower level reader for cases where the mask is specified by > non-numerical values -- i.e. there are some met files that use "MM" or some > other text, so you can't put it into a numerical array first. > >> >> Longterm wishes: beyond the scope of one file <-> one array >> but essential for larger projects: >> 1) dicts / dotdicts: >> Dotdict( A=anysizearray, N=scalar ... ) <-> a directory of little >> files >> is easy, better than np.savez >> (Haven't used hdf5, I believe Matlabv7 does.) >> >> 2) workflows: has anyone there used visTrails ? >> > > outside of the spec of this thread... > >> >> Anyway it seems to me (old grey cynic) that Numpy/scipy developers >> prefer to code first, spec and doc later. Too pessimistic ? >> >> > Well, I think many of us believe in a more agile style approach -- > incremental development. But really, as an open source project, it's really > about scratching an itch -- so there is usually a spec in mind for the itch > at hand. In this case, however, that has been a weakness -- clearly a > number of us hav written small solutions to our particular problem at hand, > but no we haven't arrived at a more general purpose solution yet. So a bit > of spec-ing ahead of time may be called for. > > On that: > > I"ve been thinking from teh botom-up -- imaging what I need for the simple > case, and how it might apply to more complex cases -- but maybe we should > think about this another way: > > What we're talking about here is really about core software engineering -- > optimization. It's easy to write a pure-python simple file parser, and > reasonable to write a complex one (genfromtxt) -- the issue is performance > -- we need some more C (or Cython) code to really speed it up, but none of > us wants to write the complex case code in C. So: > > genfromtxt is really nice for many of the complex cases. So perhaps > another approach is to look at genfromtxt, and see what > high performance lower-level functionality we could develop that could make > it fast -- then we are done. > > This actually mirrors exactly what we all usually recommend for python > development in general -- write it in Python, then, if it's really not fast > enough, write the bottle-neck in C. > > So where are the bottle necks in genfromtxt? Are there self-contained > portions that could be re-written in C/Cython? > > -Chris > > > > > > > -- > > Christopher Barker, Ph.D. > Oceanographer > > Emergency Response Division > NOAA/NOS/OR&R (206) 526-6959 voice > 7600 Sand Point Way NE (206) 526-6329 fax > Seattle, WA 98115 (206) 526-6317 main reception > > Chris.Barker(a)noaa.gov >

1 0

Apparently non-deterministic behaviour of complex array multiplication
by Karl Kappler Dec. 13, 2011

Dec. 13, 2011

Hello, I am somewhat new to scipy/numpy so please point me in the right direction if I am posting to an incorrect forum. The experience which has prompted my post is the following: I have a numpy array Y where the elements of Y are type(Y[0,0]) Out[709]: <type 'numpy.complex128'> The absolute values of the real and complex values do not far exceed say 1e-10. The shape of Y is (24, 49218). When I perform the operation: C = dot(Y,Y.conj().transpose), i.e. I form the covariance matrix by multiplying T by its conjugate transpose, I sometimes get NaN in the array C. I can imagine some reasons why this may happen, but what is truly puzzling to me is that I will be working in ipython and will execute for example: find(isnan(C)) and will be returned an list of elements of C which are NaN, fine, but then I recalculate C, and repeat the find(isnan(C)) command and I get a different answer. I type: find(isnan(dot(Y,Y.conj().transpose))) and an empty array is returned. Repeated calls of the same command however result in a non-empty array. In fact, the sequence of arrays returned from several consecutive calls varies. Sometimes there are tens of NaNs, sometimes none. I have been performing a collection of experiments for some hours and cannot get to the bottom of this; Some things I have tried: 1. Cast the array Y as a matrix X and calculate X*X.H --- in this case i get the same result in that sometimes I have NaN and sometimes I do not. 2. set A=X.H and calculate X*A --- same results* 3. set b=A.copy() and calc X*b --- same results*. 4. find(isnan(real(X*X.H))) --- same results* 5. find(isnan(real(X)*real(X.H))) - no NaN appear *N.B. "Same results" does not mean that the same indices were going NaN, simply that I was getting back a different result if I ran the command say a dozen times. So it would seem that it has something to do with the complex multiplication. I am wondering if there is too much dynamic range being used in the calculation? It absolutely amazes me that I can perform the same complex-arithmetic operation sitting at the command line and obtain different results each time. In one case I ran a for loop where I performed the multiplication 1000 times and found that 694 trials had no NaN and 306 trials had NaN. Saving X to file and then reloading it in a new ipython interpreter typically resulted in no NaN. For a fixed interpretter and instance of X or Y, the indices which go NaN (when they do) sometimes repeat many times and sometimes they vary apparently at random. Also note that I have had a similar problem with much smaller arrays, say 24 x 3076 I have also tried 'upping' the numpy array to complex256, I have like 12GB of RAM... This happens both in ipython and when I call my function from the command line. Does this sound familiar to anyone? Is my machine possessed?

6 12

NumPy book review copies
by Jitesh Y. Gawali Dec. 13, 2011

Dec. 13, 2011

Hi, As a part of our reviewing program, we are giving away limited number of copies (print & electronic) of our recent publication NumPy 1.5 Beginner's Guide <http://www.packtpub.com/numpy-1-5-using-real-world-examples-beginners-guide…> to people interested in reviewing the book. You need to publish your review/feedback on either your blog or on websites like Barnes and Noble <http://www.barnesandnoble.com/>, Slashdot <http://slashdot.org/> and so on as per your choice. We also encourage uploading the review on Amazon since it gives buyers a chance to know the book through your perspective. Please get in touch with me for more information. Thanks, Jitesh <http://www.barnesandnoble.com/> -- Jitesh Gawali *Marketing Research Executive*| Packt Publishing | www.PacktPub.com <http://www.packtpub.com/> *MSN:*_jiteshg_ <mailto:jiteshg@packtpub.com>_(a)packtpub.com_ <mailto:jiteshg@packtpub.com> Interested in becoming an author? Visit Packt's Author Website <http://authors.packtpub.com/>for all the information you need about writing for Packt.

1 0

ANN: PyCPX 0.02, a numpy/cython wrapper for CPlex
by Hoyt Koepke Dec. 13, 2011

Dec. 13, 2011

Hello, I'm pleased to announce the second release of the PyCPX wrapper for IBM's CPlex Optimizer Suite. This second release fixes several bugs and vastly optimizes the model construction stage, particularly for large models. PyCPX is a python wrapper for the CPlex optimization suite that focuses on ease of use and seamless integration with numpy. It allows one to naturally specify linear and quadratic problems over real, boolean, and integer variables. For documentation and examples, please see the website: http://www.stat.washington.edu/~hoytak/code/pycpx/index.html. Please send me any comments, questions and suggestions! I am quite open to feedback. Thanks, --Hoyt ++++++++++++++++++++++++++++++++++++++++++++++++ + Hoyt Koepke + University of Washington Department of Statistics + http://www.stat.washington.edu/~hoytak/ + hoytak(a)gmail.com ++++++++++++++++++++++++++++++++++++++++++

1 0

polynomial with negative exponents
by LASAGNA DAVIDE Dec. 12, 2011

Dec. 12, 2011

Hi, I have written a class for polynomials with negative exponents like: p(x) = a0 + a1*x**-1 + ... + an*x**-n The code is this one: class NegativeExpPolynomial( object ): def __init__ ( self, coeffs ): self.coeffs = np.array( coeffs ) def __call__( self, x ): return sum( (c*x**(-i) for i, c in enumerate( self.coeffs ) ) ) where coeffs = [a0, a1, ..., an]. I find that the way i evaluate the polynomial is kind of *slow*, especially for polynomial with order larger than ~200 and for arrays x large enough. Do you have suggestions on how to speed up this code? Regards, Davide Lasagna -- Phd Student Dipartimento di Ingegneria Aeronautica a Spaziale Politecnico di Torino, Italy tel: 011/0906871 e-mail: davide.lasagna(a)polito.it; lasagnadavide(a)gmail.com

4 4

Warning related to __multiarray_api.h
by Mads Ipsen Dec. 9, 2011

Dec. 9, 2011

Hi, I don't know if this is of importance, but when I compile code using the numpy C API, I get the warning: site-packages/numpy/core/include/numpy/__multiarray_api.h:1532: warning: 'int _import_array()' defined but not used Might be worth cleaning it up. Best regards, Mads

1 0

Problem with using PyArray_AsCArray
by Ng, Enrico Dec. 9, 2011

Dec. 9, 2011

I am trying to pass a multi-dimensional ndarray to C as a multi-dimensional C array for the purposes of passing it to mathematica. I am using PyArray_AsCArray but getting an error. ###################################################### Python Code: import Image from scipy.misc import fromimage img = Image.open("../APLS/image709_enhanced_2.tif") nimg = fromimage(img) ... mathlink_PutIntegerArray(nimg) ###################################################################### C Wrapper Code: static PyObject * mathlink_PutIntegerArray(mathlink_Link *self, PyObject *args) { npy_intp dims[3]; /* PyArray_AsCArray is for ndim <= 3 */ PyObject *o1; double *d1; long result; i=PyArray_AsCArray(&o1, (void *)&d1, dims, PyArray_NDIM(o1), PyArray_DescrFromType(PyArray_TYPE(o1))); ... } Program received signal SIGSEGV, Segmentation fault. 0x00002aaaaefa3d2e in PyArray_AsCArray (op=0x7fffffffdf68, ptr=0x7fffffffdf60, dims=0x7fffffffdf40, nd=3, typedescr=<value optimized out>) at numpy/core/src/multiarray/multiarraymodule.c:218 218 ptr3[i][j] = ap->data + i*ap->strides[0] + j*ap->strides[1] ################################################################################# I am able to read i,j,ap->data, ap->strides[0], ap->strides[1] so the error seems to be in the assignment to ptr3[i][j] This happens on the first instance i=0 j=0 ################################################################################# PyArray_AsCArray code is below. NPY_NO_EXPORT int PyArray_AsCArray(PyObject **op, void *ptr, npy_intp *dims, int nd, PyArray_Descr* typedescr) { PyArrayObject *ap; npy_intp n, m, i, j; char **ptr2; char ***ptr3; if ((nd < 1) || (nd > 3)) { PyErr_SetString(PyExc_ValueError, "C arrays of only 1-3 dimensions available"); Py_XDECREF(typedescr); return -1; } if ((ap = (PyArrayObject*)PyArray_FromAny(*op, typedescr, nd, nd, CARRAY, NULL)) == NULL) { return -1; } switch(nd) { case 1: *((char **)ptr) = ap->data; break; case 2: n = ap->dimensions[0]; ptr2 = (char **)_pya_malloc(n * sizeof(char *)); if (!ptr2) { goto fail; } for (i = 0; i < n; i++) { ptr2[i] = ap->data + i*ap->strides[0]; } *((char ***)ptr) = ptr2; break; case 3: n = ap->dimensions[0]; m = ap->dimensions[1]; ptr3 = (char ***)_pya_malloc(n*(m+1) * sizeof(char *)); if (!ptr3) { goto fail; } for (i = 0; i < n; i++) { ptr3[i] = ptr3[n + (m-1)*i]; for (j = 0; j < m; j++) { ptr3[i][j] = ap->data + i*ap->strides[0] + j*ap->strides[1]; } } *((char ****)ptr) = ptr3; } memcpy(dims, ap->dimensions, nd*sizeof(npy_intp)); *op = (PyObject *)ap; return 0; fail: PyErr_SetString(PyExc_MemoryError, "no memory"); return -1; }

2 2

Help with wrapper ndarray to c array
by Enrico Ng Dec. 9, 2011

Dec. 9, 2011

I am trying to pass a multi-dimensional ndarray to C as a multi- dimensional C array for the purposes of passing it to mathematica. They already have a wrapper for a 1-D Python list. where the list is copied to "list". Shown below: static PyObject * mathlink_PutIntegerList(mathlink_Link *self, PyObject *args) { PyObject* seq; PyObject* obj; long i, len, result; int* list; len = PyObject_Length(seq); list = PyMem_New(int, len); for(i = 0; i < len; i++) { obj = PySequence_GetItem(seq, i); list[i] = PyInt_AsLong(obj); } CheckForThreadsAndRunLink(self,result = MLPutIntegerList(self->lp, list, len)); PyMem_Free(list); CHECKNOTEQUAL(result,MLSUCCESS,self); Py_INCREF(Py_None); return Py_None; } I would like to create a similar wrapper which accepts an ndarray and provides the array laid out in memory like a C array declared explicitly as "int a[m][n]...". I also need to pass the length of the array at each level i as dim[i] and the depth. Since this is pretty much the only function I plan to wrap, I'd like to avoid using boost, swig, etc. Any help would be appreciated, looks like I should use PyArray_AsCArray but I'm having trouble finding examples that work for me.

1 0

type checking, what's recommended?
by josef.pktd＠gmail.com Dec. 8, 2011

Dec. 8, 2011

If I want to know whether something that might be an array is really a plain ndarray and not a subclass, is using `type` the safest bet? All the other forms don't discriminate against subclasses. >>> type(np.ma.zeros(3)) is np.ndarray False >>> type(np.zeros(3)) is np.ndarray True >>> isinstance(np.ma.zeros(3), np.ndarray) True >>> isinstance(np.zeros(3), np.ndarray) True >>> issubclass(np.ma.zeros(3).__class__, np.ndarray) True >>> issubclass(np.zeros(3).__class__, np.ndarray) True >>> isinstance(np.matrix(np.zeros(3)), np.ndarray) True >>> type(np.matrix(np.zeros(3))) is np.ndarray False Thanks, Josef

2 1