Mailman 3 July 2014 - NumPy-Discussion

numpy.mean still broken for large float32 arrays
by Thomas Unterthiner 25 Jul '14

25 Jul '14

Hi! The following is a known "bug" since at least 2010 [1]: import numpy as np X = np.ones((50000, 1024), np.float32) print X.mean() >>> 0.32768 I ran into this for the first time today as part of a larger program. I was very surprised by this, and spent over an hour looking for bugs in my code before noticing that the culprit was `mean` being broken for large float32 arrays. I realize that this behavior is actually documented, but it is absolutely non-intuitive. I assume most users expect `mean` to just work. This has been discussed once two years ago [2], but nothing came of that. This could be easily fixed by making `np.float64` the default dtype (as it already is for integer types), or by at least checking inside mean if the passed array was a large np.float32 array and switch the dtype to np.float64 in that case. Is there a reason why this has not been done? Cheers Thomas [1] http://mail.scipy.org/pipermail/numpy-discussion/2010-November/053697.html [2] http://numpy-discussion.10968.n7.nabble.com/Bug-in-numpy-mean-revisited-td1…

9 9

Re: [Numpy-discussion] numpy.mean still broken for large float32arrays
by Eelco Hoogendoorn 25 Jul '14

25 Jul '14

True, i suppose there is no harm in accumulating with max precision, and storing the result in the Original dtype, unless otherwise specified, although i wonder if the current nditer supports such behavior. -----Original Message----- From: "Alan G Isaac" <alan.isaac(a)gmail.com> Sent: ‎24-‎7-‎2014 18:09 To: "Discussion of Numerical Python" <numpy-discussion(a)scipy.org> Subject: Re: [Numpy-discussion] numpy.mean still broken for large float32arrays On 7/24/2014 5:59 AM, Eelco Hoogendoorn wrote to Thomas: > np.mean isn't broken; your understanding of floating point number is. This comment seems to conflate separate issues: the desirable return type, and the computational algorithm. It is certainly possible to compute a mean of float32 doing reduction in float64 and still return a float32. There is nothing implicit in the name `mean` that says we have to just add everything up and divide by the count. My own view is that `mean` would behave enough better if computed as a running mean to justify the speed loss. Naturally similar issues arise for `var` and `std`, etc. See http://www.johndcook.com/standard_deviation.html for some discussion and references. Alan Isaac _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion(a)scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

1 0

Re: [Numpy-discussion] numpy.mean still broken for large float32arrays
by Eelco Hoogendoorn 25 Jul '14

25 Jul '14

Thanks Julian, those seem like Nice improvements. The fact that it either does or doesnt work depending on the axis makes me a Little queesy; but yeah, checking that fp's do what You think they should, is unfortunately best left as the responsibility of the programmer. -----Original Message----- From: "Julian Taylor" <jtaylor.debian(a)googlemail.com> Sent: ‎24-‎7-‎2014 14:56 To: "Discussion of Numerical Python" <numpy-discussion(a)scipy.org> Subject: Re: [Numpy-discussion] numpy.mean still broken for large float32arrays On Thu, Jul 24, 2014 at 1:33 PM, Fabien <fabien.maussion(a)gmail.com> wrote: > Hi all, > > On 24.07.2014 11:59, Eelco Hoogendoorn wrote: >> np.mean isn't broken; your understanding of floating point number is. > > I am quite new to python, and this problem is discussed over and over > for other languages too. However, numpy's summation problem appears with > relatively small arrays already: > > py>import numpy as np > py>np.ones((4000,4000), np.float32).mean() > 1.0 > py>np.ones((5000,5000), np.float32).mean() > 0.67108864000000001 > > A 5000*5000 image is not unusual anymore today. > > In IDL: > IDL> mean(fltarr(5000L, 5000L)+1) > 1.0000000 > IDL> mean(fltarr(7000L, 7000L)+1) > 1.0000000 > IDL> mean(fltarr(10000L, 10000L)+1) > 0.67108864 > > I can't really explain why there are differences between the two > languages (IDL uses 32-bit, single-precision, floating-point numbers) > > Fabien > something as simple as summation is already an interesting algorithmic problem there are several ways do to with different speeds and accuracies. E.g. pythons math.fsum is always exact to one ulp but is very slow as it requires partial sorting. Then there is kahan summation that has an accuracy of O(1) ulp but its about 4 times slower than the naive sum. In practice one of the better methods is pairwise summation that is pretty much as fast as a naive summation but has an accuracy of O(logN) ulp. This is the method numpy 1.9 will use this method by default (+ its even a bit faster than our old implementation of the naive sum): https://github.com/numpy/numpy/pull/3685 but it has some limitations, it is limited to blocks fo the buffer size (8192 elements by default) and does not work along the slow axes due to limitations in the numpy iterator. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion(a)scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

1 0

change default integer from int32 to int64 on win64?
by Julian Taylor 24 Jul '14

24 Jul '14

hi, it recently came to my attention that the default integer type in numpy on windows 64 bit is a 32 bit integers [0]. This seems like a quite serious problem as it means you can't use any integers created from python integers < 32 bit to index arrays larger than 2GB. For example np.product(array.shape) which will never overflow on linux and mac, can overflow on win64. I think this is a very dangerous platform difference and a quite large inconvenience for win64 users so I think it would be good to fix this. This would be a very large change of API and probably also ABI. But as we also never officially released win64 binaries we could change it for from source compilations and give win64 binary distributors the option to keep the old ABI/API at their discretion. Any thoughts on this from win64 users? Cheers, Julian Taylor [0] https://github.com/astropy/astropy/pull/2697

5 11

__numpy_ufunc__ and 1.9 release
by Julian Taylor 24 Jul '14

24 Jul '14

hi, as you may know we want to release numpy 1.9 soon. We should have solved most indexing regressions the first beta showed. The remaining blockers are finishing the new __numpy_ufunc__ feature. This feature should allow for alternative method to overriding the behavior of ufuncs from subclasses. It is described here: https://github.com/numpy/numpy/blob/master/doc/neps/ufunc-overrides.rst The current blocker issues are: https://github.com/numpy/numpy/issues/4753 https://github.com/numpy/numpy/pull/4815 I'm not to familiar with all the complications of subclassing so I can't really say how hard this is to solve. My issue is that it there still seems to be debate on how to handle operator overriding correctly and I am opposed to releasing a numpy with yet another experimental feature that may or may not be finished sometime later. Having datetime in infinite experimental state is bad enough. I think nobody is served well if we release 1.9 with the feature prematurely based on a not representative set of users and the later after more users showed up see we have to change its behavior. So I'm wondering if we should delay the introduction of this feature to 1.10 or is it important enough to wait until there is a consensus on the remaining issues?

4 5

ANN: bcolz 0.7.0 released
by Francesc Alted 23 Jul '14

23 Jul '14

====================== Announcing bcolz 0.7.0 ====================== What's new ========== In this release, support for Python 3 has been added, Pandas and HDF5/PyTables conversion, support for different compressors via latest release of Blosc, and a new `iterblocks()` iterator. Also, intensive benchmarking has lead to an important tuning of buffer sizes parameters so that compression and evaluation goes faster than ever. Together, bcolz and the Blosc compressor, are finally fullfilling the promise of accelerating memory I/O, at least for some real scenarios: http://nbviewer.ipython.org/github/Blosc/movielens-bench/blob/master/queryi… ``bcolz`` is a renaming of the ``carray`` project. The new goals for the project are to create simple, yet flexible compressed containers, that can live either on-disk or in-memory, and with some high-performance iterators (like `iter()`, `where()`) for querying them. For more detailed info, see the release notes in: https://github.com/Blosc/bcolz/wiki/Release-Notes What it is ========== bcolz provides columnar and compressed data containers. Column storage allows for efficiently querying tables with a large number of columns. It also allows for cheap addition and removal of column. In addition, bcolz objects are compressed by default for reducing memory/disk I/O needs. The compression process is carried out internally by Blosc, a high-performance compressor that is optimized for binary data. bcolz can use numexpr internally so as to accelerate many vector and query operations (although it can use pure NumPy for doing so too). numexpr optimizes the memory usage and use several cores for doing the computations, so it is blazing fast. Moreover, the carray/ctable containers can be disk-based, and it is possible to use them for seamlessly performing out-of-memory computations. bcolz has minimal dependencies (NumPy), comes with an exhaustive test suite and fully supports both 32-bit and 64-bit platforms. Also, it is typically tested on both UNIX and Windows operating systems. Installing ========== bcolz is in the PyPI repository, so installing it is easy: $ pip install -U bcolz Resources ========= Visit the main bcolz site repository at: http://github.com/Blosc/bcolz Manual: http://bcolz.blosc.org Home of Blosc compressor: http://blosc.org User's mail list: bcolz(a)googlegroups.com http://groups.google.com/group/bcolz License is the new BSD: https://github.com/Blosc/bcolz/blob/master/LICENSES/BCOLZ.txt ---- **Enjoy data!** -- Francesc Alted

2 1

Research position in the Brazilian Research Institute for Science and Neurotechnology – BRAINN
by Thiago Franco Moraes 23 Jul '14

23 Jul '14

*Research position in the Brazilian Research Institute for Science and Neurotechnology – BRAINN Postdoc researcher to work with software development for medical imaging* The Brazilian Research Institute for Neuroscience and Neurotechnology (BRAINN) (www.brainn.org.br) focuses on the investigation of basic mechanisms leading to epilepsy and stroke, and the injury mechanisms that follow disease onset and progression. This research has important applications related to prevention, diagnosis, treatment and rehabilitation and will serve as a model for better understanding normal and abnormal brain function. The BRAINN Institute is composed of 10 institutions from Brazil and abroad and hosted by State University of Campinas (UNICAMP). Among the associated institutions is Renato Archer Information Technology Center (CTI) that has a specialized team in open-source software development for medical imaging (www.cti.gov.br/invesalius) and 3D printing applications for healthcare. CTI is located close the UNICAMP in the city of Campinas, State of São Paulo in a very technological region of Brazil and is looking for a postdoc researcher to work with software development for medical imaging related to the imaging analysis, diagnosis and treatment of brain diseases. The postdoc position is for two years with the possibility of being renovated for more two years. *Education* - PhD in computer science, computer engineering, mathematics, physics or related. *Requirements* - Digital image processing (Medical imaging) - Computer graphics (basic) * Benefits* 6.143,40 Reais per month free of taxes (about US$ 2.800,00); 15% technical reserve for conferences participation and specific materials acquisition; *Interested* Send curriculum to: jorge.silva(a)cti.gov.br with subject “Postdoc position” Applications reviews will begin August 1, 2014 and continue until the position is filled.

1 0

BLAS / LAPACK / MKL cannot be found?
by Bram Willemsen 20 Jul '14

20 Jul '14

Hi everyone, I am trying to install a package called PySparse. I have modified the paths in the template site.cfg file. It appears that during the built process it compiles numpy as well, because the error messages I get are all over the numpy mailing list (I did not find one addressing my problem exactly). My MKL libs are installed at "/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64" . There are BLAS / LAPACK libraries here (I provided the library mkl_rt.so to SuiteSparse for instance, and it resolved all LAPACK and BLAS dependencies that way). But for some reason "/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64" does not work in the built process below, in what I think is a numpy installer. The output shows that the directory is searched, but that the result is "NOT AVAILABLE". Could someone give me a pointer for why my the MKL/BLAS/LAPACK dependencies cannot be resolved? It would be very appreciated sincerely, Bram --------------------------------------------------------------------------- PART OF THE OUTPUT OF PySparse SETUP.PY --------------------------------------------------------------------------- blas_opt_info: blas_mkl_info: libraries mkl,vml,guide not found in ['/usr/local/lib', '/wgdisk/hy3300/re15/lwillemsen/local_install/lib', '/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64'] NOT AVAILABLE openblas_info: libraries not found in ['/usr/local/lib', '/wgdisk/hy3300/re15/lwillemsen/local_install/lib', '/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64'] NOT AVAILABLE atlas_blas_threads_info: Setting PTATLAS=ATLAS libraries ptf77blas,ptcblas,atlas not found in ['/usr/local/lib', '/wgdisk/hy3300/re15/lwillemsen/local_install/lib', '/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64'] NOT AVAILABLE atlas_blas_info: libraries f77blas,cblas,atlas not found in ['/usr/local/lib', '/wgdisk/hy3300/re15/lwillemsen/local_install/lib', '/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64'] NOT AVAILABLE blas_info: libraries blas not found in ['/usr/local/lib', '/wgdisk/hy3300/re15/lwillemsen/local_install/lib', '/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64'] NOT AVAILABLE blas_src_info: NOT AVAILABLE NOT AVAILABLE No blas info found Sparse:: Using BLAS info: {} Using dflt_lib_dirs = /usr/local/lib:/wgdisk/hy3300/re15/lwillemsen/local_install/lib:/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64 Using dflt_libs = [] No blas info found Eigen:: Using BLAS info: {} lapack_opt_info: lapack_mkl_info: mkl_info: libraries mkl,vml,guide not found in ['/usr/local/lib', '/wgdisk/hy3300/re15/lwillemsen/local_install/lib', '/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64'] NOT AVAILABLE NOT AVAILABLE atlas_threads_info: Setting PTATLAS=ATLAS libraries ptf77blas,ptcblas,atlas not found in /usr/local/lib libraries lapack_atlas not found in /usr/local/lib libraries ptf77blas,ptcblas,atlas not found in /wgdisk/hy3300/re15/lwillemsen/local_install/lib libraries lapack_atlas not found in /wgdisk/hy3300/re15/lwillemsen/local_install/lib libraries ptf77blas,ptcblas,atlas not found in /wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64 libraries lapack_atlas not found in /wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64 numpy.distutils.system_info.atlas_threads_info NOT AVAILABLE atlas_info: libraries f77blas,cblas,atlas not found in /usr/local/lib libraries lapack_atlas not found in /usr/local/lib libraries f77blas,cblas,atlas not found in /wgdisk/hy3300/re15/lwillemsen/local_install/lib libraries lapack_atlas not found in /wgdisk/hy3300/re15/lwillemsen/local_install/lib libraries f77blas,cblas,atlas not found in /wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64 libraries lapack_atlas not found in /wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64 numpy.distutils.system_info.atlas_info NOT AVAILABLE lapack_info: libraries lapack not found in ['/usr/local/lib', '/wgdisk/hy3300/re15/lwillemsen/local_install/lib', '/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64'] NOT AVAILABLE lapack_src_info: NOT AVAILABLE NOT AVAILABLE No lapack info found Eigen:: Using LAPACK info: {} non-existing path in 'pysparse/eigen': '/usr/local/lib:/wgdisk/hy3300/re15/lwillemsen/local_install/lib:/wgdisk/omega2dev2/env/EL5/intel/composer_xe_2013.0.079/mkl/lib/intel64' No blas info found Direct:: Using BLAS info: {}

1 1

proposal: new commit guidelines for backportable bugfixes
by Julian Taylor 19 Jul '14

19 Jul '14

hi, I have been doing a lot of backporting for the last few bugfix releases and noticed that our current approach committing to master and cherrypicking is not so good for the git history. When cherry picking a bugfix from master to a maintenance branch both branches contain a commit with the same content and git knows of no relation between them. This causes unnecessary merge conflicts when cherry picking two changes that modify the same file. The git version (1.9.1) I am using is not smart enough too figure out the changesets in both leaf commits is the same. Additionally the output of `git log maintenance/1.9.x..master` becomes very large as all already backported issues appear again in master. [0] To help with this I want to propose new best practices for pull requests of bugfixes suitable for backporting. Instead of basing the bugfix on the head commit of the master, base them on the merge base between master and the latest maintenance branch. This allows merging the PR into both master and the maintenance branch without pulling in any extra changes from either branches. Then both branches contain the same commit and gits automerging can work better and git log will only show you the commits that are only really on one branch or the other. In practice this is very simple. You can still develop your bugfix on master but before you push it you just run: git rebase --onto $(git merge-base master maintenance/1.9.x) HEAD^ In most bugfix PRs this should work without conflict as they should be relatively small. If you get a merge conflict during this operation, just do git rebase --abort and do a normal pull request, in that case the backporter should worry about the conflict. Does this sound like a reasonable procedure? Cheers, Julian [0] git cherry is supposed to help with that, but it never really worked properly for me

6 17

Re: [Numpy-discussion] `allclose` vs `assert_allclose`
by alex 19 Jul '14

19 Jul '14

On Fri, Jul 18, 2014 at 9:47 PM, Charles R Harris <charlesr.harris(a)gmail.com> wrote: > > > > On Fri, Jul 18, 2014 at 2:32 PM, Chris Barker <chris.barker(a)noaa.gov> wrote: >> >> On Fri, Jul 18, 2014 at 12:43 PM, Pauli Virtanen <pav(a)iki.fi> wrote: >>> >>> 18.07.2014 22:13, Chris Barker kirjoitti: >>> [clip] >>> > but an appropriate rtol would work there too. If only zero testing is >>> > needed, then atol=0 makes sense as a default. (or maybe atol=eps) >>> >>> There's plenty of room below eps, but finfo(float).tiny ~ 3e-308 (or >>> some big multiple) is also reasonable in the scale-freeness sense. >> >> >> right! brain blip -- eps is the difference between 1 and then next larger >> representable number, yes? So a long way away from smallest representable >> number. So yes, zero or [something]e-308 -- making zero seem like a good >> idea again.... >> >> is it totally ridiculous to have the default be dependent on dtype? >> float32 vs float64? >> > > Whatever the final decision is, if the defaults change we should start with > a FutureWarning. How we can make that work is uncertain, because I don't > know of any reliable way to detect if we are using the default value or if a > value was passed in. There are tricks like http://stackoverflow.com/questions/12265695, not that I'm suggesting to do that.

1 0

NumPy-Discussion July 2014