Mailman 3 August 2009 - NumPy-Discussion

Power distribution
by alan＠ajackson.org 08 Aug '09

08 Aug '09

Documenting my way through the statistics modules in numpy, I ran into the Power Distribution. Anyone know what that is? I Googled for it, and found a lot of stuff on electricity, but no reference for a statistical distribution of that name. Does it have a common alias? -- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan(a)ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | -----------------------------------------------------------------------

4 11

ATLAS, NumPy and Threading
by Peter Jeremy 08 Aug '09

08 Aug '09

[Apologies if anyone sees this twice - the first copy appears to have disappeared into a black hole] Should ATLAS be built with or without threading support for use with NumPy? The NumPy documentation just says that ATLAS will be used if found but gives no indication of how ATLAS should be built. I have found that system_info.py explicitly selects threaded versions of libatlas and liblapack on FreeBSD but have been unable to find any rationale behind this in the available SVN logs. -- Peter Jeremy

2 1

numpy.random.pareto, m equal zero
by josef.pktd＠gmail.com 07 Aug '09

07 Aug '09

Does it make any (statistical) sense to have numpy.random.pareto produce random numbers that start at zero? Can we change it to start at 1 which is the usual default? Notation from http://docs.scipy.org/numpy/docs/numpy.random.mtrand.RandomState.pareto/ The probability density for the Pareto distribution is .. math:: p(x) = \\frac{am^a}{x^{a+1}} where :math:`a` is the shape and :math:`m` the location constraints from Johnson, Kotz, Balakrishnan vol1 page 574 m>0, a>0, x>=m 1) as m goes to zero, the pdf goes to zero for every point, (mean, variance go to zero, essentially masspoint at zero) 2) quote from http://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/parpdf.htm (their `a` is our `m`) " Note that although the a (=m JP) parameter is typically called a location parameter (and it is in the sense that it defines the lower bound), it is not a location parameter in the technical sense that the following relation does not hold: f(x;gamma,a) = f((x-a);gamma,0) For this reason, Dataplot treats a (=m JP) as a shape parameter. In Dataplot, the a (=m JP) shape parameter is optional with a default value of 1. " my conclusion: --------------------- What numpy.random.pareto actually produces, are random numbers from a pareto distribution with lower bound m=1, but location parameter loc=-1, that shifts the distribution to the left. To actually get useful random numbers (that are correct in the usual usage http://en.wikipedia.org/wiki/Pareto_distribution), we need to add 1 to them. stats.distributions doesn't use mtrand.pareto (why?), so I never needed to check this before. rvs_pareto = 1 + numpy.random.pareto(a, size) for correction in some calculation, see the thread on the power distribution. Do we have to live with loc=-1, or can we change it, or am I misinterpreting something (which wouldn't be the first time either)? Josef

1 0

Optimized half-sizing of images?
by Christopher Barker 07 Aug '09

07 Aug '09

(second try on this message. the fist time I included a test PNG that made it too large) Hi folks, We have a need to to generate half-size version of RGB images as quickly as possible. PIL does a pretty good job, but it dawned on me that in the special case of a half-size, one might be able to do it faster with numpy, simply averaging the four pixels in the larger image to create one in the small. I'm doing tiling, and thus reducing 512x512 images to 256x256, so I imagine I'm making good use of cache (it does get pretty pokey with really large images!) What I have now is essentially this : # a is a (h, w, 3) RGB array a2 = a[0::2, 0::2, :].astype(np.uint16) a2 += a[0::2, 1::2, :] a2 += a[1::2, 0::2, :] a2 += a[1::2, 1::2, :] a2 /= 4 return a2.astype(np.uint8) time: 67.2 ms per loop I can speed it up a bit if I accumulate in a uint8 and divide as I go to prevent overflow: a2 = a[0::2, 0::2, :].astype(np.uint8) / 4 a2 += a[0::2, 1::2, :] / 4 a2 += a[1::2, 0::2, :] / 4 a2 += a[1::2, 1::2, :] / 4 return a2 time: 46.6 ms per loop That does lose a touch of accuracy, I suppose, but nothing I can see. Method 1 is about twice as slow as PIL's bilinear scaling. Can I do better? It seems it should be faster if I can avoid so many separate loops through the array. I figure there may be some way with filter or convolve or ndimage, but they all seem to return an array the same size. Any ideas? (Cython is another option, of course) -Chris Test code enclosed. -- 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

5 8

concatenate and non-contiguous arrays?
by Christopher Barker 07 Aug '09

07 Aug '09

Hi all, I just noticed that np.concatenate does not necessarily produce contiguous arrays. I has figured that it was making a copy, so would produce a C-contiguous array, but not so: In [88]: a = np.arange(60).reshape((4,5,3)) In [89]: b = np.concatenate((a, a[:, -1:, :]), axis=1) In [90]: b.flags Out[90]: C_CONTIGUOUS : False F_CONTIGUOUS : False OWNDATA : False WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False I'm also not sure why OWNDATA is false is it sharing data with somethign else? It doesn't look like it is with a. I'll toss a ascontiguous() in my code to take care of this, but I'd like to understand it. Explanation? -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

dot documentation
by T J 07 Aug '09

07 Aug '09

Hi, the documentation for dot says that a value error is raised if: If the last dimension of a is not the same size as the second-to-last dimension of b. (http://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.htm) This doesn't appear to be the case: >>> a = array([[1,2],[3,4]]) >>> b = array([1,2]) >>> dot(a,b) array([5,11]) I can see *how* 5,11 is obtained, but it seems this should have raised a ValueError since the 2 != 1. So the actual code must do something more involved. When I think about broadcasting, it seems that maybe b should have been broadcasted to: --> array([[1,2],[1,2]]) and then the multiplication done as normal (but this would give a 2x2 result). Can someone explain this to me?

2 2

Vectorize ufunc
by T J 07 Aug '09

07 Aug '09

I was wondering why vectorize doesn't make the ufunc available at the topmost level.... >>> def a(x,y): return x + y >>> b = vectorize(a) >>> b.reduce Instead, the ufunc is stored at b.ufunc. Also, b.ufunc.reduce() doesn't seem to exist until I *use* the vectorized function at least once. Can this be changed so that it exists right away (and preferably at b.reduce instead of b.ufunc.reduce)?

2 1

Test failures for rev7299
by Christopher Hanley 07 Aug '09

07 Aug '09

Hi, I receive the following test errors after building numpy rev7229 from svn: ====================================================================== FAIL: test_simple_circular (test_multiarray.TestStackedNeighborhoodIter) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Users/chanley/dev/site-packages/lib/python/numpy/core/tests/test_multia rray.py", line 1344, in test_simple_circular assert_array_equal(l, r) File "/Users/chanley/dev/site-packages/lib/python/numpy/testing/utils.py", lin e 639, in assert_array_equal verbose=verbose, header='Arrays are not equal') File "/Users/chanley/dev/site-packages/lib/python/numpy/testing/utils.py", lin e 571, in assert_array_compare raise AssertionError(msg) AssertionError: Arrays are not equal (mismatch 6.66666666667%) x: array([[ 0., 1., 2.], [ 1., 2., 3.], [ 2., 3., 0.],... y: array([[ 0., 1., 2.], [ 1., 2., 3.], [ 2., 3., 0.],... ====================================================================== FAIL: test_simple_mirror (test_multiarray.TestStackedNeighborhoodIter) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Users/chanley/dev/site-packages/lib/python/numpy/core/tests/test_multia rray.py", line 1296, in test_simple_mirror assert_array_equal(l, r) File "/Users/chanley/dev/site-packages/lib/python/numpy/testing/utils.py", lin e 639, in assert_array_equal verbose=verbose, header='Arrays are not equal') File "/Users/chanley/dev/site-packages/lib/python/numpy/testing/utils.py", lin e 571, in assert_array_compare raise AssertionError(msg) AssertionError: Arrays are not equal (mismatch 6.66666666667%) x: array([[ 0., 1., 2.], [ 1., 2., 3.], [ 2., 3., 0.],... y: array([[ 0., 1., 2.], [ 1., 2., 3.], [ 2., 3., 0.],... ---------------------------------------------------------------------- Ran 2186 tests in 10.671s FAILED (KNOWNFAIL=1, SKIP=3, failures=2) <nose.result.TextTestResult run=2186 errors=0 failures=2> I'm running on an Intel MacBook Pro running OS X 10.5.8. I am using Python 2.5.1. Chris -- Christopher Hanley Senior Systems Software Engineer Space Telescope Science Institute 3700 San Martin Drive Baltimore MD, 21218 (410) 338-4338

3 2

datetime code
by Charles R Harris 06 Aug '09

06 Aug '09

Travis, Robert, Is there any reason not to merge the c code in the datetime branch at this time? If not, I will do it. Chuck

1 0

Strange Error with NumPy
by Nanime Puloski 06 Aug '09

06 Aug '09

Can anyone explain to me why I receive an error(AtrributeError) in NumPy when I do numpy.sin(2**64), but not when I do numpy.sin(2.0**64), numpy.sin(float(2**64)) or even numpy.sin(2)? Where is the root problem in all of this?

3 2