Mailman 3 April 2002 - SciPy-Dev

RE: [SciPy-dev] RE: Remez algorithm
by Chuck Harris April 26, 2002

April 26, 2002

Hi Travis > -----Original Message----- > From: Travis Oliphant [mailto:oliphant.travis@ieee.org] > Sent: Tuesday, April 23, 2002 9:22 PM > To: scipy-dev(a)scipy.org > Subject: Re: [SciPy-dev] RE: Remez algorithm > > Have you seen the remez algorithm in the signal module? > > How does this one compare? OK, I wrote a bandpass filter driver for my routine and spent some time going through the sigtool version and comparing it to mine for somewhat realistic filter specs : taps = 90 grid_density = 16 bands = [.011,.15,.175,.239,.261,.500] desired = [0.,1.,0.] weights = [1e3,1.0,1e3] Some remarks and results: 1) the sigtools version requires a density of 32 to get accurate coefficients, I suspect an off by one bug in its peak finding routine. In any case, I have used 32 for the sigtools version and 16 for mine to get comparable results. 2) the sigtools version is significantly faster, the times for 100 iterations : sigtools -- 1.5 sec mine -- 13 sec 3) the sigtools version will use less memory, in this problem : sigtools -- about 11KiB mine -- about 259 KiB 4) the sigtools version does not return useful ancillary data, for instance, it's usefull to know just how good the fit was so that a routine can use it to make tradeoffs : number of taps, out of band rejection vs inband ripple, etc. 5) the sigtools version restricts the band edges to the range 0--.5, which makes it difficult, sometimes impossible, to design filters with complex coefficients, a practical thing to do sometimes with rad hard parts running $20K a pop. Mine goes out to 1.0. 6) the sigtools interface does not expose the function hooks available in the Parks-McClellan routine. This makes it difficult to properly design cascaded filters, or to design filters with more interesting shapes than the three basic types available. 7) the code for the sigtools version would be heck to modify or use as the basis of some other routine. It's convergence criteria are not clear; it failed to converge on some problems, but I have no idea what this meant. 8) The sigtools version is limited to filter design. Overall, I think that the sigtools version would work well for 99% of everyday filter design tasks, it's the remaining 1% that I wish to see covered, at least for my own uses. Chuck

1 0

BVP code
by Anders Lyckegaard April 26, 2002

April 26, 2002

Dear scipy-developers, This is my first post to this list. I have been working on a bit of code that implements a shooting routine for two point boundary value problems. It uses the SciPy routines odeint and fsolve, so I was wondering if the shooting routine would be of interest to the SciPy community. I am completely new to SciPy, and Scientific Computing with Python, so don't expect too much. I have only used SciPy for a few days now. I would also like to say that I think that SciPy is a great tool. It is nice to a somewhat larger collection of high quality modules for Scientific computing instead of all the small modules that are floating around the internet. Best regards Anders Lyckegaard -- Anders Lyckegaard, Ph.D student Institute of Mechanical Engineering, Aalborg University Pontoppidanstraede 101, DK-9220 Aalborg Oest, Denmark Phone: +45 9635 8080 (direct: +45 9635 9325), fax: +45 9815 1675 Email: al(a)ime.auc.dk

3 3

RE: [SciPy-dev] RE: Remez algorithm
by Chuck Harris April 24, 2002

April 24, 2002

> Have you seen the remez algorithm in the signal module? > > How does this one compare? > > -Travis > This looks to be a pretty close implementation of the original Parks-McClellan algorithm. As such it directly applies to real symmetric filters. The idea is that in the frequency domain the transfer function is a sum of terms of the form cos(2*pi*n*f) which may be rearranged into powers of n, cos(2*pi*f)**n, so that a change of variable to cos(2*pi*f) yields a Vandermonde matrix and Lagrange interpolation can be used to evaluate sums of the form a + b*cos(2*pi*1*f) + c*cos(2*pi*2*f) + ... passing thru a specified set of points without ever having to solve for a,b,c, etc. At the least, I feel that this should be generalized so as to accept transfer functions like sum( a[n]*exp(2*pi*1j*n*f)), which are also Vandermonde matrices after the change of variable x <- exp(2*pi*1j*f), that is, it would be nice to use Lagrange interpolation for points specified on the unit circle in the complex plane, which is why I sent along nevilles version of Lagrange interpolation, which should perhaps be redone in C for speed. So my first comment is that the routine should be adapted to use a general, good form of Lagrange interpolation which should be included in scipy. The routine I sent *does* solve for the coefficients directly, and then evaluates the sum by matrix multiplication. It is of slightly lower order than using Lagrange interpolation when the number of points to be interpolated at exceeds the number of coefficients, which is always, and is more general, in that it is not necessary to have matrices that could be reduced to Vandermonde form in principal. There is a generalization of Lagrange interpolation due to Muhlbach that applies to complete Chebychev systems, and such might be an interesting topic to pursue. In any case, the linear algebra routines in SciPy are highly optimized and are much faster than Lagrange interpolation in its present form. The potential drawback is that the matrices appearing in the solution may not be well conditioned, so some good sense is required on the part of the user. One could argue that the badly conditioned case reflects a certain indeterminacy in the coefficients that one wants to use, and as such reflects reality. On the other hand, Lagrange interpolation will produce pretty accurate interpolations whatever the users choice of basis, and uses less space as no matrices have to be maintained. The reason I wrote this routine is that Remez package available in MatLab was also too restrictive and I needed something better. Chuck

1 0

RE: Remez algorithm
by Chuck Harris April 24, 2002

April 24, 2002

I have attached a python version of the Remez algorithm for use in filter design. This version is for review and comment. I need to do some testing of the complex case and also include some examples. This routine has been used successfully to design complex decimating FIR filters. Any comments relative to other things are welcome. Chuck

2 1

RE: [SciPy-dev] BVP code
by Chuck Harris April 24, 2002

April 24, 2002

> -----Original Message----- > From: Jochen Küpper [mailto:jochen@unc.edu] > Sent: Tuesday, April 23, 2002 4:13 PM > To: scipy-dev(a)scipy.org > Subject: Re: [SciPy-dev] BVP code > > > On Tue, 23 Apr 2002 15:18:25 -0600 Chuck Harris wrote: > > >> ,---- > >> | @Article{11901, > >> | author = {J. W. Cooley}, > >> | title = {An Improved Eigenvalue Corrector Formula for > Solving the > >> | Schrödinger Equation for Central Fields}, > >> | journal = {Math. Comput.}, > >> | volume = {15}, > >> | pages = {363}, > >> | year = {1961}, > >> | } > >> `---- > > Chuck> Is that the one that shoots from the ends and meets in the > Chuck> middle and tries to match the slopes? > > Yep > IIRC they compute the derivative of the slope-difference with respect to the energy by integrating it along to the meeting point. There *is* an application for Newton's method with computed derivatives! It would be interesting to compare this with one of the good secant methods, brent for instance. Chuck

2 1

RE: [SciPy-dev] BVP code
by Chuck Harris April 23, 2002

April 23, 2002

> ,---- > | @Article{11901, > | author = {J. W. Cooley}, > | title = {An Improved Eigenvalue Corrector Formula for Solving the > | Schrödinger Equation for Central Fields}, > | journal = {Math. Comput.}, > | volume = {15}, > | pages = {363}, > | year = {1961}, > | } > `---- Is that the one that shoots from the ends and meets in the middle and tries to match the slopes? I believe there is also a version that includes the electron spins, so the ODE is vector instead of scalar. With the great linear algebra stuff we have, some of these are neatly solved by projecting on a finite basis set --- say sines or chebychev polynomials --- and diagonalizing the resulting matrix. Chuck

2 1

RE: [SciPy-dev] BVP code
by Chuck Harris April 23, 2002

April 23, 2002

Hi, > -----Original Message----- > From: eric [mailto:eric@scipy.org] > Sent: Tuesday, April 23, 2002 12:04 PM > To: scipy-dev(a)scipy.org > Subject: Re: [SciPy-dev] BVP code > > > Hey Anders, > > Welcome. I'm glad SciPy is proving useful to you. > > I'm not familiar with shooting routines. Can you give a > little background on > what it is and perhaps a reference. SciPy will ideally > evolve into be a > repository of "best in class" algorithms that are peer > reviewed. If your > algorithm is a general solver that fits this criteria, then > it is likely of > interest. Please send a few more details. > > thanks, > eric > I'm not in this loop, but, shooting routines solve boundary value problems for ODE's. For example, vibrating string problems where separation of variables leads to a second order ODE with one parameter which must satisfy a boundary condition at each end of the string. Satisfy the boundary condition at one end, and then vary the parameter until the integrated ODE satisfies the other boundary condition. Integrating the ODE provides the 'shooting'. Chebychev methods can also provide quick and simple solutions to these types of problems and perhaps it would be nice to include a selection of these. Chuck

1 0

Backwards compatibility for weave
by David O'Gwynn April 23, 2002

April 23, 2002

In trying to use weave on the system where I work, I found that it will not work with Python2.0. Namely because of the use of the sys._getframe method. There might be other compatibility issues as well. Is there any kind of workaround for backwards compatibility? _____________ David O'Gwynn dog(a)erc.msstate.edu

3 2

corrupted line ends
by Heiko Henkelmann April 19, 2002

April 19, 2002

Hello there, all files in scipy\optimize\Zeros have corrupted line ends. See for example zeros.h: 00000000: 2f2a 2057 7269 7474 656e 2062 7920 4368 /* Written by Ch 00000010: 6172 6c65 7320 4861 7272 6973 2063 6861 arles Harris cha 00000020: 726c 6573 2e68 6172 7269 7340 7364 6c2e rles.harris@sdl. 00000030: 7573 752e 6564 7520 2a2f 0d0d 0a0d 0d0a usu.edu */...... 00000040: 2f2a 204d 6f64 6966 6965 6420 746f 206e /* Modified to n 00000050: 6f74 2064 6570 656e 6420 6f6e 2050 7974 ot depend on Pyt 00000060: 686f 6e20 6576 6572 7977 6865 7265 2062 hon everywhere b 00000070: 7920 5472 6176 6973 204f 6c69 7068 616e y Travis Oliphan 00000080: 742e 0d0d 0a20 2a2f 0d0d 0a0d 0d0a 0d0d t.... */........ 00000090: 0a23 6966 6e64 6566 205a 4552 4f53 5f48 .#ifndef ZEROS_H 000000a0: 0d0d 0a23 6465 6669 6e65 205a 4552 4f53 ...#define ZEROS 000000b0: 5f48 0d0d 0a0d 0d0a 2364 6566 696e 6520 _H......#define 000000c0: 5a45 524f 535f 5041 5241 4d5f 4845 4144 ZEROS_PARAM_HEAD I'm checkin out with WinCVS on a Windows ME Box. All the files in the other directories seem to be OK. Heiko

1 0

import scipy fails on linux
by hoel＠germanlloyd.org April 18, 2002

April 18, 2002

Hello, I tried to install the latest scipy snapshot on a linux box. compilation went fine (other than on Solaris), but testing things failed: >python Python 2.2.1 (#1, Apr 18 2002, 17:26:54) [GCC 3.0.4] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import scipy exceptions.ImportError: /usr/local/fitools/py/linux/lib/python2.2/site-packages/scipy/linalg/flapack.so: undefined symbol: sgesdd_ exceptions.ImportError: /usr/local/fitools/py/linux/lib/python2.2/site-packages/scipy/linalg/_flinalg.so: undefined symbol: dlaswp_ Traceback (most recent call last): File "<stdin>", line 1, in ? File "/usr/local/fitools/py/linux/lib/python2.2/site-packages/scipy/__init__.py", line 42, in ? import special, io, linalg, stats, fftpack File "/usr/local/fitools/py/linux/lib/python2.2/site-packages/scipy/special/__init__.py", line 325, in ? import orthogonal File "/usr/local/fitools/py/linux/lib/python2.2/site-packages/scipy/special/orthogonal.py", line 59, in ? from scipy.linalg import eig File "/usr/local/fitools/py/linux/lib/python2.2/site-packages/scipy/linalg/__init__.py", line 40, in ? from basic import * File "/usr/local/fitools/py/linux/lib/python2.2/site-packages/scipy/linalg/basic.py", line 17, in ? import calc_lwork ImportError: /usr/local/fitools/py/linux/lib/python2.2/site-packages/scipy/linalg/calc_lwork.so: undefined symbol: ieeeck_ >>> Greetings Berthold -- Dipl.-Ing. Berthold Höllmann __ Address: hoel(a)germanlloyd.org G / \ L Germanischer Lloyd phone: +49-40-36149-7374 -+----+- Vorsetzen 32/35 P.O.Box 111606 fax : +49-40-36149-7320 \__/ D-20459 Hamburg D-20416 Hamburg

3 3