Mailman 3 September 2009 - SciPy-Dev

reimplementation of lfilter
by Sturla Molden Sept. 24, 2009

Sept. 24, 2009

Looking at this atrocious piece of C: http://svn.scipy.org/svn/scipy/trunk/scipy/signal/lfilter.c.src I got some inspiration for reusing some of my other C code to reimplement scipy.signal.lfilter. I basically just had to write this small C function + a Cython wrapper: static void lfilter_1d(double *restrict b, double *restrict a, double *x, double *y, double *zi, npy_intp K, npy_intp N) { double Z[K]; double *restrict pz = zi, *restrict z = Z; register npy_intp n, k; for (n=0; n<N; n++) { register double yn, xn = x[n]; yn = y[n] = b[0] * xn + pz[0]; for (k=0; k<K-1; k++) z[k] = b[k+1] * xn + pz[k+1] - a[k+1]*yn; z[K-1] = b[K]*xn - a[K]*yn; double *tmp = pz; pz = z; z = tmp; } if (pz != zi) memcpy(zi, pz, K*sizeof(double)); } Relative speed compared to scipy.signal.lfilter: One processor: 1D array, 1000000 floats: 150% 2D array, 1000 x 1000 floats, axis=0: 80% 2D array, 1000 x 1000 floats, axis=1: 150% It is faster than lfilter except for 2D with axis=0. For 2D with axis=1 it is probably faster than lfilter due to copy-in copy-out optimization. For arrays with more than one dimension, multithreading helps as well: Four processors: 1D array, 1000000 floats: 150% 2D array, 1000 x 1000 floats, axis=0: 160% 2D array, 1000 x 1000 floats, axis=1: 250% Multithreading obviously does not help for filtering single vector, as there is no work to share. Some other improvements over signal.lfilter, apart from readable Cython code: * GIL is released. The C code is not littered with Python C API calls that prevents this. * Filtering can be done in reverse, i.e. no need for fliplr or flipud for filtering backwards. * Arrays can be filtered in-place. * Cleaned up docstring. Still it only work as a "proof of concept" with double precision floats. Adding support for other dtypes would be easy, just replace double with any of: float float _Complex double _Complex long double long double _Complex Regards, Sturla Molden # Copyright (C) 2009 Sturla Molden import numpy as np cimport numpy as np cdef extern int _linear_filter( object a, object b, object x, object y, object z, int axis, int direction) def lfilter(object B, object A, object X, object zi=None, object axis = None, int direction=1, object inplace=False): """ Filter data along one-dimension with an IIR or FIR filter. Description =========== Filter a data sequence, x, using a digital filter. This works for many fundamental data types (including Object type). The filter is a direct form II transposed implementation of the standard difference equation (see "Algorithm"). Inputs: ======= b : vector-like The numerator coefficient vector in a 1-D sequence. a : vector-like The denominator coefficient vector in a 1-D sequence. If a[0] is not 1, then both a and b are normalized by a[0]. x : array-like An N-dimensional input array. axis : int or None The axis of the input data array along which to apply the linear filter. The filter is applied to each subarray along this axis. If axis is None, the filter is applied to x flattened. Default: axis=None. zi : array-like Initial conditions for the filter delays. It is a vector (or array of vectors for an N-dimensional input) of length max(len(a),len(b)). If zi=None or is not given then initial rest is assumed. If axis is None, zi should be a vector of length K = max(M,N), where M=len(b)-1 and N=len(a)-1. If an axis is not None, zi should either be a vector of length K, or an array with same shape as x, except for zi.shape[axis] which should be K. If zi is a vector, the same initial conditions are applied to each vector along the specified axis. Default: zi=None (initial rest). See signal.lfiltic for more information. direction : int If negative, the filter is applied in reverse direction along the axis. Default: direction=1 (filter forward). inplace : bool If True, allow x or zi to be modified inplace if possible. Default: False. Outputs: (y, {zf}) ======= y : ndarray The output of the digital filter. zf : ndarray If zi is None, this is not returned, otherwise, zf holds the final filter delay values. Algorithm: ========== The filter function is implemented as a direct II transposed structure. This means that the filter implements a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[nb]*x[n-nb] - a[1]*y[n-1] - ... - a[na]*y[n-na] using the following difference equations: y[m] = b[0]*x[m] + z[0,m-1] z[0,m] = b[1]*x[m] + z[1,m-1] - a[1]*y[m] ... z[n-3,m] = b[n-2]*x[m] + z[n-2,m-1] - a[n-2]*y[m] z[n-2,m] = b[n-1]*x[m] - a[n-1]*y[m] where m is the output sample number and n=max(len(a),len(b)) is the model order. The rational transfer function describing this filter in the z-transform domain is -1 -nb b[0] + b[1]z + ... + b[nb] z Y(z) = ---------------------------------- X(z) -1 -na a[0] + a[1]z + ... + a[na] z """ cdef np.ndarray b, a, x, y, z cdef np.npy_intp i, ierr cdef double tmp cdef object K, zshape, iaxis # get filter coeffs b, a A = np.asarray(A) B = np.asarray(B) if (A.ndim != 1) or (B.ndim != 1): raise ValueError, 'a and b must be vectors' K = max(len(A)-1,len(B)-1) b = np.zeros(K+1) a = np.zeros(K+1) a[:A.shape[0]] = A[:] b[:B.shape[0]] = B[:] # normalize by a[0] if a[0] != 1.0: tmp = a[0] for i in range(a.shape[0]): a[i] /= tmp b[i] /= tmp # set up input signal X = np.asarray(X, dtype=np.float64) if axis is None: X = X.ravel() elif type(axis) != int: raise ValueError, 'axis must be None or an integer >= 0' elif axis < 0: raise ValueError, 'axis must be None or an integer >= 0' x = X # set up output signal if inplace: y = X else: y = np.zeros(X.shape, dtype=np.float64) # set up filter delay buffer iaxis = 0 if (axis is None) else axis if zi is None: # zi not specified, assume initial rest zshape = list(X.shape) zshape[iaxis] = K z = np.zeros(zshape, dtype=np.float64) else: zi = np.asarray(zi, dtype=np.float64) # x and zi are both vectors if (x.ndim == 1) and (zi.ndim == 1): if (zi.shape[0] != K): raise ValueError, 'length of zi must be K' if inplace: z = zi else: z = zi.copy() # x is nd array, zi is vector # broadcast zi (cannot modify inplace) elif (x.ndim > 1) and (zi.ndim == 1): if (zi.shape[0] != K): raise ValueError, 'length of zi must be K' zshape = list(X.shape) zshape[iaxis] = K z = np.zeros(zshape, dtype=np.float64) zshape = [1] * X.ndim zshape[iaxis] = K z = z + zi.reshape(zshape) # x and zi are both nd arrays else: zshape = list(X.shape) zshape[iaxis] = K zshape = tuple(zshape) if (zi.shape != zshape): raise ValueError, ('bad shape of zi: expected %r, got %r' % (zshape,zi.shape)) if inplace: z = zi else: z = zi.copy() # apply the linear filter ierr = _linear_filter(a, b, x, y, z, <int> iaxis, direction) # check for error conditions on return if ierr == -1: raise MemoryError, 'malloc returned NULL' if ierr == -2: raise ValueError, 'shape of x, y, and z did not match... should never happen, debug please.' if ierr == -3: raise ValueError, 'shape of a and b did not match... should never happen, debug please.' # return y or (y, z) depending on zi # return (y if (zi is None) else (y, z)) ## Cython fails on this if (zi is None): return y else: return (y,z) /* Copyright (C) 2009 Sturla Molden */ #include <Python.h> #include <numpy/arrayobject.h> #include <string.h> #include <stdlib.h> typedef PyArrayObject ndarray; /* PyArrayObject is an ugly name */ typedef Py_ssize_t integer; /* Py_ssize_t of int for 64 bit support, npy_intp is typedef'ed incorrectly. */ /* copy data to and from temporary work arrays */ static void copy_in(double *restrict y, double *restrict x, integer n, integer s) { integer i; char *tmp; if (s == sizeof(double)) memcpy(y, x, n*sizeof(double)); else { tmp = (char *)x; for (i=0; i<n; i++) { *y++ = *x; tmp += s; x = (double *)tmp; } } } static void copy_out(double *restrict y, double *restrict x, integer n, integer s) { integer i; char *tmp; if (s == sizeof(double)) memcpy(y, x, n*sizeof(double)); else { tmp = (char *)y; for (i=0; i<n; i++) { *y = *x++; tmp += s; y = (double *)tmp; } } } /* this computes the start adress for every vector along a dimension (axis) of an ndarray */ inline char *datapointer(ndarray *a, integer *indices) { char *data = a->data; int i; integer j; for (i=0; i < a->nd; i++) { j = indices[i]; data += j * a->strides[i]; } return data; } static double ** _get_axis_pointer(integer *indices, int axis, ndarray *a, double **axis_ptr_array, int dim) { /* recursive axis pointer search for 4 dimensions or more */ integer i; double *axis_ptr; if (dim == a->nd) { /* done recursing dimensions, compute address of this axis */ axis_ptr = (double *) datapointer(a, indices); *axis_ptr_array = axis_ptr; return (axis_ptr_array + 1); } else if (dim == axis) { /* use first element only */ indices[dim] = 0; axis_ptr_array = _get_axis_pointer(indices, axis, a, axis_ptr_array, dim+1); return axis_ptr_array; } else { /* iterate range of indices */ integer length = a->dimensions[dim]; for (i=0; i < length; i++) { indices[dim] = i; axis_ptr_array = _get_axis_pointer(indices, axis, a, axis_ptr_array, dim+1); } return axis_ptr_array; } } static double **get_axis_pointers(ndarray *a, int axis, integer *count) { integer indices[a->nd]; double **out, **tmp; integer i, j; j = 1; for (i=0; i < a->nd; i++) { if (i != axis) j *= a->dimensions[i]; } *count = j; out = (double **) malloc(*count * sizeof(double*)); if (out == NULL) { *count = 0; return NULL; } tmp = out; switch (a->nd) { /* for one dimension we just return the data pointer */ case 1: *tmp = (double *)a->data; break; /* specialized for two dimensions */ case 2: switch (axis) { case 0: for (i=0; i<a->dimensions[1]; i++) *tmp++ = (double *)(a->data + i*a->strides[1]); break; case 1: for (i=0; i<a->dimensions[0]; i++) *tmp++ = (double *)(a->data + i*a->strides[0]); break; } break; /* specialized for three dimensions */ case 3: switch (axis) { case 0: for (i=0; i<a->dimensions[1]; i++) for (j=0; j<a->dimensions[2]; j++) *tmp++ = (double *)(a->data + i*a->strides[1] + j*a->strides[2]); break; case 1: for (i=0; i<a->dimensions[0]; i++) for (j=0; j<a->dimensions[2]; j++) *tmp++ = (double *)(a->data + i*a->strides[0] + j*a->strides[2]); break; case 2: for (i=0; i<a->dimensions[0]; i++) for (j=0; j<a->dimensions[1]; j++) *tmp++ = (double *)(a->data + i*a->strides[0] + j*a->strides[1]); } break; /* four or more dimensions: use recursion */ default: for (i=0; i<a->nd; i++) indices[i] = 0; _get_axis_pointer(indices, axis, a, out, 0); } done: return out; } /* applies a linear filter along one dimension */ static void lfilter_1d( double *restrict b, double *restrict a, double *x, double *y, double *zi, integer K, integer N) { double Z[K]; double *restrict pz = zi, *restrict z = Z; register integer n, k; for (n=0; n<N; n++) { register double yn, xn = x[n]; yn = y[n] = b[0] * xn + pz[0]; for (k=0; k<K-1; k++) z[k] = b[k+1] * xn + pz[k+1] - a[k+1]*yn; z[K-1] = b[K]*xn - a[K]*yn; double *tmp = pz; pz = z; z = tmp; } if (pz != zi) memcpy(zi, pz, K*sizeof(double)); } static void rev_lfilter_1d( double *restrict b, double *restrict a, double *restrict x, double *restrict y, double *zi, integer K, integer N) { double Z[K]; double *restrict pz = zi, *restrict z = Z; register integer n, k; for (n=N-1; n >= 0; n--) { register double yn, xn = x[n]; yn = y[n] = b[0] * xn + pz[0]; for (k=0; k<K-1; k++) z[k] = b[k+1]*xn + pz[k+1] - a[k+1]*yn; z[K-1] = b[K]*xn - a[K]*yn; double *tmp = pz; pz = z; z = tmp; } if (pz != zi) memcpy(zi, pz, K*sizeof(double)); } typedef void (*lfilterfun_t)(double *restrict b, double *restrict a, double *restrict x, double *restrict y, double *zi, integer K, integer N); int _linear_filter(ndarray *a, ndarray *b, ndarray *x, ndarray *y, ndarray *z, int axis, int direction) { int retval = 0; double *wx, *wy, *wz, *data_a, *data_b; integer xcount, ycount, zcount; integer xstride, ystride, zstride; double *xdata, *ydata, *zdata; double **x_axis_ptrs = NULL, **y_axis_ptrs = NULL, **z_axis_ptrs = NULL; lfilterfun_t lfilter; integer i, K, N; /* release the GIL */ Py_BEGIN_ALLOW_THREADS /* sanity check on a and b (these should never fail) */ if (a->nd != 1) goto error_ab; if (b->nd != 1) goto error_ab; if (b->dimensions[0] != a->dimensions[0]) goto error_ab; data_a = (double *)a->data; data_b = (double *)b->data; K = a->dimensions[0] - 1; /* sanity ckeck on x, y and z */ if ((x->nd != y->nd) || (y->nd != z->nd)) goto error_xyz; for (int d=0; d < a->nd; d++) { if (d == axis) continue; if ((x->dimensions[d] != y->dimensions[d]) || (y->dimensions[d] != z->dimensions[d])) goto error_xyz; } N = x->dimensions[axis]; /* extract axis pointers */ x_axis_ptrs = get_axis_pointers(x, axis, &xcount); y_axis_ptrs = get_axis_pointers(y, axis, &ycount); z_axis_ptrs = get_axis_pointers(z, axis, &zcount); /* sanity check */ if (!x_axis_ptrs) goto error_malloc; if (!y_axis_ptrs) goto error_malloc; if (!z_axis_ptrs) goto error_malloc; if ((xcount != ycount)||(ycount != zcount)) goto error_xyz; /* all ok, we can start ... */ lfilter = direction < 0 ? rev_lfilter_1d : lfilter_1d; xstride = x->strides[axis]; ystride = y->strides[axis]; zstride = z->strides[axis]; #pragma omp parallel \ firstprivate(lfilter, N, K, xcount, x_axis_ptrs, y_axis_ptrs, z_axis_ptrs, \ xstride, ystride, zstride, data_b, data_a) \ private(i, xdata, ydata, zdata, wx, wy, wz) \ shared(retval) \ default(none) { if ((xstride==sizeof(double)) && (ystride==sizeof(double)) && (xstride==sizeof(double))) { #pragma omp for schedule(static) for (i=0; i < xcount; i++) { xdata = x_axis_ptrs[i]; ydata = y_axis_ptrs[i]; zdata = z_axis_ptrs[i]; /* strictly speaking we could be aliasing xdata and ydata here, but it should not matter. */ lfilter(data_b, data_a, xdata, ydata, zdata, K, N); } } else { /* axis is strided, use copy-in copy-out */ wx = (double *)malloc((2*N+K)*sizeof(double)); if (!wx) { #pragma omp critical retval = -1; /* error_malloc */ } #pragma omp barrier #pragma omp flush(retval) if (retval < 0) /* malloc failed in some thread */ goto filtering_complete; wy = wx + N; wz = wy + N; #pragma omp for schedule(static) for (i=0; i < xcount; i++) { /* get pointe to the first element of the axis */ xdata = x_axis_ptrs[i]; ydata = y_axis_ptrs[i]; zdata = z_axis_ptrs[i]; /* copy to local contiguous buffers */ copy_in(wx, xdata, N, xstride); copy_in(wz, zdata, K, zstride); /* filter */ lfilter(data_b, data_a, wx, wy, wz, K, N); /* copy data back */ copy_out(ydata, wy, N, ystride); copy_out(zdata, wz, K, zstride); } filtering_complete: if (wx) free(wx); } } /* we are done... */ goto done; /* error conditions */ error_ab: retval = -3; goto done; error_xyz: retval = -2; goto done; error_malloc: retval = -1; /* clean up and exit */ done: if (x_axis_ptrs) free(x_axis_ptrs); if (y_axis_ptrs) free(y_axis_ptrs); if (z_axis_ptrs) free(z_axis_ptrs); /* get the GIL back and return */ Py_END_ALLOW_THREADS return retval; }

6 19

Any advice for building SciPy? LAPACK error
by Kevin Dunn Sept. 24, 2009

Sept. 24, 2009

Hi there, Yesterday I asked about building Scipy (http://article.gmane.org/gmane.comp.python.scientific.devel/11832) as I was following the various instructions * here: http://www.scipy.org/Installing_SciPy/BuildingGeneral * and here: http://www.scipy.org/Installing_SciPy/Linux Since then I discovered the SCons script, so first question: is using the SCons script the preferred way to build SciPy now? If so, I think those website should be updated. If not, please help with the next question. Using the setupscons.py script (after cloning and installing numscons from github, version 0.11.0), I encounter the error, "no lapack found, necessary for lapack module" when running "python setupscons.py build". More details in the traceback below. For reference: * Ubuntu 9.04 (Jaunty) * SciPy revision 5951 * other details about my platform: http://pastebin.com/m26f8d1d1 * and I followed the LAPACK instructions at the bottom of this page http://www.scipy.org/Installing_SciPy/BuildingGeneral (i.e. I built LAPACK from source, and I set the environment variable) Any additional pointers to solve this ? I thought ATLAS is used instead of LAPACK; can someone help me understand this? Thanks, Kevin Executing scons command (pkg is scipy.lib.lapack): /usr/bin/python "/usr/local/lib/python2.6/dist-packages/numscons/scons-local/scons.py" -f scipy/lib/lapack/SConstruct -I. scons_tool_path="" src_dir="scipy/lib/lapack" pkg_name="scipy.lib.lapack" log_level=50 distutils_libdir="../../../../../build/lib.linux-i686-2.6" distutils_clibdir="../../../../../build/temp.linux-i686-2.6" distutils_install_prefix="/usr/local/lib/python2.6/dist-packages/scipy/lib/lapack" cc_opt=gcc cc_opt_path="/usr/bin" debug=0 f77_opt=gfortran f77_opt_path="/usr/bin" cxx_opt=g++ cxx_opt_path="/usr/bin" include_bootstrap=/usr/local/lib/python2.6/dist-packages/numpy/core/include bypass=0 import_env=0 silent=0 bootstrapping=0 scons: Reading SConscript files ... Checking if gfortran needs dummy main - MAIN__. Checking gfortran name mangling - '_', '', lower-case. Checking for ATLAS ... yes Checking for F77 LAPACK ... no RuntimeError: no lapack found, necessary for lapack module: File "/home/connectmv/workspace/scipy/scipy/lib/lapack/SConstruct", line 2: GetInitEnvironment(ARGUMENTS).DistutilsSConscript('SConscript') File "/usr/local/lib/python2.6/dist-packages/numscons/core/numpyenv.py", line 133: build_dir = '$build_dir', src_dir = '$src_dir') File "/usr/local/lib/python2.6/dist-packages/numscons/scons-local/scons-local-1.2.0/SCons/Script/SConscript.py", line 553: return apply(_SConscript, [self.fs,] + files, subst_kw) File "/usr/local/lib/python2.6/dist-packages/numscons/scons-local/scons-local-1.2.0/SCons/Script/SConscript.py", line 262: exec _file_ in call_stack[-1].globals File "/home/connectmv/workspace/scipy/build/scons/scipy/lib/lapack/SConscript", line 37: raise RuntimeError("no lapack found, necessary for lapack module") error: Error while executing scons command. See above for more information. If you think it is a problem in numscons, you can also try executing the scons command with --log-level option for more detailed output of what numscons is doing, for example --log-level=0; the lowest the level is, the more detailed the output it.

2 1

compiling against custom atlas/lapack on os x
by James Kyle Sept. 23, 2009

Sept. 23, 2009

I posted this to scipy-users about a week ago, but after thinking about it the dev list seemed more appropriate. My apologies for this cross post. ---------------------------- Morning, I'm attempting to build scipy against a custom atlas/lapack library on Snow Leopard (this is for inclusion in macports). I've set my build environment according to the documentation here: > http://www.scipy.org/Installing_SciPy/BuildingGeneral Here's my environment and build command: CPATH='/Users/jkyle/Projects/macports/mports/include' LAPACK='/Users/ jkyle/Projects/macports/mports/lib' LIBRARY_PATH='/Users/jkyle/ Projects/macports/mports/lib' MACOSX_DEPLOYMENT_TARGET='10.6' ATLAS='/ Users/jkyle/Projects/macports/mports/lib' CC='/Users/jkyle/Projects/ macports/mports/bin/gcc-mp-4.3' BLAS='/Users/jkyle/Projects/macports/ mports/lib' CXX='/Users/jkyle/Projects/macports/mports/bin/g++-mp-4.3' CCFLAGS='-I/Users/jkyle/Projects/macports/mports/include -L/Users/ jkyle/Projects/macports/mports/lib' python2.6 setup.py --no-user-cfg config_fc --fcompiler gnu95 -- f77exec /Users/jkyle/Projects/macports/mports/bin/gfortran-mp-4.3 -- f90exec /Users/jkyle/Projects/macports/mports/bin/gfortran-mp-4.3 build And here is the configuration output: http://pastie.org/622884 Now, everything seems to go ok, but when I check the resulting libraries I noticed that they are still linking against the system vecLib and Accelerate frameworks (the system os x lapack/blas) when building the cblas.so and fblas.so libraries. The build also seems to create two sets of the cblas and fblas libraries. For example, it installs these linked against my libs: > /Users/jkyle/Projects/macports/mports/Library/Frameworks/ > Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/ > linalg/cblas.so: > /Users/jkyle/Projects/macports/mports/lib/liblapack.dylib > (compatibility version 0.0.0, current version 0.0.0) > /Users/jkyle/Projects/macports/mports/lib/libptf77blas.dylib > (compatibility version 0.0.0, current version 0.0.0) > /Users/jkyle/Projects/macports/mports/lib/libptcblas.dylib > (compatibility version 0.0.0, current version 0.0.0) > /Users/jkyle/Projects/macports/mports/lib/libatlas.dylib > (compatibility version 0.0.0, current version 0.0.0) > /Users/jkyle/Projects/macports/mports/lib/gcc43/libgfortran.3.dylib > (compatibility version 4.0.0, current version 4.0.0) > /Users/jkyle/Projects/macports/mports/lib/gcc43/libgcc_s.1.dylib > (compatibility version 1.0.0, current version 1.0.0) > /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current > version 124.1.1) > > /Users/jkyle/Projects/macports/mports/Library/Frameworks/ > Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/ > linalg/fblas.so: > /Users/jkyle/Projects/macports/mports/lib/liblapack.dylib > (compatibility version 0.0.0, current version 0.0.0) > /Users/jkyle/Projects/macports/mports/lib/libptf77blas.dylib > (compatibility version 0.0.0, current version 0.0.0) > /Users/jkyle/Projects/macports/mports/lib/libptcblas.dylib > (compatibility version 0.0.0, current version 0.0.0) > /Users/jkyle/Projects/macports/mports/lib/libatlas.dylib > (compatibility version 0.0.0, current version 0.0.0) > /Users/jkyle/Projects/macports/mports/lib/gcc43/libgfortran.3.dylib > (compatibility version 4.0.0, current version 4.0.0) > /Users/jkyle/Projects/macports/mports/lib/gcc43/libgcc_s.1.dylib > (compatibility version 1.0.0, current version 1.0.0) > /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current > version 124.1.1) And it installs these linked against the system libs: > > /Users/jkyle/Projects/macports/mports/Library/Frameworks/ > Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/lib/ > blas/cblas.so: > /System/Library/Frameworks/Accelerate.framework/Versions/A/ > Accelerate (compatibility version 1.0.0, current version 4.0.0) > /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current > version 124.1.1) > /Users/jkyle/Projects/macports/mports/Library/Frameworks/ > Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/lib/ > blas/fblas.so: > /Users/jkyle/Projects/macports/mports/lib/gcc43/libgfortran.3.dylib > (compatibility version 4.0.0, current version 4.0.0) > /System/Library/Frameworks/Accelerate.framework/Versions/A/ > Accelerate (compatibility version 1.0.0, current version 4.0.0) > /Users/jkyle/Projects/macports/mports/lib/gcc43/libgcc_s.1.dylib > (compatibility version 1.0.0, current version 1.0.0) > /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current > version 124.1.1) Lastly, I'm noticing that scipy is building with mixed compilers. I've set both CC and CXX in my build environment to 4.3. However, in the build logs I see that scipy is using gcc-4.2 for some things. For example: > creating build/temp.macosx-10.6-i386-2.6/scipy/cluster > creating build/temp.macosx-10.6-i386-2.6/scipy/cluster/src > compile options: '-I/Users/jkyle/Projects/macports/mports/Library/ > Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/ > numpy/core/include -I/Users/jkyle/Projects/macports/mports/Library/ > Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/ > numpy/core/include -I/Users/jkyle/Projects/macports/mports/Library/ > Frameworks/Python.framework/Versions/2.6/include/python2.6 -c' > gcc-mp-4.3: scipy/cluster/src/vq_module.c > gcc-mp-4.3: scipy/cluster/src/vq.c > /Users/jkyle/Projects/macports/mports/Library/Frameworks/ > Python.framework/Versions/2.6/lib/python2.6/site-packages/numpy/core/ > include/numpy/__multiarray_api.h:969: warning: '_import_array' > defined but not used > /usr/bin/gcc-4.2 -L/Users/jkyle/Projects/macports/mports/lib -bundle > -undefined dynamic_lookup build/temp.macosx-10.6-i386-2.6/scipy/ > cluster/src/vq_module.o build/temp.macosx-10.6-i386-2.6/scipy/ > cluster/src/vq.o -Lbuild/temp.macosx-10.6-i386-2.6 -o build/ > lib.macosx-10.6-i386-2.6/scipy/cluster/_vq.so > building 'scipy.cluster._hierarchy_wrap' extension > compiling C sources > C compiler: /Users/jkyle/Projects/macports/mports/bin/gcc-mp-4.3 - > fno-strict-aliasing -fno-common -dynamic -DNDEBUG -g -fwrapv -O3 - > Wall -Wstrict-prototypes Here, the C compiler is gcc-mp-4.3, but the one actually used is gcc-4.2. Any tips are very welcome. Cheers, -james

3 2

scipy.weave crashes wth KeyError
by Thomas Robitaille Sept. 23, 2009

Sept. 23, 2009

Hi, I'm using a recent svn revision of scipy (5925). After installing it I went to scipy/weave/examples and ran 'python array3d.py'. I get the following error message (below). Can other people reproduce this problem? If not, maybe it's some local installation issue. This problem occurs with any code that uses weave.inline. I'm using python 2.6 64-bit, but I've also seen reports of the problem on python 2.6 32-bit. I've tracked the problem down to the following file being empty .python26_compiled/darwin26compiled_catalog.dat However, it does look like there are some compiled shared object libraries in there air:.python26_compiled tom$ du -sk * 0 darwin26compiled_catalog.dat 192 m0 344 m1 1908 m2 580 m3 So it's not clear if compilation is the problem? Thanks, Thomas --- numpy: [[[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] [[12 13 14 15] [16 17 18 19] [20 21 22 23]]] Pure Inline: Traceback (most recent call last): File "array3d.py", line 105, in <module> main() File "array3d.py", line 98, in main pure_inline(arr) File "array3d.py", line 57, in pure_inline weave.inline(code, ['arr']) File "/Users/tom/Library/Python/2.6/site-packages/scipy/weave/inline_tools.py", line 324, in inline results = attempt_function_call(code,local_dict,global_dict) File "/Users/tom/Library/Python/2.6/site-packages/scipy/weave/inline_tools.py", line 392, in attempt_function_call function_list = function_catalog.get_functions(code,module_dir) File "/Users/tom/Library/Python/2.6/site-packages/scipy/weave/catalog.py", line 615, in get_functions function_list = self.get_cataloged_functions(code) File "/Users/tom/Library/Python/2.6/site-packages/scipy/weave/catalog.py", line 529, in get_cataloged_functions if cat is not None and code in cat: File "/System/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/shelve.py", line 110, in __contains__ return key in self.dict File "/Users/tom/Library/Python/2.6/site-packages/scipy/io/dumbdbm_patched.py", line 73, in __getitem__ pos, siz = self._index[key] # may raise KeyError KeyError: 0

1 0

Re: [SciPy-dev] Bug/Error with chi-squared distribution and df<1
by Kevin Dunn Sept. 23, 2009

Sept. 23, 2009

> ti, 2009-09-22 kello 14:27 -0400, josef.pktd(a)gmail.com kirjoitti: > [clip] >> Taking the hint with incomplete gamma, the following looks good. This >> would mean until Pauli fixes scipy.special if it your fix works, we >> could also use gammainc directly. I don't know the differences between >> the various implementations well enough to see whether we buy some >> other problems with this > > The C function `igam` called by `chdtr` is the same as what is exposed > as `gammainc`. Thanks again to all the helpful replies on this list. I'm using this relationship between the chi-squared distribution and the gamma function to evaluate the CDF and inverse CDF. As an example for future readers of this list, let our value be x = 3 with 0.6 degrees of freedom, Chi-square CDF ---------------- * In MATLAB: chi2cdf(3, 0.6) % 0.957890536704110 * In R: pchisq(3, 0.6) # 0.9578905 * And temporarily in SciPy: scipy.special.gammainc(0.6/2.0, 3/2.0) # note the order of parameters is flipped; result = 0.95789053670411062 * Or another temporary possibility: scipy.stats.gamma.cdf(3/2.0, 0.6/2.0) # this time the parameters are in the expected order But ideally in SciPy: scipy.stats.chi2.cdf(3, 0.6) # currently returns a NaN, but should return the same result above Inverse Chi-square CDF ---------------- * In MATLAB: chi2inv(0.95789, 0.6) % 2.999980885343446 * In R: qchisq(0.95789, 0.6) # 2.999981 * And temporarily in SciPy: scipy.special.gammaincinv(0.6/2.0, 0.95789) * 2.0 # again we reorder the parameters; result = 2.9999808853434109 * Or another temporary possibility: scipy.stats.gamma.ppf(0.95789, 0.6/2.0)*2 # 2.9999808853434109 But ideally in SciPy: scipy.stats.chi2.ppf(0.95789, 0.6) # currently returns a NaN, but should return the same result above Hope that helps in the mean time, Kevin > -- > Pauli Virtanen

2 1

Bug/Error with chi-squared distribution and df<1
by Kevin Dunn Sept. 22, 2009

Sept. 22, 2009

Thanks for the prompt reply Josef; please see my comments below. > On Tue, Sep 22, 2009 at 10:59 AM, Kevin Dunn <kgdunn(a)gmail.com> wrote: >> Hi there, >> >> I'm not an expert on distributions, but as far as I can tell, the chi2 >> distribution is defined for degrees of freedom >0. ?I'm getting "nan" >> results however when df is less than one (but still greater than 0). >> The chi2 CDF values agree between R, MATLAB and Scipy when the the >> degrees of freedom are >= 1. ?For example: >> * R: pchisq(0.95, 1) >> * MATLAB: chi2cdf(0.95, 1) >> * SciPy: scipy.stats.chi2.cdf(0.95, 1) >> > > >> However, changing the last argument to 0.5 returns a NaN in SciPy, but >> gives a result (0.8392259) in R and MATLAB. > >>>> stats.chi2.veccdf(0.95, 0.5) > array(0.83922587961194761) > >> >> I'm suspecting there is something wrong with my SciPy installation, >> because the example included with SciPy, (found using >> scipy.stats.chi2? inside iPython) calls the chi2.cdf function with >> df=0.9, yet when I run that example as follows: > > Where is this example? maybe it is an autogenerated example with wrong numbers? The example was found when typing "scipy.stats.chi2?" in ipython after importing scipy.stats. The code below was taken from the example, only I omitted the matplotlib plotting commands. >> from scipy.stats import chi2 >> import numpy as np >> numargs = chi2.numargs >> [ df ] = [0.9,]*numargs >> rv = chi2(df) >> x = np.linspace(0,np.minimum(rv.dist.b,3)) >> prb = chi2.cdf(x,df) >> >> My result for prb is as follows; which I don't think would have been >> used as an example if this is the expected output. >> array([ ?0., ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, >> ? ? ? ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, >> ? ? ? ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, >> ? ? ? ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, >> ? ? ? ?NaN, ?NaN, ?NaN, ?NaN, ?NaN, ?NaN]) >> >> Is my SciPy install messed up? ?I'm using MacOSX and downloaded >> version 0.7.1 from SourceForge this morning. ?I just tried it in >> Ubuntu Linux (version 0.7.0 though), and get the same results. >> > > chi2 uses scipy.special > > class chi2_gen > > def _cdf(self, x, df): > return special.chdtr(df, x) > > which obviously cannot handle df<1 , which I don't know if it would > ever show up in the usual statistical tests. I don't believe it is obvious why special.chdtr can't handle values df < 1. Please see my additional comments below. > > pdf doesn't seem to have any problems > >>>> stats.chi2.pdf(np.arange(5),0.99) > array([ Inf, 0.24042373, 0.10275665, 0.05078513, 0.02663761]) >>>> stats.chi2.pdf(np.arange(5),0.5) > array([ Inf, 0.14067411, 0.05073346, 0.02270277, 0.01109756]) >>>> stats.chi2.pdf(np.arange(5),0.25) > array([ Inf, 0.07382471, 0.02441481, 0.01038547, 0.00489731]) > > so numerical integration should also work > >>>> stats.chi2.veccdf(np.arange(1,5),0.25) > array([ 0.92605422, 0.96918041, 0.98539847, 0.99265551]) >>>> stats.chi2.veccdf(np.linspace(0.01,10.0,11),0.25) > array([ 0.54726537, 0.92671456, 0.96937499, 0.98547096, 0.99268483, > 0.99618434, 0.99796201, 0.99889274, 0.99939058, 0.99966116, > 0.99981006]) >>>> stats.chi2.veccdf(np.linspace(0.01,10.0,11),0.5) > array([ 0.29308089, 0.84774539, 0.93248332, 0.96674206, 0.98278044, > 0.99081467, 0.99500114, 0.99723983, 0.9984591 , 0.99913231, > 0.99950797]) > > since pdf at zero is inf, I don't know how good the numbers are if the > cdf is calculated for points close to zero > e.g. > >>>> stats.chi2.veccdf(1e-8,0.5) > array(0.0092772960765327081) > > > Since I don't think df<1 is a common case, I wouldn't want to switch > to numerical integration by default. But veccdf although only a > private function should work although it bypasses some of the checking > and broadcasting. Thanks for pointing out the numerical integration approach; 'll give it a try. As for df < 1, please see my other comments below. > > Somehow this sounds familiar, I need to check whether this or a > similar case is already on record. > > Hope that helps, > > Josef And regarding your other reply: > I didn't find anything on the scipy trac nor on the mailinglist, since > I'm subscribed to them. If you think that df<1 is an important case > you could file a ticket for scipy.special. I have no idea how > difficult it would be to extend the chi squared related functions for > this or whether there exists another workaround. I just downloaded the latest SVN code and found the code that does the work is /scipy/special/cephes/chdtr.c In that code it seems that when the degrees of freedom are less than one, the code is (artificially) forced to return a NaN: if (df < 1.0) { mtherr( "chdtrc", DOMAIN ); return(NPY_NAN); } return( igamc( df/2.0, x/2.0 ) ); This seems a somewhat artificial constraint to me, since the gamma function can accept values of 0 < df < 1. Does someone else on this list know why that constraint is there for df<1? As for a practical usage case, I'm computing CDF values for a variable by matching moments, and the value for the degrees of freedom term is a computed value. This value is always non-integer, and can sometimes be less than one. Additional justification for df < 1 can be sought by looking at plots of the chi-squared distribution (e.g. http://en.wikipedia.org/wiki/Chi-square_distribution) when plotting with a degrees of freedom parameter (they call it "k" in the Wikipedia article). There's no reason why df can't be less than one, anymore than there's reason for df being non-integer. I'll give it a go compiling SciPy with removing those constraints in chdtr.c and see how it works. I've never done this before, but I'm busy installing the NumPy/SciPy installation requirements at the moment, and I'll report how it went. If things works out, I'll file a trac report. I've also cc'd the scipy-dev list on this reply. Thanks for your help, Kevin > The only similar story was for the negative binomial for n<1, > http://projects.scipy.org/scipy/ticket/978 . In that case > scipy.special.nbdtr didn't allow the extension for n<1. > > Josef

3 2

Building SciPy: breaks when looking for get_info function
by Kevin Dunn Sept. 22, 2009

Sept. 22, 2009

Hi there, This is my first attempt at building SciPy; I'm doing this to solve a problem described in another posting here today with subject "Bug/Error with chi-squared distribution and df<1" Background: I'm using Ubuntu 9.04, I've built and installed Numpy successfully from SVN this afternoon. Other details about my platform, as requested by INSTALL.TXT, are available at http://pastebin.com/m26f8d1d1 Now I'm trying to build and install SciPy from SVN as well, but I get the following traceback regarding the `get_info` function when I call: python setup.py build When inspecting the ./scipy/special/setup.py file I notice it is trying to import get_info, however, `get_info` is from /numpy/distutils/misc_util.py where it exists only as a method of the Configuration class, not as a stand-alone function. I'm not sure how to proceed from here. Any advice? Thanks, Kevin ---- Traceback (most recent call last): File "setup.py", line 160, in <module> setup_package() File "setup.py", line 152, in setup_package configuration=configuration ) File "/usr/lib/python2.6/dist-packages/numpy/distutils/core.py", line 122, in setup return setup(**attr) File "/usr/lib/python2.6/dist-packages/numpy/distutils/core.py", line 150, in setup config = configuration() File "setup.py", line 118, in configuration config.add_subpackage('scipy') File "/usr/lib/python2.6/dist-packages/numpy/distutils/misc_util.py", line 851, in add_subpackage caller_level = 2) File "/usr/lib/python2.6/dist-packages/numpy/distutils/misc_util.py", line 834, in get_subpackage caller_level = caller_level + 1) File "/usr/lib/python2.6/dist-packages/numpy/distutils/misc_util.py", line 781, in _get_configuration_from_setup_py config = setup_module.configuration(*args) File "scipy/setup.py", line 20, in configuration config.add_subpackage('special') File "/usr/lib/python2.6/dist-packages/numpy/distutils/misc_util.py", line 851, in add_subpackage caller_level = 2) File "/usr/lib/python2.6/dist-packages/numpy/distutils/misc_util.py", line 834, in get_subpackage caller_level = caller_level + 1) File "/usr/lib/python2.6/dist-packages/numpy/distutils/misc_util.py", line 766, in _get_configuration_from_setup_py ('.py', 'U', 1)) File "scipy/special/setup.py", line 7, in <module> from numpy.distutils.misc_util import get_numpy_include_dirs, get_info ImportError: cannot import name get_info

1 0

Gauss-Lobatto points
by Ondrej Certik Sept. 22, 2009

Sept. 22, 2009

Hi, does anyone know a way to get Gauss-Lobatto points out of scipy? I know how to get Gauss integration points, but that's different: In [1]: from scipy.special.orthogonal import p_roots In [2]: p_roots(2) Out[2]: [array([-0.57735027, 0.57735027]), array([ 1., 1.])] In [3]: p_roots(3) Out[3]: [array([ -7.74596669e-01, -4.78946310e-17, 7.74596669e-01]), array([ 0.55555556, 0.88888889, 0.55555556])] In [4]: p_roots(4) Out[4]: [array([-0.86113631, -0.33998104, 0.33998104, 0.86113631]), array([ 0.34785485, 0.65214515, 0.65214515, 0.34785485])] In [5]: p_roots(5) Out[5]: [array([ -9.06179846e-01, -5.38469310e-01, -1.25197490e-16, 5.38469310e-01, 9.06179846e-01]), array([ 0.23692689, 0.47862867, 0.56888889, 0.47862867, 0.23692689])] the Gauss-Lobatto points are used for interpolation. Thanks! Ondrej

2 3

scipy trunk build issue
by Christopher Hanley Sept. 22, 2009

Sept. 22, 2009

Hi, When attempting to build scipy from the trunk I receive the following error: building extension "scipy.integrate.dop" sources error: mismatch of extension names: scipy/integrate/dop.pyf provides '_dop' but expected 'dop' I am running Python 2.5 on my fully patched MacBook Pro OS X 10.5 system. Thanks, Chris -- Christopher Hanley Senior Systems Software Engineer Space Telescope Science Institute 3700 San Martin Drive Baltimore MD, 21218 (410) 338-4338

2 2

PATCH: Update blitz++ to upstream CVS
by Norbert Nemec Sept. 21, 2009

Sept. 21, 2009

Hi there, the version of blitz++ in the repository is on the status of 2005 with just a number of patches applied. Here is a patch that updates the version to the latest CVS version, containing a large pile of bugfixes and minor extensions. I checked carefully: all changes that happened in the scipy version have been merged upstream already. Apart from the headers, this patch makes both versions identical. Is there anything I might have overlooked? If not - please, could someone with access commit this patch to the SVN? Greetings, Norbert

1 2