Mailman 3 July 2013 - NumPy-Discussion

f2py build with mkl lapack
by sunghyun Kim July 8, 2013

July 8, 2013

Hi I'm trying to use fortran wrapper f2py with intel's mkl following is my command LIB='-L/opt/intel/Compiler/11.1/064/mkl/lib/em64t/ -lguide -lpthread -lmkl_core -lmkl_intel_lp64 -lmkl_sequential' INC=-I/opt/intel/Compiler/11.1/064/mkl/include f2py --fcompiler=intelem $INC $LIB -m solveLE -c solveLE2.f solveLE2.f is simple fortran code using lapack's linear equation solver SGESV ============= CALL SGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) ============= When i use the command, compile was done. But when I use the solveLE.so, I received following error massage ================== $python test.py python: symbol lookup error: /opt/intel/Compiler/11.1/064/mkl/lib/em64t/libmkl_lapack.so: undefined symbol: mkl_lapack_sgetrf ================== I think "mkl_lapack_sgetrf" is defined in -lmkl_sequential. I don't know what should I do. Any help would be greatly appreciated! Sunghyun Kim Ph.D. Candidate Theoretical Condensed Matter Physics Group. KAIST 291 Daehak-ro(373-1 Guseong-dong), Yuseong-gu, Daejeon 305-701, Republic of Korea +10-4144-5946

3 3

Re: [Numpy-discussion] SLARTG
by Alan G Isaac July 5, 2013

July 5, 2013

On 7/5/2013 2:50 PM, Alan G Isaac wrote: > I see that CLARTG is here: > https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/arpack… > > But is there a Python interface in SciPy? > (Or any other SciPy access to Givens rotation?) Sorry, that was SLARTG, whereas CLARTG is here: https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/arpack… But the question stands: is there a Python interface to Givens rotation? Thanks, Alan Isaac

1 0

PyArray_PutTo Question
by Mark Janikas July 5, 2013

July 5, 2013

Hi All, I am a bit new to the NumPy C-API and I am having a hard time with placing results into output arrays... I am using PyArray_TakeFrom to grab an input dimension of data, then do a calculation, then I want to pack it back to the output... yet the PutTo function does not have an axis argument like the TakeFrom does... I am grabbing by column in a two-dimensional array and I would like to pack it that way. I know that I can build the result in reverse and pack the columns into rows and then reshape the output... but I am wondering why the PutTo does not behave exactly like the take-from does?... The python implementation "numpy.put" also does not have the axis... so I guess I can see the one-to-one reason for the omission. However, is building in reverse and reshaping the normal way to pack by column? Thanks much! MJ

1 0

simple select(): for anyone with too many conditions for np.select(), or scalar-valued choicelists.
by Graeme B. Bell July 5, 2013

July 5, 2013

I've made a drop-in replacement for select() which works with large numbers of conditions, and which consistently outperforms numpy.select for my use case (scalar condlist). It fixes a couple of other issues too, and (I feel) improves the internal documentation of the code. I have included benchmarks and some tests. https://github.com/gbb/numpy-simple-select The numpy dev team are welcome to include all or part of this code into the main numpy distribution or it can be kept separate or ignored if they prefer. :-) If you need more than 30 ndarrays in your 'condlist', or if you have an all-scalar choicelist, I think you will find this code particularly interesting. Formal test coverage is incomplete, but I think this is still going to be quite useful for some people. Have a nice weekend, Graeme

1 0

Covariance matrix from polyfit
by David Schaich July 4, 2013

July 4, 2013

Hi all, I recently adopted python, and am in the process of replacing my old analysis tools. For simple (e.g., linear) interpolations and extrapolations, in the past I used gnuplot. Today I set up the equivalent with polyfit in numpy v1.7.1, first running a simple test to reproduce the gnuplot result. A discussion on this list back in February alerted me that I should use 1/sigma for the weights in polyfit as opposed to 1/sigma**2. Fine -- that's not what I'm used to, but I can make a note. http://mail.scipy.org/pipermail/numpy-discussion/2013-February/065649.html Another issue mentioned in that thread is scaling the covariance matrix by fac = resids / (len(x) - order - 2.0) This wreaked havoc on the simple test I mentioned above (and include below), which fit three data points to a straight line. I spent hours trying to figure out why numpy returned negative variances, before tracking down this line. And, indeed, if I add a fake fourth data point, I end up with inf's and nan's. There is some lengthy justification for subtracting that 2 around line 590 in lib/polynomial.py. Fine -- it's nothing I recall seeing before (and I removed it from my local installation), but I'm just a new user. However, I do think it is important to fix polyfit so that it doesn't produce pathological results like those I encountered today. Here are a couple of possibilities that would let the subtraction of 2 remain: * Check whether len(x) > order + 2, and if it is not, either ** Die with an error ** Scale by resids / (len(x) - order) instead of resids / (len(x) - order - 2.0) * Don't bother with this scaling at all, leaving it to the users (who can subtract 2 if they want). This is what scipy.optimize.leastsq does, after what seems to be a good deal of discussion: "This matrix must be multiplied by the residual variance to get the covariance of the parameter estimates – see curve_fit." http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.… https://github.com/scipy/scipy/pull/448 I leave it to those of you with more numpy experience to decide what would be the best way to go. Cheers, David http://www-hep.colorado.edu/~schaich/ +++ Here's the simple example: >>> import numpy as np >>> m = np.array([0.008, 0.01, 0.015]) >>> dat = np.array([1.0822582, 1.0805417, 1.0766624]) >>> weight = np.array([1/0.000370, 1/0.000355, 1/0.000249]) >>> out, cov = np.polyfit(m, dat, 1, full=False, w=weight, cov=True) >>> print out, '\n', cov [-0.79269957 1.08854252] [[ -2.34965006e-04 2.84428412e-06] [ 2.84428412e-06 -3.66283662e-08]] >>> print np.sqrt(-1. * cov[0][0]) 0.0153285682895 >>> print np.sqrt(-1. * cov[1][1]) 0.000191385386578 +++ Gnuplot gives +++ Final set of parameters Asymptotic Standard Error ======================= ========================== A = -0.792719 +/- 0.01533 (1.934%) B = 1.08854 +/- 0.0001914 (0.01758%) +++ so up to the negative sign, all is good. For its part, scipy.optimize.leastsq needs me to do the scaling: +++ >>> import numpy as np >>> from scipy import optimize >>> m = np.array([0.008, 0.01, 0.015]) >>> dat = np.array([1.0822582, 1.0805417, 1.0766624]) >>> err = np.array([0.000370, 0.000355, 0.000249]) >>> linear = lambda p, x: p[0] * x + p[1] >>> errfunc = lambda p, x, y, err: (linear(p, x) - y) / err >>> p_in = [-1., 1.] >>> all_out = optimize.leastsq(errfunc, p_in[:], args=(m, dat, err), full_output = 1) >>> out = all_out[0] >>> cov = all_out[1] >>> print out, '\n', cov [-0.79269959 1.08854252] [[ 3.40800756e-03 -4.12544212e-05] [ -4.12544212e-05 5.31270007e-07]] >>> chiSq_dof = ((errfunc(out, m, dat, err))**2).sum() / (len(m) - len(out)) >>> cov *= chiSq_dof >>> print cov [[ 2.34964190e-04 -2.84427528e-06] [ -2.84427528e-06 3.66282716e-08]] +++

2 1

Reducing the rounding error of np.sum
by Julian Taylor July 4, 2013

July 4, 2013

hi, numpys implementation of sum is just a simple: for i in d: sum += d[i] this suffers from rather high rounding errors in the order of the d.size * epsilon. Consider: (np.ones(50000) / 10.).sum() 5000.0000000006585 There are numerous algorithms which reduce the error of this operation. E.g. python implements one in math.fsum which is accurate but slow compared to np.sum [0]. Numpy is currently lacking a precise summation, but I think it would make sense to add one. The question is whether we go with the python approach of adding a new function which is slower but more precise or if we even change the default summation algorithm (or do nothing :) ). For a new function I guess the method used in python itself makes sense, its probably well chosen by the python developers (though I did not lookup the rational for the choice yet). For replacing the default, two algorithms come to my mind, pairwise summation [1] and kahan summation (compensated sum) [2]. pairwise summation adds in pairs so usually the magnitude of the two operands is the same magnitude, this produces an error of O(log n * epsilon) for the common case. This algorithm has the advantage that it is almost as fast as the naive sum with an reasonable error. Problematic might be the buffering numpy does when reducing, this would limit the error reduction to the buffer size. kahan summation adds some extra operations to recover the rounding errors. This results in an error of o(epsilon). It is four times slower than the naive summation but this can be (partially) compensated by vectorizing it. It has the advantage of lower error, simpler implementation and buffering does not interfere. I did prototype implementations (only unit strides) of both in [3, 4]. Any thoughts on this? [0] http://docs.python.org/2/library/math.html#number-theoretic-and-representat… [1] http://en.wikipedia.org/wiki/Pairwise_summationg [2] http://en.wikipedia.org/wiki/Kahan_summation_algorithm [3] https://github.com/juliantaylor/numpy/tree/pairwise [4] https://github.com/juliantaylor/numpy/tree/kahan

3 2

Equality not working as expected with ndarray sub-class
by Thomas Robitaille July 4, 2013

July 4, 2013

Hi everyone, The following example: import numpy as np class SimpleArray(np.ndarray): __array_priority__ = 10000 def __new__(cls, input_array, info=None): return np.asarray(input_array).view(cls) def __eq__(self, other): return False a = SimpleArray(10) print (np.int64(10) == a) print (a == np.int64(10)) gives the following output $ python2.7 eq.py True False so that in the first case, SimpleArray.__eq__ is not called. Is this a bug, and if so, can anyone think of a workaround? If this is expected behavior, how do I ensure SimpleArray.__eq__ gets called in both cases? Thanks, Tom ps: cross-posting to stackoverflow

3 3

Fancy indexing oddity
by Bradley M. Froehle July 3, 2013

July 3, 2013

A colleague just showed me this indexing behavior and I was at a loss to explain what was going on. Can anybody else chime in and help me understand this indexing behavior? >>> import numpy as np >>> np.__version__ '1.7.1' >>> A = np.ones((2,3,5)) >>> mask = np.array([True]*4 + [False], dtype=bool) >>> A.shape (2, 3, 5) >>> A[:,:,mask].shape (2, 3, 4) >>> A[:,1,mask].shape (2, 4) >>> A[1,:,mask].shape (4, 3) # Why is this not (3, 4)? >>> A[1][:,mask].shape (3, 4) Thanks! Brad

2 1

trouble with numpy.float64 multiplied with cvxopt.matrix
by Jonathan Stickel July 2, 2013

July 2, 2013

I recently ran into some trouble with multiplying scalar variables of type numpy.float64 with cvxopt matrices. The cvxopt matrix is converted to a numpy array. I first reported this at the cvxopt google groups, where I give an example: https://groups.google.com/forum/#!topic/cvxopt/4suFNOY75E4 The response I got is that it would be difficult to correct this in cvxopt, and that a workaround would be to do A*s rather than s*A (where s is the scalar of type numpy.float64 and A is the cvxopt matrix). I inferred from this that the leading object takes charge of how operator works, including the type conversion. Now, I don't know much about low-level programming for operators and type conversion, but I want to ask whether this should be considered a bug or whether correcting this might be a reasonable feature request in numpy. It seems to me that scalars of any numpy type, when operated with other objects, should not change the high-level type of those objects. Thanks, Jonathan

1 0

matplotlib user survey 2013
by Michael Droettboom July 2, 2013

July 2, 2013

[Apologies for cross-posting] The matplotlib developers want to hear from you! We are conducting a user survey to determine how and where matplotlib is being used in order to focus its further development. This should only take a couple of minutes. To fill it out, visit: https://docs.google.com/spreadsheet/viewform?fromEmail=true&formkey=dHpQS25… <https://docs.google.com/spreadsheet/viewform?fromEmail=true&formkey=dHpQS25…> Please forward to your colleagues, particularly those who don't read these mailing lists. Cheers, Michael Droettboom, and the matplotlib team

1 0

np.insert with axis=-1
by Félix Hartmann July 1, 2013

July 1, 2013

Hi all, I recently upgraded from Numpy 1.6.2 to 1.7.1 on my Debian testing, and then got a bug in a program that was previously working. It turned out that the problem comes from the np.insert function when the argument `axis=-1` is given. Here is a minimal example: >>> u = np.zeros((2,3,4)) >>> ui = np.ones((2,3)) >>> u = np.insert(u, 1, ui, axis=-1) The last line should be equivalent to >>> u = np.insert(u, 1, ui, axis=2) It was indeed the case in Numpy 1.6, but in 1.7.1 it raises a ValueError exception. Note that the problem seems specific to axis=-1, and not to all negative axis values, since the following example works as expected: >>> u = np.zeros((2,3,4)) >>> ui = np.ones((2,4)) >>> u = np.insert(u, 1, ui, axis=-2) # equivalent to axis=1 I didn't check on current master, so maybe things have changed since 1.7.1. If they have not, do you think a bug report would be relevant? Cheers, Félix

2 2