Mailman 3 July 2007 - SciPy-Dev

Question about Kaiser Implementation in firwin
by Buehler, Eric (AGRE) 04 Oct '08

04 Oct '08

Hello, I have been using the signal.firwin FIR design tool and I have a question about the implementation. Scipy-0.4.8 Numpy-0.9.6 I am generating a kaiser window with the following parameters: N = 128 Cutoff = ~4e7 width = .3 All of the special kaiser generation parameters in the function follow Oppenheim and Schafer well. Even the sinc function for linear phase is follows exactly. However, the last line of the firwin function is causing me some heartburn. filter_design.py 1538 win = get_window(window,N,fftbins=1) 1539 alpha = N//2 1540 m = numpy.arange(0,N) 1541 h = win*special.sinc(cutoff*(m-alpha)) 1542 return h / sum(h) Line 1542 of filter_design.py, "return h / sum(h)", normalizes the function where it doesn't seem necessary, at least in the kaiser window case. Without the normalization, the kaiser window already returns a value of 1 at the zero frequency point. This normalization scales all of the data, making the window difficult to use in the frequency domain. Can someone point me to the rationale for this line? Looking at the code, this seems to be a pretty recent change (within the last year/year and a half). Thanks, Eric Buehler eric <dot> buehler <at> smiths-aerospace <dot> com ****************************************** The information contained in, or attached to, this e-mail, may contain confidential information and is intended solely for the use of the individual or entity to whom they are addressed and may be subject to legal privilege. If you have received this e-mail in error you should notify the sender immediately by reply e-mail, delete the message from your system and notify your system manager. Please do not copy it for any purpose, or disclose its contents to any other person. The views or opinions presented in this e-mail are solely those of the author and do not necessarily represent those of the company. The recipient should check this e-mail and any attachments for the presence of viruses. The company accepts no liability for any damage caused, directly or indirectly, by any virus transmitted in this email. ******************************************

3 2

Status of sandbox.spline and my other code.
by John Travers 07 Apr '08

07 Apr '08

Hi all, This email is intended to clarify the status of various bits of code I've committed/not committed so I don't leave too much of a mess after me. I'm about to start writing my PhD thesis, so I'll be too busy to work on scipy until September. 1. sandbox.spline stuff The module here was supposed to be a tidy up of scipy.interpolate. All of the dierkx functionality has been moved to f2py wrappers, which I think makes maintenance easier. However, it does not add any new functionality, and in retrospect appears to have been a waste of my time. So I'll leave to you to decide if you want to integrate it into scipy.interpolate or not; though it was originaly planned to be a seperate module to clear up the ambiguity between interpolation and smoothing spline functionality. However I have added 20 unit tests to the code (most of which simply check the wrapper against the pure fortran output, but still useful I think) which could easily be moved over to the current scipy.interpolate module. 2. Radial basis function module in sandbox (rbf) This code has recently had attention from Robert Hetland and is quite improved. I'm not sure if I will ever go into scipy though?? If not, I will put it into a scikit at some point. Related to this, the wiki page has been updated, however, it is at http://www.scipy.org/RadialBasisFunctions but should probably be in the Cookbook. I don't know how to move it over, so some pointers would be helpful. In addition there are four unused attachments I uploaded that could be deleted. 3. Boundary value ODE solver (BVP_SOLVER) I have fairly functional code interfacing to some very good BVP solver code (http://cs.smu.ca/~muir/BVP_SOLVER_Webpage.shtml). The license is good for scipy (as confirmed by the authors), but the code if fortran 95 - therefore cannot apparently go into scipy. I'll make a scikit when I get time (probably September). 4. An interface to the pikaia genetic algorithm routine I have a python interface to pikaia fortran code for genetic algorithm optimisation of real valued functions (http://www.hao.ucar.edu/Public/models/pikaia/pikaia.html). The fortran code itself is robust and highly optimised for numerical work, though it is nowhere near as general as the sandbox.ga module. I think it is useful to go in scipy.optimization alongside the annealing module. If people disagree I'll make it into a scikit. I hope this clears a few things up. Best regards, John Travers

4 4

ZVODE, a Double Precision Complex ODE Solver
by Nils Wagner 20 Feb '08

20 Feb '08

This might be of interest. Nils From: Alan Hindmarsh <alanh(a)llnl.gov> Date: Wed, 6 Dec 2006 18:20:32 -0500 Subject: ZVODE, a Double Precision Complex ODE Solver Of the many solvers freely available for ODE initial value problems, there seem to be very few, if any, for the case of double precision complex dependent variables, even though there have been occasional user calls for this functionality. (Of course a complex system can be treated as a real system of twice the size, but direct treatment as a complex system is preferable.) I have provided such a solver: ZVODE. ZVODE is a straightforward modification of the Fortran ODE solver DVODE, for double precision real systems, making it part of a trio, SVODE/DVODE/ZVODE. ZVODE is available in two locations: (1) Netlib, in /ode (see file zvode.f), and (2) the LLNL CASC software site, www.llnl.gov/CASC/software.html (see VODE). The latter site also includes a demonstration program, zvdemo. Alan Hindmarsh Center for Applied Scientific Computing Lawrence Livermore National Laboratory

4 5

ndimage problems
by Christopher Hanley 10 Oct '07

10 Oct '07

Hi, We have found two problems with ndimage. I have filed a ticket #455 on the scipy trac page. The first problem can be seen with this example: > import numpy as n > from scipy import ndimage as nd > a = n.ones((10,5),dtype=n.float32) * 12.3 > x = nd.rotate(a,90.0) > x Out[17]: array([[ 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 0. , 0. ], [ 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019], [ 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019], [ 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019], [ 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019, 12.30000019]], dtype=float32) }}} Notice that the last two entries of the first row are now 0. The second problem has to do with the reversing of byte order if you have big-endian data on a little endian machine. Please see the example below: >>> a = N.ones ((2,3), dtype=N.float32) * 12.3 >>> a = a.byteswap() >>> a.dtype = a.dtype.newbyteorder (">") >>> print a [[ 12.30000019 12.30000019 12.30000019] [ 12.30000019 12.30000019 12.30000019]] >>> print ndimage.rotate (a, 90.) [[ 0.00000000e+00 -4.28378144e+08] [ 0.00000000e+00 -4.28378144e+08] [ -4.28378144e+08 -4.28378144e+08]] I have taken a look at the ndimage python code and cannot find any explicit calls to byteswap. I'm guessing something is funny in one of the c-api calls. I haven't been able to track it down yet. Chris

2 3

ticket 285 (changes in optimize.brent)
by dmitrey 31 Jul '07

31 Jul '07

hi all, rev. 3209 contains changes to optimize.brent (I implemented the brent class in the way proposed by Alan Isaac + 4 tests related to brent were added). So I guess the ticket 285 should be closed. Also, I have changed some docstrings related to brent and some other funcs (): If bracket is two numbers *(a,c)* then they are assumed to be a starting interval for a downhill bracket search (see bracket);* it doesn't always mean that obtained solution will satisfy a<=x<=c*. (As for me some weeks ago I was surprised to obtain solution from outside of the (a,c) interval, MATLAB has only one func fminbound that yields solution strictly from the given interval, but some scipy.optimize line-search routines use strictly the interval and some other, like brent, use only as starting interval for bracket search, and this is not properly descriebed in documentation). Regards, D.

1 0

ticket 416 (gradient in leastsq): some questions
by dmitrey 31 Jul '07

31 Jul '07

Hi all, the 416th ticket was ascribed to me. There some problems: 1. What's the difference between these 2 funcs from __minpack.h: int jac_multipack_calling_function(int *n, double *x, double *fvec, double *fjac, int *ldfjac, int *iflag) int jac_multipack_lm_function(int *m, int *n, double *x, double *fvec, double *fjac, int *ldfjac, int *iflag) They have same description. /* This is the function called from the Fortran code it should -- use call_python_function to get a multiarrayobject result -- check for errors and return -1 if any -- otherwise place result of calculation in *fvec or *fjac. If iflag = 1 this should compute the function. If iflag = 2 this should compute the jacobian (derivative matrix) */ 2. So patch assigned to the ticket proposes to rewrite the line 152 MATRIXC2F(fjac, result_array->data, *n, *ldfjac) as MATRIXC2F(fjac, result_array->data, *ldfjac, *n) however, line 92 is the same. so maybe it needs same patch. 3. The MATRIXC2F is defined in minpack.h w/o any description: #define MATRIXC2F(jac,data,n,m) {double *p1=(double *)(jac), *p2, *p3=(double *)(data);\ int i,j;\ for (j=0;j<(m);p3++,j++) \ for (p2=p3,i=0;i<(n);p2+=(m),i++,p1++) \ *p1 = *p2; } I have no idea what does it do. So I replaced (jac,data,n,m) by (jac,data,m,n), and user's example works correctly for all cases 1-3: 1) w/o gradient info 2) with gradient info, col_deriv=0 3) with gradient info, col_deriv=1 (in the case I modified the user's gradient func so that it returns transposed gradient) scipy.test(1) also didn't yield any bugs related to leastsq. Do you agree to submit the changes to svn? Regards, D.

4 6

wiki page HowToDocument is broken
by dmitrey 30 Jul '07

30 Jul '07

hi all, the page http://projects.scipy.org/scipy/numpy/wiki/HowToDocument that is referred from http://projects.scipy.org/scipy/numpy/wiki/DocstringStandards need bugfix ("A problem occurred in a Python script..."), full output is attached below: *SubversionException* Python 2.4.3: /usr/bin/python Mon Jul 30 01:11:01 2007 A problem occurred in a Python script. Here is the sequence of function calls leading up to the error, in the order they occurred. /var/www/cgi-bin/build/bdist.linux-i686/egg/tracrst/macro.py <file:///var/www/cgi-bin/build/bdist.linux-i686/egg/tracrst/macro.py> in *render_macro*(self=<tracrst.macro.TracReSTMacro object>, req=<trac.web.api.Request object>, name=u'ReST', args=u'trunk/numpy/doc/HOWTO_DOCUMENT.txt') /var/www/cgi-bin/build/bdist.linux-i686/egg/tracrst/svn_helper.py <file:///var/www/cgi-bin/build/bdist.linux-i686/egg/tracrst/svn_helper.py> in *cat*(self=<tracrst.svn_helper.SVNHelper instance>, path='file:///home/scipy/svn/numpytrunk/numpy/doc/HOWTO_DOCUMENT.txt', rev=<libsvn.core.svn_opt_revision_t; proxy of C svn_...ion_t instance at _909de509_p_svn_opt_revision_t>) *SubversionException*: ('Unable to open an ra_local session to URL', 180001) apr_err = 180001 args = ('Unable to open an ra_local session to URL', 180001)

4 3

Launchpad hosting: an option for scipy?
by Fernando Perez 30 Jul '07

30 Jul '07

Hi all, I just wanted to mention this I ran across: http://arstechnica.com/news.ars/post/20070724-ars-at-ubuntu-live-new-launch… It sounds like it could provide an interesting, low-overhead solution for hosting .debs. I know Andrew (at least) has hosted his privately, but this appears to lower the logistical bar for others to do the same without having to know quite as much about deb packaging as Andrew does. Just a thought. f

3 2

GSoC weekly report + need discuss of ticket 285
by dmitrey 28 Jul '07

28 Jul '07

hi all, so http://projects.scipy.org/scipy/scipy/ticket/464 (optimize.fmin_powell doesn't accept a matrix input for the initial guess) rev 3199 http://projects.scipy.org/scipy/scipy/ticket/416 (MATRIXC2F transposing the wrong way in optimize.leastsq) I will commit my patch if nothing better will be proposed till some hours One more change in tnc - now return x value (optim point) is numpy.array, not Python list. Also, now it consumes x0 as numpy.array, not Python list (so http://projects.scipy.org/scipy/scipy/ticket/384 should be closed) About a day was spent for a bug related to fmin_ncg, but finally it turned out to be related to sparse matrices. So the last ticket left is bracket parameters (http://projects.scipy.org/scipy/scipy/ticket/285) Alan Isaac proposed the new format of brent (very preliminary, of course), see below So before implementing something like that I want to hear your opinions, is this the best way or something should be implemented in other way. As for me, I treat the benefits from these changes to be very suspicious. Regards, D. def brent2(func, args=(), brack=None, tol=1.48e-8, full_output=0, maxiter=500): """ copy brent docstring here """ brent = Brent(func=func, args=args, tol=tol, maxiter=maxiter) brent.set_bracket(brack=brack) brent.optimize() return brent.get_result(full_output=full_output) class Brent: #need to rethink design of __init__ def __init__(func, args=(), tol=1.48e-8, maxiter=500): self.func = func self.args = args self.tol = tol self.maxiter = maxiter self._mintol = 1.0e-11 self._cg = 0.3819660 self.xmin = None self.fval = None self.iter = 0 self.funcalls = 0 #etc...... #need to rethink design of set_bracket (new options, etc) def set_bracket(self, brack = None): self.brack = brack def get_bracket_info(self): #set up func = self.func args = self.args brack = self.brack ### BEGIN core bracket_info code ### ### carefully DOCUMENT any CHANGES in core ## if brack is None: xa,xb,xc,fa,fb,fc,funcalls = bracket(func, args=args) elif len(brack) == 2: xa,xb,xc,fa,fb,fc,funcalls = bracket(func, xa=brack[0], xb=brack[1], args=args) elif len(brack) == 3: xa,xb,xc = brack if (xa > xc): # swap so xa < xc can be assumed dum = xa; xa=xc; xc=dum assert ((xa < xb) and (xb < xc)), "Not a bracketing interval." fa = func(*((xa,)+args)) fb = func(*((xb,)+args)) fc = func(*((xc,)+args)) assert ((fb<fa) and (fb < fc)), "Not a bracketing interval." funcalls = 3 else: raise ValueError, "Bracketing interval must be length 2 or 3 sequence." ### END core bracket_info code ### return xa,xb,xc,fa,fb,fc,funcalls def optimize(self): #set up for optimization func = self.func xa,xb,xc,fa,fb,fc,funcalls = get_bracket_info(brack = self.brack, args = self.args) _mintol = self._mintol _cg = self._cg ################################# #BEGIN CORE ALGORITHM #we are making NO CHANGES in this ################################# x=w=v=xb fw=fv=fx=func(*((x,)+args)) if (xa < xc): a = xa; b = xc else: a = xc; b = xa deltax= 0.0 funcalls = 1 iter = 0 while (iter < maxiter): tol1 = tol*abs(x) + _mintol tol2 = 2.0*tol1 xmid = 0.5*(a+b) if abs(x-xmid) < (tol2-0.5*(b-a)): # check for convergence xmin=x; fval=fx break if (abs(deltax) <= tol1): if (x>=xmid): deltax=a-x # do a golden section step else: deltax=b-x rat = _cg*deltax else: # do a parabolic step tmp1 = (x-w)*(fx-fv) tmp2 = (x-v)*(fx-fw) p = (x-v)*tmp2 - (x-w)*tmp1; tmp2 = 2.0*(tmp2-tmp1) if (tmp2 > 0.0): p = -p tmp2 = abs(tmp2) dx_temp = deltax deltax= rat # check parabolic fit if ((p > tmp2*(a-x)) and (p < tmp2*(b-x)) and (abs(p) < abs(0.5*tmp2*dx_temp))): rat = p*1.0/tmp2 # if parabolic step is useful. u = x + rat if ((u-a) < tol2 or (b-u) < tol2): if xmid-x >= 0: rat = tol1 else: rat = -tol1 else: if (x>=xmid): deltax=a-x # if it's not do a golden section step else: deltax=b-x rat = _cg*deltax if (abs(rat) < tol1): # update by at least tol1 if rat >= 0: u = x + tol1 else: u = x - tol1 else: u = x + rat fu = func(*((u,)+args)) # calculate new output value funcalls += 1 if (fu > fx): # if it's bigger than current if (u<x): a=u else: b=u if (fu<=fw) or (w==x): v=w; w=u; fv=fw; fw=fu elif (fu<=fv) or (v==x) or (v==w): v=u; fv=fu else: if (u >= x): a = x else: b = x v=w; w=x; x=u fv=fw; fw=fx; fx=fu iter += 1 ################################# #END CORE ALGORITHM ################################# self.xmin = x self.fval = fx self.iter = iter self.funcalls = funcalls def get_result(self, full_output=False): if full_output: return xmin, fval, iter, funcalls else: return xmin

1 0

connecting tnc 1.3
by dmitrey 28 Jul '07

28 Jul '07

hi all, So now I'm trying to make all the tests for tnc 1.3 running ok. But there is a problem encountered: here's a code from tnc.py, lines 210-213: for i in range(n): l,u = bounds[i] if l is None: low[i] = -HUGE_VAL So if bounds are for example ([-inf, -1.5], None), as it is written in test1fg(x), it yields error "None object is not iterable". (because it tries to get lb, ub = None ) Do you think it's a sort of bug that should be fixed in svn or not? Regards, D.

4 16