Mailman 3 May 2007 - NumPy-Discussion

PyArray_ContiguousFromObject & PyDECREF
by Michiel Jan Laurens de Hoon May 17, 2011

May 17, 2011

Dear pythoneers, In Chapter 12 ("Writing a C extension") in the Numerical Python manual, in the first example PyArray_ContiguousFromObject is used together with PyDECREF, while in the last example PyArray_ContiguousFromObject is used but the arrays are not PyDECREF'ed before the function returns. From Chapter 13 ("C API Reference") it seems that a PyDECREF is needed in both cases. Or am I missing something here? Is PyDECREF not needed for some reason in the last example? Thanks in advance, Michiel de Hoon U Tokyo First example: ============== static PyObject * trace(PyObject *self, PyObject *args) { PyObject *input; PyArrayObject *array; double sum; int i, n; if (!PyArg_ParseTuple(args, "O", &input)) return NULL; array = (PyArrayObject *) PyArray_ContiguousFromObject(input, PyArray_DOUBLE, 2, 2); if (array == NULL) return NULL; n = array->dimensions[0]; if (n > array->dimensions[1]) n = array->dimensions[1]; sum = 0.; for (i = 0; i < n; i++) sum += *(double *)(array->data + i*array->strides[0] + i*array->strides[1]); Py_DECREF(array); return PyFloat_FromDouble(sum); } Last example: ============= static PyObject * matrix_vector(PyObject *self, PyObject *args) { PyObject *input1, *input2; PyArrayObject *matrix, *vector, *result; int dimensions[1]; double factor[1]; double real_zero[1] = {0.}; long int_one[1] = {1}; long dim0[1], dim1[1]; extern dgemv_(char *trans, long *m, long *n, double *alpha, double *a, long *lda, double *x, long *incx, double *beta, double *Y, long *incy); if (!PyArg_ParseTuple(args, "dOO", factor, &input1, &input2)) return NULL; matrix = (PyArrayObject *) PyArray_ContiguousFromObject(input1, PyArray_DOUBLE, 2, 2); if (matrix == NULL) return NULL; vector = (PyArrayObject *) PyArray_ContiguousFromObject(input2, PyArray_DOUBLE, 1, 1); if (vector == NULL) return NULL; if (matrix->dimensions[1] != vector->dimensions[0]) { PyErr_SetString(PyExc_ValueError, "array dimensions are not compatible"); return NULL; } dimensions[0] = matrix->dimensions[0]; result = (PyArrayObject *)PyArray_FromDims(1, dimensions, PyArray_DOUBLE); if (result == NULL) return NULL; dim0[0] = (long)matrix->dimensions[0]; dim1[0] = (long)matrix->dimensions[1]; dgemv_("T", dim1, dim0, factor, (double *)matrix->data, dim1, (double *)vector->data, int_one, real_zero, (double *)result->data, int_one); return PyArray_Return(result); } -- Michiel de Hoon, Assistant Professor University of Tokyo, Institute of Medical Science Human Genome Center 4-6-1 Shirokane-dai, Minato-ku Tokyo 108-8639 Japan http://bonsai.ims.u-tokyo.ac.jp/~mdehoon

3 2

Compact way of performing array math with specified result type?
by Russell E. Owen Feb. 8, 2010

Feb. 8, 2010

I often find myself doing simple math on sequences of numbers (which might or might not be numpy arrays) where I want the result (and thus the inputs) coerced to a particular data type. I'd like to be able to say: numpy.divide(seq1, seq2, dtype=float) but ufuncs don't allow on to specify a result type. So I do this instead: numpy.array(seq1, dtype=float) / numpy.array(seq2, dtype=float) Is there a more compact solution (without having to create the result array first and supply it as an argument)? -- Russell

4 3

Re: [Numpy-discussion] Floating point exception with numpy and embedded python interpreter
by Andrew Straw Sept. 20, 2007

Sept. 20, 2007

Arkaitz Bitorika wrote: > > On 19 Apr 2006, at 18:37, Andrew Straw wrote: > >> >>> I've been embedding Python for ages on the same program and other >>> modules work fine, only numpy fails. >> >> >> Most other modules don't use the SSE units, so wouldn't get hit by such >> a bug. > > > Is there a way of not using those units from numpy, to check if > that's what's going on? I think that numpy only accesses the SSE units through ATLAS or other external library. So, build numpy without ATLAS. But I'm not 100% sure anymore if there aren't any optimizations that directly use SSE if it's available. > Or alternatively, how would I check if my program is messing with the > SSE bits? Hmm, I think that's a bit hairy. I'd suggest simply asking the C++ library's mailing list if they alter the error bits on the control registers of the SSE unit. (Out of curiousity, what library is it?) If you want hairy, though, I think you'd have to check from C with the appropriate calls -- I'd start with the source code in that bug report. It looks like they're inlining an assembly statement to query a SSE control register.

5 7

rant against from numpy import * / from pylab import *
by Sebastian Haase Aug. 9, 2007

Aug. 9, 2007

Hi! I use the wxPython PyShell. I like especially the feature that when typing a module and then the dot "." I get a popup list of all available functions (names) inside that module. Secondly, I think it really makes code clearer when one can see where a function comes from. I have a default import numpy as N executed before my shell even starts. In fact I have a bunch of my "standard" modules imported as <some single capital letter>. This - I think - is a good compromise to the commonly used "extra typing" and "unreadable" argument. a = sin(b) * arange(10,50, .1) * cos(d) vs. a = N.sin(b) * N.arange(10,50, .1) * N.cos(d) I would like to hear some comments by others. On a different note: I just started using pylab, so I did added an automatic "from matplotlib import pylab as P" -- but now P contains everything that I already have in N. It makes it really hard to *find* (as in *see* n the popup-list) the pylab-only functions. -- what can I do about this ? Thanks, Sebastian

2 3

Aligning an array on Windows
by Francesc Altet July 7, 2007

July 7, 2007

Hi, Some time ago I made an improvement in speed on the numexpr version of PyTables so as to accelerate the operations with unaligned arrays (objects that can appear quite commonly when dealing with columns of recarrays, as PyTables does). This improvement has demostrated to work correctly and flawlessly in Linux machines (using GCC 4.x and in both 32-bit and 64-bit Linux boxes) for several weeks of intensive testing. Moreover, its speed-up is ranging from a 40% on modern processors and up to a 70% in older ones, so I'd like to keep it. The surprise came today when I tried to compile the same code on a Windows box (Win XP Service Pack 2) using MSVC 7.1, through the free (as in beer) Toolkit 2003. The compilation process went fine, but the problem is that I'm getting crashes from time to time when running the numexpr test suite. After some in-depth investigation, I'm pretty sure that the problem is in a concrete part of the code that I'd modified for this improvement. IMO, the affected code is in numexpr/interp_body.c and reads like: case OP_COPY_II: VEC_ARG1(memcpy(dest, x1, sizeof(int)); dest += sizeof(int); x1 += stride1); case OP_COPY_LL: VEC_ARG1(memcpy(dest, x1, sizeof(long long)); dest += sizeof(long long); x1 += stride1); case OP_COPY_FF: VEC_ARG1(memcpy(dest, x1, sizeof(double)); dest += sizeof(double); x1 += stride1); case OP_COPY_CC: VEC_ARG1(memcpy(dest, x1, sizeof(double)*2); dest += sizeof(double)*2; x1 += stride1); This might seem complicated, but it is not. Each of the OP_COPY_XX is a function that has to copy source (x1) to destination (dest) for int, long long, double and complex data types (this function will be called in a loop for copying all the data in array). The code for the original numexpr reads as: case OP_COPY_BB: VEC_ARG1(b_dest = b1); case OP_COPY_II: VEC_ARG1(i_dest = i1); case OP_COPY_FF: VEC_ARG1(f_dest = f1); case OP_COPY_CC: VEC_ARG1(cr_dest = c1r; ci_dest = c1i); i.e. the copy is done through direct assignment. This can be done because, in the original numexpr, an array is always guaranteed to reach this part of the code (the computing kernel) in the aligned form. But in my code, this is not guaranteed (the copy is made precisely for alignment purposes), so this is why I need to make use of memcpy/memmove calls. The thing I don't see is why my version of the code can create problems on Windows platforms and work perfectly on Linux ones. I've tried to use memmove instead of memcpy, but the problem persists. I've had a look at how numpy makes an 'aligned' copy of an unaligned array, and it seems to me that it uses memcpy/memmove (not sure when you should use one or another) just as I use it above. It might be possible that the problem is in another place, but my tests reaffirm me in the possibility that something is wrong with my copy code above (but again, I can't see where). Of course, we can get rid of this optimization but it is a bit depressing to have renounce to it just because it doesn't work on Windows :( Thanks in advance for any hint you may provide! -- Francesc Altet | Be careful about using the following code -- Carabos Coop. V. | I've only proven that it works, www.carabos.com | I haven't tested it. -- Donald Knuth

2 4

SciPy Journal
by Travis Oliphant June 7, 2007

June 7, 2007

Hi everybody, I'm sorry for the cross posting, but I wanted to reach a wide audience and I know not everybody subscribes to all the lists. I've been thinking more about the "SciPy Journal" that we discussed before and I have some thoughts. 1) I'd like to get it going so that we can push out an electronic issue after the SciPy conference (in September) 2) I think it's scope should be limited to papers that describe algorithms and code that are in NumPy / SciPy / SciKits. Perhaps we could also accept papers that describe code that depends on NumPy / SciPy that is also easily available. 3) I'd like to make a requirement for inclusion of new code in SciPy that it have an associated journal article describing the algorithms, design approach, etc. I don't see this journal article as being user-interface documentation for the code. I see this is as a place to describe why the code is organized as it is and to detail any algorithms that are used. 4) The purpose of the journal as I see it is to a) provide someplace to document what is actually done in SciPy and related software. b) provide a teaching tool of numerical methods with actual "people use-it" code that would be useful to researchers, students, and professionals. c) hopefully clever new algorithms will be developed for SciPy by people using Python that could be show-cased here d) provide a peer-review publication opportunity for people who contribute to open-source software 5) We obviously need associate editors and people willing to review submitted articles as well as people willing to submit articles. I have two articles that can be submitted within the next two months. What do other people have? As an example of the kind of thing a SciPy Journal would be useful for. I have recently over-hauled the interpolation.py file for SciPy by incorporating the B-spline stuff that is partly in fitpack. In the process I noticed two things: 1) I have (what seems to me) a different recursive algorithm for calculating derivatives of B-splines than I could find in fitpack. 2) I have developed a different way to determine the K-1 extra degrees of freedom for Kth-order spline fitting than I have seen before. The SciPy Journal would be a great place to document both of these things while describing the spline interpolation design of scipy.interpolate It is true that I could submit this stuff to other journals, but it seems like that doing that makes the information harder to find in the future and not easier. I'm also dissatisfied with how information exclusionary academic journals seem to be. They are catching up, but they are still not as accessible as other things available on the internet. Given the open nature of most scientific research, it is remarkable that getting access to the information is not as easy as it should be with modern search engines (if your internet domain does not subscribe to the e-journal). Comments and feedback is welcome. -Travis

10 12

Re: [Numpy-discussion] Vista installer?
by Ryan Krauss June 4, 2007

June 4, 2007

Success. I think this is officially a viable alternative for running Scipy/Numpy/Matplotlib/IPython in Vista. I turned off my firewall on my host (XP) and I was able to browse the internet in Ubuntu. Andrew was right that I still couldn't ping, but I seem to have everything else. And with shared folders enabled I seem to have some default SMB stuff set up. I couldn't successfully browse the shares of my partitions (C$, E$, ...), but I did specifically share a folder (from windows) and editted a Python file in that folder from within VMware. The available shared folders showed up under Places > Network as I should have expected for Ubuntu SMB shares. Another of my students supposedly has everything working in Vista with no problems. I will have to check that out and report back. Ryan On 5/29/07, Andrew Straw <strawman(a)astraw.com> wrote: > Ryan Krauss wrote: > > Thanks for your help in getting me up and going with VMWare. Andrew, > > your virtual appiance seems to work quite well. > > > > I have two hurdles left (I hope there are only 2). The first is still > > networking. I have NAT chosen in the VMWare player and that seems to > > give me limited connectivity. I seem to be able to use apt but can't > > ping or browse the internet. Any thoughts on how to trouble shoot > > this? > I'm guessing that the firewall on the host computer is preventing a lot > of this. I'm not sure I'd expect ping to work ever over NAT, as that's a > UDP thing, which is stateless and thus my guess is that NAT doesn't know > where to direct the incoming packet. However, I think you should be able > to browse, and the fact that apt works but browsing doesn't is certainly > odd, since the standard apt repositories are simply web servers. > > You can try bridged mode. That's usually what I do, but I figured it was > less likely to work than NAT (because it would require 2 external IP > addresses, and I'm not sure how your network is configured), so although > it allows full connectivity, including UDP, I chose NAT for that image. > > > > My second hurdle is how to move files between the host and guest OS. > > I see a shared folders option in the player, but don't seem to be able > > to add to the list of folders. My students need to be able to import > > python modules I give them and then submit files they would be > > developing within the guest OS. I googled for vmware shared folders > > and it seems like a non-free option. > Hmm, I haven't done any sharing through VMware server. On my mac, I have > using VMware fusion, but generally found that using SSH was more > reliable. So, maybe you can have them SSH (well, scp, perhaps using > WinSCP on the Windows side) between their two virtual computers? > > > > I haven't found a good "VMWare for dummies" or getting started guide, > > so feel free to redirect me. > VMware's docs are pretty good but necessarily get indexed by Google -- > make sure you browse their website. > > Good luck, > Andrew > > > > Thanks, > > > > Ryan > > > > On 5/29/07, David Cournapeau <david(a)ar.media.kyoto-u.ac.jp> wrote: > >> Ryan Krauss wrote: > >> > Hey David, > >> > > >> > Thanks for all your help. I run Ubuntu stand alone on my dual boot > >> > machines and have no network issues (actually, the computer running > >> > VMWare is dual boot). So, I think my network issues are VMWare > >> > specific. I don't know how it handles the hardware. > >> In this case, it should be even easier. In ethernet options of vmware, > >> you should choose NAT, and then vmware will have a dhcp server on the > >> HOST which will provide an address by dhcp to your guest. > >> > >> David > >> > >

2 2

corrcoef of masked array
by Jesper Larsen June 4, 2007

June 4, 2007

Hi numpy users, I have a masked array of dimension (nvariables, nobservations) that contain missing values at arbitrary points. Is it safe to rely on numpy.corrcoeff to calculate the correlation coefficients of a masked array (it seems to give reasonable results)? Cheers, Jesper

3 5

flatindex
by Tobias Knopp June 1, 2007

June 1, 2007

Hi! I was looking for a method to find the indices of the smallest element of an 3-dimensional array a. Therefore i used a.argmax() The problem was, that argmax gives me a flat index. My question is, if there is a build-in function to convert the flat index back to a multidimensional one. I know how to write such a procedure but was curious if one exists in numpy. Thanks for response and the really impressive software package. Tobi

4 4

ATLAS,LAPACK compilation - help!
by sittner＠lkb.ens.fr June 1, 2007

June 1, 2007

Hello there, I'm new here, so excuse me if the solution is trivial: i have installed ATLAS and LAPACK on my ubuntu 7 dual core intel machine. now, when i try to install numpy, it tells me it doesn't find these libraries: " $ python setup.py install Running from numpy source directory. F2PY Version 2_3816 blas_opt_info: blas_mkl_info: libraries mkl,vml,guide not found in /usr/local/lib libraries mkl,vml,guide not found in /usr/lib NOT AVAILABLE atlas_blas_threads_info: Setting PTATLAS=ATLAS libraries lapack,blas not found in /usr/local/lib/ATLAS/src/ libraries lapack,blas not found in /usr/local/lib/ATLAS libraries lapack,blas not found in /usr/local/lib libraries lapack,blas not found in /usr/lib NOT AVAILABLE atlas_blas_info: libraries lapack,blas not found in /usr/local/lib/ATLAS/src/ libraries lapack,blas not found in /usr/local/lib/ATLAS libraries lapack,blas not found in /usr/local/lib libraries lapack,blas not found in /usr/lib NOT AVAILABLE ......" I have installed ATLAS and lapack with no errors. ATLAS is in usr/local/lib/ATLAS/: $ ls /usr/local/lib/ATLAS bin doc interfaces Make.Linux_UNKNOWNSSE2_2 README tune CONFIG include lib makes src config.c INSTALL.txt Makefile Make.top tst.o so, what seems to be the problem? thanks, t

5 7