Mailman 3 February 2010 - SciPy-User

Lagrange Multipliers in optimize.slsqp
by Guilherme P. de Freitas 10 Mar '10

10 Mar '10

Hi all, Would it be possible for optimize.slsqp to return the Lagrange multipliers of the constraints? I don't know Fortran, but it seems that they are available in the original Fortran code: http://projects.scipy.org/scipy/browser/trunk/scipy/optimize/slsqp/slsqp_op… Thanks, Guilherme -- Guilherme P. de Freitas http://www.gpfreitas.com

4 7

Signal to noise ratio
by Nils Wagner 04 Mar '10

04 Mar '10

Hi all, I have two questions concerning signal processing I have used scipy.stats.signaltonoise to compute the signal-to-noise ratio. The value is 0.0447. How can I judge it ? How can I filter out high frequencies using scipy ? How can I eliminate noise from the signal ? Nils

5 13

scipy.interpolate.rbf sensitive to input noise ?
by denis 03 Mar '10

03 Mar '10

A followup to the 15feb thread "creating a 3D surface plot from collected data": the doc for scipy.interpolate.rbf says that 3 of the 7 radial functions use a scale paramater `epsilon`, 4 do not. Odd -- how can some be scale-free ? Running rbf on 100 1d random.uniform input points (after changing the linalg.solve in rbf.py to lstsq) shows that all but linear and gaussian are noisy, giving interpolants way outside the input range: N: 100 ngrid: 50 sigma: 0.1 # min max max |delta| Rbf -1.0 1.0 0.5 gaussian -1.0 1.2 0.3 linear -1.4 2.1 2.5 thin-plate -1.0 2.2 3.0 inverse multiquadric -1.0 2.4 3.2 multiquadric -2.1 12.8 7.4 quintic -10.6 7.0 11.7 cubic Should Rbf be moved off to noisymethods/... until experts can look it over, or have I done something stupid ? (A back-of-the-envelope method for choosing epsilon aka r0 would be nice, and a rough noise amplification model would be nice too.) Also, the doc should have a link to http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data and should say that the initial rbf() is N^2. cheers -- denis <pre><code> # from http://scipy.org/Cookbook/RadialBasisFunctions import sys import numpy as np from rbf import Rbf # lstsq # from scipy.interpolate import Rbf Rbfuncs = ( 'gaussian', 'linear', 'thin-plate', 'inverse multiquadric', 'multiquadric', 'quintic', 'cubic', ) N = 100 ngrid = 50 sigma = .1 exec "\n".join( sys.argv[1:] ) # N= np.random.seed(1) np.set_printoptions( 2, threshold=50, edgeitems=5, suppress=True ) # .2f x = np.sort( np.random.uniform( 0, 10, N+1 )) y = np.sin(x) + np.random.normal( 0, sigma, x.shape ) xi = np.linspace( 0, 10, ngrid+1 ) print "N: %d ngrid: %d sigma: %.2g" % (N, ngrid, sigma) print "# min max max |delta| Rbf" for func in Rbfuncs: rbf = Rbf( x, y, function=func ) fi = rbf(xi) maxdiff = np.amax( np.abs( np.diff( fi ))) print "%5.1f %4.1f %4.1f %s" % ( fi.min(), fi.max(), maxdiff, func ) </pre></code>

3 16

frequency components of a signal buried in a noisy time domain signal
by Nils Wagner 27 Feb '10

27 Feb '10

Hi all, A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal. I found a Matlab template at http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml Matlab has a function nextpow2(L) Is there a similar build-in function in numpy/scipy ? I tried to convert the m-file into a pythonic form. What is needed to obtain a similar figure of the single-sided amplitude spectrum using numpy/scipy/matplotlib ? from numpy import sin, linspace, pi from numpy.random import randn from pylab import plot, show, title, xlabel, ylabel from scipy.fft import fft Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = arange(0,L)*T x = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t) y = x + 2*randn(len(t)) plot(Fs*t[:50],y[:50]) title('Signal corrupted with zero-mean random noise') xlabel('Time (milliseconds)') # # # #NFFT = 2^nextpow2(L); # Next power of 2 from length of y Y = fft(y,NFFT)/L f = Fs/2*linspace(0,1,NFFT/2+1) plot(f,2*abs(Y[:NFFT/2+1])) title('Single-sided amplitude spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') show() What can be done in case of nonequispaced data ? http://dx.doi.org/10.1137/0914081 Thanks in advance Nils

8 34

Difference between SQP and SLSQP
by jagan prabhu 27 Feb '10

27 Feb '10

Hi, I have a basic problem in understanding SLSQP method, what is the difference between SQP (sequential quadratic programming) and SLSQP(Sequential Least SQuares Programming)? please can any one explain the working of SLSQP. Your Mail works best with the New Yahoo Optimized IE8. Get it NOW! http://downloads.yahoo.com/in/internetexplorer/

2 1

Pixel Offset example no longer working
by mikey 26 Feb '10

26 Feb '10

Hi there, Just wondering if anyone knows how to do a pixel offset in matplotlib? I've tried running the example in http://www.scipy.org/Cookbook/Matplotlib/Transformations#line-58 but it appears to have been trashed by changes to the api. This is a very useful function if you're doing discrete data plots like I am and you want uniform clusters of points. Thanks for your time. Regards, Mikey

2 2

Moving median with O(log n) updates
by Wes McKinney 25 Feb '10

25 Feb '10

At Pycon this past weekend Raymond Hettinger gave a nice talk about leveraging Python's built-in data structures to implement various algorithms in a clear and elegant way. http://us.pycon.org/2010/conference/schedule/event/86/ video: http://pycon.blip.tv/file/3261295/ He posted a full implementation of a the moving median using a skiplist here (MIT license): http://code.activestate.com/recipes/576930-efficient-running-median-using-a… Some time ago for pandas I had implemented moving median is a fairly naive way (O(n * window)) using the linear-time wirth method (quickselect would provide equivalent performance). I decided to translate the implementation above into Cython and tailor it for ndarray input and float data. In most cases it achieves 5-10x speedup over the pure Python O(n * log(window)) implementation, and clearly trouncing other implementations. Here it is: http://code.google.com/p/pandas/source/browse/trunk/pandas/lib/src/skiplist… (the actual rolling_median is in moments.pyx) and some timings: In [3]: import pandas In [4]: arr = randn(10000) In [5]: %timeit pandas.rolling_median(arr, 1000) 1 loops, best of 3: 401 ms per loop of course this is roughly linear in the window size: In [11]: %time pandas.rolling_median(arr, 5000) CPU times: user 1.07 s, sys: 0.00 s, total: 1.07 s Wall time: 1.06 s new version (not yet ready for primetime with NaN-handling etc.) In [6]: import pandas.lib.tseries as ts In [7]: %timeit ts.rolling_median(arr, 1000) 10 loops, best of 3: 120 ms per loop Larger window size in this case has minimal effect: In [13]: %timeit ts.rolling_median(arr, 5000) 10 loops, best of 3: 120 ms per loop I'm pretty sure scikits.timeseries is doing something O(n * window^2) In [8]: import scikits.timeseries.lib as l In [9]: %time l.mov_median(arr, 1000) CPU times: user 86.48 s, sys: 0.00 s, total: 86.48 s Wall time: 86.48 s I had a fun time with Cython making the implementation fairly good, but I have the feeling that there's too much Python list access (which is faster than accessing an object ndarray using the buffer interface I found). If someone with better C / Cython skills wants to take a crack at speeding it up even more please be my guest. - Wes

1 0

typo in odeint docstring
by Tony S Yu 25 Feb '10

25 Feb '10

I think the shape of the solutions returned by scipy.integrate.odeint is incorrectly documented: the number of rows matches the number of time steps (in the docstring, N_rows = len(y0)) Example: >>> from scipy import integrate >>> import numpy as np >>> >>> x0 = np.ones(5) >>> dxdt = lambda x, t: x + 1 >>> times = np.linspace(0, 1, 10) >>> print integrate.odeint(dxdt, x0, times).shape (10, 5) Cheers, -Tony PS: I know there's some sort of doc editor for scipy. Should I be making the change there instead of emailing the list? Patch: Index: scipy/integrate/odepack.py =================================================================== --- scipy/integrate/odepack.py (revision 6250) +++ scipy/integrate/odepack.py (working copy) @@ -57,7 +57,7 @@ Returns ------- - y : array, shape (len(y0), len(t)) + y : array, shape (len(t), len(y0)) Array containing the value of y for each desired time in t, with the initial value y0 in the first row.

3 3

Vector comparison
by eric＠depagne.org 25 Feb '10

25 Feb '10

Hi all, I have two sets of vectors, one made of 3-dim vectors and the other one of 4- dims vectors. I want to find for each of my 3D vectors (let's call them the input vector), which one, among my 4D vectors (let's call them my models), has the 3 first coordinates that are the closest of my input vector, so that i can assign to my input vector the 4th dimension taken from the models. I am currently maximizing the cosine of two vectors, but I was wondering if scipy had some more clever routines, that could for instance take into account measurement errors for my vectors. Thanks a lot. Éric. -- Un clavier azerty en vaut deux ----------------------------------------------------------

3 2

stats.distributions skew, kurtosis bugs
by josef.pktd＠gmail.com 25 Feb '10

25 Feb '10

I found an old script of mine to check the skew and kurtosis of the distributions in scipy.stats with bootstrapping, and I started to produce a nicer printout. The attached textfile shows all distributions where the stats() method gives results that deviate by at least by 3 bootstrap standard deviations from the bootstrap mean. This can be because the moment is infinite or not defined, or because there is a bug. There are bugs in some skew and kurtosis for sure, but I don't know which. Also, this is a random test, so some false results might show up because of the random sampling. Bootstrap were generated by drawing an initial sample with rvs(size=1000) from the distribution and then run a bootstrap sampling until the bootstrap standard deviation settles down or a maximum of 30000 samples has been reached. It's again a warning not to rely on higher moment in stats.distributions without checking, and if anyone finds an incorrect calculation or can verify false positives, then I would appreciate any help. Josef

1 0