Mailman 3 September 2010 - SciPy-Dev

SLSQP Constrained Optimizer Status
by Rob Falck Nov. 24, 2013

Nov. 24, 2013

I'm currently implementing the Sequential Least Squares Quadratic Programming (SLSQP) optimizer, by Dieter Kraft, for use in Scipy. The Fortran code being wrapped with F2PY is here: http://www.netlib.org/toms/733 (its use within Scipy has been cleared) SLSQP provides for bounds on the independent variables, as well as equality and inequality constraint functions, which is a capability that doesn't exist in scipy.optimize. Currently the code works, although the constraint normals are being generated by approximation. I'm working on a way to pass these in. I think the most elegant way will be a single function that returns the matrix of constraint normals. For a demonstration of what the code can do, here is an optimization of f(x,y) = 2xy + 2x - x**2 - 2y**2 Example 14.2 in Chapra & Canale gives the maximum as x=2.0, y=1.0. The unbounded optimization tests find this solution. As expected, its faster when derivatives are provided rather than approximated. Unbounded optimization. Derivatives approximated. Elapsed time: 1.45792961121 ms Results [[1.9999999515712266, 0.99999996181577444], -1.9999999999999984, 4, 0, 'Optimization terminated successfully.'] Unbounded optimization. Derivatives provided. Elapsed time: 1.03211402893 ms Results [[1.9999999515712266, 0.99999996181577444], -1.9999999999999984, 4, 0, 'Optimization terminated successfully.'] The following example uses an equality constraint to find the optimal when x=y. Bound optimization. Derivatives approximated. Elapsed time: 1.384973526 ms Results [[0.99999996845920858, 0.99999996845920858], -0.99999999999999889, 4, 0, 'Optimization terminated successfully.'] I've tried to conform to the syntax used by the other optimizers in scipy.optimize. The function definition and doc string are below. If anyone is interested in testing it out, let me know. def fmin_slsqp( func, x0 , eqcons=[], f_eqcons=None, ieqcons=[], f_ieqcons=None, bounds = [], fprime = None, fprime_cons=None,args = (), iter = 100, acc = 1.0E-6, iprint = 1, full_output = 0, epsilon = _epsilon ): """ Minimize a function using Sequential Least SQuares Programming Description: Python interface function for the SLSQP Optimization subroutine originally implemented by Dieter Kraft. Inputs: func - Objective function (in the form func(x, *args)) x0 - Initial guess for the independent variable(s). eqcons - A list of functions of length n such that eqcons[j](x0,*args) == 0.0 in a successfully optimized problem f_eqcons - A function of the form f_eqcons(x, *args) that returns an array in which each element must equal 0.0 in a successfully optimized problem. If f_eqcons is specified, eqcons is ignored. ieqcons - A list of functions of length n such that ieqcons[j](x0,*args) >= 0.0 in a successfully optimized problem f_ieqcons - A function of the form f_ieqcons(x0, *args) that returns an array in which each element must be greater or equal to 0.0 in a successfully optimized problem. If f_ieqcons is specified, ieqcons is ignored. bounds - A list of tuples specifying the lower and upper bound for each independent variable [ (xl0, xu0), (xl1, xu1), ...] fprime - A function that evaluates the partial derivatives of func fprime_cons - A function of the form f(x, *args) that returns the m by n array of constraint normals. If not provided, the normals will be approximated. Equality constraint normals precede inequality constraint normals. args - A sequence of additional arguments passed to func and fprime iter - The maximum number of iterations (int) acc - Requested accuracy (float) iprint - The verbosity of fmin_slsqp. iprint <= 0 : Silent operation iprint == 1 : Print summary upon completion (default) iprint >= 2 : Print status of each iterate and summary full_output - If 0, return only the minimizer of func (default). Otherwise, output final objective function and summary information. epsilon - The step size for finite-difference derivative estimates. Outputs: ( x, { fx, gnorm, its, imode, smode }) x - The final minimizer of func. fx - The final value of the objective function. its - The number of iterations. imode - The exit mode from the optimizer, as an integer. smode - A string describing the exit mode from the optimizer. Exit modes are defined as follows: -1 : Gradient evaluation required (g & a) 0 : Optimization terminated successfully. 1 : Function evaluation required (f & c) 2 : Number of equality constraints larger than number of independent variables 3 : More than 3*n iterations in LSQ subproblem 4 : Inequality constraints incompatible 5 : Singular matrix E in LSQ subproblem 6 : Singular matrix C in LSQ subproblem 7 : Rank-deficient equality constraint subproblem HFTI 8 : Positive directional derivative for linesearch 9 : Iteration limit exceeded -- - Rob Falck

4 16

Subversion scipy.stats irregular problem with source code example
by James Phillips Dec. 9, 2010

Dec. 9, 2010

Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered. Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like: Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line: Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/optimize/minpack.py", line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func and program flow is stopped. In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code. James Phillips ######################################################## import numpy, scipy, scipy.stats # test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data) # test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data) print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print

7 42

Improvements to SciPy.org Server
by Aaron River Sept. 29, 2010

Sept. 29, 2010

Greetings Gentle Folk, I'm happy to announce that the server hosting SciPy.org, as well as other OSS projects, has been migrated into a cluster assigned to our community support efforts. The base instance itself remains relatively unchanged, but is now running atop a more capable set of servers. (Improvements to the services themselves are next in line.) Yesterday afternoon, I allowed the system to run with the configurations carried over from the old hardware, while watching the services for warning signs of trouble. Then, that evening, I loosened things up and adjusted configurations to better utilize the available resources. Though I always expect there to be "gotchas", I've thus far observed much improved response times and throughput. If you do experience any problems, please do not hesitate to let us know. We hope you enjoy the improvements to SciPy.org. Thanks, -- Aaron River Sr IT Administrator Enthought, Inc. --

1 0

Re: [SciPy-Dev] spatial.test() failure
by Warren Weckesser Sept. 28, 2010

Sept. 28, 2010

> > Nils Wagner wrote: > >>>> from scipy import spatial > >>>> spatial.test() > >>>> > > Running unit tests for scipy.spatial > > NumPy version 2.0.0.dev8714 > > NumPy is installed in > > /home/nwagner/local/lib64/python2.6/site-packages/numpy > > SciPy version 0.9.0.dev6803 > > SciPy is installed in > > /home/nwagner/local/lib64/python2.6/site-packages/scipy > > Python version 2.6.5 (r265:79063, Jul 5 2010, 11:46:13) > > [GCC 4.5.0 20100604 [gcc-4_5-branch revision 160292]] > > nose version 0.11.2 > > > ............................................................................................................................F.................................................................................................................................................. > > ====================================================================== > > FAIL: Tests pdist(X, 'minkowski') on iris data. (float32) > > ---------------------------------------------------------------------- > > Traceback (most recent call last): > > File > > > "/home/nwagner/local/lib64/python2.6/site-packages/scipy/spatial/tests/test_distance.py", > > line 837, in test_pdist_minkowski_3_2_iris_float32 > > self.failUnless(within_tol(Y_test1, Y_right, eps)) > > AssertionError > > > > ---------------------------------------------------------------------- > > Ran 271 tests in 38.005s > > > > FAILED (failures=1) > > <nose.result.TextTestResult run=271 errors=0 failures=1> > > > > > In r6799, the names of some tests were changed so there are no longer > duplicated test function names. This means tests are now being run that > had previously been shadowed by the duplicate names. I also get the > failure of the test test_pdist_minkowski_3_2_iris_float32(). I created > a ticket for the problem: http://projects.scipy.org/scipy/ticket/1278 > > Warren > If the tolerance in the test is changed from eps=1e-7 to 1e-6, the test passes on my computer. Several other float32 tests using the iris data also use 1e-6. I added a note about the tolerance of similar tests in the ticket. It does not seem unreasonable to simply loosen the tolerance a bit. Any objections? Warren

3 4

scipy.optimize.curve_fit does not work with single parameter function
by Daniele Nicolodi Sept. 28, 2010

Sept. 28, 2010

Hello, I copied the code of scipy.optimize.curve_fit into my data analysis code to use it on a not so recent versions of scipy available on my systems. However it does not work when I try to fit a simple function that has only one parameter, for example: def func(x, p1): return x * p1 The problem is that scipy.optimize.leastsq() returns a list of parameters values for cases where the function has more than one parameter it returns a tuple of parameters values, while in the case of one single parameters it returns a numpy scalar value. This causes an error in the curve_fit _general_function() and _weighted_general_function() where the target function is called: function(xdata, *params) This works only if params is a python list or tuple. Also the return value suffers from the same problem. If the target function has more than one parameter you get a parameters list, otherwise you get a scalar. It would be more consistent to return always a list, also because the covariance is always returned as a matrix. Please ignore this if that has been resolved in more recent version of scipy, I have acces to 0.7.2 (last available on debian sid). Cheers, -- Daniele

2 1

problems configuring scipy (0.8.0) from numpy (1.5.0) in a itanium cluster
by Miguel Pan Fidalgo Sept. 28, 2010

Sept. 28, 2010

Technical information ++++++++++++++++++++++ $ python -c 'import sys;print sys.version' 2.4.2 (#1, Apr 13 2007, 16:14:04) [GCC 4.1.0 (SUSE Linux)] $ module load numpy/1.5.0 $ module list Currently Loaded Modulefiles: 1) icc/11.1.056 3) mkl/10.2.2 5) nose/0.11.2 7) numpy/1.5.0 2) ifort/11.1.056 4) fftw/3.2.2 6) SuiteSparse/3.4.0 $ python -c 'import numpy;print numpy.__version__' 1.5.0 $ uname -a Linux cn002 2.6.16.53-0.8_SFS2.3_0-default #1 SMP Mon Nov 5 14:15:47 CET 2007 ia64 ia64 ia64 GNU/Linux numpy-1.5.0 +++++++++++++++++ I have installed numpy-1.5.0 in our Itanium cluster. This is my site.cfg file: $ cat /opt/cesga/numpy-1.5.0/src_numpy-1.5.0/site.cfg [DEFAULT] library_dirs = /opt/cesga/fftw-3.2.2/lib:/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64:/opt/cesga/SuiteSparse-3.4.0/lib include_dirs = /opt/cesga/fftw-3.2.2/include:/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/include:/opt/cesga/SuiteSparse-3.4.0/include # UMFPACK # ------- [amd] amd_libs = amd [umfpack] umfpack_libs = umfpack # FFT libraries # ------------- [fftw] libraries = fftw3 # MKL # ---- [mkl] library_dirs = /opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64 lapack_libs = mkl_lapack mkl_libs = mkl_solver_lp64, mkl_intel_lp64, mkl_intel_thread, mkl_i2p, mkl_core These are the hacks: $ diff system_info.py.original system_info.py 818c818 < _lib_mkl = ['mkl','vml','guide'] --- > _lib_mkl = ['mkl','mkl_lapack','guide'] $ diff intelccompiler.py.original intelccompiler.py 28a29 > cc_exe=cc_exe + ' -liomp5 -lpthread -fPIC' scipy-0.8.0 +++++++++++++++ When I try to install scipy-0.8.0, the configuration process fails: $ tar zxvf scipy-0.8.0.tar.gz $ mv scipy-0.8.0 src $ cp -va run.sh src/ $ cd src/ $ cat run.sh #!/bin/bash module load numpy/1.5.0 python setup.py config_fc --fcompiler=intele #install --prefix=/sfs/home/cesga/mpan/scipy-0.8.0 $ ./run.sh Warning: No configuration returned, assuming unavailable.blas_opt_info: blas_mkl_info: FOUND: libraries = ['mkl_solver_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_i2p', 'mkl_core', 'pthread'] library_dirs = ['/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64'] define_macros = [('SCIPY_MKL_H', None)] include_dirs = ['/opt/cesga/fftw-3.2.2/include', '/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/include', '/opt/cesga/SuiteSparse-3.4.0/include'] FOUND: libraries = ['mkl_solver_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_i2p', 'mkl_core', 'pthread'] library_dirs = ['/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64'] define_macros = [('SCIPY_MKL_H', None)] include_dirs = ['/opt/cesga/fftw-3.2.2/include', '/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/include', '/opt/cesga/SuiteSparse-3.4.0/include'] lapack_opt_info: lapack_mkl_info: mkl_info: FOUND: libraries = ['mkl_solver_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_i2p', 'mkl_core', 'pthread'] library_dirs = ['/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64'] define_macros = [('SCIPY_MKL_H', None)] include_dirs = ['/opt/cesga/fftw-3.2.2/include', '/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/include', '/opt/cesga/SuiteSparse-3.4.0/include'] FOUND: libraries = ['mkl_lapack', 'mkl_solver_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_i2p', 'mkl_core', 'pthread'] library_dirs = ['/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64'] define_macros = [('SCIPY_MKL_H', None)] include_dirs = ['/opt/cesga/fftw-3.2.2/include', '/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/include', '/opt/cesga/SuiteSparse-3.4.0/include'] FOUND: libraries = ['mkl_lapack', 'mkl_solver_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_i2p', 'mkl_core', 'pthread'] library_dirs = ['/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64'] define_macros = [('SCIPY_MKL_H', None)] include_dirs = ['/opt/cesga/fftw-3.2.2/include', '/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/include', '/opt/cesga/SuiteSparse-3.4.0/include'] umfpack_info: libraries umfpack not found in /opt/cesga/fftw-3.2.2/lib libraries umfpack not found in /opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64 amd_info: libraries amd not found in /opt/cesga/fftw-3.2.2/lib libraries amd not found in /opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64 FOUND: libraries = ['amd'] library_dirs = ['/opt/cesga/SuiteSparse-3.4.0/lib'] swig_opts = ['-I/opt/cesga/SuiteSparse-3.4.0/include'] define_macros = [('SCIPY_AMD_H', None)] include_dirs = ['/opt/cesga/SuiteSparse-3.4.0/include'] FOUND: libraries = ['umfpack', 'amd'] library_dirs = ['/opt/cesga/SuiteSparse-3.4.0/lib'] swig_opts = ['-I/opt/cesga/SuiteSparse-3.4.0/include', '- I/opt/cesga/SuiteSparse-3.4.0/include'] define_macros = [('SCIPY_UMFPACK_H', None), ('SCIPY_AMD_H', None)] include_dirs = ['/opt/cesga/SuiteSparse-3.4.0/include'] running config_fc unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options the warning (first line) implies that BLAS primitives have to be recompiled. At the end the installation ends properly... this is not the point, but scipy doesn't use MKL BLAS primitives. BUT according to numpy... blas_mkl and blas_opt should be the same (like lapack_mkl and lapack_opt): $ module load numpy/1.5.0 $ python Python 2.4.2 (#1, Apr 13 2007, 16:14:04) [GCC 4.1.0 (SUSE Linux)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> from numpy.distutils.system_info import get_info >>> get_info('blas_mkl') {'libraries': ['mkl_solver_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_i2p', 'mkl_core', 'pthread'], 'library_dirs': ['/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64'], 'define_macros': [('SCIPY_MKL_H', None)], 'include_dirs': ['/opt/cesga/fftw-3.2.2/include', '/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/include', '/opt/cesga/SuiteSparse-3.4.0/include']} >>> get_info('blas_opt') {'libraries': ['mkl_solver_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_i2p', 'mkl_core', 'pthread'], 'library_dirs': ['/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/lib/64'], 'define_macros': [('SCIPY_MKL_H', None)], 'include_dirs': ['/opt/cesga/fftw-3.2.2/include', '/opt/cesga/intel/intel11.0/Compiler/11.1/056/mkl/include', '/opt/cesga/SuiteSparse-3.4.0/include']} >>> There is no need to recompile BLAS!!!! any idea or suggestion ? WHAT I'm doing wrong ? -- Atentamente, Miguel Pan Fidalgo (Applications Technician) mail: mpan at cesga.es web: http://www.cesga.es Avda. de Vigo s/n 15705, Santiago de Compostela Telf.: +34 981 569810 - Fax: 981 594616 -------------------------------------------------------------------------

1 0

scipy error compiling csr_wrap under Python 2.7
by Bruce Southey Sept. 27, 2010

Sept. 27, 2010

Hi, I failed to get the scipy 0.8 rc1 and SVN to build under Python2.7 with numpy '2.0.0.dev8469' and gcc version 4.4.4 20100503 (Red Hat 4.4.4-2) (GCC) . But scipy does compile with Python 2.6. building 'scipy.sparse.sparsetools._csr' extension compiling C++ sources C compiler: g++ -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -fPIC compile options: '-I/usr/local/lib/python2.7/site-packages/numpy/core/include -I/usr/local/include/python2.7 -c' g++: scipy/sparse/sparsetools/csr_wrap.cxx scipy/sparse/sparsetools/csr_wrap.cxx: In function ‘int require_size(PyArrayObject*, npy_intp*, int)’: scipy/sparse/sparsetools/csr_wrap.cxx:2910: error: expected ‘)’ before ‘PRIdPTR’ scipy/sparse/sparsetools/csr_wrap.cxx:2910: warning: spurious trailing ‘%’ in format scipy/sparse/sparsetools/csr_wrap.cxx:2910: warning: too many arguments for format scipy/sparse/sparsetools/csr_wrap.cxx:2917: error: expected ‘)’ before ‘PRIdPTR’ scipy/sparse/sparsetools/csr_wrap.cxx:2917: warning: spurious trailing ‘%’ in format scipy/sparse/sparsetools/csr_wrap.cxx:2917: warning: too many arguments for format scipy/sparse/sparsetools/csr_wrap.cxx: In function ‘int require_size(PyArrayObject*, npy_intp*, int)’: scipy/sparse/sparsetools/csr_wrap.cxx:2910: error: expected ‘)’ before ‘PRIdPTR’ scipy/sparse/sparsetools/csr_wrap.cxx:2910: warning: spurious trailing ‘%’ in format scipy/sparse/sparsetools/csr_wrap.cxx:2910: warning: too many arguments for format scipy/sparse/sparsetools/csr_wrap.cxx:2917: error: expected ‘)’ before ‘PRIdPTR’ scipy/sparse/sparsetools/csr_wrap.cxx:2917: warning: spurious trailing ‘%’ in format scipy/sparse/sparsetools/csr_wrap.cxx:2917: warning: too many arguments for format error: Command "g++ -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -fPIC -I/usr/local/lib/python2.7/site-packages/numpy/core/include -I/usr/local/include/python2.7 -c scipy/sparse/sparsetools/csr_wrap.cxx -o build/temp.linux-x86_64-2.7/scipy/sparse/sparsetools/csr_wrap.o" failed with exit status 1 The error is possibly related the variable 'NPY_INTP_FMT' (defined in numpy/core/include/numpy/ndarraytypes.h). A grep shows these files: scipy/sparse/sparsetools/coo_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/coo_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/bsr_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/bsr_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/numpy.i: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/numpy.i: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/csr_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/csr_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/csc_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/csc_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/dia_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/dia_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); Note that the scipy/sparse/sparsetools/numpy.i contains tabs in some places instead of spaces. The ndarraytypes.h is present in /usr/local/lib/python2.7/site-packages/numpy/core/include/numpy/ Bruce

5 10

script for trailing whitespace cleanup
by josef.pktd＠gmail.com Sept. 24, 2010

Sept. 24, 2010

Is the script for trailing whitespace cleanup that is used for numpy, scipy available somewhere? Thanks, Josef

3 3

Problem with N-dimensional interpolation using a new griddata function for N>=3
by Adam Machnik Sept. 22, 2010

Sept. 22, 2010

Hello , I have downloaded from SVN the latest `griddata` implementation for working in N-D with irregular grids (scipy v0.9.0dev6812). It works great for 2D problems. However for 3D (and more) I have a problem: the NaN values appears always in the same position, and it does not depend on the scalar field being interpolated. Here is a small example of simple 3x3x3 regular grid and a constant scalar field of 1.0 values, that is interpolated at the input data points. Example: ****************************** from numpy import * from scipy.interpolate import * grid_x, grid_y, grid_z = mgrid[0:1:3j, 0:1:3j, 0:1:3j] pts = [[ 0.0 , 0.0 , 0.0 ], [ 0.0 , 0.0 , 0.5 ], [ 0.0 , 0.0 , 1.0 ], [ 0.0 , 0.5 , 0.0 ], [ 0.0 , 0.5 , 0.5 ], [ 0.0 , 0.5 , 1.0 ], [ 0.0 , 1.0 , 0.0 ], [ 0.0 , 1.0 , 0.5 ], [ 0.0 , 1.0 , 1.0 ], [ 0.5 , 0.0 , 0.0 ], [ 0.5 , 0.0 , 0.5 ], [ 0.5 , 0.0 , 1.0 ], [ 0.5 , 0.5 , 0.0 ], [ 0.5 , 0.5 , 0.5 ], [ 0.5 , 0.5 , 1.0 ], [ 0.5 , 1.0 , 0.0 ], [ 0.5 , 1.0 , 0.5 ], [ 0.5 , 1.0 , 1.0 ], [ 1.0 , 0.0 , 0.0 ], [ 1.0 , 0.0 , 0.5 ], [ 1.0 , 0.0 , 1.0 ], [ 1.0 , 0.5 , 0.0 ], [ 1.0 , 0.5 , 0.5 ], [ 1.0 , 0.5 , 1.0 ], [ 1.0 , 1.0 , 0.0 ], [ 1.0 , 1.0 , 0.5 ], [ 1.0 , 1.0 , 1.0 ]] vls = ones(27) points = array(pts).astype('float64') values = array(vls).astype('float64') grid_z1 = griddata(points, values, (grid_x, grid_y, grid_z), method='linear', fill_value=0) print 'grid_z1=', grid_z1 *********************** I obtain the following result: grid_z1= [[[ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.]] [[ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.]] [[ 1. 1. 1.] [ 1. 1. nan] [ 1. 1. 1.]]] Thanks, Adam

4 11

Problem with manipulating dates in the time series object
by gharras Sept. 22, 2010

Sept. 22, 2010

Hi there, I really like the time series object from scipy but I have a problem with manipulating dates. Here is a copy from my ipython session: ----------------------------------------------------------------------------- In [748]: myts=ts.time_series(data=n.arange(20).reshape(2,10).T, dates=n.arange(10)+720000, freq='D') In [749]: myts Out[749]: timeseries( [[ 0 10] [ 1 11] [ 2 12] [ 3 13] [ 4 14] [ 5 15] [ 6 16] [ 7 17] [ 8 18] [ 9 19]], dates = [17-Apr-1972 ... 26-Apr-1972], freq = D) In [750]: myts.dates[myts.weekday==1]+=1 In [751]: myts.dates Out[751]: DateArray([17-Apr-1972, 19-Apr-1972, 19-Apr-1972, 20-Apr-1972, 21-Apr-1972, 22-Apr-1972, 23-Apr-1972, 24-Apr-1972, 26-Apr-1972, 26-Apr-1972], freq='D') In [752]: myts.has_duplicated_dates() Out[752]: False ----------------------------------------------------------------------------- In words: I can change the dates by hand (in line 750), but the object is not aware of it, i.e. it still thinks that it does not contain any duplicates but the dates were however affected by the manipulation in line 750, cf. line 751. Is there a work-around? all the best. Cheers, Georges

2 1

[ANN] SciPy India 2009 Call for Presentations
by Jarrod Millman Sept. 21, 2010

Sept. 21, 2010

========================== SciPy 2010 Call for Papers ========================== The second `SciPy India Conference <http://scipy.in>`_ will be held from December 13th to 18th, 2010 at `IIIT-Hyderabad <http://www.iiit.ac.in/>`_. At this conference, novel applications and breakthroughs made in the pursuit of science using Python are presented. Attended by leading figures from both academia and industry, it is an excellent opportunity to experience the cutting edge of scientific software development. The conference is followed by two days of tutorials and a code sprint, during which community experts provide training on several scientific Python packages. We invite you to take part by submitting a talk abstract on the conference website at: http://scipy.in Talk/Paper Submission ========================== We solicit talks and accompanying papers (either formal academic or magazine-style articles) that discuss topics regarding scientific computing using Python, including applications, teaching, development and research. Papers are included in the peer-reviewed conference proceedings, published online. Please note that submissions primarily aimed at the promotion of a commercial product or service will not be considered. Important Dates ========================== Monday, Oct. 11: Abstracts Due Saturday, Oct. 30: Schedule announced Tuesday, Nov. 30: Proceedings paper submission due Monday-Tuesday, Dec. 13-14: Conference Wednesday-Friday, Dec. 15-17: Tutorials/Sprints Saturday, Dec. 18: Sprints Organizers ========================== * Jarrod Millman, Neuroscience Institute, UC Berkeley, USA (Conference Co-Chair) * Prabhu Ramachandran, Department of Aerospace Engineering, IIT Bombay, India (Conference Co-Chair) * FOSSEE Team

1 0

utility for composing sparse matrices
by Robert Cimrman Sept. 21, 2010

Sept. 21, 2010

FYI: I have created a function that creates a sparse matrix from sparse (or dense) blocks, somewhat in the spirit of np.r_ and friends, see the example below. Shall we add it to scipy.sparse? It may help to go around some of the fancy-indexing limitations. r. In [1]: from compose_sparse import compose_sparse In [2]: %doctest_mode *** Pasting of code with ">>>" or "..." has been enabled. Exception reporting mode: Plain Doctest mode is: ON >>> >>> import scipy.sparse as sps >>> >>> A = sps.csr_matrix([[1, 0], [0, 1]]) >>> >>> B = sps.coo_matrix([[1, 1]]) >>> >>> K = compose_sparse([[A, B.T], [B, 0]]) >>> >>> print K.todense() --> print(K.todense()) [[1 0 1] [0 1 1] [1 1 0]]

3 4

PATCH: enhancing fmin_bfgs to handle functions that are not defined everywhere
by Enrico Avventi Sept. 20, 2010

Sept. 20, 2010

hello, i would like to propose a patch for enhancing the optimization module. the svn diff output is attached. let me know what do you think about it. *Rationale:* Convergence theory for a large class of uncostrained optimization algorithms (descent methods) does not require the function to be defined everywhere. In fact you just need a function that is defined in an open convex set and whose sublevel set: {x: f(x) <= f(x0) } is compact where x0 is the starting point. It would be helpful to let the optimization routines to handle the most general case for which they converge. *Assumptions:* I assumed that the function tends to (plus) infinity as it approaches its domain boundary. This is reasonable in the sense that if not true the function either tends to minus infinity or can be extended to a larger set. The function is assumed to be infinite ouside its domain with derivative, if any, non defined (None). *Implementation:* The line search methods is the only block of the optimization routine(s) that needs to be changed significantly in order to handle such functions. In particular i modified the line search method (line_search_wolfe2) to call a new zoom algorith (_barrier_zoom) whenever one of the bracketing points leads to an unfeasible point. _barrier_zoom uses linear-logarithmic and quadratic-logarithmic interpolants instead of quadratic and cubic ones in order to generate new candidates steps. Note that the interpolation strategy may not be ideal and could be certainly improved. All descent algorithms can be easily adapted to use such a line search, here i started with fmin_bfgs. In order to adapt it i added a new exception (DomainError) and modified phi and derphi in line_search_wolfe1 to raise it when they are evaluated outside the domain. Whenever that happens fmin_bfgs catches the exception and calls line_search_wolfe2. Additionally line_search_wolfe2 will call the original _zoom method if the bracketing phase leads to finite-valued points. this ensures that for problems that are defined everywhere the old code-path is always taken. * Testing:* I added two functions for testing the scalar_search_wolfe2 phi1(x) = -100 x - log(1-x) and phi2(x) = -100 x + (1-x)^-k. All old and new tests pass. regards, Enrico

2 1

Error in scipy.ndimage documentation
by Matthieu Brucher Sept. 20, 2010

Sept. 20, 2010

Hi, The following function is referenced on the scipy.ndimage page (http://docs.scipy.org/doc/scipy/reference/tutorial/ndimage.html#ndimage-cca…) : static int _shift_function(int *output_coordinates, double* input_coordinates, int output_rank, int input_rank, void *callback_data) { int ii; /* get the shift from the callback data pointer: */ double shift = *(double*)callback_data; /* calculate the coordinates: */ for(ii = 0; ii < irank; ii++) icoor[ii] = ocoor[ii] - shift; /* return OK status: */ return 1; } It obviously cannot be compiled (irank -> input_rank, icoor -> input_coordinates, ocoor -> output_coordinates) Matthieu -- Information System Engineer, Ph.D. Blog: http://matt.eifelle.com LinkedIn: http://www.linkedin.com/in/matthieubrucher

2 2

dopri ODE integrator solout integration constraints
by Janick Martinez Esturo Sept. 17, 2010

Sept. 17, 2010

Hi, I need to be able to observe and possibly stop an ODE integration using the scipy.integrate.ode dopri5 and dopri853 integrators. However, the supporting call of the _solout method is completely disabled as the fortran iout parameter is set to zero inside the wrapper code. I therefore attached a patch that does not break the existing interface and enables the usage of a costum _solout method by overwriting the method in a derived class as follows: from scipy.integrate.ode import ode, dopri5, IntegratorBase class mydopri5(dopri5): name = 'mydopri5' def _solout(self, nr,xold,x,y,nd,icomp,con): print("Costum solout with variable access ", x, y) return(1) #Note: return -1 to break def reset(self,n,has_jac): super(mydopri5, self).reset(n, has_jac) self._iout = 1 self.call_args = [self.rtol,self.atol,self._solout,self._iout,self.work,self.iwork] if mydopri5.runner: IntegratorBase.integrator_classes.append(mydopri5) y0, t0 = [1.0, 2.0], 0 def f(t, y): return [y[0] + y[1], -y[1]**2] r = ode(f).set_integrator('mydopri5') r.set_initial_value(y0, t0) r.integrate(10) print r.t, r.y Here I derived from dopri5 overwrote the reset method to be called with iout=1 and implemented a simple _solout callback. I'm using this code in a productive environment and just wanted to know if there is a chance that the patch might get included in a future release of scipy or what would be necessary to to so? I'm not experienced with Fortran libs, but it might be interesting to make the Fortran CONTD5 function also avaible to be called from the _solout method in order to interpolate intermediate values correctly, but that would require a little bit more efford i guess. Any comments about this would be nice. Greets Janick

2 1

Scipy on Python 3
by Pauli Virtanen Sept. 16, 2010

Sept. 16, 2010

Hi, I flushed the Python 3 branch containing work from me and David to SVN trunk. Scipy now builds with Python 3, and all tests pass, except for scipy.weave which still needs to be ported. More testing is welcome. I suspect Scipy's test suite does not cover all of the code, so there might be some work left to do. You'll probably need the latest 1.5.x branch Numpy to build, due to some fixes in Numpy's distutils that are not in 1.5.0. Pauli

8 16

Re: [SciPy-Dev] Scipy.org down, again
by Jonathan Stickel Sept. 15, 2010

Sept. 15, 2010

I notice this happens a lot, too. Interestingly, and fortunately for me, the docs seems to be much less affected: http://docs.scipy.org/ So I at least can reach the documentation. Jonathan On 9/15/10 11:00 , scipy-dev-request(a)scipy.org wrote: > Date: Tue, 14 Sep 2010 16:40:31 -0700 > From: Fernando Perez<fperez.net(a)gmail.com> > Subject: [SciPy-Dev] Scipy.org down, again > To: SciPy Developers List<scipy-dev(a)scipy.net>

1 0

Scipy.org down, again
by Fernando Perez Sept. 14, 2010

Sept. 14, 2010

f

1 0

spatial.test() failure
by Nils Wagner Sept. 13, 2010

Sept. 13, 2010

>>> from scipy import spatial >>> spatial.test() Running unit tests for scipy.spatial NumPy version 2.0.0.dev8714 NumPy is installed in /home/nwagner/local/lib64/python2.6/site-packages/numpy SciPy version 0.9.0.dev6803 SciPy is installed in /home/nwagner/local/lib64/python2.6/site-packages/scipy Python version 2.6.5 (r265:79063, Jul 5 2010, 11:46:13) [GCC 4.5.0 20100604 [gcc-4_5-branch revision 160292]] nose version 0.11.2 ............................................................................................................................F.................................................................................................................................................. ====================================================================== FAIL: Tests pdist(X, 'minkowski') on iris data. (float32) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/spatial/tests/test_distance.py", line 837, in test_pdist_minkowski_3_2_iris_float32 self.failUnless(within_tol(Y_test1, Y_right, eps)) AssertionError ---------------------------------------------------------------------- Ran 271 tests in 38.005s FAILED (failures=1) <nose.result.TextTestResult run=271 errors=0 failures=1>

2 1

problem building scipy with qhull
by Martin Hofsäß Sept. 9, 2010

Sept. 9, 2010

Hi all together, I had the newest scipy version 0.9 from the svn trunk and I get one error by building scipy under windows with the intel mkl and the MSVC & IFORT compiler. Who can me help? Executing scons command (pkg is scipy.spatial): C:\Python26\python.exe "C:\Python26\lib\site-packages\numscons\scons-local\scons.py" -f scipy\spatial\SConstruct -I. scons_tool_path="" src_dir="scipy\spatial" pkg_path="scipy\spatial" pkg_name="scipy.spatial" log_level=0 distutils_libdir="..\..\..\..\build\lib.win32-2.6" distutils_clibdir="..\..\..\..\build\temp.win32-2.6" distutils_install_prefix="C:\Python26\Lib\site-packages\scipy\spatial" cc_opt=msvc cc_opt=msvc debug=0 f77_opt=ifort cxx_opt=msvc include_bootstrap=C:\Python26\lib\site-packages\numpy\core\include bypass=1 import_env=0 silent=0 bootstrapping=0 scons: Reading SConscript files ... Mkdir("build\scons\scipy\spatial") C:\Programme\Intel\Compiler\11.0\074\fortran\bin\ifortvars.bat scons: done reading SConscript files. scons: Building targets ... cl /Fobuild\scons\scipy\spatial\ckdtree.obj /c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo /IC:\Python26\lib\site-packages\numpy\core\include /IC:\Python26\include build\scons\scipy\spatial\ckdtree.c ckdtree.c link /nologo /dll /OUT:build\scons\scipy\spatial\ckdtree.pyd /LIBPATH:C:\Python26\libs build\scons\scipy\spatial\ckdtree.obj Bibliothek "build\scons\scipy\spatial\ckdtree.lib" und Objekt "build\scons\scipy\spatial\ckdtree.exp" werden erstellt. Install file: "build\scons\scipy\spatial\ckdtree.pyd" as "build\lib.win32-2.6\scipy\spatial\ckdtree.pyd" cl /Fobuild\scons\scipy\spatial\qhull.obj /c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo /IC:\Python26\lib\site-packages\numpy\core\include /IC:\Python26\include build\scons\scipy\spatial\qhull.c qhull.c build\scons\scipy\spatial\qhull.c(1143) : warning C4133: '=': Inkompatible Typen - von 'PyObject *' zu 'PyArrayObject *' build\scons\scipy\spatial\qhull.c(1768) : warning C4018: '>=': Konflikt zwischen 'signed' und 'unsigned' build\scons\scipy\spatial\qhull.c(2057) : warning C4018: '>=': Konflikt zwischen 'signed' und 'unsigned' build\scons\scipy\spatial\qhull.c(2777) : warning C4013: 'dgesv_' undefiniert; Annahme: extern mit Rckgabetyp int cl /Fobuild\scons\scipy\spatial\qhull\src\geom.obj /c build\scons\scipy\spatial\qhull\src\geom.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo geom.c cl /Fobuild\scons\scipy\spatial\qhull\src\geom2.obj /c build\scons\scipy\spatial\qhull\src\geom2.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo geom2.c cl /Fobuild\scons\scipy\spatial\qhull\src\global.obj /c build\scons\scipy\spatial\qhull\src\global.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo global.c cl /Fobuild\scons\scipy\spatial\qhull\src\io.obj /c build\scons\scipy\spatial\qhull\src\io.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo io.c cl /Fobuild\scons\scipy\spatial\qhull\src\mem.obj /c build\scons\scipy\spatial\qhull\src\mem.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo mem.c cl /Fobuild\scons\scipy\spatial\qhull\src\merge.obj /c build\scons\scipy\spatial\qhull\src\merge.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo merge.c cl /Fobuild\scons\scipy\spatial\qhull\src\poly.obj /c build\scons\scipy\spatial\qhull\src\poly.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo poly.c cl /Fobuild\scons\scipy\spatial\qhull\src\poly2.obj /c build\scons\scipy\spatial\qhull\src\poly2.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo poly2.c cl /Fobuild\scons\scipy\spatial\qhull\src\qset.obj /c build\scons\scipy\spatial\qhull\src\qset.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo qset.c cl /Fobuild\scons\scipy\spatial\qhull\src\user.obj /c build\scons\scipy\spatial\qhull\src\user.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo user.c cl /Fobuild\scons\scipy\spatial\qhull\src\stat.obj /c build\scons\scipy\spatial\qhull\src\stat.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo stat.c cl /Fobuild\scons\scipy\spatial\qhull\src\qhull.obj /c build\scons\scipy\spatial\qhull\src\qhull.c /MD /nologo /GS- /D_CRT_SECURE_NO_WARNINGS /openmp /Ox /DNDEBUG /W3 /nologo qhull.c lib /nologo /OUT:build\scons\scipy\spatial\qhullsrc.lib build\scons\scipy\spatial\qhull\src\geom.obj build\scons\scipy\spatial\qhull\src\geom2.obj build\scons\scipy\spatial\qhull\src\global.obj build\scons\scipy\spatial\qhull\src\io.obj build\scons\scipy\spatial\qhull\src\mem.obj build\scons\scipy\spatial\qhull\src\merge.obj build\scons\scipy\spatial\qhull\src\poly.obj build\scons\scipy\spatial\qhull\src\poly2.obj build\scons\scipy\spatial\qhull\src\qset.obj build\scons\scipy\spatial\qhull\src\user.obj build\scons\scipy\spatial\qhull\src\stat.obj build\scons\scipy\spatial\qhull\src\qhull.obj Install file: "build\scons\scipy\spatial\qhullsrc.lib" as "build\temp.win32-2.6\qhullsrc.lib" link /nologo /dll /OUT:build\scons\scipy\spatial\qhull.pyd /LIBPATH:C:\Python26\libs build\scons\scipy\spatial\qhullsrc.lib build\temp.win32-2.6\qhullsrc.lib build\scons\scipy\spatial\qhull.obj Bibliothek "build\scons\scipy\spatial\qhull.lib" und Objekt "build\scons\scipy\spatial\qhull.exp" werden erstellt. qhull.obj : error LNK2019: Verweis auf nicht aufgelöstes externes Symbol "_dgesv_" in Funktion "___pyx_pf_5scipy_7spatial_5qhull__get_barycentric_transforms". build\scons\scipy\spatial\qhull.pyd : fatal error LNK1120: 1 nicht aufgelöste externe Verweise. scons: building terminated because of errors.

2 1

Unexpected behaviour from timeseries.filled (in Python 2.7)
by Åsmund Hjulstad Sept. 6, 2010

Sept. 6, 2010

It seems something in Python 2.7 / Numpy 1.5 (or my Python 2.7 installation) is breaking scikits.timeseries.filled import scikits.timeseries from numpy import array t3 = scikits.timeseries.time_series(array([0.]), mask=[True],dates=['2010-Aug'],freq='M') t3.fill_value= 0. print(t3) print(t3.data) print(t3.fill_value) print(t3.filled()) # It seems it ignores the fill_value parameter print(t3.filled(0.)) In Python 2.6 this outputs c:\python26\python scikits_bug.py [--] [ 0.] 0.0 [ 0.] [ 0.] , but in Python 2.7 I get: c:\python27\python scikits_bug.py [--] [ 0.] 0.0 [ 1.00000000e+20] # This must be wrong! [ 0.] Python 2.6: scikits.timeseries version 0.91.3, numpy 1.3.0 Python 2.7: scikits.timeseries version 0.91.3, numpy 1.5.0 -- Åsmund Hjulstad, asmund(a)hjulstad.com

2 4

N-dimensional interpolation
by Pauli Virtanen Sept. 6, 2010

Sept. 6, 2010

Hi all, I took the Qhull by the horns, and wrote a straightforward `griddata` implementation for working in N-D: http://github.com/pv/scipy-work/tree/qhull http://github.com/pv/scipy-work/blob/qhull/scipy/interpolate/qhull.pyx Sample (4-D): ------------------------------------------------------------- import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import griddata def f(x, y, z, u): return np.cos(x + np.cos(y + np.cos(z + np.cos(u)))) points = np.random.rand(2000, 4) values = f(*points.T) p = np.linspace(0, 1, 1500) z = griddata(points, values, np.c_[p, p, p, p], method='linear') plt.plot(p, z, '-', p, f(p, p, p, p), '-') plt.show() ------------------------------------------------------------- It performs the N-D Delaunay tesselation with Qhull, and after that does linear barycentric interpolation on the simplex containing each point. (The nearest-neighbor "interpolation" is, on the other hand, implemented using scipy.spatial.KDTree.) I ended up writing some custom simplex-walking code for locating the simplex containing the points -- this is much faster than a brute-force search. However, it is slightly different from what `qh_findbestfacet` does. [If you're a computational geometry expert, see below...] Speed-wise, Qhull appears to lose in 2-D to `scikits.delaunay` by a factor of about 10 in triangulation in my tests; this, however, includes the time taken by computing the barycentric coordinate transforms. Interpolation, however, seems to be faster. *** I think this would be a nice addition to `scipy.optimize`. I'd like to get it in for 0.9. Later on, we can specialize the `griddata` function to perform better in low dimensions (ndim=2, 3) and add more interpolation methods there; for example natural neighbors, interpolating splines, etc. In ndim > 3, the `griddata` is however already feature-equivalent to Octave and MATLAB. Comments? Pauli PS. A question for computational geometry experts: `qh_findbestfacet` finds the facet on the lower convex hull whose oriented hyperplane distance to the point lifted on the paraboloid is maximal --- however, this does not always give the simplex containing the point in the projected Delaunay tesselation. There's a counterexample in qhull.pyx:_find_simplex docstring. On the other hand, Qhull documentation claims that this is the way to locate the Delaunay simplex containing a point. What gives? It's clear, however, that if a simplex contains the point, then the hyperplane distance is positive. So currently, I just walk around positive-distance simplices checking the inside condition for each and hoping for the best. If nothing is found, I just fall back to brute force search. This seems to work well in practice. However, if a point is outside the triangulation (but inside the bounding box), I have to fall back to a brute-force search. Is there a better option here? *** Another sample (2-D, object-oriented interface): ------------------------------------------------------------- import time import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import LinearNDInterpolator # some random data x = np.array([(0,0), (0,1), (1, 1), (1, 0)], dtype=np.double) np.random.seed(0) xr = np.random.rand(200*200, 2).astype(np.double) x = np.r_[xr, x] y = np.arange(x.shape[0], dtype=np.double) # evaluate on a grid nn = 500 xx = np.linspace(-0.1, 1.1, nn) yy = np.linspace(-0.1, 1.1, nn) xx, yy = np.broadcast_arrays(xx[:,None], yy[None,:]) xi = np.array([xx, yy]).transpose(1,2,0).copy() start = time.time() ip = LinearNDInterpolator(x, y) print "Triangulation", time.time() - start start = time.time() zi = ip(xi) print "Interpolation", time.time() - start from scikits.delaunay import Triangulation start = time.time() tri2 = Triangulation(x[:,0], x[:,1]) print "scikits.delaunay triangulation", time.time() - start start = time.time() ip2 = tri2.linear_interpolator(y) zi2 = ip2[-0.1:1.1:500j,-0.1:1.1:500j].T print "scikits.delaunay interpolation", time.time() - start print "rel-difference", np.nanmax(abs(zi - zi2))/np.nanmax(np.abs(zi)) plt.figure() plt.imshow(zi) plt.clim(y.min(), y.max()) plt.colorbar() plt.figure() plt.imshow(zi2) plt.clim(y.min(), y.max()) plt.colorbar() plt.show() ------------------------------------------------------------- -- Pauli Virtanen

5 8

Nonlinear large-scale root finding
by Pauli Virtanen Sept. 5, 2010

Sept. 5, 2010

Hi, I finally committed the scipy.optimize.nonlin rewrite that was pending for a long time. What is now there for large-scale root finding is: * Newton-Krylov method (``scipy.optimize.newton_krylov``) * Generalized secant methods: - Limited-memory Broyden (``scipy.optimize.broyden1``, ``scipy.optimize.broyden2``) - Anderson (``scipy.optimize.anderson``) * Simple iterations (``scipy.optimize.diagbroyden``, ``scipy.optimize.excitingmixing``, ``scipy.optimize.linearmixing``) Here's a simple demo: http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#large-prob… -- Pauli Virtanen

1 0

Re: [SciPy-Dev] RFR: N-dimensional interpolation
by Eric Bruning Sept. 5, 2010

Sept. 5, 2010

I made a trial run with some 3D data from a weather radar this evening. I'm interpolating to a horizontal slice through a (distorted) spherical grid, as shown here: http://bitbucket.org/deeplycloudy/pyrsl/src/tip/example.py Both linear and nearest-neighbor work great. Image of linear interpolation at http://yfrog.com/msspop Many thanks for cooking this up. I've been hoping for something like this for a long time. -Eric

1 0

Scipy.org down (again)
by Fernando Perez Sept. 1, 2010

Sept. 1, 2010

Howdy, the usual... Cheers, f

1 0