Mailman 3 April 2014 - SciPy-User

Question on spectral derivative
by جاك الفضة April 11, 2014

April 11, 2014

I have data that form a scalar field on a 2D grid, evenly spaced. The grid has a finite size. There is no particular periodicity patern in my data. I want to calculate the value of the gradient at each grid point with a spectral method. Using simple FFT does not seem to work, especially near the borders (probably because the data are not periodic). I have ears about DCT, but I am not sure. I think that I could extend my grid with one row/column at each of the 4 boundaries, duplicate the real boundary on these new boundary points to get a zero derivative boundary for the DCT algorithm. However, I am not sure how to get it working. Another possibility would be to extend my grid in some way that FFT will not create spurious oscillation near the boundaries, but again, I my attempts failed so far. I would appreciate advice on the best way to do this.

2 2

Optimization problem using fmin_slsqp
by Qing Yu April 11, 2014

April 11, 2014

Dear SciPy experts, I'm using the fmin_slsqp function in scipy.optimize to fit my data, I first defined the function f and the residual of fitting, and then used actual data x and y to fit the curve. However, I encountered an error that I could not explain. Here's a short version of my script: def f(x,a): return a[0]*exp((a[2]**2)*cos(x-a[1])-1)+a[3] def residual(a,x,y): return y-f(x,a) x=[0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6] y=[0.17,0.36,0.61,0.38,0.17,0.16] a0=[1,pi/3,pi/2,0] out, fx, its, imode, smode = fmin_slsqp(residual,a0,args=(x,y),bounds=[(None,None),(0,pi),(0,None),(None,None)], full_output=True) I ended up with the error: "failed in converting 8th argument `g' of _slsqp.slsqp to C/Fortran array." Does anyone have any idea of why this happened? Thanks so much in advance for any help! Best, Qing

3 11

sparse distance matrix ?
by josef.pktd＠gmail.com April 10, 2014

April 10, 2014

(sunday morning puzzle) Is there a way to get thresholded distance matrices as sparse matrices without going through the full dense distance array or explicit looping? I didn't see anything in scipy.spatial or in scikit-learn. usecase is linear algebra not nearest neighbor search example converts 2d score levels into 1d score levels. The example could be done with kronecker product, but not in the general case. >>> ng, ni = 5, 4 >>> x = 100 * np.repeat(np.arange(1,ng+1), ni) + np.tile(np.arange(ni), ng) >>> x array([100, 101, 102, 103, 200, 201, 202, 203, 300, 301, 302, 303, 400, 401, 402, 403, 500, 501, 502, 503]) dense absolute distance >>> dx1 = np.abs(x[:,None] - x) >>> dx1[:10, :10] array([[ 0, 1, 2, 3, 100, 101, 102, 103, 200, 201], [ 1, 0, 1, 2, 99, 100, 101, 102, 199, 200], [ 2, 1, 0, 1, 98, 99, 100, 101, 198, 199], [ 3, 2, 1, 0, 97, 98, 99, 100, 197, 198], [100, 99, 98, 97, 0, 1, 2, 3, 100, 101], [101, 100, 99, 98, 1, 0, 1, 2, 99, 100], [102, 101, 100, 99, 2, 1, 0, 1, 98, 99], [103, 102, 101, 100, 3, 2, 1, 0, 97, 98], [200, 199, 198, 197, 100, 99, 98, 97, 0, 1], [201, 200, 199, 198, 101, 100, 99, 98, 1, 0]]) I'd like to get the following as scipy.sparse matrices: close distances >>> (dx1) * (dx1 < 3) array([[0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [2, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 2, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 2, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 2], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 1], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0]]) reversed >>> (5 - dx1) * (dx1 < 3) array([[5, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 5, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [3, 4, 5, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 5, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 4, 5, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 3, 4, 5, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 5, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 4, 5, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 3, 4, 5, 4, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 4, 3, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 5, 4, 3, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 4, 5, 4, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 4, 5, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 4, 3, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 5, 4, 3], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 4, 5, 4], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 4, 5]]) boolean >>> (dx1 < 3).astype(int) array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1]]) >>> Josef

3 4

EuroScipy Reminder: call for abstracts closes in 4 days
by Nelle Varoquaux April 10, 2014

April 10, 2014

Hello everyone, Just a quick reminder that the EuroScipy call for abstracts closes on the 14th: don't forget to submit your talk proposal! It is in four days only! In short, EuroScipy is a cross-disciplinary gathering focused on the use and development of the Python language in scientific research. This event strives to bring together both users and developers of scientific tools, as well as academic research and state of the art industry. EuroSciPy 2014, the Seventh Annual Conference on Python in Science, takes place in *Cambridge, UK on 27 - 30 August 2014*. The conference features two days of tutorials followed by two days of scientific talks. The day after the main conference, developer sprints will be organized on projects of interest to attendees. The topics presented at EuroSciPy are very diverse, with a focus on advanced software engineering and original uses of Python and its scientific libraries, either in theoretical or experimental research, from both academia and the industry. The program includes keynotes, contributed talks and posters. Submissions for talks and posters are welcome on our website ( http://www.euroscipy.org/2014/). In your abstract, please provide details on what Python tools are being employed, and how. The deadline for submission is 14 April 2014. Also until 14 April 2014, you can apply for a sprint session on 31 August 2014. See https://www.euroscipy.org/2014/calls/sprints/ for details. Thanks, N

1 0

[OT] Job / Post doc position
by Erickson, Bradley J., M.D., Ph.D. April 9, 2014

April 9, 2014

Hi Everyone, I want to let anyone interested in a post-doctoral position be aware that I am currently hiring. The position is for a PhD with experience in image processing preferably using python. The exact field of your PhD is flexible. There is more than 1 position open, and the focus can range from computation of features in medical images (mostly brain MRI) to developing and applying machine learning methods to find important diagnostic and therapeutic features in the images. Our main development environment is python. The salary is competitive, and Mayo is very respected in medical image processing. If interested in learning more, please send me a note indicating your area of interest and a current CV to bje(a)mayo.edu<mailto:bje@mayo.edu>. This is a 100% FTE position with initial appointment of 1 year, but with option of renewal for more years. Mayo Clinic is an equal opportunity employer. Thank you and pardon the interruption if you are not looking for a job.

1 0

[OT, FFT] FFT of signal larger than processor
by OlyDLG April 8, 2014

April 8, 2014

Hi! This Q came up at my work and I said I'd post here to see if anyone had any reference suggestions. We need to FFT a series that uses up to 64K RAM to process, but we're downsizing the processing core to 8K RAM; we have 2 kinds of supplementary flash memory available: "on-chip" (limited supply, faster read/writes) & "off-chip" (ample supply, slower read/writes). Are there known algorithms for FFT-ing a series piecemeal? Such that minimize the number of data exchanges between storage and processor? References? Thanks! OlyDLG

4 4

Segment fault when run svd
by Wenlei Xie April 7, 2014

April 7, 2014

Hi, I got segment fault when I run svd on a huge matrix. The matrix contains n lines, each line is a vector with 2 million dimension vector. I use the following command for SVD: s = linalg.svdvals(mat) It works well when for n = 100, 300, 500. But it got segment fault when n = 1000. Any suggestions? Thank you! Best, Wenlei

2 2

simple optimization but very complex constraints
by Ben Orr April 6, 2014

April 6, 2014

I am trying to find the optimal student loan payoff schedule. I believe the scipy.optimize.minimize function will do what I need but I'm not sure if it's the best solution. I represent the loan payment schedule as a 2D array; each column is a loan, each row is a month. I start off by determining the initial guess (x0). There are 'number-of-loans' columns. The number of rows is the max of the number of months to pay off each loan if making the minimum payments. I then fill in each column with the minimum payment for each loan until it is payed off, followed by zeros. I would like to pay "y" amount each month (must be greater than the sum of the minimum payments). Therefore, I would like to adjust the amount payed each month (aka the sum of each row) to be equal to or less than my target monthly payment such that the total payoff amount (sum of the array) is minimized. Thus, my function to minimize is simply "sum(x)". The trick is all the constraints. Each cell must be greater-than-or-equal-to the minimum payment, unless the remaining balance is less than the minimum payment, in which case the minimum of that cell is the remaining balance (0 if paid off). The maximum for any given cell is the target monthly payment "y". The maximum of the sum of each row is also the target monthly payment "y". Each column must pay-off the loan that column represents. Note that the sum(column)!=initialBalance. "Balance" is only the principal. The following formula calculates the remaining balance after a monthly payment: i = interest rate (ex: 6.7%, i=6.7 not 0.067) b = balance (after a payment) bp = balance previous (before a payment) p = payment b = bp*(1+(i/1200))-((p/(i/1200))*((1+(i/1200))-1)) I hoped to simplify the problem by having an array the same shape as "x" but with the remaining balance (formula above) after each month's payment. This array would be used in for both constraints pertaining to remaining balance (remaining balance < min payment and each column must pay off each loan). See attached spreadsheet for calculating both x0 and the remaining balance. Given all that (sorry for the long description), any help choosing and implementing the appropriate optimization method is greatly appreciated. Thank you. Ben

2 1

scipy.sparse parallelism?
by Christopher Jordan-Squire April 6, 2014

April 6, 2014

Does parallelism in the scipy.sparse solvers only come from multi-threaded BLAS? I.e., the solvers themselves aren't inherently parallel other than through a multi-threaded BLAS? Follow-up question: Numpy/Scipy are using the intel mkl, will scipy.sparse matrices use intel mkl's sparse blas? -Chris

2 1

ANN: IEP 3.4 - the Interactive Editor for Python (now with IPython)
by Almar Klein April 4, 2014

April 4, 2014

Dear all, We are pleased to announce that IEP version 3.4 is now available for download. IEP is a cross-platform Python IDE focused on interactivity and introspection, which makes it very suitable for scientific computing. Its practical design is aimed at simplicity and efficiency. This version is marked by many improvements to the shell. Most notably, we now have *IPython integration*! Further, there are more ways to customize the shell (e.g. define command line and environment args), there is support for coloured text, you can click on filenames in tracebacks to open them, and it just behaves better. Also, many improvements have been made to the debug system. Website: http://iep-project.org Release notes: https://bitbucket.org/iep-project/iep/wiki/Release_notes To install: - Binaries are available for Windows, OSX and Linux. - conda install -c pyzo iep (in a p33 env) - pip install iep (IEP depends on Python3 and PySide/PyQt4.) Happy coding, Almar

2 2