Mailman 3 January 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

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

Re: [SciPy-dev] Callback feature for iterative solvers
by Nils Wagner 04 May '07

04 May '07

Ed Schofield wrote: > > On 1/24/07, *Nils Wagner* <nwagner(a)iam.uni-stuttgart.de > <mailto:nwagner@iam.uni-stuttgart.de>> wrote: > > Hi Ed, > > Great ! > Thank you very much ! > It seems to work well with the exception of GMRES. > Please find attached a short test. > Am I missing something ? > > > No, I agree. It seems that there's something wrong with GMRES -- the > iter_ return parameter from gmresrevcom() is never incremented for > this example. I'm not familiar with the algorithm -- perhaps it really > does converge in one iteration? If not, it's just a bug. Could you > file a bug against gmres() in Tracker? > > -- Ed > Hi Ed, Done. Maybe you can add some comments. Again thank you very much for the new functionality !! http://projects.scipy.org/scipy/scipy/ticket/360 Cheers, Nils

4 3

sparse revisited
by Nils Wagner 01 Feb '07

01 Feb '07

Hi, A program of my colleague is broken due to the recent changes wrt sparse. So I turned back to r2450 and it works fine. This is the result with the latest svn version Solving frequency 450.0 Traceback (most recent call last): File "main_ce_v3.py", line 49, in ? x = spsolve(k_dyn, fq) File "/usr/local/lib64/python2.4/site-packages/scipy/linsolve/linsolve.py", line 86, in spsolve autoTranspose = True ) File "/usr/local/lib64/python2.4/site-packages/scipy/linsolve/umfpack/umfpack.py", line 581, in linsolve self.numeric( mtx ) File "/usr/local/lib64/python2.4/site-packages/scipy/linsolve/umfpack/umfpack.py", line 429, in numeric raise RuntimeError, '%s failed with %s' % (self.funs.numeric, RuntimeError: <function umfpack_zi_numeric at 0x2aff5e4412a8> failed with UMFPACK_ERROR_out_of_memory The matrices are very large and I am not sure if I should file a ticket. Any comments how to proceed (Tim, Nathan) ? Nils THANKS.txt might need some update. I cannot find some developers e.g. Tim Leslie, Nathan Bell.

3 5

Patch: Add 'insert' method to fitpack.py (in scipy.interpolate)
by Zachary Pincus 30 Jan '07

30 Jan '07

Hello, In the fortran fitpack library (aka ddierckx) that underlies scipy.interpolate, an insert.f function is provided to insert new knots into a given spline. Attached is a patch which exposes this function to python via the low- level fitpack.py interface. This patch properly handles both regular and parametric splines, and allows a 'knot multiplicity' to be specified so that a given knot can be inserted multiple times (which can be handled more efficiently in the c-level code than in python- level code). The inclusion of this function will make it easy to, among other things, convert a b-spline into a Bezier curve (useful for plotting purposes). (This can be accomplished by inserting new knots at the current knot points until each knot has a multiplicity equal to the spline degree.) I hope the patch will prove useful, Zach Pincus Program in Biomedical Informatics and Department of Biochemistry Stanford University School of Medicine

2 1

mlabwrap and scipy -- calling matlab from python
by Alexander Schmolck 30 Jan '07

30 Jan '07

Hi, I'm the author of mlabwrap a very high-level python to matlab bridge that's currently hosted at <http://mlabwrap.sourceforge.net>. There has been talk off-list about integrating mlabwrap into scipy and I'm following up on Jarrod Millman's suggestion to move this on-list. Robert Kern has suggested that mlabwrap could be the first scikit package hosted under scikits.scipy.org and I think it would be great to have mlabwrap integrated into the greater scipy-cosmos -- whilst mlabwrap is clearly too specific to go into scipy proper, matlab is very widely used and as python is getting increasingly popular for scientific computing I think there should be considerably demand for painlessly integrating matlab and python and being somehow part of scipy should increase exposure and help maintenance. The only caveat is that I myself have limited resources (I'm currently writing up a thesis and am also looking for a job [1]), so I can't pour much time into helping to set up the scikits infrastructure. I have however spent a good chunk of this weekend to bring mlabwrap up to scratch and just uploaded a new (alpha) version of mlabwrap to above sourceforge page; numpy and 64bit machines ought to work now (Numeric is also still supported) and thanks to Matthew Brett installation should be even easier (under linux that is; I have only been able to test the new version under 32 and 64-bit linux, confirmation that windows and OS-X installs also work would be very welcome). Mlabwrap has been pretty stable so far (no major bug reports as of yet and it has been around for a few years with users on all major platforms), but unfortunately I haven't had as much time to test this release as I'd like to, hence the alpha -- the unittests pass fine but there might be unspotted memory violations etc. There are 2 changes related to 64-bit readiness that might have introduced issues * I've changed the handle for the matlab session (lHandle) from int to a Py(C)Object * I've changed int to npy_intp at various places; I'm still effectively using ints and PyArray_FromDims in the matlab-array->ndarray direction -- as far as I can tell The Mathworks haven't yet properly nailed 64 bit support yet; there are mwSize and mwIndex types analogous to npy_intp, only that the docs in most places claim they are the same as int (somewhat implausibly, since it appears that from 7.3 you need to use them for 64-bit sparse arrays: <http://www.mathworks.com/access/helpdesk/help/techdoc/rn/index.html?/access…>). What I've done is that I've used mwSize and mwIndex in place of int where appropriate (and #define'd them as int if they are undefined), but cast to int (or int*) on calling PyArray_FromDims. I'm not sure what a better solution would currently be; this appears to work on 32 and 64-bit linux. I'd be grateful if people could give it a spin and even more for any code reviews or improvements (the code isn't very long -- about 500 non-comment lines of C++ and python each and only the C++ is problematic); valgrinding would also be very welcome (I haven't got around to it so far with the newest version). If you have experience porting Numeric code to numpy and with 64-bit issues and you'd like to use existing matlab code from python, now would be a great time to try mlabwrap out and/or have a look at the code. The sooner a shiny new mlabwrap 1.0 version ends up wherever it's meant to end up, the better. cheers, alexander schmolck Footnotes: [1] So should you happen to be currently looking for someone with scientific computing, pattern recognition and cognitive science expertise and pretty strong skills in python, matlab as well as a number of other languages...

1 0

Fwd: lapack license
by Nils Wagner 27 Jan '07

27 Jan '07

This might be of interest. Nils --- the forwarded message follows ---

1 0

Lapack tests too stringent?
by Ravikiran Rajagopal 24 Jan '07

24 Jan '07

Hi, Running scipy.test(10,10) on revision 2591 causes a failure in check_syevr_irange (scipy.lib.tests.test_lapack.test_flapack_float): AssertionError: Arrays are not almost equal (mismatch 33.3333333333%) x: array([-0.66992444, 0.48769474, 9.18222523], dtype=float32) y: array([-0.66992434, 0.48769389, 9.18223045]) The problem with this test is that the double precision version of this test passes on my machine but not the single precision version. To eliminate SciPy and ATLAS from the chain of causation, I ran the attached program with the reference BLAS and the reference LAPACK from netlib. I have similar failures on 32-bit and on 64-bit machines, both running up to date Fedora Core 6. With LAPACK 3.1.0, here are the eigenvalue results for double precision for the array from the test: -0.669924337185138 0.487693886153334 9.1822304510318 The above agrees with the reference values in the scipy test above. However, with single precision, the values are: -0.669923424720764 0.487694501876831 9.18223190307617 This does not meet the requirement for arrays being almost equal; is this because eps is dependent on double precision and is too low for single precision? To confirm that this is indeed an issue with LAPACK, I ran it on old RHEL3 32-bit machine with gcc/g77 (as opposed to gcc/gfortran on my FC6 machines) and with lapack 3.0 (up to date with all patches). In this case too, there is a difference between single precision and double precision values: Single Double -0.669924139976501 -0.669924337185137 0.487694054841995 0.487693886153333 9.18223094940186 9.1822304510318 Do others see similar numbers? The attached program can be compiled using gcc -o sytest syevr.c -llapack -lblas -lg2c where -lg2c is not required if you are using gfortran. To switch between single and double precision tests, change the typedef at the top of the file, and the first letter of the call on line 66. In my tests, changing the value of abstol did not make any difference. I wonder why I have not encountered this before since the tests have been in SciPy for a long time according to the SVN log. Regards, Ravi

1 0

Docstring standards for NumPy and SciPy
by Travis Oliphant 24 Jan '07

24 Jan '07

There was a lively discussion on the SciPy List before Christmas regarding establishing a standard for documentation strings for NumPy / SciPy. I am very interested in establishing such a standard. A hearty thanks goes to William Stein for encouraging the discussion. I hope it is very clear that the developers of NumPy / SciPy are quite interested in good documentation strings but recognize that producing them can be fairly tedious and un-interesting work. This is the best explanation I can come up with for the relative paucity of documentation rather than some underlying agenda *not* to produce them. I suspect a standard has not been established largely because of all the discussions taking place within the documentation communities of epydoc, docutils, etc. and a relative unclarity on what to do about Math in docstrings. I'd like to get something done within the next few days (like placing the standard on a wiki and/or placing a HOWTO_DOCUMENT file in the distribution of NumPy). My preference is to use our own basic format for documentation with something that will translate the result into something that the epydoc package can process (like epytext or reStructuredText). The reason, I'd like to use our own simple syntax, is that I'm not fully happy with either epytext or reStructuredText. In general, I don't like a lot of line-noise and "formatting" extras. Unfortuntately both epytext and reStructuredText seem to have their fair share of such things. Robert wrote some great documentation for a few functions (apparently following a reStructuredText format). While I liked that he did this, it showed me that I don't very much like all the line-noise needed for structured text. I've looked through a large number of documentation strings that I've written over the years and believe that the following format suffices. I would like all documentation to follow this format. This format attempts to be a combination of epytext and restructured text with additions for latex-math. The purpose is to make a docstring readable but also allowing for some structured text directives. At some point we will have a sub-routine that will translate docstrings in this format to pure epytext or pure restructured text. """ one-line summary not using variable names or the function name A few sentences giving an extended description. Inputs: var1 -- Explanation variable2 -- Explanation Outputs: named, list, of, outputs named -- explanation list -- more detail of -- blah, blah. outputs -- even more blah Additional Inputs: kwdarg1 -- A little-used input not always needed. kwdarg2 -- Some keyword arguments can and should be given in Inputs Section. This is just for "little-used" inputs. Algorithm: Notes about the implemenation algorithm (if needed). This can have multiple paragraphs as can all sections. Notes: Additional notes if needed Authors: name (date): notes about what was done name (date): major documentation updates can be included here also. See also: * func1 -- any notes about the relationship * func2 -- * func3 -- (or this can be a comma separated list) func1, func2, func3 (For NumPy functions, these do not need to have numpy. namespace in front of them) (For SciPy they don't need the scipy. namespace in front of them). (Use np and sp for abbreviations to numpy and scipy if you need to reference the other package). Examples: examples in doctest format Comments: This section should include anything that should not be displayed in a help or other hard-copy output. Such things as substitution-directed directives should go here. """ Additional Information: For paragraphs, indentation is significant and indicates indentation in the output. New paragraphs are marked with blank line. Text-emphasis: Use *italics*, **bold**, and `courier` if needed in any explanations (but not for variable names and doctest code or multi-line code) Math: Use \[...\] or $...$ for math in latex format (remember to use the raw-format for your text string in such cases). Place it in a new-paragraph for displaystyle or in-line for inline style. References: Use L{reference-link} for any code links (except in the see-also section). The reference-link should contain the full path-name (unless the function is in the same name-space as this one is. Use http:// for any URL's Lists: * item1 - subitem + subsubitem * item2 * item3 or 1. item1 a. subitem i. subsubitem1 ii. subsubitem2 2. item2 3. item3 for lists. Definitions: descripition This is my description for any definitions needed. Addtional Code-blocks: {{{ for multi-line code-blocks that are not examples to be run but should be formatted as code. }}} Tables: Tables should be constructed as either: +------------------------+------------+----------+ | Header row, column 1 | Header 2 | Header 3 | +========================+============+==========+ | body row 1, column 1 | column 2 | column 3 | +------------------------+------------+----------+ | body row 2 | Cells may span | +------------------------+-----------------------+ or || Header row, column 1 || Header 2 || Header 3 || ------------------------------------------------------- || body row, column 1 || column 2 || column 3 || || body row 2 |||| Cells may span columns || Footnotes: [1] or [CITATION3] for Footnotes which are placed at the bottom of the docstring as [1] Footnote [CITATION3] Additional note. Substitution: Use |somename|{optional text} with (the next line is placed at the bottom of the docstring in the Comments: section) .. |somename| image::myfile.png or .. |somename| somedirective:: {optional text} for placing restructured text directives in the main text. Please address comments to this proposal, very soon. I'd like to finalize it within a few days. -Travis

5 11

A patch for weave.catalog
by Fernando Perez 19 Jan '07

19 Jan '07

Hi all, does anyone object to the patch at the end? Rationale: when running using scipy in parallel on a cluster whose members share an NFS filesystem, the current code can blow up because the test file is hardcoded to dir+'/dummy', and all the processes race to manipulate the same file. My patch tries to fix this using portable calls (as far as the docs say, both socket.gethostname() and os.getpid() are fully portable). If nobody objects in 24 hours or I get an explicit OK before that from a core dev, I'll be happy to commit this. I've tested it and it solves my recurrent deadlocks on my system. I suspect scipy hadn't seen too much parallel use before, so I'm sure we'll begin finding other similar little gremlins as time goes by. I'll be happy to fix the ones I can. Cheers, f ############ Index: Lib/weave/catalog.py =================================================================== --- Lib/weave/catalog.py (revision 2579) +++ Lib/weave/catalog.py (working copy) @@ -33,6 +33,7 @@ import os,sys,string import pickle +import socket import tempfile try: @@ -127,13 +128,30 @@ os.mkdir(p) def is_writable(dir): - dummy = os.path.join(dir, "dummy") + """Determine whether a given directory is writable in a portable manner. + + :Parameters: + - dir: string + A string represeting a path to a directory on the filesystem. + + :Returns: + True or False. + """ + + # Do NOT use a hardcoded name here due to the danger from race conditions + # on NFS when multiple processes are accessing the same base directory in + # parallel. We use both hostname and pocess id for the prefix in an + # attempt to ensure that there can really be no name collisions (tempfile + # appends 6 random chars to this prefix). + prefix = 'dummy_%s_%s_' % (socket.gethostname(),os.getpid()) try: - open(dummy, 'w') - except IOError: - return 0 - os.unlink(dummy) - return 1 + tmp = tempfile.TemporaryFile(prefix=prefix,dir=dir) + except OSError: + return False + # The underlying file is destroyed upon closing the file object (under + # *nix, it was unlinked at creation time) + tmp.close() + return True def whoami(): """return a string identifying the user."""

2 2