Mailman 3 June 2013 - NumPy-Discussion

low level optimization in NumPy and minivect
by Frédéric Bastien 27 Jun '13

27 Jun '13

Hi, I saw that recently Julian Taylor is doing many low level optimization like using SSE instruction. I think it is great. Last year, Mark Florisson released the minivect[1] project that he worked on during is master thesis. minivect is a compiler for element-wise expression that do some of the same low level optimization that Julian is doing in NumPy right now. Mark did minivect in a way that allow it to be reused by other project. It is used now by Cython and Numba I think. I had plan to reuse it in Theano, but I didn't got the time to integrate it up to now. What about reusing it in NumPy? I think that some of Julian optimization aren't in minivect (I didn't check to confirm). But from I heard, minivect don't implement reduction and there is a pull request to optimize this in NumPy. The advantage to concentrate some functionality in one common package is that more project benefit from optimization done to it. (after the work to use it first!) How this could be done in NumPy? NumPy have its own code generator for many dtype. We could call minivect code generator to replace some of it. What do you think of this idea? Sadly, I won't be able to spent time on the code for this, but I wanted to raise the idea while people are working on that, in case it is helpful. Frédéric [1] https://github.com/markflorisson88/minivect

7 13

slicing with boolean in numpy master
by josef.pktd＠gmail.com 26 Jun '13

26 Jun '13

Is there a change in the behavior of boolean slicing in current master? If not I have to find another candidate in numpy master. (py27d) E:\Josef\testing\tox\py27d\Scripts>python Python 2.7.1 (r271:86832, Nov 27 2010, 18:30:46) [MSC v.1500 32 bit (Intel)] on win32 Type "help", "copyright", "credits" or "license" for more information. >>> import numpy as np >>> np.__version__ '1.7.1' >>> x = np.ones((5,3)) >>> mask = np.arange(5) < 4 >>> x[mask, :] array([[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]]) We get errors like the following when running the statsmodels tests with a current or recent numpy master, but not with numpy 1.7.1 ====================================================================== ERROR: Failure: IndexError (too many indices) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Users/tom/python2.7/lib/python2.7/site-packages/nose/loader.py", line 518, in makeTest return self._makeTest(obj, parent) File "/Users/tom/python2.7/lib/python2.7/site-packages/nose/loader.py", line 577, in _makeTest return MethodTestCase(obj) File "/Users/tom/python2.7/lib/python2.7/site-packages/nose/case.py", line 345, in __init__ self.inst = self.cls() File "/Users/tom/python2.7/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-macosx-10.8-x86_64.egg/statsmodels/emplike/tests/test_aft.py", line 19, in __init__ super(Test_AFTModel, self).__init__() File "/Users/tom/python2.7/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-macosx-10.8-x86_64.egg/statsmodels/emplike/tests/test_aft.py", line 12, in __init__ self.mod1 = sm.emplike.emplikeAFT(endog, exog, data.censors) File "/Users/tom/python2.7/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-macosx-10.8-x86_64.egg/statsmodels/emplike/aft_el.py", line 248, in __init__ self.uncens_endog = self.endog[np.bool_(self.censors), :].\ IndexError: too many indices Thanks, Josef

2 5

PEP 208 and upstreaming numpy
by Todd 26 Jun '13

26 Jun '13

I see that the plan to merge old Numeric into the python standard library, PEP 208, is listed as withdrawn, although no reasons are given as far as I could see. Considering how mature Numpy has become, and how common it is as a dependency for python packages, I was wondering if there were still plans on the table to ultimately merge at least the ndarray class into the python standard library, and if so where those plans stand.

2 3

[ANN] Using some MATLAB optimization solvers from Python (OpenOpt/FuncDesigner)
by Dmitrey 25 Jun '13

25 Jun '13

Hi all, > > FYI some MATLAB solvers now can be involved with OpenOpt or FuncDesigner : > * LP linprog * QP quadprog * LLSP lsqlin * MILP bintprog > > Sparsity handling is supported. > > You should have * MATLAB (or MATLAB Component Runtime) * mlabwrap > Unfortunately, it will hardly work out-of-the-box, you have to adjust some paths and some environment variables. > > As for nonlinear solvers, e.g. fmincon, probably they could be connected via involving C MEX files, but it is not possible with current state of mlabwrap yet. > > Read MATLAB entry for details. > > Regards, D. >

1 0

Proposed enhancement to "Returns" section for numpy doc strings
by Rob Ruana 25 Jun '13

25 Jun '13

I am an admirer of the numpy doc string style. I'd like to propose a simple enhancement to the Returns section that should make the common case more readable, while maintaining backward compatibility with existing documentation. The current numpy documentation spec [1] requires that return values are named. Most of the time I only care about the return *type* of a function, and I can't be bothered to come up with a name. Of course, I can use the default "out" return name, but that adds nothing of value to the documentation, and takes away from the readability. I'd really like to just write: Returns ------- bool True if some condition, False otherwise. Of course there are times when names can *enhance* the documentation, for example numpy.meshgrid: Returns ------- X1, X2,..., XN : ndarray For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` , return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij' or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy' with the elements of `xi` repeated to fill the matrix along the first dimension for `x1`, the second for `x2` and so on. I suggest that the *name* of each return value should be optional, while the *type* is required. This would allow for the common case, while still supporting the documentation as it stands. I've implemented this change in the numpydoc_returns_unnamed feature branch on my github fork. The numpydoc sphinx extension has been updated, as well as the reStructuredText example.py and the guide to numpy documentation. Please review my pull request (https://github.com/numpy/numpy/pull/3468) and let me know what you think. [1]: https://github.com/numpy/numpy/blob/7af95b0f779f587dfe56c44c5091731b122ceb1…

1 0

[ANN] Multidimensional Array - MDArray 0.5.3
by Rodrigo Botafogo 24 Jun '13

24 Jun '13

*Announcement* MDArray version 0.5.3 has Just been released. MDArray is a multi dimensional array implemented for JRuby inspired by NumPy (www.numpy.org) and Masahiro Tanaka´s Narray (narray.rubyforge.org). MDArray stands on the shoulders of Java-NetCDF and Parallel Colt. At this point MDArray has libraries for mathematical, trigonometric and descriptive statistics methods. NetCDF-Java Library is a Java interface to NetCDF files<http://www.unidata.ucar.edu/software/netcdf/index.html>, as well as to many other types of scientific data formats. It is developed and distributed by Unidata (http://www.unidata.ucar.edu). Parallel Colt ( http://grepcode.com/snapshot/repo1.maven.org/maven2/net.sourceforge.paralle…) is a multithreaded<http://en.wikipedia.org/wiki/Thread_%28computer_science%29> version of Colt <http://dsd.lbl.gov/~hoschek/colt/> ( http://acs.lbl.gov/software/colt/). Colt provides a set of Open Source Libraries for High Performance Scientific and Technical Computing in Java. Scientific and technical computing is characterized by demanding problem sizes and a need for high performance at reasonably small memory footprint. *What´s new**:* *Performance Improvement* On previous versions, array operations were done by passing a Ruby Proc to a loop for all elements of the given arrays. For instance, adding two MDArrays was done by passing Proc.new { |a, b| a + b } and looping through all elements of the arrays. Procs are very flexible in Ruby; however, from my experience with MDArray, also very slow. On this version, when available, instead of passing a Proc to the loop, we pass a native Java method. Available Java methods are those extracted from Parallel Colt and listed below. Note that Parallel Colt has native methods for the following types only: “double”, “float”, “long” and “int”. With this change, there was a performance improvement of over 90%, and using MDArray operations is close to native Java operations. We expect (but have not yet benchmarking data) that this brings MDArray performance close to similar solutions such as NArray, NMatrix and NumPy (please try it, and if this assertion is false, I´ll be glad to change it in future announcements). Methods not available in Parallel Colt but supported by Ruby, such as “sinh”, “tanh”, and “add” for “byte” type, etc. are still supported by MDArray. Again, to improve performance, instead of passing a Proc we now create a class as follows class Add * *def self.apply(a, b) a + b end end This change brought performance improvement of over 60% for MDArray operations with Ruby methods. *Experimental Lazy Evaluation* Usual MDArray operations are done eagerly, i.e., if @a, @b, @c are three MDArrays then the following: @d = @a + @b + @c will be evaluated as follows: first @a + @b is performed and stored in a temporary variable, then this temporary variable is added to @c. For large expressions, temporary variables can have significant performance impact. This version of MDArray introduces lazy evaluation of expressions. Thus, when in lazy mode: @lazy_d = @a + @b + @c will not evaluate immediately. Rather, the expression is preprocessed and only executed when required. Since at execution time the whole expression is known, there is no need for temporary variables as the whole expression is executed at once. To put MDArray in lazy mode we only need to set its mode to lazy with the following command “MDArray.lazy = true”. All expressions after that are by default lazy. In lazy mode, MDArray resembles Numexpr, however, there is no need to write the expression as a string, and there is no compilation involved. MDArray does not implement broadcasting rules as NumPy. As a result, trying to operate on arrays of different shape raises an exception. On lazy mode, this exception is raise only at evaluation time, so it is possible to have an invalid lazy array. To evaluate a lazy array one should use the “[]” method as follows: @d = lazy_d[] @d is now a normal MDArray. Lazy MDArrays are really lazy, so let’s assume that @a = [1, 2, 3, 4] and @b = [5, 6, 7, 8]. Lets also have @l_c = @a + @b. Now doing @c = @l_c[], will evaluate @c to [6, 8, 10, 12]. Now, let´s do @a[1] = 20 and then @d = @l_c[]. Now @d evaluates to [25, 8, 10, 12] as the new value of @a is used. Lazy arrays can be evaluated inside expressions: @l_c = (@a + @b)[] + @c In this example, @l_c is a lazy array, but (@a + @b) is evaluated when the “[]” method is called and then added to @c. If now the value of @a or @b is changed, the evaluation of @l_c will not be changed as in the previous example. Finally, laziness is contagious. So, let´s assume that we have @l_c as above, a lazy array and we do MDArray.lazy = false. From this point on in the code, operations will be done eagerly. Now doing: @e = @d + @l_c, @e is a lazy array as its construction involves a lazy array. One should be careful when in eager mode mixing lazy and eager arrays: @c = @l_a + (@b + @c) then, with parenthesis, first (@b + @c) is evaluated eagerly and then added lazily to @l_a, giving a lazy array. In this version, Lazy evaluation is around 40% *less* efficient in one machine I tested up to approximately the same performance in another equipment than eager evaluation when only native Java methods (Parallel Colt methods described below) are used in the expression. If expression involves any Ruby method, evaluation of lazy expressions becomes much slower than eager evaluation. In order to improve performance, I believe that compilation of expression will be necessary. *MDArray and SciRuby**:* MDArray subscribes fully to the SciRuby Manifesto (http://sciruby.com/). *“**Ruby* <http://www.ruby-lang.org/>* has for some time had no equivalent to the beautifully constructed **NumPy* <http://numpy.scipy.org/>*, SciPy<http://www.scipy.org/>, and **matplotlib* <http://matplotlib.sourceforge.net/>* libraries for ** Pytho <http://www.python.org/>n.* *We believe that the time for a Ruby science and visualization package has come. Sometimes when a solution of sugar and water becomes super-saturated, from it precipitates a pure, delicious, and diabetes-inducing crystal of sweetness, induced by no more than the tap of a finger. So is occurring now, we believe, with numeric and visualization libraries for Ruby.”* *MDArray main properties are**:* · Homogeneous multidimensional array, a table of elements (usually numbers), all of the same type, indexed by a tuple of positive integers; · Easy calculation for large numerical multi dimensional arrays; · Basic types are: boolean, byte, short, int, long, float, double, string, structure; · Based on JRuby, which allows importing Java libraries; · Operator: +,-,*,/,%,**, >, >=, etc.; · Functions: abs, ceil, floor, truncate, is_zero, square, cube, fourth; · Binary Operators: &, |, ^, ~ (binary_ones_complement), <<, >>; · Ruby Math functions: acos, acosh, asin, asinh, atan, atan2, atanh, cbrt, cos, erf, exp, gamma, hypot, ldexp, log, log10, log2, sin, sinh, sqrt, tan, tanh, neg; · Boolean operations on boolean arrays: and, or, not; · Fast descriptive statistics from Parallel Colt (complete list found bellow); · Easy manipulation of arrays: reshape, reduce dimension, permute, section, slice, etc.; · Reading of two dimensional arrays from CSV files (mainly for debugging and simple testing purposes); · StatList: a list that can grow/shrink and that can compute Parallel Colt descriptive statistics; · Experimental lazy evaluation (still slower than eager evaluation). *Descriptive statistics methods imported from Parallel Colt**:* auto_correlation, correlation, covariance, durbin_watson, frequencies, geometric_mean, harmonic_mean, kurtosis, lag1, max, mean, mean_deviation, median, min, moment, moment3, moment4, pooled_mean, pooled_variance, product, quantile, quantile_inverse, rank_interpolated, rms, sample_covariance, sample_kurtosis, sample_kurtosis_standard_error, sample_skew, sample_skew_standard_error, sample_standard_deviation, sample_variance, sample_weighted_variance, skew, split, standard_deviation, standard_error, sum, sum_of_inversions, sum_of_logarithms, sum_of_powers, sum_of_power_deviations, sum_of_squares, sum_of_squared_deviations, trimmed_mean, variance, weighted_mean, weighted_rms, weighted_sums, winsorized_mean. *Double and Float methods from Parallel Colt*: acos, asin, atan, atan2, ceil, cos, exp, floor, greater, IEEEremainder, inv, less, lg, log, log2, rint, sin, sqrt, tan. *Double, Float, Long and Int methods from Parallel Colt*: abs, compare, div, divNeg, equals, isEqual (is_equal), isGreater (is_greater), isles (is_less), max, min, minus, mod, mult, multNeg (mult_neg), multSquare (mult_square), neg, plus (add), plusAbs (plus_abs), pow (power), sign, square. *Long and Int methods from Parallel Colt* and, dec, factorial, inc, not, or, shiftLeft (shift_left), shiftRightSigned (shift_right_signed), shiftRightUnsigned (shift_right_unsigned), xor. *MDArray installation and download**:* · Install Jruby · jruby –S gem install mdarray *MDArray Homepages**:* · http://rubygems.org/gems/mdarray · https://github.com/rbotafogo/mdarray/wiki *Contributors**:* Contributors are welcome. *MDArray History**:* · 24/05/2013: Version 0.5.0 – Over 90% Performance improvements for methods imported from Parallel Colt and over 40% performance improvements for all other methods (implemented in Ruby); · 16/05/2013: Version 0.5.0 - All loops transferred to Java with over 50% performance improvements. Descriptive statistics from Parallel Colt; · 19/04/2013: Version 0.4.3 - Fixes a simple, but fatal bug in 0.4.2. No new features; · 17/04/2013: Version 0.4.2 - Adds simple statistics and boolean operators; · 05/04/2013: Version 0.4.0 – Initial release. -- Rodrigo Botafogo Integrando TI ao seu negócio

1 0

Computing least squares fit with masked array
by Thomas Goebel 24 Jun '13

24 Jun '13

Hi list, i want to calculate the least squares fit of some data with missing values. So i masked all values using numpys masked array module. Unfortunately using linalg.lstsq i only get nan data back. Can somebody help me solve this problem? The output of np.linalg.lstsq ist: (masked_array(data = [ nan nan nan nan], mask = False, fill_value = 1e+20) , masked_array(data = [ nan], mask = False, fill_value = 1e+20) , 4, array([ 2.08406074e+03, 1.73066556e+01, 1.32393450e+01, 1.09943429e-01])) Thanks for your help, Tom The python code i use is: #!/usr/bin/env python2 # -*- coding: utf-8 -*- import numpy as np x = np.array([-30, -25, -20]) y = np.array([25, 20, 35]) z = np.array([ [1285, 1613, 2190], [np.nan, 1283, 1897], [np.nan, np.nan, 1619] ]) N, M = len(x), len(y) x_new = np.hstack([x] * N) y_new = np.repeat(y, M) z_new = np.repeat(z, 1) z_new = np.ma.masked_array(z_new, np.isnan(z_new)) A = np.array([x_new, y_new, x_new*y_new, np.ones(N*M)]) out = np.linalg.lstsq(A.T, z_new) print out

1 0

ANN: XDress v0.2
by Anthony Scopatz 24 Jun '13

24 Jun '13

Hello All, I am pleased to announce the latest version of xdress, in preparation for SciPy 2013. For more information please visit the website: http://xdress.org Be Well Anthony ======================== XDress 0.2 Release Notes ======================== XDress is an automatic wrapper generator for C/C++ written in pure Python. Currently, xdress may generate Python bindings (via Cython) for C++ classes, functions, and certain variable types. It also contains idiomatic wrappers for C++ standard library containers (sets, vectors, maps). In the future, other tools and bindings will be supported. The main enabling feature of xdress is a dynamic type system that was designed with the purpose of API generation in mind. Release highlights: - First class support for C via pycparser. - Python 3 support, by popular demand! - A plethora of awesome type system updates, including: - type matching - lambda-valued converters - function pointer support - Easy to use and implement plugin architecture. Please visit the website for more information: http://xdress.org/ Ask questions on the mailing list: xdress(a)googlegroups.com Download the code from GitHub: http://github.com/scopatz/xdress XDress is free & open source (BSD 2-clause license) and requires Python 2.7, NumPy 1.5+, Cython 0.19+, and optionally GCC-XML, pycparser, and lxml. New Features ============ First Class C Support --------------------- Wrapping C code is now fully handled through with the optional pycparser. This means that you don't have to worry about whether or not the GCC-XML parser will work on your particular bit of code. Furthermore, C structs and their members are now handled idiomatically. (C++ structs are actually treated as C++ classes, which means that their are allowed to have constructors and other C++ concepts not present in C.) Python 3 Support ---------------- The entire code base is built and tested under Python 3.3. This uses the single code base model, as most development takes pace in Python 2.7. The Cython code that is generated may be used by both Python 2 & 3. Type System ----------- The type system has been expanded and hardened to handle additional use cases, largely motivated by the desire for realistic C support. A new function pointer (``fucntion_pointer``) refinement type has been added. When converting from C to Python, a new function object is created that wraps the underlying call. For converting from Python to C, a virtual table of callback functions is constructed that have the same signature of the pointer but hold a reference to a Python function. The size of the table and thus how many callbacks you can have before overwriting previous ones, is set by the ``max_callbacks`` key in the extra dictionary in class descriptions. This defaults to 8. A new enum refinement type now also comes stock. These may be exposed to Python in the ``rc.variables`` list. The type system now also comes with basic type matching tools. There is a new ``TypeMatcher`` class, a ``matches()`` function, and a singleton ``MatchAny`` that may be used for determining whether a type adheres to a pattern. The TypeMatcher class itself is immutable and hashable and therefore may be used anywhere other type elements (tuples, str, int) may be used including as dict keys! This is helpful for specifying conversions for large groups of types. Finally, the type system conversion dictionaries now accept callable objects as values. This was put in to handle templated types where the number of argument is not in general know beforehand, e.g. enums and function pointers. The values must be callable with only a single argument -- the type itself. For example, ``lambda t: rtnvalue`` is valid. Plugins! -------- XDress is a suite of tools written on top of a type system. Thus the entire core has been refactored to implement a very nimble and easy plugin system. The plugin mechanism enables external projects with their own code generators to easily hook into the xdress tools. Additionally, this allows users to decide at run time which plugins they want to use. Mini-FAQ ======== * Why not use an existing solution (eg, SWIG)? Their type systems don't support run-time, user provided refinement types, and thus are unsuited for verification & validation use cases that often arise in computational science. Furthermore, they tend to not handle C++ dependent types well (i.e. vector<T> does not come back as a np.view(..., dtype=T)). * Why GCC-XML and not Clang's AST? I tried using Clang's AST (and the remnants of a broken visitor class remain in the code base). However, the official Clang AST Python bindings lack support for template argument types. This is a really big deal. Other C++ ASTs may be supported in the future -- including Clang's. * I run xdress and it creates these files, now what?! It is your job to integrate the files created by xdress into your build system. Join in the Fun! ================ If you are interested in using xdress on your project (and need help), contributing back to xdress, starting up a development team, or writing your own code generation plugin tool on top of the type system and autodescriber, please let us know. Participation is very welcome! Authors ======= - `Anthony Scopatz <http://scopatz.com/>`_ - Spencer Lyon * - Gerald Dalley * - Alexander Eisenhuth * An * indicates a first time contributor. Links ===== 1. Homepage - http://xdress.org/ 2. Mailing List - xdress(a)googlegroups.com 3. GitHub Organization - https://github.com/xdress 4. Pycparser - https://pypi.python.org/pypi/pycparser

1 0

request for SWIG numpy.i users
by Ralf Gommers 20 Jun '13

20 Jun '13

Hi, If you're using or are very familiar with SWIG and the numpy.i interface to it, please help to test and/or review https://github.com/numpy/numpy/pull/3148. It's a fairly major update to numpy.i by Bill Spotz, containing the following: - support for 4D arrays and memory managed output arguments - rework for the deprecated API's in numpy 1.6 and 1.7 - a bug fix in a 3D typemap - documentation improvements It would be good to have this merged before branching 1.8.x. Not many of the regular reviewers of numpy PRs are familiar with numpy.i, therefore help would be much appreciated. Thanks, Ralf

3 6

Profiling (was GSoC : Performance parity between numpy arrays and Python scalars)
by David Cournapeau 20 Jun '13

20 Jun '13

> * Re: the profiling, I wrote a full oprofile->callgrind format script > years ago: http://vorpus.org/~njs/op2calltree.py > Haven't used it in years either but neither oprofile nor kcachegrind > are terribly fast-moving projects so it's probably still working, or > could be made so without much work. > Or easier is to use the gperftools CPU profiler: > https://gperftools.googlecode.com/svn/trunk/doc/cpuprofile.html I don't have experience with gperftools, but on recent linux kernels, you can also use perf, which can't be made easier to use (no runtime support needed), but you need a 'recent' kernel http://indico.cern.ch/getFile.py/access?contribId=20&sessionId=4&resId=0&ma… I am hoping to talk a bit about those for our diving into numpy c code tutorial in June, what's the + of gperf in your opinion ? David

5 17