Christoph Sawade | 16 Apr 14:16 2014

Fwd: Binomial proportion confidence interval

Hi,

I would like to add approximate intervals [1] to the scipy project. Approximate intervals often have higher statistical power than exact intervals for binomial random variables.
Example usage can be found in the "binom R" package [2]. In analogy to [3], I see two possible ways to integrate these methods into the scipy API:

(a,b) = scipy.stats.binom.approx_interval(alpha, loc, scale, method='wilson')
(a,b) = scipy.stats.binom.wilson_interval(alpha, loc, scale), 

where (a,b) are the end-points of range that contain 100*alpha % of the rv's possible values.

What do you think?

Best, Christoph

[1] http://www-stat.wharton.upenn.edu/~tcai/paper/Binomial-StatSci.pdf
[2] http://cran.r-project.org/web/packages/binom/index.html
[3] 
http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_discrete.interval.html#scipy.stats.rv_discrete.interval




_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
Manoj Kumar | 15 Apr 15:24 2014
Picon

Mp3 support for scipy.io and some help

Hello Scipy-Devs,

I was working on an application that involves analyzing sound files. Basically, I would like to extract information like frequencies, amplitudes so that I can plot a frequency time transform. Since mp3 is the most widely used format, I would have to write my own function to transform the mp3 format into a numpy array (similar to what .wav
does). Would scipy be interested in such a thing?

However, Lars tells me there a lot of patent issues. Is there any workaround to this?
Sorry if this was off topic and any help would be greatly appreciated.

--
Regards,
Manoj Kumar,
Mech Undergrad
http://manojbits.wordpress.com
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
Benny Malengier | 15 Apr 08:32 2014
Picon

libflame integrated LAPACK

I think the announcement of libflame is of interest for scipy, http://www.cs.utexas.edu/~flame/web/.
It seems all is under a license scipy can  use:

From: Field G. Van Zee field <at> cs.utexas.edu
Date: April 07, 2014
Subject: The union of libflame and LAPACK

Sponsored by an NSF Software Infrastructure for Sustained
Innovation grant, we have been developing a new, vertically
integrated dense linear algebra software stack. At the bottom of
this software stack is the BLAS-like Library Instantiation Software
(BLIS). Above this, targeting sequential and multithreaded
architectures is libflame. At the top of the stack is Elemental
for distributed memory architectures.

libflame targets roughly the same layer as does LAPACK, and now we
have incorporated the LAPACK code base into libflame. For those
operations where libflame has the native functionality, the LAPACK
code becomes an interface. For all other operations, the netlib
implementation provides that functionality. We affectionately call
this new union "flapack", which offers the following benefits:

1) The libflame implementation of LAPACK is entirely coded in C.
No Fortran libraries or compilers are required.
2) The libflame library builds upon the BLIS interface. This
interface, unlike the BLAS, allows for arbitrary row and column
stride. While some applications may benefit from this (e.g., those
that perform computation with slices of tensors), from a
development and maintainability point of view it allows more
functionality to be supported with less code.
3) The union of the two libraries allows users to benefit from both
the LAPACK and libflame code base, within one package.
4) "flapack" passes the LAPACK test suite on platforms where we
have tested this. (There is one exception of a test case that
involves packed matrices that we believe is not in general use.)

The library is available under a 3-clause BSD license at:
https://github.com/flame/libflame

_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
Charles R Harris | 13 Apr 07:42 2014
Picon

Errors with numpy-devel

Hi All,

I get 75 errors and 3 failures when testing against current numpy on my machine. Most of the errors are due to either the deprecation of the binary '-' operator for booleans or to the deprecation of double ellipsis for indexing, i.e., '(..., ...)' . The remainder look like two numerical precision problems and one I can't immediately identify.

The main question I have is what is the best way to deal with the deprecations?


FAIL: test_lsmr.TestLSMR.testBidiagonalA
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/usr/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/home/charris/.local/lib/python2.7/site-packages/scipy/sparse/linalg/isolve/tests/test_lsmr.py", line 60, in testBidiagonalA
    self.assertCompatibleSystem(A,xtrue)
  File "/home/charris/.local/lib/python2.7/site-packages/scipy/sparse/linalg/isolve/tests/test_lsmr.py", line 40, in assertCompatibleSystem
    assert_almost_equal(norm(x - xtrue), 0, 6)
  File "/home/charris/.local/lib/python2.7/site-packages/numpy/testing/utils.py", line 486, in assert_almost_equal
    raise AssertionError(_build_err_msg())
AssertionError:
Arrays are not almost equal to 6 decimals
 ACTUAL: 6.048630163037888e-07
 DESIRED: 0

======================================================================
FAIL: test_qhull.TestUtilities.test_degenerate_barycentric_transforms
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/usr/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/home/charris/.local/lib/python2.7/site-packages/numpy/testing/decorators.py", line 146, in skipper_func
    return f(*args, **kwargs)
  File "/home/charris/.local/lib/python2.7/site-packages/scipy/spatial/tests/test_qhull.py", line 296, in test_degenerate_barycentric_transforms
    assert_(bad_count < 20, bad_count)
  File "/home/charris/.local/lib/python2.7/site-packages/numpy/testing/utils.py", line 50, in assert_
    raise AssertionError(smsg)
AssertionError: 20

======================================================================
FAIL: test_trim (test_mstats_basic.TestTrimming)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/charris/.local/lib/python2.7/site-packages/scipy/stats/tests/test_mstats_basic.py", line 270, in test_trim
    assert_equal(trimx._mask.ravel(),[1]*20+[0]*70+[1]*20)
  File "/home/charris/.local/lib/python2.7/site-packages/numpy/ma/testutils.py", line 123, in assert_equal
    return assert_array_equal(actual, desired, err_msg)
  File "/home/charris/.local/lib/python2.7/site-packages/numpy/ma/testutils.py", line 196, in assert_array_equal
    header='Arrays are not equal')
  File "/home/charris/.local/lib/python2.7/site-packages/numpy/ma/testutils.py", line 189, in assert_array_compare
    verbose=verbose, header=header)
  File "/home/charris/.local/lib/python2.7/site-packages/numpy/testing/utils.py", line 660, in assert_array_compare
    raise AssertionError(msg)
AssertionError:
Arrays are not equal

(mismatch 9.09090909091%)
 x: array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True, False, False, False, False, False, False, False,...
 y: array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...

Chuck


_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
Nathaniel Smith | 10 Apr 21:45 2014
Picon

Fwd: [Python-Dev] PEP 465: A dedicated infix operator for matrix multiplication

Hey all,

Given the sometimes rocky history of collaboration between numerical
Python and core Python, I thought it might be helpful to flag this
posting for broader distribution -- it gives one perspective on how
the core devs see things. (And is certainly consistent with my
experience around PEP 465.)

(Nick is, among other things, core Python's packaging czar, and
previously on record [1] as wishing for more feedback on how the
current energy around python packaging could take our needs into
account.)

-n

[1] https://ncoghlan_devs-python-notes.readthedocs.org/en/latest/pep_ideas/core_packaging_api.html#a-long-caveat-on-this-essay

---------- Forwarded message ----------
From: Nick Coghlan <ncoghlan <at> gmail.com>
Date: Tue, Apr 8, 2014 at 1:32 PM
Subject: Re: [Python-Dev] PEP 465: A dedicated infix operator for
matrix multiplication
To: Björn Lindqvist <bjourne <at> gmail.com>
Cc: Sturla Molden <sturla.molden <at> gmail.com>, "python-dev <at> python.org"
<python-dev <at> python.org>

On 8 April 2014 21:24, Björn Lindqvist <bjourne <at> gmail.com> wrote:
> 2014-04-08 12:23 GMT+02:00 Sturla Molden <sturla.molden <at> gmail.com>:
>> Björn Lindqvist <bjourne <at> gmail.com> wrote:
>>
>>> import numpy as np
>>> from numpy.linalg import inv, solve
>>>
>>> # Using dot function:
>>> S = np.dot((np.dot(H, beta) - r).T,
>>>            np.dot(inv(np.dot(np.dot(H, V), H.T)), np.dot(H, beta) - r))
>>>
>>> # Using dot method:
>>> S = (H.dot(beta) - r).T.dot(inv(H.dot(V).dot(H.T))).dot(H.dot(beta) - r)
>>>
>>> Don't keep your reader hanging! Tell us what the magical variables H,
>>> beta, r and V are. And why import solve when you aren't using it?
>>> Curious readers that aren't very good at matrix math, like me, should
>>> still be able to follow your logic. Even if it is just random data,
>>> it's better than nothing!
>>
>> Perhaps. But you don't need to know matrix multiplication to see that those
>> expressions are not readable. And by extension, you can still imagine that
>> bugs can easily hide in unreadable code.
>>
>> Matrix multiplications are used extensively in anything from engineering to
>> statistics to computer graphics (2D and 3D). This operator will be a good
>> thing for a lot of us.
>
> All I ask for is to be able to see that with my own eyes. Maybe there
> is some drastic improvement I can invent to make the algorithm much
> more readable? Then the PEP:s motivation falls down. I don't want to
> have to believe that the code the pep author came up with is the most
> optimal. I want to prove that for myself.

Note that the relationship between the CPython core development team
and the Numeric Python community is based heavily on trust - we don't
expect them to teach us to become numeric experts, we just expect them
to explain themselves well enough to persuade us that a core language
or interpreter change is the only realistic way to achieve a
particular goal. This does occasionally result in quirky patterns of
feature adoption, as things like extended slicing, full rich
comparison support, ellipsis support, rich buffer API support, and now
matrix multiplication support, were added for the numeric community's
benefit without necessarily offering any immediately obvious benefit
for those not using the numeric Python stack - it was only later that
they became pervasively adopted throughout the standard library (with
matmul, for example, a follow on proposal to allow tuples, lists and
arrays to handle vector dot products may make sense).

This particular problem has been kicking around long enough, and is
sufficiently familiar to several of us, that what's in the PEP already
presents a compelling rationale for the *folks that need to be
convinced* (which is primarily Guido, but if enough of the other core
devs think something is a questionable idea, we can often talk him out
of it - that's not the case here though).

Attempting to condense that preceding 10+ years of history *into the
PEP itself* wouldn't be a good return on investment - the links to the
earlier PEPs are there, as are the links to these discussion threads.

Cheers,
Nick.

P.S. We've been relatively successful in getting a similar trust based
dynamic off the ground for the packaging and distribution community
over the past year or so. The next big challenge in trust based
delegation for the core development team will likely be around a
Python 3.5+ only WSGI 2.0 (that can assume the Python 3 text model,
the restoration of binary interpolation, the availability of asyncio,
etc), but most of the likely principals there are still burned out
from the WSGI 1.1 debate and the Python 3 transition in general :(

>
>
> --
> mvh/best regards Björn Lindqvist
> _______________________________________________
> Python-Dev mailing list
> Python-Dev <at> python.org
> https://mail.python.org/mailman/listinfo/python-dev
> Unsubscribe: https://mail.python.org/mailman/options/python-dev/ncoghlan%40gmail.com

--
Nick Coghlan   |   ncoghlan <at> gmail.com   |   Brisbane, Australia
_______________________________________________
Python-Dev mailing list
Python-Dev <at> python.org
https://mail.python.org/mailman/listinfo/python-dev
Unsubscribe: https://mail.python.org/mailman/options/python-dev/njs%40pobox.com

--

-- 
Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
Warren Weckesser | 8 Apr 14:00 2014
Picon

wiki.scipy.org is down.

The wiki site, wiki.scipy.org, appears to be down.

Warren
Severin Holzer-Graf | 6 Apr 13:18 2014
Picon

sparse.block_diag improvement

Hi,

i had to generate very large sparse block diagonal matrices.
I found out that the existing sparse.block_diag is way to slow for
this and took over 30 secs to generate the matrix i needed.
The reason why the current sparse.block_diag implementation does that
very slow is because it generates everything in a row-wise manner,
which involves looping and a lot of memory.

I wrote my own method that creates the indptr and indices directly.
This works for equally shaped matrices that should be arranged in a
sparse block diagonal matrix.
A benchmark can be found here: http://de.tinypic.com/r/30ue59h/8
It's about a x600 speedup.

My question is, is there interest in integrating this in Scipy, and if
yes, how can this be done best?
I tried extending the current block_diag in a fork, but the indexing
trick only works for equally sized matrices. So either add a check in
the beginning of the existing block_diag and than executing the fast
way, or writing another function, especially for equally sized
matrices.

Best
Severin
Jaime Fernández del Río | 3 Apr 23:35 2014
Picon

Algorithm used in minimum_filter1d

Hi,

A PR of mine was merged yesterday that uses an improved algorithm for the computation of 1D minimum filters in ndimage, see here:


Looking a little more into this, I have found out that pandas' rolling_min function uses a different algorithm. They attribute it to bottleneck, which attributes it to Richard Harter:


Diving a little deeper in google scholar, the origin of the algorithm seems to be this 1996 paper by Scott C. Douglas:


In this last paper I found out that I had reinvented the MINLINE algorithm, and that the one used in pandas and bottleneck is known as MINLIST.

If the sequence has n items and the filter is of size m, MINLIST has guaranteed O(n) performance, but needs O(m) extra storage, while MINLINE only has expected O(n) performance, and O(n.m) worst case, but uses O(1) extra storage. Comparing pandas rolling_min with the new scipy's minimum_filter1d, scipy runs about 3-4x faster for random data of any length, although it is also several times slower for large windows on strictly increasing sequences.

Would implementing MINLIST be a better option than MINLINE? It will likely make the function slower than the old scipy implementation for small windows, and the extra storage is annoying. But it does have a more consistent performance over all possible inputs. In a way it is like choosing between quicksort and mergesort. Perhaps implementing both and having a user selectable switch is the way to go.

But before putting it all together, I'd like to hear some opinions on whether this is worth it at all, and what the preferred option is.

Jaime

--
(\__/)
( O.o)
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
Daniel Farrell | 2 Apr 07:46 2014
Picon

Is odeint re-entrant?

Dear list,

The docs warn that the lsoda, vode, zvode solvers in scipy.integrate.ode are not re-entrant. How about the the implementation for scipy.integrate.odeint? The docs state that is used the lsoda solver, so I guess that it is not re-entrant.

To check my understanding, re-entrant prevents me from running multiple instances of the above solvers on different processes?

Best wishes,

Dan

_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
Matt Newville | 1 Apr 20:11 2014

Re: SciPy-Dev Digest, Vol 126, Issue 1

> Date: Tue, 1 Apr 2014 13:57:01 +0200
> From: Johann cohen-tanugi <johann.cohentanugi <at> gmail.com>
> Subject: Re: [SciPy-Dev] Levenberg-Marquardt Implementation
> To: SciPy Developers List <scipy-dev <at> scipy.org>
> Message-ID:
>         <CAJ_f96c6-U3r5QopnHeP_SkOrn1k-Vkt=ht=ifD3bPrF=Us6sw <at> mail.gmail.com>
> Content-Type: text/plain; charset="iso-8859-1"
>
> FWIW I came across a python LM implementation of the minpack LM routine in
> https://code.google.com/p/astrolibpy/source/browse/mpfit/mpfit.py
> The license is GPLv3 though.
> Johann

The code (originally in IDL by Craig Markwardt at U Wisconsin and
available at http://www.physics.wisc.edu/~craigm/idl/) and the
translation into Python by Mark Rivers (with whom I work) shown in the
link you gave (and I could post separately) is actually not restricted
by GPLv3 and is under a more permissive BSD-like license:

 Translated from MINPACK-1 in FORTRAN, Apr-Jul 1998, CM
 Copyright (C) 1997-2002, Craig Markwardt
 This software is provided as is without any warranty whatsoever.
 Permission to use, copy, modify, and distribute modified or
 unmodified copies is granted, provided this copyright and disclaimer
 are included unchanged.

   Translated from MPFIT (Craig Markwardt's IDL package) to Python,
   August, 2002.  Mark Rivers

Astrolib.py may be distributing this under the GPLv3, but this code is
more permissively available.

Whether one wants to use this version (as opposed to the version in
scipy.optimize.leastsq) is another matter -- mpfit.py is in pure
python, and considerably slower than leastsq.   It might be worth
considering wrapping Markwardt's C version of MPFIT (which came after
the 2002 translation by Mark Rivers), but I have not looked into doing
this.

MPFIT does handle constraints to some extent (box constraints are
handled  tied values between parameters rare handled in a klunky way).
  Lmfit-py (which we're now using in place of MPFIT), uses a different
bounds implementations (that as described by MINUIT which makes
getting uncertainties for values near the boundaries more reliable)
and has a more flexible means for constraining parameter values to one
another.

I don't know of any implementation that handles "sparse Jacobians"
well.... I'm not sure I know of a problem that clearly fits that
category, though I can believe they exist.

--Matt Newville
Brian Lee Newsom | 27 Mar 17:02 2014
Picon

Re: Multivariate ctypes integration PR

Hi,

I believe I have corrected most of the issues brought up by Pauli with the multivariate code.  I have not merged the python callback machinery as of now, but that can still be done.

I ran my c code through indent using the suggested command to fix formatting issues, removed printfs, fixed comments, added some documentation, and check the return value of c_array_from_tuple.

I added a few test cases showing the accuracy of solutions and increase in speed.  Unfortunately the unusual function signature means that I cannot use built-in c library functions, so I created the test_multivariate.c file that must be compiled to run the tests.  I know this can be built into the process, but I was unsure of how to do that, so I currently only run these tests on a linux system.

Please let me know what else I can do to get this merged as soon as possible.

Thanks,
Brian
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev

Gmane