Mailman 3 June 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

chi-square test for a contingency (R x C) table
by Warren Weckesser July 12, 2010

July 12, 2010

I've been digging into some basic statistics recently, and developed the following function for applying the chi-square test to a contingency table. Does something like this already exist in scipy.stats? If not, any objects to adding it? (Tests are already written :) Warren ----- def chisquare_contingency(table): """Chi-square calculation for a contingency (R x C) table. This function computes the chi-square statistic and p-value of the data in the table. The expected frequencies are computed based on the relative frequencies in the table. Parameters ---------- table : array_like, 2D The contingency table, also known as the R x C table. Returns ------- chisquare statistic : float The chisquare test statistic p : float The p-value of the test. """ table = np.asarray(table) if table.ndim != 2: raise ValueError("table must be a 2D array.") # Create the table of expected frequencies. total = table.sum() row_sum = table.sum(axis=1).reshape(-1,1) col_sum = table.sum(axis=0) expected = row_sum * col_sum / float(total) # Since we are passing in 1D arrays of length table.size, the default # number of degrees of freedom is table.size-1. # For a contingency table, the actual number degrees of freedom is # (nr - 1)*(nc-1). We use the ddof argument # of the chisquare function to adjust the default. nr, nc = table.shape dof = (nr - 1) * (nc - 1) dof_adjust = (table.size - 1) - dof chi2, p = chisquare(np.ravel(table), np.ravel(expected), ddof=dof_adjust) return chi2, p -----

4 39

Is imresize() defined twice in scipy.misc.pilutil.py
by Vincent Davis July 3, 2010

July 3, 2010

Is imresize defined twice in scipy.misc.pilutil.py or am i missing something? Thanks Vincent Below is a separate issue maybe related to the above. According to the docs size : int, float or tuple * int - Percentage of current size. * float - Fraction of current size. * tuple - Size of the output image. an integer should work. Just wonder if someone can confirm this >>> from scipy.misc.pilutil import imresize >>> x = np.random.randn(3,3) >>> imresize(x, 2) Traceback (most recent call last): File "/Applications/WingIDE.app/Contents/MacOS/src/debug/tserver/_sandbox.py", line 1, in <module> # Used internally for debug sandbox under external interpreter File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/misc/pilutil.py", line 316, in imresize imnew = im.resize(size) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/PIL/Image.py", line 1302, in resize im = self.im.resize(size, resample) TypeError: must be 2-item sequence, not float >>> imresize(x, 2.0) array([[255, 255, 252, 252, 171, 171], [255, 255, 252, 252, 171, 171], [163, 163, 0, 0, 185, 185], [163, 163, 0, 0, 185, 185], [255, 255, 223, 223, 167, 167], [255, 255, 223, 223, 167, 167]], dtype=uint8) >>> imresize(x, (5,4)) array([[255, 255, 252, 171], [255, 255, 252, 171], [163, 163, 0, 185], [163, 163, 0, 185], [255, 255, 223, 167]], dtype=uint8)

3 4

function that are duplicates of numpy, such as scipy.special.basic.log
by Vincent Davis July 3, 2010

July 3, 2010

Is there a reason to keep function such as scipy.special.basic.log() that appears to duplicate numpy.log() Vincent

4 5

Warnings raised (from fit in scipy.stats)
by Skipper Seabold July 2, 2010

July 2, 2010

Since the raising of warning behavior has been changed (I believe), I have been running into a lot of warnings in my code when say I do something like In [120]: from scipy import stats In [121]: y = [-45, -3, 1, 0, 1, 3] In [122]: v = stats.norm.pdf(y)/stats.norm.cdf(y) Warning: invalid value encountered in divide Sometimes, this is useful to know. Sometimes, though, it's very disturbing when it's encountered in some kind of iteration or optimization. I have been using numpy.clip to get around this in my own code, but when it's buried a bit deeper, it's not quite so simple. Take this example. In [123]: import numpy as np In [124]: np.random.seed(12345) In [125]: B = 6.0 In [126]: x = np.random.exponential(scale=B, size=5000) In [127]: from scipy.stats import expon In [128]: expon.fit(x) <dozens of warnings clipped> Out[128]: (0.21874043533906118, 5.7122829778172939) The fit is achieved by fmin (as far as I know, since disp=0 in the rv_continuous.fit...), but there are a number of warnings emitted. Is there any middle ground to be had in these type of situations via context management perhaps? Should I file a ticket? Skipper

3 6

Contribute an evolutionary optimization algorithm
by bowie_22 June 30, 2010

June 30, 2010

Hello together, at the moment I am evaluationg the optimization features of scipy. During that evaluation I missed a evolutionary optimization algorithm in the"Global" section (there is anneal and brute). Just for fun I have started to implement an evo algorithm that uses (as far as possible) the same signature as the other algorithms in scipy.optimize. My plan is to use it as if it would be a part of scipy. In the meantime I ask my self, if it would make sense to contribute the evo stuff to scipy. But I have no idea how to do this. Who is the person who decides if this is possible or not? Are there means to ensure a minimum of qualitiy standards (checklists, reviews ...)? Can someone please give me some hints. Regs Marcus

1 0

Add to SearchBar firefox plugin.
by Vincent Davis June 29, 2010

June 29, 2010

Just thought I would share this. Sorry if you consider it noise :-) Allows you to add almost any search box as a firefox searchbar option, like scipy bug tracker search ...... https://addons.mozilla.org/en-US/firefox/addon/3682/ Vincent

1 0

Understanding source of docs, source file ".py" or ".rst" file
by Vincent Davis June 29, 2010

June 29, 2010

There are two sets of docs for scipy.constants.constants, that in the source constants.py and that in that in scipy/doc/source/constants.rst. The rst file is what is used on the website but I am not sure how to (can you) view it from an interactive py session. Is it necessary to maintain both? Or more generally whats the protocol for how and why this is done. Thanks Vincent

4 6

Status of trunk
by Gael Varoquaux June 29, 2010

June 29, 2010

Sorry for the silly question, but I have been unable to follow properly the mailing lists lately. What is the current status of trunk? I am sitting with Robert Cimrman and we would like to push a few new feature in (sparse matrix stuff), but it does not seem right to be pushing in what is going to be released soon. I am just asking to plan a bit how Robert and I are going to work on this. Cheers, Gaël PS: Basically, we want to tackle ticket #1057

2 2

FAILED (KNOWNFAIL=11, SKIP=39, errors=6, failures=1), py27 64bit osx
by Vincent Davis June 29, 2010

June 29, 2010

Not sure if these are know, I looked for the failure in the tickets but found only old ones. Vincent ERROR: test_decomp.test_lapack_misaligned(<function solve at 0x10253e140>, (array([[ 1.734e-255, 8.189e-217, 4.025e-178, 1.903e-139, 9.344e-101, 4.422e-062, 2.169e-023, 1.028e+016, 5.035e+054, 2.388e+093], [ 1.169e+132, 5.549e+170, 2.713e+209, 1.290e+248, 6.299e+286, -9.272e-292, -4.524e-253, -2.155e-214, -1.050e-175, -5.007e-137], [ -2.438e-098, -1.163e-059, -5.659e-021, -2.703e+018, -1.314e+057, -6.281e+095, -3.049e+134, -1.459e+173, -7.078e+211, -3.391e+250], [ -1.643e+289, 3.524e-294, 1.734e-255, 8.189e-217, 4.025e-178, 1.903e-139, 9.344e-101, 4.422e-062, 2.169e-023, 1.028e+016], [ 5.035e+054, 2.388e+093, 1.169e+132, 5.549e+170, 2.713e+209, 1.290e+248, 6.299e+286, -9.272e-292, -4.524e-253, -2.155e-214], [ -1.050e-175, -5.007e-137, -2.438e-098, -1.163e-059, -5.659e-021, -2.703e+018, -1.314e+057, -6.281e+095, -3.049e+134, -1.459e+173], [ -7.078e+211, -3.391e+250, -1.643e+289, 3.524e-294, 1.734e-255, 8.189e-217, 4.025e-178, 1.903e-139, 9.344e-101, 4.422e-062], [ 2.169e-023, 1.028e+016, 5.035e+054, 2.388e+093, 1.169e+132, 5.549e+170, 2.713e+209, 1.290e+248, 6.299e+286, -9.272e-292], [ -4.524e-253, -2.155e-214, -1.050e-175, -5.007e-137, -2.438e-098, -1.163e-059, -5.659e-021, -2.703e+018, -1.314e+057, -6.281e+095], [ -3.049e+134, -1.459e+173, -7.078e+211, -3.391e+250, -1.643e+289, 3.524e-294, 1.734e-255, 8.189e-217, 4.025e-178, 1.903e-139]]), array([ NaN, NaN, NaN, NaN, NaN, -3.264e+62, 6.682e+23, 9.710e-15, -1.988e-53, -2.520e-94])), {'overwrite_a': True, 'overwrite_b': True}) test_decomp.test_lapack_misaligned(<function solve at 0x10253e140>, (array([[ 1.734e-255, 8.189e-217, 4.025e-178, 1.903e-139, 9.344e-101, ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/nose/case.py", line 186, in runTest self.test(*self.arg) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/linalg/tests/test_decomp.py", line 1074, in check_lapack_misaligned func(*a,**kwargs) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/linalg/basic.py", line 49, in solve a1, b1 = map(asarray_chkfinite,(a,b)) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/lib/function_base.py", line 528, in asarray_chkfinite "array must not contain infs or NaNs") ValueError: array must not contain infs or NaNs ====================================================================== ERROR: test_none (test_signaltools.TestMedFilt) Ticket #1124. ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/signal/tests/test_signaltools.py", line 289, in test_none signal.medfilt(None) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/signal/signaltools.py", line 317, in medfilt return sigtools._order_filterND(volume,domain,order) ValueError: order_filterND not available for this type ====================================================================== ERROR: test_complex_nonsymmetric_modes (test_arpack.TestEigenComplexNonSymmetric) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py", line 267, in test_complex_nonsymmetric_modes self.eval_evec(m,typ,k,which) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py", line 248, in eval_evec eval,evec=eigen(a,k,which=which) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 397, in eigen params.iterate() File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 222, in iterate raise RuntimeError("Error info=%d in arpack" % self.info) RuntimeError: Error info=-8 in arpack ====================================================================== ERROR: test_nonsymmetric_modes (test_arpack.TestEigenNonSymmetric) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py", line 204, in test_nonsymmetric_modes self.eval_evec(m,typ,k,which) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py", line 186, in eval_evec eval,evec=eigen(a,k,which=which,**kwds) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 397, in eigen params.iterate() File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 222, in iterate raise RuntimeError("Error info=%d in arpack" % self.info) RuntimeError: Error info=-8 in arpack ====================================================================== ERROR: test_starting_vector (test_arpack.TestEigenNonSymmetric) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py", line 214, in test_starting_vector self.eval_evec(self.symmetric[0],typ,k,which='LM',v0=v0) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py", line 186, in eval_evec eval,evec=eigen(a,k,which=which,**kwds) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 397, in eigen params.iterate() File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 222, in iterate raise RuntimeError("Error info=%d in arpack" % self.info) RuntimeError: Error info=-8 in arpack ====================================================================== FAIL: test_complex_symmetric_modes (test_arpack.TestEigenComplexSymmetric) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py", line 156, in test_complex_symmetric_modes self.eval_evec(self.symmetric[0],typ,k,which) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py", line 145, in eval_evec assert_array_almost_equal(eval,exact_eval,decimal=_ndigits[typ]) File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/testing/utils.py", line 774, in assert_array_almost_equal header='Arrays are not almost equal') File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/testing/utils.py", line 618, in assert_array_compare raise AssertionError(msg) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([ 1.07188725 +6.23436023e-08j, 4.91291142 -3.25412906e-08j], dtype=complex64) y: array([ 5.+0.j, 6.+0.j], dtype=complex64) >> raise AssertionError('\nArrays are not almost equal\n\n(mismatch 100.0%)\n x: array([ 1.07188725 +6.23436023e-08j, 4.91291142 -3.25412906e-08j], dtype=complex64)\n y: array([ 5.+0.j, 6.+0.j], dtype=complex64)') ---------------------------------------------------------------------- Ran 4608 tests in 177.171s FAILED (KNOWNFAIL=11, SKIP=39, errors=6, failures=1) <nose.result.TextTestResult run=4608 errors=6 failures=1>

2 6