Mailman 3 September 2004 - SciPy-Dev

MCMC, Kalman Filtering, AI for SciPy?
by chris＠fisher.forestry.uga.edu 01 Oct '04

01 Oct '04

For the past couple of years, I have been maintaining a package called PyMC (http://pymc.sourceforge.net) which provides a general Markov chain Monte Carlo sampler (Metropolis-Hastings), linear and non-linear Kalman filtering algorithms and several reinforcement learning algorithms for AI applications. To date, it has been used mostly for ecological modeling applications, but is far more general than that. This package currently relies on SciPy, particularly for such tasks as plotting, and uses f2py extensively. I am wondering if there is interest in integrating such functionality into SciPy. If so, I would be happy to contribute, integrate and maintain the code under the auspices of the SciPy project. If not, I am equally content to continue maintaining PyMC separately. I encourage those interested to check out PyMC and let me know if it would be worthwhile to add any of this to SciPy. Cheers, Chris Fonnesbeck chris (at) fisher.forestry.uga.edu

12 62

RE: [SciPy-dev] MCMC, Kalman Filtering, AI for SciPy?
by Charles Harris 30 Sep '04

30 Sep '04

From: scipy-dev-bounces(a)scipy.net on behalf of eric jones > >Robert Kern wrote: >> Charles Harris wrote: >> >> [snip] >>> >>> I agree that search and indexing are the best ways to find stuff, but >>> I am mostly concerned as to where to commit stuff. Clustering, where >>> does that go? >> >> >> scipy.cluster I would imagine. ;-) >> >>> Lattice methods, where do they go? How about useful data structures >>> or combinatorics? So on and so forth. I think the upper level GAMS >>> categories cover sufficient range that most things can be put into a >>> directory without embarrassment. As to the detailed breakdown in the >>> GAMS sub-classifications, I am not so sure. >> >> >> To make the discussion a bit more concrete, here is an example >> directory structure corresponding to the top-level GAMS >> classifications. The names are all my own, so feel free to pretend >> they are something more to your liking. >> >> <snip hierarchy> >> Now that I see it, it is somewhat appealing. I would probably want to >> break up some of those into two or more top-level groups. I definitely >> don't want to see too many subpackages under each of the top-level >> groups ("Flat is better than nested."). >> >>Here is where the current SciPy modules would likely get lumped in the >>GAMS hierarchy. >scipy/ > analysis/ > numbertheory/ > functions/ special > linalg/ linalg, sparse > interpolation/ interpolate > rootfinding/ > optimization/ optimize, ga > calculus/ integrate > diffeq/ > integraltransforms/ fftpack > approximation/ > probstat/ stats > simulation/ > datahandling/ io > symbolic/ > geometry/ > graphics/ xplt, gplt, plt > service/ gui_thread > develop/ > other/ cow, cluster ??, signal ?? >(Cluster and signal didn't fit anywhere obvious to me) Gams has clustering under probstat. The union/find algorithm could also go under data where they have trees, etc. Lattice stuff goes into numbertheory. Service also contains the machine parameters, eps and such. Remez algorithm goes into approximation (L_inf). Weave into develop. Hmm, signal processing is missing somewhere. Markov chain into simulation. Grey codes, permutations into other, although GAMS has permutations under data. We could split diffeq into ode, pde, but it is probably ok as is. Control is under optimization (for optimal control) but could be brought to the top. The addition of rootfinding is good. Convolutions and such are under integraltransforms, Haar transforms would go there also. Filters are under probstats in time series analysis, so it might make sense to create a signal (time series?) directory, probstats seems to be an overloaded GAMS category that could use some upperlevel subdivision. >The naming conventions are often quite similar. The SciPy names are >generally shorter which is nice for typing. Where SciPy has multiple >packages [(linalg, sparse), (optimize, ga), etc.], it is likely a good >idea. Like you, I don't want to see a deep nesting in the package >structure. >Looking at this, I don't see any real reason to reorganize top level >package names. Are any of them that bad or misleading? On the other >hand, I do think we should reorganize the functions within them some to >fix the places where they are organized based on "build" convenience >instead of actual function. This will probably necessitate the addition >of new top level groups and maybe the pruning of one of the current >ones. I've made a Wiki page to keep suggestions that people have: > http://www.scipy.org/wikis/featurerequests/PackageReorganization >If you update the page, you might also post to python-dev so that people >know to go check on the Wiki (that is so painful...). We can obviously >also just discuss it here and then transfer to the Wiki later. [side >note: this using a wiki and a mailing list for communication is also a >little painful]. >> Fernando, could you give an example or two where you would want to >> replicate a function across sub-packages? I'm wary of doing so as >> there is already the enormous amount of replication with respect to, >> at least, the base Numeric functions. Try scipy.special.<tab> in >> IPython. I realize what you're proposing doesn't even come close to >> that, but I'd like an example in any case. >I don't like the replication idea very well. I think things should live >in one place. Otherwise people will wonder if two functions that are >actually the same have different purposes, implementation, etc. >> >> And since we are talking about re-organization, is there anything we >> can do about the problem I just mentioned? It wreaks havoc with not >> only tab-completion but also automatic documentation generation [1]. >> Is it practical to be careful about what we import into __init__.py? >> By which I mean not doing "from foo import *" in __init__.py where >> foo.py does "from scipy_base import *". On the other hand, explicitly >> listing all of the names in special is gonna be a major pain and >> fragile to boot. <snip> _______________________________________________ Scipy-dev mailing list Scipy-dev(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-dev

2 1

Status of iterative solvers
by Nils Wagner 30 Sep '04

30 Sep '04

Hi all, BTW, when are the iterative solvers up for general testing in scipy ? Nils

1 0

Re: Slicot library and scipy
by Travis Oliphant 28 Sep '04

28 Sep '04

Nils Wagner wrote: > Dear Travis, > > As you encouraged me to learn f2py I would like to ask you, > if you can help me with an interface for continuous and discrete Lyapunov > equations as a starting point. Once I know how to manage it in general, > I would be glad to submit more slicot routines to scipy. It looks like your problems are related to making sure SB03MD has the libraries it needs to run. It looks like it needs a dummy SELECT function for DGEES (it actually never calls it but the interface requires it). You also need to be sure that SB03MY and SB03MX and MB01RD are compiled and linked along with SB03MD (it calls them). I'm attaching a modified .pyf file that you could use and a util.f file with the SELECT function defined. You need to compile like this (assuming you have a libslicot.a placed in <path-to-slicot>) f2py -L<path-to-slicot> -lslicot -L<path-to-atlas> -llapack -lf77blas -latlas -c slicot.pyf utils.f This is the better way to do it and will allow you to grow slicot.pyf until you've made available all of the routines. You can start the process of wrapping all of the slicot library by letting f2py construct the basic .pyf file by going to the slicot src directory and doing f2py -h slicot.pyf -m slicot *.f This will give you a basic, raw interface to all of the slicot routines. You will then want to edit the slicot.pyf file (like I did here for sb03md) to make the routine more useful --- hide some variables, show that some are input and some are output variables --- and so forth. You will find along the way that f2py will fail to parse some of the SLICOT fortran routines. Just don't include them in the list of routines to make an interface for. I needed to move IB03AD.f and IB03BD.f out of the way before running the command to get it to complete. The resulting slicot.pyf file is huge and contains a very raw interface to the entire slicot library (minus those two routines --- I'm submitting a bug report to Pearu now). The slicot.pyf file can be compiled with f2py as before. If you want to write a setup.py file you can (but that is another story). > > I am quite sure that many scipy users are interested in features > provided by > the slicot library. > > What do you think ? > SLICOT would be a nice addition to SciPy. But... I noticed that the SLICOT folks don't want people distributing SLICOT with their library, so it looks like only the interface could be distributed with SciPy and you would have to go and download the SLICOT tools from them, which is too bad unless some kind of agreement and separate license for SciPy could be obtained --- it looks like they mainly want credit, so in SciPy we could put in the help for the tools that if you use them you must credit them in a paper. At any rate, perhaps the interface could be distributed and then a check made at compile time (by default off) to see if the library is already there then the interface could be built. -Travis > I look forward to heraring from you. > > Cheers, > Nils > > > LOGICAL FUNCTION SELECT (X,Y) DOUBLE PRECISION X, Y SELECT = .TRUE. RETURN END ! -*- f90 -*- python module slicot ! in interface ! in :slicot subroutine sb03md(dico,job,fact,trana,n,a,lda,u,ldu,c,ldc,scale,sep,ferr,wr,wi,iwork,dwork,ldwork,info) ! in :slicot:SB03MD.f character :: dico character :: job character :: fact character :: trana integer intent(hide), depend(a) :: n=shape(a,1) double precision intent(in,out), dimension(lda,n) :: a integer intent(hide), depend(a) :: lda=shape(a,0) double precision intent(in,out),dimension(ldu,n) :: u integer intent(hide),depend(u) :: ldu=shape(u,0) double precision intent(in,out),dimension(ldc,n) :: c integer intent(hide),depend(c) :: ldc=shape(c,0) double precision intent(out) :: scale double precision intent(out) :: sep double precision intent(out) :: ferr double precision intent(out),dimension(n) :: wr double precision intent(out),dimension(n) :: wi integer intent(hide), dimension(n*n) :: iwork double precision intent(out), dimension(ldwork) :: dwork integer intent(hide), depend(n) :: ldwork = (2*n+2)*n integer intent(out) :: info end subroutine sb03md end interface end python module slicot ! This file was auto-generated with f2py (version:2.39.235_1703). ! See http://cens.ioc.ee/projects/f2py2e/

3 2

matrixmultiply!=dot when dotblas is present
by Fernando Perez 28 Sep '04

28 Sep '04

Hi all, I just reported a Numpy bug to the numpy list, but scipy can work around it trivially. The problem is the following (I'm copy/pasting my original msg): I found something today a bit unpleasant: if you install numeric without any BLAS support, 'matrixmultiply is dot==True', so they are fully interchangeable. However, to my surprise, if you build numeric with the blas optimizations, they are NOT identical. The reason is a bug in Numeric.py. After they define dot, they do this: #This is obsolete, don't use in new code matrixmultiply = dot and at the very end of the file, they do: # try to import blas optimized dot, innerproduct and vdot, if available try: from dotblas import dot, innerproduct, vdot except ImportError: pass Obviously this means that matrixmultiply is stuck with the _old_ definition of dot, and does not benefit from the blas optimizations. This is BAD, as for a 1024x1024 matrix the difference is staggering: planck[Numeric]> pylab In [1]: a=rand(1024,1024) In [2]: b=rand(1024,1024) In [3]: from IPython.genutils import timing In [4]: timing 1,dot,a,b ------> timing(1,dot,a,b) Out[4]: 0.55591500000000005 In [5]: timing 1,matrixmultiply,a,b ------> timing(1,matrixmultiply,a,b) Out[5]: 68.142640999999998 In [6]: _/__ Out[6]: 122.57744619231356 Pretty significant difference, eh? :) Scipy can work around this problem (which is still there in Numpy 23.4) by manually doing 'matrixmultiply=dot' AFTER dot has been imported from Numeric. This will guarantee that users of scipy.matrixmultiply() actually get the fast version. Cheers, f

1 0

erf
by Charles Harris 27 Sep '04

27 Sep '04

Which erf is scipy using? There is one in cephes/ndtr.c and another in cdflib/erf.f . The two don't give the same results for smallish x, and I suspect that erf.f is the more accurate because ndtr computes it by 1.0 - erfc(x) .

1 0

Scipy for Numarray
by Charles Harris 22 Sep '04

22 Sep '04

Hi All, I've been thinking about getting started on a version of scipy for numarray. It's going to be a big job, so I have questions to put to the community. 1. How should it be set up in cvs? There is currently a scipy module to be checked out, should we have a numarray_scipy module at the same level? 2. Where should it install? To avoid name conflicts, I think it might be reasonable to install it, or its components, in the numarray packages directory. 3. There is stuff like Pycrust and ipython that is not directly related to scipy, should it be taken over also, or perhaps broken out into a scipy addons module. 4. There is a lot of infrastructure for the installation and testing of scipy. I am not very familiar with that, so how does it need to be modified to install a numarray version of scipy. 5. How can we keep both versions updated? Chuck

4 9

Weave - strange behaviour
by Manuel Metz 22 Sep '04

22 Sep '04

Hi, I was really excited about weave. But I found a (to me) somewhat strange behaviour, which I do not understand. I used the following code to test weave: [...] class test(object): def wtest(self): a = numarray.array( [1,2,3] ) assert( type(a)==type( numarray.array([]) ) ) ccode2 = """ int val; using namespace std; for(int i=0; i<a.length(); i++) { val = py_to_int( a[i], "val" ); cout << val << endl; }""" weave.inline( ccode2, ["a",] ) [...] Okay, the code worked fine, output of the program was: >>>>> file changed None 1 2 3 <<<<< Now if I only change the last line of the C-code from }""" to } """ so nothing more than just a line-break, I always (even after the C-code is already compiled) get the following output: >>>>> repairing catalog by removing key file changed None 1 2 3 <<<<< And it seams that EACH time I call the function the C-code gets compiled. This behavior also showed up depending on the first line of the C-code (line break or not; additionen spaces). Has anyone else seen this? Any idea what is going wrong? How can I prevent this strange behaviour (except try and error)? My system: - Debian linux (testing) - Package: python2.3-scipy Version: 0.3.0+266.4239-1 Manuel -- ------------------------------------- Manuel Metz Sternwarte der Universitaet Bonn Auf dem Huegel 71 (room 3.06) D - 53121 Bonn E-Mail: mmetz(a)astro.uni-bonn.de Phone: (+49) 228 / 73-3660 Fax: (+49) 228 / 73-3672 -------------------------------------

3 2

Performance Python with Weave article updated
by Prabhu Ramachandran 22 Sep '04

22 Sep '04

Hi folks, Just to let you know that I've updated the article on weave compared with other options here: http://www.scipy.org/documentation/weave/weaveperformance.html The examples now work with weave from SciPy-0.3. I've also added a comparison with Psyco and Pyrex. The tarball has also been updated and includes all the sources along with a trivial setup.py to build the f2py module and the Pyrex module. cheers, prabhu

2 2

Re: Python & Scientific computing links?
by Alexey Goldin 16 Sep '04

16 Sep '04

I just found a great resource --- "social bookmarking tool" at http://del.icio.us. You can use it as a bookmark manager like bookmarks.yahoo.com, however this is much more then that. You can tag your bookmarks with one or more tag and browse other people bookmarks if they use the same tag. The database is accessible through RSS and some kind of RPC as well. Posting a bookmark via bookmarklets or Mozilla extension is trivial. For example, here is RSS feed for python related resources: http://del.icio.us/rss/hedgehog/python and here is one for scientific python: http://del.icio.us/rss/hedgehog/scipy (I think for now I am the only one using scipy tag for now) Using RSS feed is a nice way to discover resources other people found....

1 0