Mailman 3 June 2007 - SciPy-User

object too deep??
by Emanuele Zattin June 29, 2007

June 29, 2007

I have this optimization problem: this function returns the sum of some gaussians given their parameters in arrays: def gaussian(height, center_x, center_y, width): """Returns a gaussian function with the given parameters""" width = float(width) return lambda x,y: sum(height*exp(-(((center_x-x)/width)**2+((center_y-y)/width)**2)/2)) this function tries to fit given a starting image: def fitgaussian(data, obj_x, obj_y, obj_v): """Returns (height, x, y, width) the gaussian parameters of a 2D distribution found by a fit""" #params = moments(data) params = obj_v, obj_x-obj_x[0]+2, obj_y-obj_y[0]+2, ones(len(obj_x)) errorfunction = lambda p: ravel(gaussian(*p)(*indices(data.shape)) - data) p, success = leastsq(errorfunction, params) return p and i use them with: # how many maxima here? max_list = [i] for j in range(len(obj_x)): if obj_x[j] >= x1 and obj_x[j] < x2 and obj_y[j] >= y1 and obj_y[j] < y2 and j != i: max_list.append(j) #for indices in max_list: ml = array(max_list) params = fitgaussian(neigh, obj_x[ml], obj_y[ml], obj_v[ml]) print len(max_list), params but i get an error like: In [9]: run cutoff --------------------------------------------------------------------------- <type 'exceptions.ValueError'> Traceback (most recent call last) /home/emanuelez/Tesi/Code/cutoff.py in <module>() 174 # FIND OBJECTS PROPERTIES 175 # ----------------------- --> 176 get_objects_info(blurred, 2, obj_x, obj_y, obj_v) 177 178 /home/emanuelez/Tesi/Code/cutoff.py in get_objects_info(image, size, obj_x, obj_y, obj_v) 143 #for indices in max_list: 144 ml = array(max_list) --> 145 params = fitgaussian(neigh, obj_x[ml], obj_y[ml], obj_v[ml]) 146 print len(max_list), params 147 /home/emanuelez/Tesi/Code/cutoff.py in fitgaussian(data, obj_x, obj_y, obj_v) 124 params = obj_v, obj_x-obj_x[0]+2, obj_y-obj_y[0]+2, ones(len(obj_x)) 125 errorfunction = lambda p: ravel(gaussian(*p)(*indices(data.shape)) - data) --> 126 p, success = leastsq(errorfunction, params) 127 return p 128 /usr/lib/python2.5/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag) 264 if (maxfev == 0): 265 maxfev = 200*(n+1) --> 266 retval = _minpack._lmdif(func,x0,args,full_output,ftol,xtol,gtol,maxfev,epsfcn,factor,diag) 267 else: 268 if col_deriv: <type 'exceptions.ValueError'>: object too deep for desired array WARNING: Failure executing file: <cutoff.py> What does "object too deep for desired array" mean? I'm really puzzled about this. Thanks for any help or suggestion! Emanuele

2 3

http://www.scipy.org/Installing_SciPy is dead
by Christian K June 29, 2007

June 29, 2007

Hi, I can't access http://www.scipy.org/Installing_SciPy. Could somebody please check that, fix it? Christian

2 1

[scikits] Updated generic optimizer (and egg download link)
by Matthieu Brucher June 29, 2007

June 29, 2007

Hi, I did not improve much of the code, but here is a link to the scikit : http://download.gna.org/pypeline/ I will add other step and perhaps line search modules, as well as helper functions for fitting purpose (least squares or more robust fits), but I can't give you a dead line. The code works very well, I use it for my PhD searches, and it really saves me a lot of time/troubles. I'll add in another scikit a neighboorhood search (Kd-tree in biopython) as well once I've found a solution for the matrix lib I use.. Matthieu P.S. : pypeline is an "old" project similar in many aspects to what David wanted to do for his SoC, but much more general, for those who might be intrigued. It's being developped slowly, but steadily.

4 9

verify whether a matrix is positive definite or not
by Nils Wagner June 29, 2007

June 29, 2007

Hi all, I have a parameter-dependent matrix B(x) = B_0 + x B_1, 0 \le x \le 1 where B_0 and B_1 are symmetric. How can I determine critical values x* (if any) such that B(x*) is not positive definite ? from scipy import * def B(x): return array(([[11.,8.],[8.,7.]])) - x*array(([[20.,1.],[1.,26]])) X = linspace(0,1,100) for x in X: print x L=linalg.cholesky(B(x),lower=1) I mean it would be nice if cholesky could return info=1 if the matrix is not spd. The current behaviour is Traceback (most recent call last): File "test_spd.py", line 11, in ? L=linalg.cholesky(B(x),lower=1) File "/usr/lib64/python2.4/site-packages/scipy/linalg/decomp.py", line 552, in cholesky if info>0: raise LinAlgError, "matrix not positive definite" numpy.linalg.linalg.LinAlgError: matrix not positive definite Helpful suggestions would be appreciated. Nils

3 3

Re: [SciPy-user] [SciPy-dev] question about scipy.optimize.line_search
by Alan G Isaac June 29, 2007

June 29, 2007

On Thu, 28 Jun 2007, Dmitrey apparently wrote: > help(line_search) yields > -------------------------------------------------------------------- > line_search(f, myfprime, xk, pk, gfk, old_fval, old_old_fval, args=(), > c1=0.0001, c2=0.90000000000000002, amax=50) > Find alpha that satisfies strong Wolfe conditions. > Uses the line search algorithm to enforce strong Wolfe conditions > Wright and Nocedal, 'Numerical Optimization', 1999, pg. 59-60 > For the zoom phase it uses an algorithm by > Outputs: (alpha0, gc, fc) > -------------------------------------------------------------------- > So I need to know what are other args, especially gfk (is it a gradient > in point xk?), old_fval, old_old_fval (I guess I know what do c1 & c2 mean) This is certainly lacking documentation! A little is here: http://docs.neuroinf.de/api/scipy/scipy.optimize.optimize-pysrc.html#line_s… Can anyone help Dmitrey more? Thank you, Alan Isaac

3 2

3D density calculation
by Chris Lee June 28, 2007

June 28, 2007

Hi everyone, I was hoping this list could point me in the direction of a more efficient solution to a problem I have. I have 4 vectors: x, y, z, and t that are about 1 million in length that describe the positions of photons. As my simulation progresses it updates the positions so x, y, z, and t change by an unknown (and unknowable) amount every update. This worked very well for its original purpose but now I need to calculate the photon density change over time. Currently after each update, I iterate over time slices, x slices, and y slices and then make an histogram of z which I then stitch together to create a density. However, this becomes very slow as the photons spread out in space and time. Does anyone know how to take such a large vector set and return a density efficiently?

3 7

lambda forms?
by Emanuele Zattin June 28, 2007

June 28, 2007

Hello, i'm trying to define a function to calculate the sum of some gaussian given their paramers. Something like: def gaussian(height, center_x, center_y, width): """Returns a gaussian function with the given parameters""" width = float(width) return lambda x,y: sum(height*exp(-(((center_x-x)/width)**2+((center_y-y)/width)**2)/2)) where all the parameters are arrays of the same length (the length is variable). This version of course does not work, but i wonder if it is still possible using lambda forms or if i have to write the function in an explicit way. Emanuele

2 2

signal.butter questions
by Ryan Krauss June 27, 2007

June 27, 2007

I have two questions about signal.butter. First, I don't see it in the docstring, but I am assuming the Wn must be normalized by dividing by the Nyquist frequency. Is that true? It seems like that is what buttord is saying and I assume they use the same conventions. Second, what is meant by the analog keyword? Thanks, Ryan

2 2

Another leastsq Jacobian bug
by Lin Shao June 26, 2007

June 26, 2007

Hi there, It seems like when calling leastsq(), if the Jacobian matrix is organized as rows of partial derivatives (w.r.t. each variable), then no optimization is done at all -- the return value is the same as the initial guess. It only works when the matrix is columns of derivatives and set the parameter col_deriv to 1. --lin

3 7

How to Enable Delaunay Package
by Lorenzo Isella June 26, 2007

June 26, 2007

Dear All, I am posting this after a discussion originated on the matplotlib mailing list. Fundamentally, I need to plot data on irregular (i.e. non equi-spaced rectangular) grids. I finally was recommended to look at the Delaunay package (see approach2 at the link: http://scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data ). The problem is that the approach: from scipy.sandbox.delaunay import * does not work (the system does not find the requested module). Now, I am running Debian testing on my box and I have Python2.3,2.4,2.5 installed beside SciPy as taken from the standard repositories. Under /usr/lib/python2.4/site-packages/scipy/sandbox I have the file setup.py which I copy and paste at the end of the email. I try uncommenting the line dealing with Delaunay, but that did not help me out (probably it is useful only if I am rebuilding SciPy, which I would like to avoid). Anyone has experienced the same problem or has any suggestions? I am really in need to get this working in order to be able to perform some non-trivial data plotting with matplotlib. Many thanks Lorenzo import os def configuration(parent_package='',top_path=None): from numpy.distutils.misc_util import Configuration config = Configuration('sandbox',parent_package,top_path) sandbox_packages = [] try: sandbox_file = open(os.path.join(config.package_path, 'enabled_packages.txt'), 'rU') except IOError: pass else: for line in sandbox_file: p = line.strip() if line.startswith('#'): continue sandbox_packages.append(p) sandbox_file.close() for p in sandbox_packages: config.add_subpackage(p) # All subpackages should be commented out in the version # committed to the repository. This prevents build problems # for people who are not actively working with these # potentially unstable packages. # You can put a list of modules you want to always enable in the # file 'enabled_packages.txt' in this directory (you'll have to create it). # Since this isn't under version control, it's less likely you'll # check it in and screw other people up :-) # An example package: #config.add_subpackage('exmplpackage') # Monte Carlo package #config.add_subpackage('montecarlo') # PySparse fork with NumPy compatibility #config.add_subpackage('pysparse') # Robert Kern's corner: #config.add_subpackage('rkern') # ODRPACK #config.add_subpackage('odr') # Delaunay triangulation and Natural Neighbor interpolation config.add_subpackage('delaunay') # Gist-based plotting library for X11 #config.add_subpackage('xplt') # elementwise numerical expressions #config.add_subpackage('numexpr') # Statistical models #config.add_subpackage('models') # Adaptation of Scientific.IO (2.4.9) to use NumPy #config.add_subpackage('netcdf') # Finite Difference Formulae package #config.add_subpackage('fdfpack') # Package with useful constants and unit-conversions defined #config.add_subpackage('constants') # Interpolating between sparse samples #config.add_subpackage('buildgrid') # Package for Support Vector Machine #config.add_subpackage('svm') # Package for Gaussian Mixture Models #config.add_subpackage('pyem') # David Cournapeau's corner: autocorrelation, lpc, lpc residual #config.add_subpackage('cdavid') # New spline package (based on scipy.interpolate) #config.add_subpackage('spline') return config if __name__ == '__main__': from numpy.distutils.core import setup setup(**configuration(top_path='').todict())

5 7