Mailman 3 October 2013 - SciPy-Dev

Linalg: function for diagonalize hermitian commuting matrices simultaneously
by Lars-Hendrik Frahm 08 Oct '14

08 Oct '14

Dear SciPy Developers, I like to contribute a function for simultaneously diagonalize a set of hermitian commuting matrices. Especially for matrices with degenerate eigenvalues. Is there need for something like that in SciPy ? I think it could be helpful for someone else and I didn't find it in recent SciPy versions or on the mailing list. Hope to hear from you.

3 3

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

1.0 roadmap: weave
by Arnd Baecker 06 May '14

06 May '14

Hi, > On Sat, Sep 21, 2013 at 8:54 PM, Ralf Gommers <ralf.gommers(a)gmail.com> wrote: [...] > weave > ````` > This is the only module that was not ported to Python 3. Effectively it's > deprecated (not recommended to use for new code). In the future it should be > removed from scipy (can be made into a separate module). Let me just mention that I am worried about the removal of weave: a substantial amount of code developed in our group over the years uses weave. Porting this to cython seems not trivial, in particular, because many of the original authors (diploma students, PhD students, ...) have left academia by now. Well, as porting weave to Python 3 is way out of my capabilities, this is just a comment in the hope that someone could give this a try (if at some point testing would be needed, I could try to help, of course). Best, Arnd

9 12

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

Re: [SciPy-Dev] [SciPy-User] Linear Programming via Simplex Algorithm
by Rob Falck 03 Nov '13

03 Nov '13

Hello Ralf, First, my apologies that this went to SciPy Users, it was intended for the Developers list and I'm redirecting this message there. On Tue, Oct 22, 2013 at 2:08 PM, Ralf Gommers <ralf.gommers(a)gmail.com>wrote: > > > On Mon, Oct 21, 2013 at 3:59 AM, Rob Falck <robfalck(a)gmail.com> wrote: > >> Hello, >> >> I've spent some time recently polishing a simplex-based linear >> programming function. I've seen traffic from time to time about including >> such functionality in scipy.optimize but it always seems to have been >> closed without inclusion. >> > > Hi Rob, I assume you have seen https://github.com/scipy/scipy/pull/218then? Looks like it's functionality that we'd like to see in scipy.optimize > but needs some more effort than anyone has been willing to spend so far. > I've seen pull request 218 and parts of it are very good. The simplex method itself is lean, but it doesn't support equality or lower-bound constraints. It also differs in that it returns a new class LPResult. In my approach I had tried to stay with the general scipy.optimize.Result setup and use only those fields which are applicable. Do you have any guidance on which way would be preferable? > My intent was to develop a linear-programming routine that, given a >> non-standard linear programming problem, generates a canonical tableau and >> solves it using the simplex algorithm. By nonstandard I mean that variables >> may have negative values, and three types of constraints are supported >> (lower-bound, upper-bound, and equality). >> >> I've named a top-level routine "linprog". I suspect in the future that >> scipy may include other algorithms besides the simplex to solve linear >> programming problems, in which case linprog would serve as the main >> function (similarly to the way minimize serves as the interfact to all nlp >> routines). >> > > A generic `linprog` interface sounds like a good idea. It looks like the > API and the current implementation are fairly specific to your lpsimplex > algorithm however. Compare how little minimize() does to the 400 LoC in > linprog(). If you want to make linprog() generic maybe it would help if you > would consider how you could integrate the algorithm of > https://github.com/scipy/scipy/pull/218 into it without changing the API > besides adding a "method" keyword. > I'd be happy to incorporate both, although theres a good bit of overlap between the two contributions. I'd hate for that author's work to go unused, so perhaps I can replace the core simplex solving routine of my code with that one. It's leaner and honestly I like a few things about it better than mine. That way my top-level simplex routine would rework a problem with lower-bound and equality constraints into standard form, to be solved by that underlying simplex routine. Also, my implementation of linprog is long because it contains a great deal of error checking. I will rework a top-level linprog interface which leaves most of the error-checking to the underlying linear programming routines. > For example, consider a relatively simple problem: >> >> Maximize: f = 3*x[0] + 2*x[1] >> Subject to: 2*x[0] + 1*x[1] <= 10 >> 1*x[0] + 1*x[1] <= 8 >> 1*x[0] + 0*x[1] <= 4 >> >> where: 0 <= x[0] <= inf >> 0 <= x[1] <= inf >> >> The inputs are such that the objective is specified by array c where f = >> dot(c,x). >> We have three upper-bound constraints, which we define using dot(A_ub,x) >> <= b_ub >> >> >>>c = [3,2] >> >>>b_ub = [10,8,4] >> >>>A_ub = [[2,1], >> >>> [1,1], >> >>> [1,0]] >> >>>res=linprog(c=c,A_ub=A_ub,b_ub=b_ub,objtype='max',disp=True) >> Optimization terminated successfully. >> Current function value: 18.000000 >> Iterations: 3 >> >> I've written a suite of unit tests and have documented all functions. >> It's somewhat non-standard but I've also included two callback functions. >> linprog_verbose_callback prints the simplex tableau, pivot element, and >> basic variables at each iteration of the simplex method. >> linprog_terse_callback simply prints the value of the solution vector x at >> each iteration. >> > > Do the callbacks need to be separate public functions? My impression is > no. And what about lpsimplex()? If linprog() is the generic interface maybe > lpsimplex can be private? > > I felt like a verbose callback that displays the simplex tableau at each iteration would be useful for people who want it. I can remove it and include it as an example callback in documentation or something like that. lpsimplex will be made private (or more likely replaced by the implementation from pull request 218) The only top-level function will be linprog, which like minimize, will support a "method" argument. For now, the only valid method will be "simplex". > The code is currently in the linprog branch at >> https://github.com/robfalck/scipy/tree/linprog/scipy >> (relevant files are scipy/optimize/linprog.py, >> scipy/optimize/__init__.py, and scipy/optimize/tests/test_linprog.py) >> > > This looks pretty good from a quick browse (disclaimer: I haven't tried to > understand your algorithm in detail). > > >> I welcome constructive criticism of the algorithm. It's been designed >> to detect degenerate cases due to unbounded problems, and it has a >> relatively basic way to avoid cycling in the simplex solution. >> > > Have you benchmarked your algorithm against other implementations? And/or > checked the efficiency (typical should be 2N to 3N iterations, with N the > size of the constraint matrix)? > I've checked my algorithm against some textbook examples and a few randomly generated cases from http://demonstrations.wolfram.com/TwoPhaseSimplexMethod/ Basically there are two important things I had to get right. First, was the linprog routine (to become _linprog_simplex) creating a correct tableau given the problem at hand? Secondly, given that the tableau was correct was the algorithm choosing the correct pivots to get to the solution as quickly as possible. In both cases, from all my testing thus far, the answer has been yes. I've checked a basic cycling case and was able to successfully break out of the cycling to converge to the correct solution, although I wouldn't say my method of avoiding cycling is advanced. If someone has a relatively complex example to test this against, I'd be happy to try it. I doubt this implementaiton will be fast with dozens of variables and constraints, but I don't have a good feeling for where the maximum appropriate problem size should be. > If there's interest in including this I'd be happy to proceed with >> submitting a PR. I still have a bit of cleanup to perform but it's >> relatively close to being ready (pep8 compliance, etc) >> > > Sounds good to me! > > In terms of completing the cleanup, I think fixing line length to 80 chars > and breaking up linprog() into smaller logical units would be good to > before submitting imho. Also, linprog.py should be renamed to _linprog.py. > I'll will work on this and submit a pull request when I'm satisfied with the state of things. I will keep you posted. > Cheers, > Ralf > > > _______________________________________________ > SciPy-User mailing list > SciPy-User(a)scipy.org > http://mail.scipy.org/mailman/listinfo/scipy-user > > -- - Rob Falck

6 12

ANN: NumPy 1.8.0 release.
by Charles R Harris 31 Oct '13

31 Oct '13

I am pleased to announce the availability of NumPy 1.8.0. This release is the culmination of over a years worth of work by the NumPy team and contains many fixes, enhancements, and new features. Highlights are: - New, no 2to3, Python 2 and Python 3 are supported by a common code base. - New, gufuncs for linear algebra, enabling operations on stacked arrays. - New, inplace fancy indexing for ufuncs with the ``.at`` method. - New, ``partition`` function, partial sorting via selection for fast median. - New, ``nanmean``, ``nanvar``, and ``nanstd`` functions skipping NaNs. - New, ``full`` and ``full_like`` functions to create value initialized arrays. - New, ``PyUFunc_RegisterLoopForDescr``, better ufunc support for user dtypes. - Numerous performance improvements in many areas. This release requires Python 2.6, 2.7 or 3.2-3.3, support for Python 2.4 and 2.5 has been dropped. Sources and binaries can be found at http://sourceforge.net/projects/numpy/files/NumPy/1.8.0/. Some 119 people contributed to this release. This is a remarkable increase and shows that there is still life in this venerable code that had its beginning in Numeric some 18 years ago. Many thanks to you all. Enjoy, Chuck NumPy 1.8.0 Release Notes ************************* This release supports Python 2.6 -2.7 and 3.2 - 3.3. Highlights ========== * New, no 2to3, Python 2 and Python 3 are supported by a common code base. * New, gufuncs for linear algebra, enabling operations on stacked arrays. * New, inplace fancy indexing for ufuncs with the ``.at`` method. * New, ``partition`` function, partial sorting via selection for fast median. * New, ``nanmean``, ``nanvar``, and ``nanstd`` functions skipping NaNs. * New, ``full`` and ``full_like`` functions to create value initialized arrays. * New, ``PyUFunc_RegisterLoopForDescr``, better ufunc support for user dtypes. * Numerous performance improvements in many areas. Dropped Support =============== Support for Python versions 2.4 and 2.5 has been dropped, Support for SCons has been removed. Future Changes ============== The Datetime64 type remains experimental in this release. In 1.9 there will probably be some changes to make it more useable. The diagonal method currently returns a new array and raises a FutureWarning. In 1.9 it will return a readonly view. Multiple field selection from a array of structured type currently returns a new array and raises a FutureWarning. In 1.9 it will return a readonly view. The numpy/oldnumeric and numpy/numarray compatibility modules will be removed in 1.9. Compatibility notes =================== The doc/sphinxext content has been moved into its own github repository, and is included in numpy as a submodule. See the instructions in doc/HOWTO_BUILD_DOCS.rst.txt for how to access the content. .. _numpydoc: https://github.com/numpy/numpydoc The hash function of numpy.void scalars has been changed. Previously the pointer to the data was hashed as an integer. Now, the hash function uses the tuple-hash algorithm to combine the hash functions of the elements of the scalar, but only if the scalar is read-only. Numpy has switched its build system to using 'separate compilation' by default. In previous releases this was supported, but not default. This should produce the same results as the old system, but if you're trying to do something complicated like link numpy statically or using an unusual compiler, then it's possible you will encounter problems. If so, please file a bug and as a temporary workaround you can re-enable the old build system by exporting the shell variable NPY_SEPARATE_COMPILATION=0. For the AdvancedNew iterator the ``oa_ndim`` flag should now be -1 to indicate that no ``op_axes`` and ``itershape`` are passed in. The ``oa_ndim == 0`` case, now indicates a 0-D iteration and ``op_axes`` being NULL and the old usage is deprecated. This does not effect the ``NpyIter_New`` or ``NpyIter_MultiNew`` functions. The functions nanargmin and nanargmax now return np.iinfo['intp'].min for the index in all-NaN slices. Previously the functions would raise a ValueError for array returns and NaN for scalar returns. NPY_RELAXED_STRIDES_CHECKING ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ There is a new compile time environment variable ``NPY_RELAXED_STRIDES_CHECKING``. If this variable is set to 1, then numpy will consider more arrays to be C- or F-contiguous -- for example, it becomes possible to have a column vector which is considered both C- and F-contiguous simultaneously. The new definition is more accurate, allows for faster code that makes fewer unnecessary copies, and simplifies numpy's code internally. However, it may also break third-party libraries that make too-strong assumptions about the stride values of C- and F-contiguous arrays. (It is also currently known that this breaks Cython code using memoryviews, which will be fixed in Cython.) THIS WILL BECOME THE DEFAULT IN A FUTURE RELEASE, SO PLEASE TEST YOUR CODE NOW AGAINST NUMPY BUILT WITH:: NPY_RELAXED_STRIDES_CHECKING=1 python setup.py install You can check whether NPY_RELAXED_STRIDES_CHECKING is in effect by running:: np.ones((10, 1), order="C").flags.f_contiguous This will be ``True`` if relaxed strides checking is enabled, and ``False`` otherwise. The typical problem we've seen so far is C code that works with C-contiguous arrays, and assumes that the itemsize can be accessed by looking at the last element in the ``PyArray_STRIDES(arr)`` array. When relaxed strides are in effect, this is not true (and in fact, it never was true in some corner cases). Instead, use ``PyArray_ITEMSIZE(arr)``. For more information check the "Internal memory layout of an ndarray" section in the documentation. Binary operations with non-arrays as second argument ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Binary operations of the form ``<array-or-subclass> * <non-array-subclass>`` where ``<non-array-subclass>`` declares an ``__array_priority__`` higher than that of ``<array-or-subclass>`` will now unconditionally return *NotImplemented*, giving ``<non-array-subclass>`` a chance to handle the operation. Previously, `NotImplemented` would only be returned if ``<non-array-subclass>`` actually implemented the reversed operation, and after a (potentially expensive) array conversion of ``<non-array-subclass>`` had been attempted. (`bug <https://github.com/numpy/numpy/issues/3375>`_, `pull request <https://github.com/numpy/numpy/pull/3501>`_) Function `median` used with `overwrite_input` only partially sorts array ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If `median` is used with `overwrite_input` option the input array will now only be partially sorted instead of fully sorted. Fix to financial.npv ~~~~~~~~~~~~~~~~~~~~ The npv function had a bug. Contrary to what the documentation stated, it summed from indexes ``1`` to ``M`` instead of from ``0`` to ``M - 1``. The fix changes the returned value. The mirr function called the npv function, but worked around the problem, so that was also fixed and the return value of the mirr function remains unchanged. Runtime warnings when comparing NaN numbers ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Comparing ``NaN`` floating point numbers now raises the ``invalid`` runtime warning. If a ``NaN`` is expected the warning can be ignored using np.errstate. E.g.:: with np.errstate(invalid='ignore'): operation() New Features ============ Support for linear algebra on stacked arrays ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The gufunc machinery is now used for np.linalg, allowing operations on stacked arrays and vectors. For example:: >>> a array([[[ 1., 1.], [ 0., 1.]], [[ 1., 1.], [ 0., 1.]]]) >>> np.linalg.inv(a) array([[[ 1., -1.], [ 0., 1.]], [[ 1., -1.], [ 0., 1.]]]) In place fancy indexing for ufuncs ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The function ``at`` has been added to ufunc objects to allow in place ufuncs with no buffering when fancy indexing is used. For example, the following will increment the first and second items in the array, and will increment the third item twice: ``numpy.add.at(arr, [0, 1, 2, 2], 1)`` This is what many have mistakenly thought ``arr[[0, 1, 2, 2]] += 1`` would do, but that does not work as the incremented value of ``arr[2]`` is simply copied into the third slot in ``arr`` twice, not incremented twice. New functions `partition` and `argpartition` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ New functions to partially sort arrays via a selection algorithm. A ``partition`` by index ``k`` moves the ``k`` smallest element to the front of an array. All elements before ``k`` are then smaller or equal than the value in position ``k`` and all elements following ``k`` are then greater or equal than the value in position ``k``. The ordering of the values within these bounds is undefined. A sequence of indices can be provided to sort all of them into their sorted position at once iterative partitioning. This can be used to efficiently obtain order statistics like median or percentiles of samples. ``partition`` has a linear time complexity of ``O(n)`` while a full sort has ``O(n log(n))``. New functions `nanmean`, `nanvar` and `nanstd` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ New nan aware statistical functions are added. In these functions the results are what would be obtained if nan values were ommited from all computations. New functions `full` and `full_like` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ New convenience functions to create arrays filled with a specific value; complementary to the existing `zeros` and `zeros_like` functions. IO compatibility with large files ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Large NPZ files >2GB can be loaded on 64-bit systems. Building against OpenBLAS ~~~~~~~~~~~~~~~~~~~~~~~~~ It is now possible to build numpy against OpenBLAS by editing site.cfg. New constant ~~~~~~~~~~~~ Euler's constant is now exposed in numpy as euler_gamma. New modes for qr ~~~~~~~~~~~~~~~~ New modes 'complete', 'reduced', and 'raw' have been added to the qr factorization and the old 'full' and 'economic' modes are deprecated. The 'reduced' mode replaces the old 'full' mode and is the default as was the 'full' mode, so backward compatibility can be maintained by not specifying the mode. The 'complete' mode returns a full dimensional factorization, which can be useful for obtaining a basis for the orthogonal complement of the range space. The 'raw' mode returns arrays that contain the Householder reflectors and scaling factors that can be used in the future to apply q without needing to convert to a matrix. The 'economic' mode is simply deprecated, there isn't much use for it and it isn't any more efficient than the 'raw' mode. New `invert` argument to `in1d` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The function `in1d` now accepts a `invert` argument which, when `True`, causes the returned array to be inverted. Advanced indexing using `np.newaxis` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ It is now possible to use `np.newaxis`/`None` together with index arrays instead of only in simple indices. This means that ``array[np.newaxis, [0, 1]]`` will now work as expected and select the first two rows while prepending a new axis to the array. C-API ~~~~~ New ufuncs can now be registered with builtin input types and a custom output type. Before this change, NumPy wouldn't be able to find the right ufunc loop function when the ufunc was called from Python, because the ufunc loop signature matching logic wasn't looking at the output operand type. Now the correct ufunc loop is found, as long as the user provides an output argument with the correct output type. runtests.py ~~~~~~~~~~~ A simple test runner script ``runtests.py`` was added. It also builds Numpy via ``setup.py build`` and can be used to run tests easily during development. Improvements ============ IO performance improvements ~~~~~~~~~~~~~~~~~~~~~~~~~~~ Performance in reading large files was improved by chunking (see also IO compatibility). Performance improvements to `pad` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The `pad` function has a new implementation, greatly improving performance for all inputs except `mode=<function>` (retained for backwards compatibility). Scaling with dimensionality is dramatically improved for rank >= 4. Performance improvements to `isnan`, `isinf`, `isfinite` and `byteswap` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ `isnan`, `isinf`, `isfinite` and `byteswap` have been improved to take advantage of compiler builtins to avoid expensive calls to libc. This improves performance of these operations by about a factor of two on gnu libc systems. Performance improvements via SSE2 vectorization ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Several functions have been optimized to make use of SSE2 CPU SIMD instructions. * Float32 and float64: * base math (`add`, `subtract`, `divide`, `multiply`) * `sqrt` * `minimum/maximum` * `absolute` * Bool: * `logical_or` * `logical_and` * `logical_not` This improves performance of these operations up to 4x/2x for float32/float64 and up to 10x for bool depending on the location of the data in the CPU caches. The performance gain is greatest for in-place operations. In order to use the improved functions the SSE2 instruction set must be enabled at compile time. It is enabled by default on x86_64 systems. On x86_32 with a capable CPU it must be enabled by passing the appropriate flag to the CFLAGS build variable (-msse2 with gcc). Performance improvements to `median` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ `median` is now implemented in terms of `partition` instead of `sort` which reduces its time complexity from O(n log(n)) to O(n). If used with the `overwrite_input` option the array will now only be partially sorted instead of fully sorted. Overrideable operand flags in ufunc C-API ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When creating a ufunc, the default ufunc operand flags can be overridden via the new op_flags attribute of the ufunc object. For example, to set the operand flag for the first input to read/write: PyObject \*ufunc = PyUFunc_FromFuncAndData(...); ufunc->op_flags[0] = NPY_ITER_READWRITE; This allows a ufunc to perform an operation in place. Also, global nditer flags can be overridden via the new iter_flags attribute of the ufunc object. For example, to set the reduce flag for a ufunc: ufunc->iter_flags = NPY_ITER_REDUCE_OK; Changes ======= General ~~~~~~~ The function np.take now allows 0-d arrays as indices. The separate compilation mode is now enabled by default. Several changes to np.insert and np.delete: * Previously, negative indices and indices that pointed past the end of the array were simply ignored. Now, this will raise a Future or Deprecation Warning. In the future they will be treated like normal indexing treats them -- negative indices will wrap around, and out-of-bound indices will generate an error. * Previously, boolean indices were treated as if they were integers (always referring to either the 0th or 1st item in the array). In the future, they will be treated as masks. In this release, they raise a FutureWarning warning of this coming change. * In Numpy 1.7. np.insert already allowed the syntax `np.insert(arr, 3, [1,2,3])` to insert multiple items at a single position. In Numpy 1.8. this is also possible for `np.insert(arr, [3], [1, 2, 3])`. Padded regions from np.pad are now correctly rounded, not truncated. C-API Array Additions ~~~~~~~~~~~~~~~~~~~~~ Four new functions have been added to the array C-API. * PyArray_Partition * PyArray_ArgPartition * PyArray_SelectkindConverter * PyDataMem_NEW_ZEROED C-API Ufunc Additions ~~~~~~~~~~~~~~~~~~~~~ One new function has been added to the ufunc C-API that allows to register an inner loop for user types using the descr. * PyUFunc_RegisterLoopForDescr Deprecations ============ The 'full' and 'economic' modes of qr factorization are deprecated. General ~~~~~~~ The use of non-integer for indices and most integer arguments has been deprecated. Previously float indices and function arguments such as axes or shapes were truncated to integers without warning. For example `arr.reshape(3., -1)` or `arr[0.]` will trigger a deprecation warning in NumPy 1.8., and in some future version of NumPy they will raise an error. Authors ======= This release contains work by the following people who contributed at least one patch to this release. The names are in alphabetical order by first name: * 87 * Adam Ginsburg + * Adam Griffiths + * Alexander Belopolsky + * Alex Barth + * Alex Ford + * Andreas Hilboll + * Andreas Kloeckner + * Andreas Schwab + * Andrew Horton + * argriffing + * Arink Verma + * Bago Amirbekian + * Bartosz Telenczuk + * bebert218 + * Benjamin Root + * Bill Spotz + * Bradley M. Froehle * Carwyn Pelley + * Charles Harris * Chris * Christian Brueffer + * Christoph Dann + * Christoph Gohlke * Dan Hipschman + * Daniel + * Dan Miller + * daveydave400 + * David Cournapeau * David Warde-Farley * Denis Laxalde * dmuellner + * Edward Catmur + * Egor Zindy + * endolith * Eric Firing * Eric Fode * Eric Moore + * Eric Price + * Fazlul Shahriar + * Félix Hartmann + * Fernando Perez * Frank B + * Frank Breitling + * Frederic * Gabriel * GaelVaroquaux * Guillaume Gay + * Han Genuit * HaroldMills + * hklemm + * jamestwebber + * Jason Madden + * Jay Bourque * jeromekelleher + * Jesús Gómez + * jmozmoz + * jnothman + * Johannes Schönberger + * John Benediktsson + * John Salvatier + * John Stechschulte + * Jonathan Waltman + * Joon Ro + * Jos de Kloe + * Joseph Martinot-Lagarde + * Josh Warner (Mac) + * Jostein Bø Fløystad + * Juan Luis Cano Rodríguez + * Julian Taylor + * Julien Phalip + * K.-Michael Aye + * Kumar Appaiah + * Lars Buitinck * Leon Weber + * Luis Pedro Coelho * Marcin Juszkiewicz * Mark Wiebe * Marten van Kerkwijk + * Martin Baeuml + * Martin Spacek * Martin Teichmann + * Matt Davis + * Matthew Brett * Maximilian Albert + * m-d-w + * Michael Droettboom * mwtoews + * Nathaniel J. Smith * Nicolas Scheffer + * Nils Werner + * ochoadavid + * Ondřej Čertík * ovillellas + * Paul Ivanov * Pauli Virtanen * peterjc * Ralf Gommers * Raul Cota + * Richard Hattersley + * Robert Costa + * Robert Kern * Rob Ruana + * Ronan Lamy * Sandro Tosi * Sascha Peilicke + * Sebastian Berg * Skipper Seabold * Stefan van der Walt * Steve + * Takafumi Arakaki + * Thomas Robitaille + * Tomas Tomecek + * Travis E. Oliphant * Valentin Haenel * Vladimir Rutsky + * Warren Weckesser * Yaroslav Halchenko * Yury V. Zaytsev + A total of 119 people contributed to this release. People with a "+" by their names contributed a patch for the first time.

1 0

deprecating stats.randwppf and randwcdf
by Ralf Gommers 26 Oct '13

26 Oct '13

Hi, https://github.com/scipy/scipy/pull/3007 proposes to deprecate stats.randwppf and stats.randwcdf. That makes a lot of sense, so we'll merge that PR in a few days unless someone objects. Cheers, Ralf

1 0

ANN: SciPy 0.13.0 release
by Ralf Gommers 25 Oct '13

25 Oct '13

On behalf of the SciPy development team I'm pleased to announce the availability of SciPy 0.13.0. This release contains some interesting new features (see highlights below) and half a year's worth of maintenance work. 65 people contributed to this release. Some of the highlights are: - support for fancy indexing and boolean comparisons with sparse matrices - interpolative decompositions and matrix functions in the linalg module - two new trust-region solvers for unconstrained minimization This release requires Python 2.6, 2.7 or 3.1-3.3 and NumPy 1.5.1 or greater. Support for Python 2.4 and 2.5 has been dropped as of this release. Sources and binaries can be found at http://sourceforge.net/projects/scipy/files/scipy/0.13.0/, release notes are copied below. Enjoy, Ralf ========================== SciPy 0.13.0 Release Notes ========================== .. contents:: SciPy 0.13.0 is the culmination of 7 months of hard work. It contains many new features, numerous bug-fixes, improved test coverage and better documentation. There have been a number of deprecations and API changes in this release, which are documented below. All users are encouraged to upgrade to this release, as there are a large number of bug-fixes and optimizations. Moreover, our development attention will now shift to bug-fix releases on the 0.13.x branch, and on adding new features on the master branch. This release requires Python 2.6, 2.7 or 3.1-3.3 and NumPy 1.5.1 or greater. Highlights of this release are: - support for fancy indexing and boolean comparisons with sparse matrices - interpolative decompositions and matrix functions in the linalg module - two new trust-region solvers for unconstrained minimization New features ============ ``scipy.integrate`` improvements -------------------------------- N-dimensional numerical integration ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A new function `scipy.integrate.nquad`, which provides N-dimensional integration functionality with a more flexible interface than ``dblquad`` and ``tplquad``, has been added. ``dopri*`` improvements ^^^^^^^^^^^^^^^^^^^^^^^ The intermediate results from the ``dopri`` family of ODE solvers can now be accessed by a *solout* callback function. ``scipy.linalg`` improvements ----------------------------- Interpolative decompositions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Scipy now includes a new module `scipy.linalg.interpolative` containing routines for computing interpolative matrix decompositions (ID). This feature is based on the ID software package by P.G. Martinsson, V. Rokhlin, Y. Shkolnisky, and M. Tygert, previously adapted for Python in the PymatrixId package by K.L. Ho. Polar decomposition ^^^^^^^^^^^^^^^^^^^ A new function `scipy.linalg.polar`, to compute the polar decomposition of a matrix, was added. BLAS level 3 functions ^^^^^^^^^^^^^^^^^^^^^^ The BLAS functions ``symm``, ``syrk``, ``syr2k``, ``hemm``, ``herk`` and ``her2k`` are now wrapped in `scipy.linalg`. Matrix functions ^^^^^^^^^^^^^^^^ Several matrix function algorithms have been implemented or updated following detailed descriptions in recent papers of Nick Higham and his co-authors. These include the matrix square root (``sqrtm``), the matrix logarithm (``logm``), the matrix exponential (``expm``) and its Frechet derivative (``expm_frechet``), and fractional matrix powers (``fractional_matrix_power``). ``scipy.optimize`` improvements ------------------------------- Trust-region unconstrained minimization algorithms ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The ``minimize`` function gained two trust-region solvers for unconstrained minimization: ``dogleg`` and ``trust-ncg``. ``scipy.sparse`` improvements ----------------------------- Boolean comparisons and sparse matrices ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ All sparse matrix types now support boolean data, and boolean operations. Two sparse matrices `A` and `B` can be compared in all the expected ways `A < B`, `A >= B`, `A != B`, producing similar results as dense Numpy arrays. Comparisons with dense matrices and scalars are also supported. CSR and CSC fancy indexing ^^^^^^^^^^^^^^^^^^^^^^^^^^ Compressed sparse row and column sparse matrix types now support fancy indexing with boolean matrices, slices, and lists. So where A is a (CSC or CSR) sparse matrix, you can do things like:: >>> A[A > 0.5] = 1 # since Boolean sparse matrices work >>> A[:2, :3] = 2 >>> A[[1,2], 2] = 3 ``scipy.sparse.linalg`` improvements ------------------------------------ The new function ``onenormest`` provides a lower bound of the 1-norm of a linear operator and has been implemented according to Higham and Tisseur (2000). This function is not only useful for sparse matrices, but can also be used to estimate the norm of products or powers of dense matrices without explictly building the intermediate matrix. The multiplicative action of the matrix exponential of a linear operator (``expm_multiply``) has been implemented following the description in Al-Mohy and Higham (2011). Abstract linear operators (`scipy.sparse.linalg.LinearOperator`) can now be multiplied, added to each other, and exponentiated, producing new linear operators. This enables easier construction of composite linear operations. ``scipy.spatial`` improvements ------------------------------ The vertices of a `ConvexHull` can now be accessed via the `vertices` attribute, which gives proper orientation in 2-D. ``scipy.signal`` improvements ----------------------------- The cosine window function `scipy.signal.cosine` was added. ``scipy.special`` improvements ------------------------------ New functions `scipy.special.xlogy` and `scipy.special.xlog1py` were added. These functions can simplify and speed up code that has to calculate ``x * log(y)`` and give 0 when ``x == 0``. ``scipy.io`` improvements ------------------------- Unformatted Fortran file reader ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The new class `scipy.io.FortranFile` facilitates reading unformatted sequential files written by Fortran code. ``scipy.io.wavfile`` enhancements ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ `scipy.io.wavfile.write` now accepts a file buffer. Previously it only accepted a filename. `scipy.io.wavfile.read` and `scipy.io.wavfile.write` can now handle floating point WAV files. ``scipy.interpolate`` improvements ---------------------------------- B-spline derivatives and antiderivatives ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ `scipy.interpolate.splder` and `scipy.interpolate.splantider` functions for computing B-splines that represent derivatives and antiderivatives of B-splines were added. These functions are also available in the class-based FITPACK interface as ``UnivariateSpline.derivative`` and ``UnivariateSpline.antiderivative``. ``scipy.stats`` improvements ---------------------------- Distributions now allow using keyword parameters in addition to positional parameters in all methods. The function `scipy.stats.power_divergence` has been added for the Cressie-Read power divergence statistic and goodness of fit test. Included in this family of statistics is the "G-test" (http://en.wikipedia.org/wiki/G-test). `scipy.stats.mood` now accepts multidimensional input. An option was added to `scipy.stats.wilcoxon` for continuity correction. `scipy.stats.chisquare` now has an `axis` argument. `scipy.stats.mstats.chisquare` now has `axis` and `ddof` arguments. Deprecated features =================== ``expm2`` and ``expm3`` ----------------------- The matrix exponential functions `scipy.linalg.expm2` and `scipy.linalg.expm3` are deprecated. All users should use the numerically more robust `scipy.linalg.expm` function instead. ``scipy.stats`` functions ------------------------- `scipy.stats.oneway` is deprecated; `scipy.stats.f_oneway` should be used instead. `scipy.stats.glm` is deprecated. `scipy.stats.ttest_ind` is an equivalent function; more full-featured general (and generalized) linear model implementations can be found in statsmodels. `scipy.stats.cmedian` is deprecated; ``numpy.median`` should be used instead. Backwards incompatible changes ============================== LIL matrix assignment --------------------- Assigning values to LIL matrices with two index arrays now works similarly as assigning into ndarrays:: >>> x = lil_matrix((3, 3)) >>> x[[0,1,2],[0,1,2]]=[0,1,2] >>> x.todense() matrix([[ 0., 0., 0.], [ 0., 1., 0.], [ 0., 0., 2.]]) rather than giving the result:: >>> x.todense() matrix([[ 0., 1., 2.], [ 0., 1., 2.], [ 0., 1., 2.]]) Users relying on the previous behavior will need to revisit their code. The previous behavior is obtained by ``x[numpy.ix_([0,1,2],[0,1,2])] = ...`. Deprecated ``radon`` function removed ------------------------------------- The ``misc.radon`` function, which was deprecated in scipy 0.11.0, has been removed. Users can find a more full-featured ``radon`` function in scikit-image. Removed deprecated keywords ``xa`` and ``xb`` from ``stats.distributions`` -------------------------------------------------------------------------- The keywords ``xa`` and ``xb``, which were deprecated since 0.11.0, have been removed from the distributions in ``scipy.stats``. Changes to MATLAB file readers / writers ---------------------------------------- The major change is that 1D arrays in numpy now become row vectors (shape 1, N) when saved to a MATLAB 5 format file. Previously 1D arrays saved as column vectors (N, 1). This is to harmonize the behavior of writing MATLAB 4 and 5 formats, and adapt to the defaults of numpy and MATLAB - for example ``np.atleast_2d`` returns 1D arrays as row vectors. Trying to save arrays of greater than 2 dimensions in MATLAB 4 format now raises an error instead of silently reshaping the array as 2D. ``scipy.io.loadmat('afile')`` used to look for `afile` on the Python system path (``sys.path``); now ``loadmat`` only looks in the current directory for a relative path filename. Other changes ============= Security fix: ``scipy.weave`` previously used temporary directories in an insecure manner under certain circumstances. Cython is now required to build *unreleased* versions of scipy. The C files generated from Cython sources are not included in the git repo anymore. They are however still shipped in source releases. The code base received a fairly large PEP8 cleanup. A ``tox pep8`` command has been added; new code should pass this test command. Scipy cannot be compiled with gfortran 4.1 anymore (at least on RH5), likely due to that compiler version not supporting entry constructs well. Authors ======= This release contains work by the following people (contributed at least one patch to this release, names in alphabetical order): * Jorge Cañardo Alastuey + * Tom Aldcroft + * Max Bolingbroke + * Joseph Jon Booker + * François Boulogne * Matthew Brett * Christian Brodbeck + * Per Brodtkorb + * Christian Brueffer + * Lars Buitinck * Evgeni Burovski + * Tim Cera * Lawrence Chan + * David Cournapeau * Dražen Lučanin + * Alexander J. Dunlap + * endolith * André Gaul + * Christoph Gohlke * Ralf Gommers * Alex Griffing + * Blake Griffith + * Charles Harris * Bob Helmbold + * Andreas Hilboll * Kat Huang + * Oleksandr (Sasha) Huziy + * Gert-Ludwig Ingold + * Thouis (Ray) Jones * Juan Luis Cano Rodríguez + * Robert Kern * Andreas Kloeckner + * Sytse Knypstra + * Gustav Larsson + * Denis Laxalde * Christopher Lee * Tim Leslie * Wendy Liu + * Clemens Novak + * Takuya Oshima + * Josef Perktold * Illia Polosukhin + * Przemek Porebski + * Steve Richardson + * Branden Rolston + * Skipper Seabold * Fazlul Shahriar * Leo Singer + * Rohit Sivaprasad + * Daniel B. Smith + * Julian Taylor * Louis Thibault + * Tomas Tomecek + * John Travers * Richard Tsai + * Jacob Vanderplas * Patrick Varilly * Pauli Virtanen * Stefan van der Walt * Warren Weckesser * Pedro Werneck + * Nils Werner + * Michael Wimmer + * Nathan Woods + * Tony S. Yu + A total of 65 people contributed to this release. People with a "+" by their names contributed a patch for the first time.

3 2

Bug in stats.fisher_exact
by Warren Weckesser 24 Oct '13

24 Oct '13

A problem with stats.fisher_exact was reported on stackoverflow: http://stackoverflow.com/questions/19548854/python-scipy-stats-module-value… The culprit appears to be this line: https://github.com/scipy/scipy/blame/master/scipy/stats/stats.py#L2633 I assume that `np.max(pexact, pmode)` should be `np.maximum(pexact, pmode)`. Could someone more familiar with this calculation confirm this? Warren

1 1