Mailman 3 May 2012 - SciPy-Dev

SLSQP Constrained Optimizer Status
by Rob Falck 25 Nov '13

25 Nov '13

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

scipy.signal.normalize problems?
by Josh Lawrence 29 Jun '12

29 Jun '12

Hey all, I've been having some problems designing a Chebyshev filter and I think I have narrowed down the hang-up to scipy.signal.normalize. I think what's going on in my case is that line 286 of filter_design.py (the first allclose call in the normalize function) is producing a false positive. Here's the function definition: def normalize(b, a): """Normalize polynomial representation of a transfer function. If values of b are too close to 0, they are removed. In that case, a BadCoefficients warning is emitted. """ b, a = map(atleast_1d, (b, a)) if len(a.shape) != 1: raise ValueError("Denominator polynomial must be rank-1 array.") if len(b.shape) > 2: raise ValueError("Numerator polynomial must be rank-1 or" " rank-2 array.") if len(b.shape) == 1: b = asarray([b], b.dtype.char) while a[0] == 0.0 and len(a) > 1: a = a[1:] outb = b * (1.0) / a[0] outa = a * (1.0) / a[0] if allclose(outb[:, 0], 0, rtol=1e-14): <------------------ Line 286 warnings.warn("Badly conditioned filter coefficients (numerator): the " "results may be meaningless", BadCoefficients) while allclose(outb[:, 0], 0, rtol=1e-14) and (outb.shape[-1] > 1): outb = outb[:, 1:] if outb.shape[0] == 1: outb = outb[0] return outb, outa I marked line 286. If I reproduce all the steps carried out by scipy.signal.iirdesign, I end up with a (b, a) pair which results of scipy.signal.lp2lp and looks like this: In [106]: b_lp2 Out[106]: array([ 1.55431359e-06+0.j]) In [107]: a_lp2 Out[107]: array([ 1.00000000e+00 +0.00000000e+00j, 3.46306104e-01 -2.01282794e-16j, 2.42572185e-01 -6.08207573e-17j, 5.92946943e-02 +0.00000000e+00j, 1.82069156e-02 +5.55318531e-18j, 2.89328123e-03 +0.00000000e+00j, 4.36566281e-04 -2.95766719e-19j, 3.50842810e-05 -3.19180568e-20j, 1.64641246e-06 -1.00966301e-21j]) scipy.signal.iirdesign takes b_lp2, a_lp2 (my local variable names to keep track of what's going on) and runs them through scipy.signal.bilinear (in filter_design.py bilinear is called on line 624 within iirfilter. iirdesign calls iirfilter which calls bilinear). Inside bilinear, normalize is called on line 445. I've made my own class with bilinear copied and pasted from filter_design.py to test things. In bilinear, the input to normalize is given by b = [ 1.55431359e-06 1.24345087e-05 4.35207804e-05 8.70415608e-05 1.08801951e-04 8.70415608e-05 4.35207804e-05 1.24345087e-05 1.55431359e-06] a = [ 72269.02590913 -562426.61430468 1918276.19173089 -3745112.83646825 4577612.13937628 -3586970.61385926 1759651.18184723 -494097.93515708 60799.46134722] In normalize, right before the allclose() call, outb is defined by outb = [[ 2.15073272e-11 1.72058618e-10 6.02205162e-10 1.20441032e-09 1.50551290e-09 1.20441032e-09 6.02205162e-10 1.72058618e-10 2.15073272e-11]] From what I read of the function normalize, it should only evaluate true if all of the coefficients in outb are smaller than 1e-14. However, that's not what is going on. I have access to MATLAB and if I go through the same design criteria to design a chebyshev filter, I get the following: b = 1.0e-08 * Columns 1 through 5 0.002150733144728 0.017205865157826 0.060220528052392 0.120441056104784 0.150551320130980 Columns 6 through 9 0.120441056104784 0.060220528052392 0.017205865157826 0.002150733144728 which matches up rather well for several significant figures. I apologize if this is not clearly explained, but I'm not sure what to do. I tried messing around with the arguments to allclose (switching it to be allclose(0, outb[:,0], ...), or changing the keyword from rtol to atol). I am also not sure why normalize is setup to run on rank-2 arrays. I looked through filter_design.py and none of the functions contained in it send a rank-2 array to normalize from what I can tell. Any thoughts? Cheers, Josh Lawrence

2 8

scipy 0.11.0 release schedule
by Ralf Gommers 05 Jun '12

05 Jun '12

Hi all, It's time for a new release I think; 0.10.0 was released just over 6 months ago. There are a lot of PRs to merge, but no release blockers anymore as far as I know. Does anyone still have important fixes or other things that should go in? For the release schedule I propose: June 7: beta 1 June 17: rc 1 June 30: rc 2 July 7: final release Ralf

4 8

test failures with python 2.4
by Denis Laxalde 05 Jun '12

05 Jun '12

Hi, Here is a list of test failures with python 2.4 (there are also failures with 2.7, see below). Interpolate: test_interpnd.TestCloughTocher2DInterpolator.test_dense ... Exception exceptions.TypeError: <exceptions.TypeError instance at 0x7fd4d67fecf8> in 'scipy.spatial.qhull._get_delaunay_info' ignored Segmentation fault Spatial: test_qhull.TestUtilities.test_convex_hull ... ERROR test_qhull.TestUtilities.test_degenerate_barycentric_transforms ... ERROR test_qhull.TestUtilities.test_find_simplex ... Exception exceptions.TypeError: <exceptions.TypeError instance at 0x7f9e857963b0> in 'scipy.spatial.qhull._get_delaunay_info' ignored FAIL test_qhull.TestUtilities.test_more_barycentric_transforms ... ERROR test_qhull.TestUtilities.test_plane_distance ... ok ====================================================================== ERROR: test_qhull.TestUtilities.test_convex_hull ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/nose/case.py", line 197, in runTest self.test(*self.arg) File "/home/denis/.local/lib/python2.4/site-packages/scipy/spatial/tests/test_qhull.py", line 75, in test_convex_hull assert_equal(tri.convex_hull, [[1, 2], [3, 2], [1, 0], [3, 0]]) File "qhull.pyx", line 1109, in scipy.spatial.qhull.Delaunay.convex_hull (scipy/spatial/qhull.c:8406) ValueError: cannot resize an array references or is referenced by another array in this way. Use the resize function ====================================================================== ERROR: test_qhull.TestUtilities.test_degenerate_barycentric_transforms ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/nose/case.py", line 197, in runTest self.test(*self.arg) File "/home/denis/.local/lib/python2.4/site-packages/scipy/spatial/tests/test_qhull.py", line 133, in test_degenerate_barycentric_transforms bad_count = np.isnan(tri.transform[:,0,0]).sum() File "qhull.pyx", line 1028, in scipy.spatial.qhull.Delaunay.transform (scipy/spatial/qhull.c:7449) File "qhull.pyx", line 398, in scipy.spatial.qhull._get_barycentric_transforms (scipy/spatial/qhull.c:4300) File "/home/denis/.local/lib/python2.4/site-packages/numpy/core/numeric.py", line 235, in asarray return array(a, dtype, copy=False, order=order) File "stringsource", line 366, in View.MemoryView.memoryview.__getitem__ (scipy/spatial/qhull.c:14704) File "stringsource", line 650, in View.MemoryView._unellipsify (scipy/spatial/qhull.c:17965) TypeError: Cannot index with type '<type 'int'>' ====================================================================== ERROR: test_qhull.TestUtilities.test_more_barycentric_transforms ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/nose/case.py", line 197, in runTest self.test(*self.arg) File "/home/denis/.local/lib/python2.4/site-packages/scipy/spatial/tests/test_qhull.py", line 157, in test_more_barycentric_transforms unit_cube=True) File "/home/denis/.local/lib/python2.4/site-packages/scipy/spatial/tests/test_qhull.py", line 96, in _check_barycentric_transforms c = barycentric_transform(tri.transform, centroids) File "qhull.pyx", line 1028, in scipy.spatial.qhull.Delaunay.transform (scipy/spatial/qhull.c:7449) File "qhull.pyx", line 398, in scipy.spatial.qhull._get_barycentric_transforms (scipy/spatial/qhull.c:4300) File "/home/denis/.local/lib/python2.4/site-packages/numpy/core/numeric.py", line 235, in asarray return array(a, dtype, copy=False, order=order) File "stringsource", line 366, in View.MemoryView.memoryview.__getitem__ (scipy/spatial/qhull.c:14704) File "stringsource", line 650, in View.MemoryView._unellipsify (scipy/spatial/qhull.c:17965) TypeError: Cannot index with type '<type 'int'>' ====================================================================== FAIL: test_qhull.TestUtilities.test_find_simplex ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/nose/case.py", line 197, in runTest self.test(*self.arg) File "/home/denis/.local/lib/python2.4/site-packages/scipy/spatial/tests/test_qhull.py", line 33, in test_find_simplex assert_equal(i, p[2], err_msg='%r' % (p,)) File "/home/denis/.local/lib/python2.4/site-packages/numpy/testing/utils.py", line 256, in assert_equal return assert_array_equal(actual, desired, err_msg, verbose) File "/home/denis/.local/lib/python2.4/site-packages/numpy/testing/utils.py", line 707, in assert_array_equal verbose=verbose, header='Arrays are not equal') File "/home/denis/.local/lib/python2.4/site-packages/numpy/testing/utils.py", line 636, in assert_array_compare raise AssertionError(msg) AssertionError: Arrays are not equal (0.25, 0.25, 1) (mismatch 100.0%) x: array(-1, dtype=int32) y: array(1) ---------------------------------------------------------------------- Ran 291 tests in 13.335s FAILED (errors=3, failures=1) Weave: FAIL: test_set_complex (test_scxx_object.TestObjectSetItemOpKey) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/scipy/weave/tests/test_scxx_object.py", line 943, in test_set_complex assert_equal(sys.getrefcount(key),4) # should be 3 File "/home/denis/.local/lib/python2.4/site-packages/numpy/testing/utils.py", line 313, in assert_equal raise AssertionError(msg) AssertionError: Items are not equal: ACTUAL: 3 DESIRED: 4 Failures not specific to python 2.4. Sparse: ====================================================================== ERROR: adding a dense matrix to a sparse matrix ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/tests/test_base.py", line 527, in test_add_dense sum1 = self.dat + self.datsp File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/dok.py", line 133, in __getitem__ raise TypeError('index must be a pair of integers or slices') TypeError: index must be a pair of integers or slices ====================================================================== ERROR: test_matmat_sparse (test_base.TestDOK) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/tests/test_base.py", line 425, in test_matmat_sparse assert_array_almost_equal( a2*bsp, a*b) File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/dok.py", line 133, in __getitem__ raise TypeError('index must be a pair of integers or slices') TypeError: index must be a pair of integers or slices ====================================================================== ERROR: test_radd (test_base.TestDOK) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/tests/test_base.py", line 279, in test_radd c = a + b File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/dok.py", line 133, in __getitem__ raise TypeError('index must be a pair of integers or slices') TypeError: index must be a pair of integers or slices ====================================================================== ERROR: test_rsub (test_base.TestDOK) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/tests/test_base.py", line 290, in test_rsub assert_array_equal((self.dat - self.datsp),[[0,0,0,0],[0,0,0,0],[0,0,0,0]]) File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/dok.py", line 133, in __getitem__ raise TypeError('index must be a pair of integers or slices') TypeError: index must be a pair of integers or slices ====================================================================== ERROR: subtracting a dense matrix to/from a sparse matrix ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/tests/test_base.py", line 535, in test_sub_dense sum1 = 3*self.dat - self.datsp File "/home/denis/.local/lib/python2.4/site-packages/scipy/sparse/dok.py", line 133, in __getitem__ raise TypeError('index must be a pair of integers or slices') TypeError: index must be a pair of integers or slices ---------------------------------------------------------------------- Ran 1210 tests in 12.543s FAILED (KNOWNFAIL=4, errors=5) Special: ................../home/denis/.local/lib/python2.7/site-packages/scipy/special/tests/test_basic.py:1645: RuntimeWarning: divide by zero encountered in iv c1 = special.iv(v, x) /home/denis/.local/lib/python2.7/site-packages/scipy/special/tests/test_basic.py:1645: RuntimeWarning: overflow encountered in iv c1 = special.iv(v, x) /home/denis/.local/lib/python2.7/site-packages/scipy/special/tests/test_basic.py:1646: RuntimeWarning: invalid value encountered in iv c2 = special.iv(v, x+0j) /home/denis/.local/lib/python2.7/site-packages/scipy/special/tests/test_basic.py:1652: RuntimeWarning: divide by zero encountered in divide dc = abs(c1/c2 - 1) /home/denis/.local/lib/python2.7/site-packages/scipy/special/tests/test_basic.py:1652: RuntimeWarning: invalid value encountered in divide dc = abs(c1/c2 - 1) /home/denis/.local/lib/python2.7/site-packages/scipy/special/tests/test_basic.py:1659: RuntimeWarning: overflow encountered in iv assert_(dc[k] < 1e-9, (v[k], x[k], special.iv(v[k], x[k]), special.iv(v[k], x[k]+0j))) F...................................................K.K.............................................................................................................................................................................................................................................................................................................................................................................................K........K.................................................. ====================================================================== FAIL: test_iv_cephes_vs_amos_mass_test (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/denis/.local/lib/python2.7/site-packages/scipy/special/tests/test_basic.py", line 1659, in test_iv_cephes_vs_amos_mass_test assert_(dc[k] < 1e-9, (v[k], x[k], special.iv(v[k], x[k]), special.iv(v[k], x[k]+0j))) File "/usr/lib/pymodules/python2.7/numpy/testing/utils.py", line 34, in assert_ raise AssertionError(msg) AssertionError: (189.2947429454936, 3.0238805556481037, 4.089165443940765e-317, 0j) ---------------------------------------------------------------------- Ran 514 tests in 6.908s FAILED (KNOWNFAIL=4, failures=1) Building/testing environment for python 2.4 (Debian wheezy/sid x86_64): NumPy version 1.6.2 NumPy is installed in /home/denis/.local/lib/python2.4/site-packages/numpy SciPy version 0.11.0.dev-127acfd SciPy is installed in /home/denis/.local/lib/python2.4/site-packages/scipy Python version 2.4.6 (#2, Sep 25 2009, 22:22:06) [GCC 4.3.4] nose version 1.1.2 -- Denis

2 5

New contribution: bode() function for LTI systems
by Bjørn Forsman 04 Jun '12

04 Jun '12

Hi all, I'm new on this list. I have been wanting to have a bode() function for LTI systems in Python for some time now. Today I got tired of waiting and wrote one myself :-) I have sent pull request on github: https://github.com/scipy/scipy/pull/224 - ENH: ltisys: new bode() function While at it I also fixed a small bug in the lti class init function: https://github.com/scipy/scipy/pull/223 - ENH: ltisys: make lti zpk initialization work with plain lists Is this OK? Should I post the patches here too? Best regards, Bjørn Forsman

3 6

Adding t-test with unequal variances to stats.py
by Junkshops 03 Jun '12

03 Jun '12

Hi all, I've issued a pull request (http://github.com/scipy/scipy/pull/227) for a version of scipy/stats/stats.py with the following changes: 1) Adds a method for running a t-test with unequal or unknown population variances. ttest_ind assumes that population variances are equal. 2) Refactored common code in the 4 t-test methods into shared methods. 3) This section of code, which has variations in multiple methods, looks buggy to me: d = np.mean(a,axis) - np.mean(b,axis) svar = ((n1-1)*v1+(n2-1)*v2) / float(df) t = d/np.sqrt(svar*(1.0/n1 + 1.0/n2)) t = np.where((d==0)*(svar==0), 1.0, t) #define t=0/0 = 0, identical means Surely if d=0, regardless of svar, t should be set to 0, not 1. Similarly, if svar = 0 then both variances are zero (assuming that each data set has at least 2 points - perhaps there should be a check for this?). In that case, if d==0 t should be zero. Otherwise, t should be +/-inf. Hence, (svar==0) is redundant. Accordingly, I've changed the lines in all functions to be the equivalent of t = np.where((d==0), 0.0, t) This handles the case where both d and svar are 0. The respective tests have also been changed. If I'm missing something here, please let me know. Thanks, Gavin

2 2

How to get "correct" values for unit tests
by Andreas Hilboll 31 May '12

31 May '12

Hi, while working on the tests for the new scipy.interpolate.{Smooth,LSQ}SphereBiavariateSpline classes, I'm wondering how to come up with sensible TRUE example values to test against. In the case mentioned (see https://github.com/scipy/scipy/pull/192), I simply wrapped a routine (sphere.f) from FITPACK. So I could write a direct FORTRAN program using sphere.f to calculate some "TRUE" values. However, that would just check that the wrapping actually works. Is this considered enough? Ultimately, I would like a test to assure that the results are correct. But for that, wouldn't it be "better" (whatever that means) to use a different library to calculate the TRUE results? Sorry, this might be a confusing email. Cheers, Andreas.

3 2

sum_angle() and sum_polar() functions
by Robert Jördens 30 May '12

30 May '12

Hello, I submitted two functions to numpy that sum 2d matrices along angled cartesian or polar coordinates. https://github.com/numpy/numpy/pull/230 The two functions certainly have their main application in image processing and might be better suited for scipy of scikits-image. sum_angle() is not much more than the old scipy.misc.pilutils.radon() transform. But the later is deprecated and has several problems (floats, non-conserved sum(), interpolation, speed) as discussed in the pull request. The new scikits-image.transform.radon() appears to be more generic but a bit complicated and potentially even slower than the imrotate()-based version in scipy. Could sum_angle() and sum_polar() find a place in scipy or scikits-image or are they simple enough to be useful for numpy? -- Robert Jordens.

4 6

ANN: SfePy 2012.2
by Robert Cimrman 29 May '12

29 May '12

I am pleased to announce release 2012.2 of SfePy. Description ----------- SfePy (simple finite elements in Python) is a software for solving systems of coupled partial differential equations by the finite element method. The code is based on NumPy and SciPy packages. It is distributed under the new BSD license. Home page: http://sfepy.org Downloads, mailing list, wiki: http://code.google.com/p/sfepy/ Git (source) repository, issue tracker: http://github.com/sfepy Highlights of this release -------------------------- - reimplemented acoustic band gaps code using the homogenization engine - high order quadrature rules - unified dot product and mass terms, lots of other term updates/fixes, - updated the PDE solver application For full release notes see http://docs.sfepy.org/doc/release_notes.html#id1 (rather long and technical). Best regards, Robert Cimrman and Contributors (*) (*) Contributors to this release (alphabetical order): Vladimír Lukeš, Andre Smit

1 0

Scipy Docs permissions
by Karmel Allison 29 May '12

29 May '12

Hi all, After many years of my living off the work you do, I thought I'd start helping out. Can I get edit permissions for the docs to start? I just registered with the username karmel. Thanks! Karmel

4 6