Mailman 3 October 2013 - NumPy-Discussion

Efficient square distance computation
by Ke Sun Oct. 10, 2013

Oct. 10, 2013

Dear all, I have written the following function to compute the square distances of a large matrix (each sample a row). It compute row by row and print the overall progress. The progress output is important and I didn't use matrix multiplication. I give as input a 70,000x800 matrix. The output should be a 70,000x70,000 matrix. The program runs really slow (16 hours for 1/3 progress). And it eats 36G memory (fortunately I have enough). Could you give some insights on how to modify the code to be efficient and to eat less memory? thanks, Ke Sun def dist2_large( data ): import time if data.ndim != 2: raise RuntimeError( "data should be a matrix" ) N,D = data.shape print 'using the sample-wise implementation' print '%d samples, %d dimensions' % (N,D) start_t = time.time() d2 = np.zeros( [N,N] ) for i in range( N ): print "\r%5d/%d" % (i+1, N), for j in range( N ): d2[i,j] = ((data[i] - data[j])**2).sum() total_t = time.time() - start_t hours = (total_t / 3600) minutes = (total_t % 3600) / 60 print "\nfinished in %2dh%2dm" % (hours, minutes) return d2

5 11

Reading from binary file with memmap, with offset
by Andreas Hilboll Oct. 10, 2013

Oct. 10, 2013

Hi, I have a problem using memmap correctly. I need to read a data file which consists of an ASCII header and appended binary single precision floating point values. memmap complains that the "Size of available data is not a multiple of the data-type size." But as far as I can tell, the size *doas* match the data-type size. The problem is illustrated by the following code: ---8<------- url = "http://www.iup.uni-bremen.de/~hilboll/download/20120204.XD4_N2" localfile = "np_memmap.dat" import os import urllib import numpy as np # download data file if not os.access(localfile, os.R_OK): urllib.urlretrieve(url, localfile) with open(localfile, "r") as fd: # read offset from first line of file offset = int(fd.readlines()[0].split()[-2]) # jump to begin of data block fd.seek(offset) # read until EOF blob = fd.read() print("Size of data blob [bytes]: {}".format(len(blob))) print("This can actually be divided by 4: {} / 4.0 = {}".format( len(blob), len(blob) / 4.0)) # go back to begin of data block fd.seek(offset) print("But it cannot be loaded as np.memmap with dtype float32:") np.memmap(fd, dtype="float32", mode="r") ---8<------- Any help is greatly appreciated :) -- Andreas.

3 4

Bug in numpy.correlate documentation
by Bernhard Spinnler Oct. 9, 2013

Oct. 9, 2013

The numpy.correlate documentation says: correlate(a, v) = z[k] = sum_n a[n] * conj(v[n+k]) In [1]: a = [1, 2] In [2]: v = [2, 1j] In [3]: z = correlate(a, v, 'full') In [4]: z Out[4]: array([ 0.-1.j, 2.-2.j, 4.+0.j]) However, according to the documentation, z should be z[-1] = a[1] * conj(v[0]) = 4.+0.j z[0] = a[0] * conj(v[0]) + a[1] * conj(v[1]) = 2.-2.j z[1] = a[0] * conj(v[1]) = 0.-1.j which is the time reversed version of what correlate() calculates. IMHO, the correlate() code is correct. The correct formula in the docs (which is also the correlation formula in standard text books) should be z[k] = sum_n a[n+k] * conj(v[n]) Cheers, Bernhard

3 3

Segfault with QR Decomposition
by Charanpal Dhanjal Oct. 9, 2013

Oct. 9, 2013

Oops, the second line should have been export OPENBLAS_NUM_THREADS=8

2 2

Segfault with QR Decomposition
by Charanpal Dhanjal Oct. 9, 2013

Oct. 9, 2013

I get a segmentation fault upon running the following: import numpy A = numpy.ones((7000000, 8)) Q, R = numpy.linalg.qr(A) on Python 2.7.3, Linux 64-bit using numpy 1.9.0.dev-ec3603f linked against OpenBLAS. If A is a smaller matrix then the QR decomposition works (for example A has shape (4000000, 8)). I haven't quite narrowed down the exact threshold where the crash occurs, however I know that the above A is 448 MB (Q and R are no bigger), and the machine in question has 32GB of RAM. I also tested scipy.linalg.qr (version 0.14.0.dev-ced994c) with the same results. I don't get the same problem on my laptop which is Python 2.7.3, Linux 64-bit but with numpy 1.8.0rc1 linked to OpenBLAS. Both machines have OpenBLAS 0.2.6. Does anyone have some insight into why this problem is occuring? Thanks very much for any help, Charanpal

2 2

Re: [Numpy-discussion] select column based on another column value
by Roelf Schreurs Oct. 9, 2013

Oct. 9, 2013

Hi I have the following array and want to calculate the average per day. Import numpy as np from datetime import date today = date(today.year, today.month, 9) yesterday = date(today.year, today.month, 8) twodaysago = date(today.year, today.month, 7) testdata = np.array([[today, "r", 3.2],[today, "r", 4.3],[today, "r", 5.3],[yesterday, "r", 6.3],[yesterday, "r", 9.3],[yesterday, "r", 1.3],[twodaysago, "r", 3.3],[twodaysago, "r", 1.2],[twodaysago, "r", 4.3]]) which produces array([[datetime.date(2013, 10, 9), 'r', 3.2], [datetime.date(2013, 10, 9), 'r', 4.3], [datetime.date(2013, 10, 9), 'r', 5.3], [datetime.date(2013, 10, 8), 'r', 6.3], [datetime.date(2013, 10, 8), 'r', 9.3], [datetime.date(2013, 10, 8), 'r', 1.3], [datetime.date(2013, 10, 7), 'r', 3.3], [datetime.date(2013, 10, 7), 'r', 1.2], [datetime.date(2013, 10, 7), 'r', 4.3]], dtype=object) And I want the output as (numbers aren't correct below): 2013-10-09, 4.54 2013-10-08, 5.43 2013-10-07, 2.76 But not sure how. Thanks Roelf Tibra Trading Europe Limited is Authorised and Regulated by the FCA (Company No. 06061713) Tibra Global Services Pty Limited London Branch (Company No. BR014370) Tibra Equities Europe Limited (Company No. 07283912) Registered in England & Wales - Level 12, 5 Aldermanbury Square London EC2V 7HR Unless otherwise stated, this communication is neither an offer nor solicitation of an offer to sell or purchase any investment The contents of this email including any attachments are confidential. If you have received this email in error, please advise the sender by return email and delete this email. Any unauthorised use of the contents of the email is prohibited and you must not disseminate, copy or distribute the message or use the information contained in the email or its attachments in any way. The views or opinions expressed are the author's own and may not reflect the views or opinions of Tibra. Tibra does not guarantee the integrity of any emails or attached files. E-mails may be interfered with, may contain computer viruses or other defects. Under no circumstances do we accept liability for any loss or damage which may result from your receipt of this message or any attachments.

2 2

Question about typenum
by Valentin Haenel Oct. 8, 2013

Oct. 8, 2013

Hi, I have a quick question about typenum. Certain functions, like `PyArray_SimpleNewFromData` `PyArray_SimpleNew` take a typeenum argument. Is there any way to go from typeenum to something that can be passed to the dtype constructor, like mapping 12 -> '<f8'? thanks, V-

2 4

Re: [Numpy-discussion] NumPy-Discussion Digest, Vol 85, Issue 19
by Siegfried Gonzi Oct. 7, 2013

Oct. 7, 2013

On 7 Oct 2013, at 21:16, numpy-discussion-request(a)scipy.org wrote: > > > Message: 5 > Date: Mon, 7 Oct 2013 13:38:53 -0500 > From: Aronne Merrelli <aronne.merrelli(a)gmail.com> > Subject: Re: [Numpy-discussion] Equivalent to IDL's help function > > > There isn't anything quite the same. (I think what you are really asking > for is a way to print the essential info about one variable name, at least > that is how I would use the IDL "help" procedure). In IPython, I use the > whos magic to do something similar, although it just prints this info for > all variables or all variables of one class, rather than just one variable. > I do not think there is a way to do it for just one variable. > > Here are some examples - you can see this works quite well but it will > become unwieldy if your interactive namespace becomes large: > > In [1]: x = 1; y = 2; z = 3.3; d = {'one':1, 'two':2, 'three':3} > > In [2]: whos > Variable Type Data/Info > ----------------------------- > d dict n=3 > x int 1 > y int 2 > z float 3.3 > > In [3]: whos dict > Variable Type Data/Info > ---------------------------- > d dict n=3 > > In [4]: xa = np.arange(111); ya = np.ones((22,4)) > > In [5]: whos ndarray > Variable Type Data/Info > ------------------------------- > xa ndarray 111: 111 elems, type `int64`, 888 bytes > ya ndarray 22x4: 88 elems, type `float64`, 704 bytes > > Hi [I hope I am not screwing up the digest reply function here]. I am after a "whos" which would work in a script. It is not very often that I develop code at the command line. I am definitely not one of the best programmers out there but I used "help" a lot in my IDL scripts and code. Our research group is migrating away from IDL towards Python. I think Python's help is not the same than IDL's help. I know copying things from other languages is not always a good idea but one cannot argue that IDL's help comes in rather handy while developing and testing code. > > > On Mon, Oct 7, 2013 at 12:15 PM, Siegfried Gonzi > <sgonzi(a)staffmail.ed.ac.uk>wrote: > >> Hi all >> >> What is the equivalent to IDL its help function, e.g. >> >> == >> IDL> a = make_array(23,23,) >> >> IDL> help,a >> >> will result in: >> >> A FLOAT = Array[23, 23] >> >> or >> >> IDL> a = create_struct('idl',23) >> >> IDL> help,a >> >> gives: >> >> A STRUCT = -> <Anonymous> Array[1] >> >> == >> >> I have been looking for it ever since using numpy. It would make my life >> so much easier. >> >> >> Thanks, Siegfried >> >> >> -- >> The University of Edinburgh is a charitable body, registered in >> Scotland, with registration number SC005336. >> >> _______________________________________________ >> NumPy-Discussion mailing list >> NumPy-Discussion(a)scipy.org >> http://mail.scipy.org/mailman/listinfo/numpy-discussion >> > The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

2 1

Re: [Numpy-discussion] [ANN] MATLAB ODE solvers - now available in Python (Dmitrey)
by David Goldsmith Oct. 7, 2013

Oct. 7, 2013

On Sun, Oct 6, 2013 at 10:00 AM, <numpy-discussion-request(a)scipy.org> wrote: > Message: 2 > Date: Sat, 05 Oct 2013 21:36:48 +0300 > From: Dmitrey <tmp50(a)ukr.net> > Subject: Re: [Numpy-discussion] [ANN] MATLAB ODE solvers - now > available in Python > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > Cc: numpy-discussion(a)scipy.org > Message-ID: <1380997576.559804301.aoynahlj(a)frv43.ukr.net> > Content-Type: text/plain; charset="utf-8" > > Seems like using the MATLAB solvers with MCR requires my wrappers > containing in several files to be compiled with MATLAB Compiler before. I > have no license for MATLAB thus I may have problems if I'll make it done > and will spread it with OpenOpt suite code, also, binary files are > incompatible with BSD license. > Darn, knew it was too good to be true. > On the other hand, IIRC a little bit obsolete MATLAB versions (I don't > think difference is essential) have more liberal licenses. > As for MATLAB solvers examples, I have already mentioned them in the mail > list, you could see them in http://openopt.org/ODE (just replace solver > name from scipy_lsoda to ode23s or any other), http://openopt.org/NLP , > http://openopt.org/SNLE Oooops, so sorry. :-o DG > > ---------------------- > Regards, D. http://openopt.org/Dmitrey >

2 1

Problem with numpy's array reference assignment?
by Bernhard Spinnler Oct. 6, 2013

Oct. 6, 2013

I have problems to get a piece of code to work with a new numpy/scipy version. The code essentially sets up a matrix Ryy and a vector Rya and solves the system of linear equations Ryy*c = Rya for c. Then it checks whether the resulting vector c satisfies the equation: Ryy*c must be equal to Rya. While this code runs fine on - python-2.7.5.p1, numpy-1.7.0, scipy-0.12.0.p0 - python-2.7.3, numpy-1.7.1, scipy-0.12.0 it fails on - python 2.7.2, numpy 1.9.0.dev-fde3dee, scipy 0.14.0.dev-4938da3 with error: AssertionError: Arrays are not almost equal to 6 decimals (mismatch 100.0%) x: array([ 9.+0.j, 8.+0.j, 7.+0.j]) y: array([ 7.+0.j, 6.+0.j, 5.+0.j]) The piece of code is: ------------------------------------ import numpy as np from numpy.testing import assert_array_almost_equal lag0 = 5 N = 3 D = 2 corr_ya = np.array([0,1,2,3,4,5,6,7,8,9+0j]) corr_yy = np.array([0,1,2,3,4,5,4,3,2,1]) Rya = corr_ya[lag0+D:lag0+D-N:-1] Ryy = np.zeros((N, N), complex) for i in np.arange(N): Ryy[i,:] = corr_yy[lag0-i:lag0-i+N] c = np.linalg.solve(Ryy, Rya) # check result Ryy_x_c = np.dot(Ryy, c) print "Ryy*c =", repr(Ryy_x_c) print "Rya =", repr(Rya) assert_array_almost_equal(Ryy_x_c, Rya) ------------------------------------ I guess that it has something to do with numpy assigning arrays by reference and not copying them, since the problem goes away when vector Rya is copied before passing it to solve, i.e. doing Rya_copy = copy(Rya) and passing Rya_copy to solve. The problem also does not occur when the initial array corr_ya is an integer array, e.g. by deleting the "+0j" from the last element of corr_ya above. (Normally, my initial arrays are complex. I just used artificially simplified arrays above to show the problem.) I could imagine that this is due to numpy copying the integers into a new complex array instead of creating a reference that is passed to solve. Could it be that a bug has been introduced in recent numpy/scipy version? Or am I misunderstanding something? Thanks, Bernhard

2 2