Mailman 3 April 2007 - SciPy-Dev

Status of the doc format for scipy code ?
by David Cournapeau April 26, 2007

April 26, 2007

Hi, I wanted to know what is the status of the doc format for scipy (http://projects.scipy.org/scipy/numpy/wiki/DocstringStandards) ? Is anything settled down ? cheers, David

4 10

Which matrix library in C++ for scipy
by Matthieu Brucher April 26, 2007

April 26, 2007

Hello, I'm wondering what C++ matrix library I should for future inclusion of code in scipy. What I want to do is to port a class that can do Parzen window or K-neighboors, and it uses my own matrix library. Should I use Blitz ++ ? Matthieu

4 15

Updated generic optimizers proposal
by Matthieu Brucher April 25, 2007

April 25, 2007

Hi, I'm lauching a new thread, the last was pretty big, and as I almost put every advice in this proposal, I thought it would be better. First, I used scipy coding standard, I hope I didn't forget something. I do not know where it would be put at the moment on my scipy tree, and the tests are visual for the moment, I have to make them automatic, but I do not know the framework used by scipy, I have to check it first. So, the proposal : - combining several objects to make an optimizer - a function should be an object defining the __call__ method and graient, hessian, ... if needed. It can be passed as several separate functions as Alan suggested it, a new object is then created - an optimizer is a combination of a function, a step_kind, a line_search, a criterion and a starting point x0. - the result of the optimization is return after a call to the optimize() method - every object (step or line_search) saves its modification in a state variable in the optimizer. This variable can be accessed if needed after the optimization. - after each iteration, a record function is called with this state variable - it is a dict, BTW -, if you want to save the whole dict, don't forget to copy it, as it is modified during the optimization For the moment are implemented : - a standard algorithm, only calls step_kind then line_search for a new candidate - the next optimizer would be one that calls a modifying function on the computed result, that can be useful in some cases - - criteria : - monotony criterion : the cost is decreasing - a factor can be used to allow an error - - relative value criterion : the relative value error is higher than a fixed error - absolute value criterion : the same with the absolute error - step : - gradient step - Newton step - Fletcher-Reeves conjugate gradient step - other conjugate gradient will be available - - line search : - no line search, just take the step - damped search, it's an inexact line search, that searches in the step direction a set of parameters than decreases the cost by dividing by two the step size while the cost is not decreasing - Golden section search - Fibonacci search I'm not pulling other criterion, step or line search, as my time is finite when doing a structural change. There are 3 classic optimization test functions in the package, Rosenbrock, Powell and a quadratic function, feel free to try them. Sometimes, the optimizer converges to the true minimum, sometimes it does not, I tried to propose several solutions to show that every combinaison does not manage to find the minimum. Matthieu

4 7

Dataset for examples and license
by David Cournapeau April 25, 2007

April 25, 2007

Hi, I would like to know what should be done when including some dataset in scipy ? For example, during the development of my project pymachine, I would like to include some famous data like iris/old faithful data, etc... for demo of classic machine learning algorithms. R has some intereseting data, but is licensed under the GPL, and I am not quite sure what the status of the data are wrt the license ? Does GPL also cover raw data ? cheers, David

5 9

Possible bug in scipy.stats.kurtosistest()
by Rudolph van der Merwe April 24, 2007

April 24, 2007

I think I found a bug in the kurtosistest function in the scipy.stats module. The relevant file is scipy/stats/stats.py The kurtosistest function seems to be a direct implementation of the algorithm described by D'Agostino et all in the following paper: R. B. D'Agostina, A. Belanger, and R. B. D'Agostino, Jr. "A Suggestion for Using Powerful and Informative Tests of Normality", The American Statistician, Vol. 44, No. 4. (Nov., 1990), pp. 316-321 One of the first steps in the algorithm is calculating the kurtosis of the distribution. This is done using the stats.kurtosis function. This function can either calculate the Fisher or the Pearson form of the kurtosis. The scipy code in kurtosistest makes use of the Fisher form which is the default for the kurtosis function. The algorithm in D'Agostino's paper, however, makes use of the Pearson form of the kurtosis. I verified this issue on the numerical example given in the paper. The results calculated by scipy.stats.kurtosistest only agrees with the paper if the call to kurtosis in kurtosistest is changed from b2 = kurtosis(a, axis) to b2 = kurtosis(a, axis, fisher=False) Since kurtosistest is used by the scipy.stats.normaltest function (which does a D'Agostino-Pearson normality test), this bug affects that function's correctness as well. Is there a more formal way of filing this bug than posting a message on this forum? Rudolph P.S. I'm using Scipy version 0.5.2 -- Rudolph van der Merwe KAT (Karoo Array Telescope) / www.kat.ac.za

3 2

SciPy improvements
by Ondrej Certik April 24, 2007

April 24, 2007

Hi, I am studying theoretical physics and I have collected a lot of useful python code, that I believe could go to SciPy. So first I'll review what I have and if you find it interesting, I would like to discuss how I could implement it in SciPy. 1) Optimization http://chemev.googlecode.com/svn/trunk/chemev/optimization/ I have a Differential Evolution optimizer, Simplex optimizer, mcmc (not well tested yet), I took a code from someone else, but adapted the interface to the SciPy's one: def fmin_de(f,x0,callback=None,iter=None): Those are unconstrained optimizers. Then I have a constrains code, that applies a logistic function to the fitting variable and allows me to do constrained optimization. For example the L-BFGS with my constrains converges 7x faster, than the original L-BFGS-B on my problem. 2) Nonlinear solvers I have written these nonlinear solvers for the problem R(x) = 0, where x and R has a dimension "n": broyden1 - Broyden's first method - is a quasi-Newton-Raphson method for updating an approximate Jacobian and then inverting it broyden2 - Broyden's second method - the same as broyden1, but updates the inverse Jacobian directly broyden3 - Broyden's second method - the same as broyden2, but instead of directly computing the inverse Jacobian, it remembers how to construct it using vectors, and when computing inv(J)*F, it uses those vectors to compute this product, thus avoding the expensive NxN matrix multiplication. broyden_generalized - Generalized Broyden's method, the same as broyden2, but instead of approximating the full NxN Jacobian, it construct it at every iteration in a way that avoids the NxN matrix multiplication. This is not as precise as broyden3. anderson - extended Anderson method, the same as the broyden_generalized, but added w_0^2*I to before taking inversion to improve the stability anderson2 - the Anderson method, the same as anderson, but formulated differently linear_mixing exciting_mixing I use them in the self-consistent cycle of the Density Functional Theory (so I use a terminology of DFT literature in the names of the methods). Also I am writing a BFGS solver with linesearch, that should behave even better than the Broyden scheme. Of course I am trying to use SciPy's code (like linesearch) wherever possible. 3) PETSC bindings I found these nice petsc bindings: http://cheeseshop.python.org/pypi/petsc4py/0.7.2 I believe this could also be an optional package in SciPy. Because if SciPy has some sparse matrices code, then it should definitely also has this. 4) Finite element code? I have my own code, that uses libmesh: http://libmesh.sourceforge.net/ and calls tetgen and parses input from gmsh etc. Can convert the mesh, can refine it, can solve it, etc. Webpages are here: http://code.google.com/p/femgeom/ http://code.google.com/p/libmeshpetscpackage/ http://code.google.com/p/grainmodel/ I am not sure here, if it should belong to SciPy. Probably not. 5) Symbolic manipulation in Python http://code.google.com/p/sympy/ We'll have some google summer of code students working on this and also I am not sure if it belongs to SciPy. However, this project looks very promising. ----------------- So that's it. I have some comments to SciPy: 1) Documentation Virtually none, I just use the source code to understand, what SciPy can do and how. But the docstrings are good though. I would suggest to update the http://www.scipy.org/doc/api_docs/ more often (for example I didn't find there the new l-bfgs code). 2) What is the official NumPy's page? I believe it should be http://numpy.org/ however, it points to a sourceforge page. I believe this homepage should contain all the relevant information about numpy and contain links to the fee based documentation and possibly some tutorials. The documentation of NumPy is scattered across many pages and I find it confusing. I know you have some list here: http://www.scipy.org/MigratingFromPlone But I am quite confused from the whole SciPy page. I think less but unconfusing information is better, but that's just my opinion. I think the front page of both SciPy and NumPy should be clean and simple with a clear link to a documentation. Like http://www.pytables.org/moin/PyTables http://matplotlib.sourceforge.net/ However, the documentation: http://www.scipy.org/Documentation is confusing, because except the fee based guide to NumPy, it's not clear, what is the official SciPy doc and what the best way of learning SciPy is. So I am interested in your opinions and then I'll just integrate my code into scipy.optimize and send you a patch to this list? Ondrej

8 27

Disappearing websites and unlicensed code
by Jack Twilley April 24, 2007

April 24, 2007

I was looking for some examples on using Python's wave module to read and write WAV files and came across the tutorial formerly found at http://scipy.mit.edu/tutorials/wave.pdf -- I say formerly found because I haven't been able to get to that website since I first found the link via Google. I was able to access the file via Google's cache and I have been using the code for private work, but the code in the tutorial is unlicensed as it sits. I checked the web for scipy and found the website which points to this mailing list and thought I'd ask here if anyone knows how the aforementioned code is licensed and what happened to the MIT site. Jack.

3 2

BVP solver and dependence on Fortran 90/95
by John Travers April 23, 2007

April 23, 2007

Hi all, I have a number of modules I'm working on for my own use which depend on Fortran 90/95 language features. One is a python wrapper to a boundary value problem solver (http://cs.smu.ca/~muir/BVP_SOLVER_Webpage.shtml) which would be a very good addition to scipy I think (the authors have explicitly agreed that it can be used in scipy). My question is: what is the policy on scipy modules depending on Fortran 90/95 code. There currently doesn't appear to be any in the repository. However with the availability of gfortran now broadening and the fact that most (all?) commercial compilers now have full support, is it any longer a problem (putting aside f2py support)? Cheers, John

6 8

Solution for problem with f90 allocatable arrays on 64-bit systems
by Neilen Marais April 21, 2007

April 21, 2007

Hi, I've found the problem in the generated code that results in the bug described here: http://projects.scipy.org/scipy/numpy/ticket/147 The problem is an index which should be specified as pointer size but is specified as integer in the fortran part of the module code. When wrapping the following code: (geomwrap.f90) MODULE mesh ! The 4 node indices per element that define all the mesh elements INTEGER, DIMENSION(:,:), ALLOCATABLE :: element_nodes CONTAINS SUBROUTINE init_mesh() ALLOCATE(element_nodes(4,2)) element_nodes(:,1) = (/1, 2, 3, 4/) element_nodes(:,2) = (/1, 3, 4, 5/) END SUBROUTINE init_mesh END MODULE mesh and generating wrappers: f2py --fcompiler=gnu95 -m geo --build-dir test -c geomwrap.f90 in the file geo-f2pywrappers2.f90 ! -*- f90 -*- ! This file is autogenerated with f2py (version:2_3718) ! It contains Fortran 90 wrappers to fortran functions. subroutine f2py_mesh_getdims_element_nodes(r,s,f2pysetdata,flag) use mesh, only: d => element_nodes integer flag external f2pysetdata logical ns integer s(*),r,i,j the definition for s(*) should be interger(8) (at least on 64-bit machines), i.e. ! -*- f90 -*- ! This file is autogenerated with f2py (version:2_3718) ! It contains Fortran 90 wrappers to fortran functions. subroutine f2py_mesh_getdims_element_nodes(r,s,f2pysetdata,flag) use mesh, only: d => element_nodes integer flag external f2pysetdata logical ns integer(8) s(*) integer r,i,j is correct. I'm going to patch f2py on my machine to generate integer(8) sized s(*), but this will of course be wrong on 32-bit machines. What is the recommended way of getting a pointer-sized int in fortran? Thanks Neilen

1 1

Genetic algorithm optimization
by John Travers April 20, 2007

April 20, 2007

Hi all, I have recently needed a genetic algorithm code for some work. There appear to be various python implementations, which I didn't try and the package in the sandbox looked promising initially but I think a large amount of work needs to be done to get it going again. As I was pushed for time a took an easy option and wrapped what appears to be a fairly robust implementation in Fortran 77 called pikaia: http://www.hao.ucar.edu/Public/models/pikaia/pikaia.html It has apparently been used in a number of scientific research projects and appears to be quite mature. I have a very simple f2py interface (attached along with a patch to the pikaia source which I slightly modified) and python calling function (also attached). While it doesn't provide a fully customizable, all bells and whistles GA library which some people may want (see attached ga.py for docs on what features it does have), it does work very well for standard global optimization problems (see example, also attached). It would fit very neatly in the optimization package of scipy as it is a small one function routine, see attached test script for an example. The GA code itself is explicitly placed in the public domain. The file you can download from the above website contains a quicksort and uniform random number routine which are under the ACM license which is incompatible with scipy. But these parts are easily replaced with compatibly licensed versions (suggestions for such routines would be welcome). What do people think about including it in the optimization module? Cheers, John

1 0