Mailman 3 August 2008 - SciPy-Dev

SLSQP Constrained Optimizer Status
by Rob Falck 24 Nov '13

24 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

PyArray_FromDims and friends
by Zachary Pincus 06 Oct '08

06 Oct '08

Hi all, In testing out svn scipy and numpy, I noticed some run-time errors from scipy.interpolate because the _fitpack module's c sources use PyArray_FromDims and PyArray_FromDimsAndData, which are now deprecated in numpy svn. I opened a ticket and made a patch for this particular case, but I'm not sure if there's an overall strategy (someone writes a very good regexp, say, and lets it loose), or if these will be fixed piecemeal. If the latter, here's the patch: http://scipy.org/scipy/scipy/ticket/723 Zach

7 9

Question about Kaiser Implementation in firwin
by Buehler, Eric (AGRE) 04 Oct '08

04 Oct '08

Hello, I have been using the signal.firwin FIR design tool and I have a question about the implementation. Scipy-0.4.8 Numpy-0.9.6 I am generating a kaiser window with the following parameters: N = 128 Cutoff = ~4e7 width = .3 All of the special kaiser generation parameters in the function follow Oppenheim and Schafer well. Even the sinc function for linear phase is follows exactly. However, the last line of the firwin function is causing me some heartburn. filter_design.py 1538 win = get_window(window,N,fftbins=1) 1539 alpha = N//2 1540 m = numpy.arange(0,N) 1541 h = win*special.sinc(cutoff*(m-alpha)) 1542 return h / sum(h) Line 1542 of filter_design.py, "return h / sum(h)", normalizes the function where it doesn't seem necessary, at least in the kaiser window case. Without the normalization, the kaiser window already returns a value of 1 at the zero frequency point. This normalization scales all of the data, making the window difficult to use in the frequency domain. Can someone point me to the rationale for this line? Looking at the code, this seems to be a pretty recent change (within the last year/year and a half). Thanks, Eric Buehler eric <dot> buehler <at> smiths-aerospace <dot> com ****************************************** The information contained in, or attached to, this e-mail, may contain confidential information and is intended solely for the use of the individual or entity to whom they are addressed and may be subject to legal privilege. If you have received this e-mail in error you should notify the sender immediately by reply e-mail, delete the message from your system and notify your system manager. Please do not copy it for any purpose, or disclose its contents to any other person. The views or opinions presented in this e-mail are solely those of the author and do not necessarily represent those of the company. The recipient should check this e-mail and any attachments for the presence of viruses. The company accepts no liability for any damage caused, directly or indirectly, by any virus transmitted in this email. ******************************************

3 2

Roadmap scipy 0.7
by Nils Wagner 30 Aug '08

30 Aug '08

Hi all, I had a short look at http://projects.scipy.org/scipy/scipy/roadmap Is there a new deadline for scipy 0.7 and numpy 1.2 ? Cheers, Nils

4 8

hypergeom or me?
by josef.pktd＠gmail.com 29 Aug '08

29 Aug '08

Is this correct? scipy 0.6.0: >>> stats.hypergeom.rvs(33, 18, 6, size=10) array([3, 4, 4, 4, 3, 3, 2, 2, 4, 3]) >>> stats.hypergeom.rvs(33, 6, 18, size=10) array([24, 24, 25, 22, 24, 23, 24, 25, 25, 23]) >>> stats.hypergeom.pmf(range(20),33, 18, 6) array([ 0.00451891, 0.04880423, 0.18856179, 0.33522095, 0.29009506, 0.11603802, 0.01676105, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ]) >>> stats.hypergeom.pmf(range(20),33, 6, 18) array([ 0.00451891, 0.04880423, 0.18856179, 0.33522095, 0.29009506, 0.11603802, 0.01676105, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ]) >>> stats.hypergeom.cdf(range(20),33, 18, 6) array([ 0.00451891, 0.05332314, 0.24188492, 0.57710588, 0.86720093, 0.98323895, 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. ]) >>> stats.hypergeom.cdf(range(20),33, 6, 18) array([ 0.00451891, 0.05332314, 0.24188492, 0.57710588, 0.86720093, 0.98323895, 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. ]) Josef

1 1

[mailinglist] Re: NNLS
by Uwe Schmitt 29 Aug '08

29 Aug '08

Alan Isaac worte: >Uwe Schmitt wrote: >>/ The NNLS code is now via SVN at />>/ http://public.procoders.net/nnls/nnls_with_f2py/ />>/ How can I contribute this code now ? / >Did you get the needed information for this? No. But I had some problems receiving mails from this list, maybe my spamfilter ate them. >>/ Is there further any interest in code for />>/ * ICA (Independent componenent ananlysis) ? />>/ I wrapped existing C-Code with f2py. />>/ * NMF/NNMA (nonnegative matrix factorization / - approximation) ? />>/ which is pure Python/numpy code. / >I'd like to see the latter find its way into SciPy. Nice. I'm just playing with it and try to implement some variations as Sparse NNMA... >(I'm not familiar with ICA; what's the application >area?) ICA means "Independent Component Analysis" which can be used for solving the cocktailparty problem or for blind deconvolution. Look at http://www.procoders.net/?p=30 for further information. Greetings, Uwe >Cheers, >Alan Isaac -- Dr. rer. nat. Uwe Schmitt F&E Mathematik mineway GmbH Science Park 2 D-66123 Saarbrücken Telefon: +49 (0)681 8390 5334 Telefax: +49 (0)681 830 4376 uschmitt(a)mineway.de www.mineway.de Geschäftsführung: Dr.-Ing. Mathias Bauer Amtsgericht Saarbrücken HRB 12339

2 1

dae solvers
by Benny Malengier 29 Aug '08

29 Aug '08

Hello, I need a differential algebraic equation solver, and did not see this already in scipy. I started a new module, dae.py, based on ode.py to interface with ddaspk.f and possibly also lsodi.f from netlib.org (less good, but there is a patch in the tracker) I think this approach is better than the patch given in http://projects.scipy.org/scipy/scipy/ticket/615 for lsodi, as DAE essentially mean to solve for G(y', y, t) =0 and does not fit in the y' = f(y,t) scheme needed for ode.py Is there broader interest for this? That is, should I take the extra effort to create a nice patch for scipy. Time is money and all that. Benny Malengier

3 7

ticket 422 applies for other methods: vectorize and *args bugs
by josef.pktd＠gmail.com 29 Aug '08

29 Aug '08

I think I found more errors with vectorize (sdf) and *args in scipy.stats.distributions: When I try with scipy 0.6.0 print 'rice.cdf', stats.rice.cdf(5,1.5) raise ValueError, "mismatch between python function inputs"\ ValueError: mismatch between python function inputs and received arguments if I add in class rv_continuous def _cdf(self, x, *args): + self.veccdf.nin = self.numargs+1 #JP return self.veccdf(x,*args) I also needed to add def _cdf_single_call(self, x, *args): + import scipy.integrate #JP return scipy.integrate.quad(self._pdf, self.a, x, args=args)[0] then the error is gone: >>> print 'rice.cdf', stats.rice.cdf(5,1.5) rice.cdf 0.999557350428 in order to get rice.rvs(1.5) to work, I did also def _ppf(self, q, *args): + self.vecfunc.nin = self.numargs+1 #JP return self.vecfunc(q,*args) def _isf(self, q, *args): + self.vecfunc.nin = self.numargs+1 #JP return self.vecfunc(1.0-q,*args) now, print 'rice.cdf', stats.rice.cdf(5,1.5) print 'rice.rvs',stats.rice.rvs([1.5]) works but this does not print 'rice.rvs',stats.rice.rvs([1.5]) print 'rice.cdf', stats.rice.cdf(5,1.5) cdf needs to be called before rvs. So my changes are not in the right position, the cdf needs to be initialized before, I can call rvs. But somewhere along these lines, it should be possible to make the generic methods function correctly. This was with scipy 0.60, but I haven't seen any relevant changes in the current trunk Josef for discrete distribution such as zipf, see my other email

2 1

Re: [SciPy-dev] Problem with F distribution, or with me? - ticket 422 might apply for discrete distributions
by josef.pktd＠gmail.com 28 Aug '08

28 Aug '08

On Tue, Aug 26, 2008 at 8:18 PM, <josef.pktd(a)gmail.com> wrote: > Can somebody run: > > stats.zipf.pmf(range(10),1.5) > stats.zipf.cdf(10,1.5) > stats.zipf.cdf(range(10),1.5) > > I get >>>> stats.zipf.cdf(10,1.5) > Traceback (most recent call last): > File "<stdin>", line 1, in ? > File "C:\Josef\_progs\Subversion\scipy_trunk\dist\Programs\Python24\Lib\site-p > ackages\scipy\stats\distributions.py", line 3549, in cdf > place(output,cond,self._cdf(*goodargs)) > File "C:\Josef\_progs\Subversion\scipy_trunk\dist\Programs\Python24\Lib\site-p > ackages\scipy\stats\distributions.py", line 3458, in _cdf > return self._cdfvec(k,*args) > File "C:\Programs\Python24\Lib\site-packages\numpy\lib\function_base.py", line > 1092, in __call__ > raise ValueError, "mismatch between python function inputs"\ > ValueError: mismatch between python function inputs and received arguments > > but this might be, because I am running either an older version of > scipy or from my own faulty (?) build from subversion. > If the current version still has this error, then I think it is > related to the generic calculation of the cdf, > similar to ticket 422, changeset 3797 for the moment calculation. > > In my tests I get many value errors, but again I don't know whether it > is my version/setup, whether the parameters are not ruled out but > don't make sense, or whether there is actually a bug somewhere. > > Josef > in changeset 3797 this was added to the continuous random variable: 315 self.generic_moment.nin = self.numargs+1 # Because of the *args argument 316 # of _mom0_sc, vectorize cannot count the number of arguments correctly. in current trunk these are lines 318 and 319 I don't understand the details of numpy and vectorize, but I think that the same problem with the number of arguments also applies to the generic calculation for the discrete distribution, i.e. lines 3401 self._ppf = new.instancemethod(sgf(_drv_ppf,otypes='d'), 3402 self, rv_discrete) 3403 self._pmf = new.instancemethod(sgf(_drv_pmf,otypes='d'), 3404 self, rv_discrete) 3405 self._cdf = new.instancemethod(sgf(_drv_cdf,otypes='d'), 3406 self, rv_discrete) since all _drv_??? methods also use the *args argument and the exception message I get looks the same as in ticket 422. Josef

1 2

Re: [SciPy-dev] Problem with F distribution, or with me? - error in stats.fatiguelife.rvs
by josef.pktd＠gmail.com 27 Aug '08

27 Aug '08

It looks like that there is an error in stats.fatiguelife.rvs Kolmogorov test fails >>> stats.kstest('fatiguelife','',args=(5,),N=1000) (0.093216666807115545, array(2.5853082230575808e-008)) Mean of sample >>> stats.fatiguelife.stats(5,moments='m') array(13.5) >>> np.mean(stats.fatiguelife.rvs(5,size=1000)) 26.683858360164475 >>> np.mean(stats.fatiguelife.rvs(5,size=10000)) 26.841525716395847 >>> np.mean(stats.fatiguelife.rvs(5,size=100000)) 26.730694604009678 >>> np.mean(stats.fatiguelife.rvs(5,size=100000)/2) 13.469823793800416 >>> stats.fatiguelife.stats(3,moments='m') array(5.5) >>> np.mean(stats.fatiguelife.rvs(3,size=100000)) 10.922712537094393 >>> np.mean(stats.fatiguelife.rvs(3,size=100000)/2) 5.5340854278553246 Variance of sample >>> stats.fatiguelife.stats(3,moments='v') array(110.25) >>> np.var(stats.fatiguelife.rvs(3,size=1000000)) 440.1793356094052 >>> np.var(stats.fatiguelife.rvs(3,size=1000000)/2) 110.445022957997 >>> np.var(stats.fatiguelife.rvs(3,size=10000000)/2) 110.03364894832275 >>> stats.fatiguelife.stats(5,moments='v') array(806.25) >>> np.var(stats.fatiguelife.rvs(5,size=1000000)) 3222.4271388000293 >>> np.var(stats.fatiguelife.rvs(5,size=1000000)/2) 809.29193071702855 theoretical mean and cdf look correct, according to http://www.itl.nist.gov/div898/handbook/eda/section3/eda366a.htm but random number generator, is wrong by approximately the scale of 1/2 Josef

3 5