Mailman 3 April 2013 - SciPy-Dev

Bootstrap confidence limits code
by Constantine Evans 05 Oct '14

05 Oct '14

Hello everyone, A few years ago I implemented a scikit for bootstrap confidence limits (https://github.com/cgevans/scikits-bootstrap). I didn’t think much about it after that until recently, when I realized that some people are actually using it, and that there’s apparently been some talk about implementing this functionality in either scipy.stats or statsmodels (I should thank Randal Olson for discussing this and bringing it to my attention). As such I’ve rewritten most of the code, and written up some docstrings. The current code can do confidence intervals with basic percentile interval, bias-corrected accelerated, and approximate bootstrap confidence methods, and can also provide bootstrap and jackknife indexes. Most of it is implemented from the descriptions in Efron and Tibshirani’s Introduction to the Bootstrap, but the ABC code at the moment is a port from the modified-BSD-licensed bootstrap package for R (not the boot package) as I’m not entirely confident in my understanding of the method. And so, I have a few questions for everyone: * Is there any interest in including this sort of code in either scipy.stats or statsmodels? If so, where do people think would be the better place? The code is relatively small; at the moment it is less than 200 lines, with docstrings probably making up 100 of those lines. * Also, if so, what would need to be changed, added, and improved beyond what is mentioned in the Contributing to Scipy part of the reference guide? I’m never a fan of my own code, and imagine quite a bit would need to be fixed; I know tests will need to be added too. In addition, I have a few questions about what would be better practice for the API, and I haven’t really found a guide on best practices for Scipy: * When I started writing the code, I wrote a single function ci for confidence intervals, with a method argument to choose the method. This is easy for users, especially so that they don’t have to look through documentation to realize that BCA is the most generally useful method (at least from everything I’ve read) and that there really isn’t any reason to use many of the simpler methods. However, ABC takes different paramenters, and needs a statistic function that takes weights, which makes this single-function organization trickier. At the moment, I have a separate function for ABC. Would it be better to split up all the methods to their own functions? * ABC requires a statistic function that takes weights. I’ve noticed that things like np.average takes a weights= argument. Would it be better to require input of a stat(data,weights) function, or input of a stat(data,weights=) with weights as a named argument? The latter would be nice in terms of allowing the same function to be used for all methods, but would make it impossible to use a lambda for the function. Is there some other method of doing this entirely? * Are there any missing features that anyone thinks should be added? I apologize if much of this is answered elsewhere, I just haven’t found any of it; I also apologize if this is far too long-winded and confusing! Regards, Constantine Evans

4 3

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

Vector Strength function
by Todd 22 Nov '13

22 Nov '13

I am interested in implementing a function for scipy. The function is called "vector strength". It is basically a measure of how reliably a set of events occur at a particular phase. It was originally developed for neuroscience research, to determine how well a set of neural events sync up with a periodic stimulus like a sound waveform. However, it is useful for determining how periodic a supposedly periodic set of events really are, for example: 1. Determining whether crime is really more common during a full moon and by how much 2. Determining how concentrated visitors to a coffee shop are during rush hour 3. Determining exactly how concentrated hurricanes are during hurricane season My thinking is that this could be implemented in stages: First would be a Numpy function that would add a set of vectors in polar coordinates. Given a number of magnitude/angle pairs it would provide a summed magnitude/angle pair. This would probably be combined with a cartesian<->polar conversion functions. Making use of this function would be a scipy function that would actually implement the vector strength calculation. This is done by treating each event as a unit vector with a phase, then taking the average of the vectors. If all events have the same phase, the result will have an amplitude of 1. If they all have a different phases, the result will have an amplitude of 0. It may even be worth having a dedicated polar dtype, although that may be too much. What does everyone think of this proposal?

4 12

issues trac migration review
by Skipper Seabold 15 May '13

15 May '13

Many of the issues have been migrated. Please have a look and let me know of any problems. https://github.com/jseabold/scipy-trac-migration/issues?state=open I got rate limited by github, so I will have to continue later tonight to make sure there are no problems. You'll notice that I used the same @ -> atmention: substitution for now so that people don't get notifications, as I'm working out the kinks. I used the same trac -> github mapping as the numpy migration and added several more by hand of names I recognized as being a duplicate e-mail address of a regular contributor. If we want different issue labels, let me know. It's probably easier to format them now than to do it later. Anything else? Skipper

9 60

Multilinear interpolation
by Pablo Winant 10 May '13

10 May '13

Hi, Continuing a series of thread about interpolation, I have an almost ready implementation of multilinear implementation, that is linear interpolation in each dimension. I have been toying with several implementations - in pure numpy (quite memory hungry) - in cython - using pycuda Obviously, if there is some interest for it in scipy, it would be a good motivation to finish it. There are two main differences with respect to the existing interp routine : - it is independent of fitpack - it extrapolates linearly while the existing routine returns a constant (or NaN) Now I have some questions about it: - in this case I several implementations, in the same library that can potentially be enabled on different systems (with a compiler, with a gpu). I made a python file that tries to import one version after all shadowing slower ones by faster ones. It is meant to replicate the .m/.mex file behavior of Matlab. Is there a better approach ? - the cython code is meant to operate on double* inputs, and is limited to 4 dimensions. It is however very easy to generate single or higher order code. What would you use for that purpose ? Best, Pablo

4 13

GSoC Proposal draft -- Improvements to the sparse package of Scipy: support for bool dtype and better interaction with Numpy
by Blake Griffith 08 May '13

08 May '13

Seeking feedback, comments, and criticism Link: https://github.com/cowlicks/GSoC-proposal/blob/master/proposal.markdown

3 4

Request for extension to scipy.integrate
by Nathan Woods 06 May '13

06 May '13

SciPy's multidimensional integration capabilities are somewhat limited, as mentioned previously: https://github.com/scipy/scipy/issues/2098. Although dblquad and tplquad attempt to address this problem, they do not allow any detailed control over the underlying quad algorithm and there is no option for 4-dimensional integration at all. One obvious instance where this functionality would be desired is in the evaluation of volume integrals in 4-dimensional space-time. Another instance is the integration of discontinuous functions, where access to the "points" option of quad is critical. As it stands, SciPy does not provide any good functionality for either of these situations. I've written a recursive wrapper function for quad that resolves both of these problems, allowing n-dimensional integration with complete specification of options for quad. I've included a test problem that demonstrates the resolution of both of the above problems (just execute the script). I've done my best to document everything, and I've verified the result of the test problem against Mathematica. I would like this code (or equivalent functionality) to be included in the integrate package. I'm open to any feedback or suggestions, particularly with the interface. Nathan Woods """ Recursive integration using SciPy's quad function. Contains the following: class Error - module wrapper for Exception. For future expansion only. function mul_quad - Evaluate multiple integrals using recursive calls to quad. """ from scipy.integrate import quad import numpy class Error(Exception): pass def mul_quad(func,ranges,args=(),opts=(),_depth=0,_int_vars=()): """ Evaluate multiple integrals through recursive calls to scipy.integrate.quad. mul_quad takes care of the programming required to perform nested integration using the 1-dimensional integration routines provided by SciPy's adaptive quadrature function. It extends the capabilities of dblquad and tplquad by allowing for more levels of integration, and allowing the user to specify nearly the full range of options allowed by quad, for each level of integration. Users are cautioned that nested Gaussian integration of this kind is computationally intensive, and may be unsuitable for many nested integrals. Usage: mul_quad(func,ranges,args=(),opts=(),_depth=0,_int_vars=()) Inputs: func - callable, acceptable to SciPy's quad, returning a number. Should accept a float, followed by the contents of _int_vars and args, e.g. if x is a float, args=(1.4,string), and _int_vars = (1,1,3), then func should be of the form func(x,1,1,3,1.4,string). ranges - sequence describing the ranges of integration. Integrals are performed in order, so ranges[0] corresponds to the first argument of func, ranges[1] to the second, and so on. Each element of ranges may be either a constant sequence of length 2 or else a function that returns such a sequence. If a function, then it will be called with all of the integration arguments available to that point. e.g. for func = f(x0,x1,x2,x3), the range of integration for x0 may be defined as either a constant such as (0,1) or as a function range0(x1,x2,x3). The functional range of integration for x1 will be range1(x2,x3), x2 will be range2(x3), and so on. args - optional sequence of arguments. Contains only arguments of func beyond those over which the integration is being performed. opts - optional sequence of options for SciPy's quad. Each element of opts may be specified as either a dictionary or as a function that returns a dictionary similarly to ranges. opts must either be left empty (), or it must be the same length as ranges. Options are passed in the same order as ranges, so opts[0] corresponds to integration over the first argument of func, and so on. The full_output option from quad is not available, due to the difficulty of consolidating the large number of additional outputs. For reference, the default options from quad are: - epsabs = 1.49e-08 - epsrel = 1.49e-08 - limit = 50 - points = None - weight = None - wvar = None - wopts = None (As of Apr 2013) _depth - used to determine level of integration. Should be omitted by the user, except for debugging purposes. _int_vars - contains values of integration variables in inner integration loops. Should not be used manually except for debugging. Returns: out - value of multiple integral in the specified range. abserr - estimate of the absolute error in the result. The maximum value of abserr among all the SciPy quad evaluations. """ global abserr if _depth == 0: abserr = None if not (len(opts) in [0,len(ranges)]): raise Error('opts must be given for all integration levels or none!') total_args = _int_vars+args # Select the range and opts for the given depth of integration. ind = -_depth-1 if callable(ranges[ind]): current_range = ranges[ind](*total_args) else: current_range = ranges[ind] if len(opts) != 0: if callable(opts[ind]): current_opts = opts[ind](*total_args) else: current_opts = opts[ind] else: current_opts = {} try: if current_opts["full_output"] != 0: raise Error('full_output option is disabled!') except(KeyError): pass # Check to make sure that all points lie within the specified range. try: for point in current_opts["points"]: if point < current_range[0] or point > current_range[1]: current_opts["points"].remove(point) except(KeyError): pass if current_range is ranges[0]:# Check to see if down to 1-D integral. func_new = func else: # Define a new integrand. def func_new(*_int_vars): return mul_quad(func,ranges,args=args,opts=opts, _depth=_depth+1,_int_vars=_int_vars) out = quad(func_new,*current_range,args=_int_vars+args,**current_opts) # Track the maximum value of abserr from quad if abserr is None: abserr = out[1] if out[1] > abserr: abserr = out[1] if _depth == 0: return out[0],abserr else: return out[0] if __name__=='__main__': func = lambda x0,x1,x2,x3 : x0**2+x1*x2-x3**3+numpy.sin(x0)+( 1 if (x0-.2*x3-.5-.25*x1>0) else 0) points=[[lambda (x1,x2,x3) : .2*x3+.5+.25*x1],[],[],[]] def opts0(*args,**kwargs): return {'points':[.2*args[2]+.5+.25*args[0]]} print mul_quad(func,[[0,1],[-1,1],[.13,.8],[-.15,1]], opts=[opts0,{},{},{}])

2 3

Re: [SciPy-Dev] GSoC Proposal draft -- Improvements to the sparse package of Scipy: support for bool dtype and better interaction with Numpy
by Daniel Smith 01 May '13

01 May '13

>> Seeking feedback, comments, and criticism >> >> Link: >> >> https://github.com/cowlicks/GSoC-proposal/blob/master/proposal.markdown >> >> >> > Do any of the sparse folks have some feedback here? > > Chuck I don't know much about the details/scope of the GSoC or their proposal process, but I can give some thoughts about the proposed changes. Adding the bool data type to would be cool. I think pv has already discussed some of the issues with that on here. For the ufuncs improvements, I would just recommend that you be careful about potential ndarray/matrix incompatibilities. The SciPy sparse class has a bunch of backends and an unpleasant amount of repeated code to implement those interfaces. Consolidating interface code while adding the ufunc functionality would be very helpful. Even if you only consolidate some CSC/CSR repeated code, that could make longer term upkeep a lot easier. NumPy has moved towards the matrix class being a wrapper that simply changes the interface to the existing ndarray functionality. It seems that they are also trying to deprecate the use of matrices. Moving the sparse class closer to a universal ndarray interface would be an excellent goal, but just adding another interface could cause more problems than help. Daniel

2 2

SciPy Central in GSoC 2013
by Surya Kasturi 30 Apr '13

30 Apr '13

Hi everyone, How about SciPy Central on GSoC 2013? I was thinking to participate in GSoC this year and thought of contributing to this project as GSoC student. Firstly, I was luck enough to find SPC which closely matched my interests. Since, I am already slightly familiar with the existing code, made some small contributions, I thought it would be very nice to continue my work as GSoC student during summer. Is it possible? Surya

6 35

GSoC proposal -- Help!
by Blake Griffith 30 Apr '13

30 Apr '13

Hello, sorry to send this out again, but no one responded last time. I'm proposing to work on scipy sparse for this GSoC. My proposal is titled: Improvements to the sparse package of Scipy: support for bool dtype and better interaction with NumPy Please let me know if you think there is anything I can do to improve it! Link: https://github.com/cowlicks/GSoC-proposal/blob/master/proposal.markdown

1 0