Re: [Numpy-discussion] Silent Broadcasting considered harmful
I personally use Octave and/or Numpy for several years now and never ever needed braodcasting. But since it is still there there will be many users who need it, there will be some use for it. Uhm, yeah, there is some use for it. Im all for explicit over implicit, but personally current broadcasting rules have never bothered me, certainly not to the extent of justifying massive backwards compatibility violations. Take It from someone who relies on broadcasting for every other line of code. On Sun, Feb 8, 2015 at 10:03 PM, Stefan Reiterer dom...@gmx.net wrote: Hi! As shortly discussed on github: https://github.com/numpy/numpy/issues/5541 I personally think that silent Broadcasting is not a good thing. I had recently a lot of trouble with row and column vectors which got bradcastet toghether altough it was more annoying than useful, especially since I had to search deep down into the code to find out that the problem was nothing else than Broadcasting... I personally use Octave and/or Numpy for several years now and never ever needed braodcasting. But since it is still there there will be many users who need it, there will be some use for it. So I suggest that the best would be to throw warnings when arrays get Broadcasted like Octave do. Python warnings can be catched and handled, that would be a great benefit. Another idea would to provide warning levels for braodcasting, e.g 0 = Never, 1=Warn once, 2=Warn always, 3 = Forbid aka throw exception, with 0 as default. This would avoid breaking other code, and give the user some control over braodcasting. Kind regards, Stefan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Optimizing numpy's einsum expression (again)
Thanks for taking the time to think about this; good work. Personally, I don't think a factor 5 memory overhead is much to sweat over. The most complex einsum I have ever needed in a production environment was 5/6 terms, and for what this anecdote is worth, speed was a far bigger concern to me than memory. On Fri, Jan 16, 2015 at 12:30 AM, Daniel Smith dgasm...@icloud.com wrote: Hello everyone, I originally brought an optimized einsum routine forward a few weeks back that attempts to contract numpy arrays together in an optimal way. This can greatly reduce the scaling and overall cost of the einsum expression for the cost of a few intermediate arrays. The current version (and more details) can be found here: https://github.com/dgasmith/opt_einsum I think this routine is close to a finalized version, but there are two problems which I would like the community to weigh in on: Memory usage- Currently the routine only considers the maximum size of intermediates created rather than cumulative memory usage. If we only use np.einsum it is straightforward to calculate cumulative memory usage as einsum does not require any copies of inputs to be made; however, if we attempt to use a vendor BLAS the memory usage can becomes quite complex. For example, the np.tensordot routine always forces a copy for ndarrays because it uses ndarray.transpose(…).reshape(…). A more memory-conscious way to do this is to first try and do ndarray.reshape(…).T which does not copy the data and numpy can just pass a transpose flag to the vendor BLAS. The caveat here is that the summed indices must be in the correct order- if not a copy is required. Maximum array size is usually close to the total overhead of the opt_einsum function, but can occasionally be 2-5 times this size. I see the following ways forward: - Ignore cumulative memory and stick with maximum array size as the limiting memory factor. - Implement logic to figure out if the input arrays needs to be copied to use BLAS, compute the extra memory required, and add an extra dimension to the optimization algorithms (use BLAS or do not use BLAS at each step). Some of this is already done, but may get fairly complex. - Build an in-place nd-transpose algorithm. - Cut out BLAS entirely. Keeping in mind that vendor BLAS can be orders of magnitude faster than a pure einsum implementation, especially if the BLAS threading is used. Path algorithms- There are two algorithms “optimal” (a brute force algorithm, scales like N!) and “opportunistic” (a greedy algorithm, scales like N^3). The optimal path can take seconds to minutes to calculate for a 7-10 term expression while the opportunistic path takes microseconds even for 20+ term expressions. The opportunistic algorithm works surprisingly well and appears to obtain the correct scaling in all test cases that I can think of. Both algorithms use the maximum array size as a sieve- this is very beneficial from several aspects. The problem occurs when a much needed intermediate cannot be created due to this limit- on occasions not making this intermediate can have slowdowns of orders of magnitude even for small systems. This leads to certain (and sometimes unexpected) performance walls. Possible solutions: - Warn the user if the ratio between an unlimited memory solution and a limited memory solution becomes large. - Do not worry about it. Thank you for your time, -Daniel Smith ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Sorting refactor
I don't know if there is a general consensus or guideline on these matters, but I am personally not entirely charmed by the use of behind-the-scenes parallelism, unless explicitly requested. Perhaps an algorithm can be made faster, but often these multicore algorithms are also less efficient, and a less data-dependent way of putting my cores to good use would have been preferable. Potentially, other code could slow down due to cache trashing if too many parallel tasks run in parallel. Id rather be in charge of such matters myself; but I imagine adding a keyword arg for these matters would not be much trouble? On Fri, Jan 16, 2015 at 12:43 PM, Julian Taylor jtaylor.deb...@googlemail.com wrote: On 16.01.2015 12:33, Lars Buitinck wrote: 2015-01-16 11:55 GMT+01:00 numpy-discussion-requ...@scipy.org: Message: 2 Date: Thu, 15 Jan 2015 21:24:00 -0800 From: Jaime Fern?ndez del R?o jaime.f...@gmail.com Subject: [Numpy-discussion] Sorting refactor To: Discussion of Numerical Python numpy-discussion@scipy.org Message-ID: capowhwkf6rnwcrgmcwsmq_lo3hshjgbvlsrn19z-mdpe25e...@mail.gmail.com Content-Type: text/plain; charset=utf-8 This changes will make it easier for me to add a Timsort generic type function to numpy's arsenal of sorting routines. And I think they have value by cleaning the source code on their own. Yes, they do. I've been looking at the sorting functions as well and I've found the following: * The code is generally hard to read because it prefers pointer over indices. I'm wondering if it would get slower using indices. The closer these algorithms are to the textbook, the easier to insert fancy optimizations. * The heap sort exploits undefined behavior by using a pointer that points before the start of the array. However, rewriting it to always point within the array made it slower. I haven't tried rewriting it using indices. * Quicksort has a quadratic time worst case. I think it should be turned into an introsort [1] for O(n log n) worst case; we have the heapsort needed to do that. This probably rarely happens in numeric data, and we do have guaranteed nlog runtime algorithms available. But it also is not costly to do, e.g. the selection code is a introselect instead of a normal quickselect. I'd say not high priority, but if someone wants to do it I don't see why not. * Quicksort is robust to repeated elements, but doesn't exploit them. It can be made to run in linear time if the input array has only O(1) distinct elements [2]. This may come at the expense of some performance on arrays with no repeated elements. * Using optimal sorting networks instead of insertion sort as the base case can speed up quicksort on float arrays by 5-10%, but only if NaNs are moved out of the way first so that comparisons become cheaper [3]. This has consequences for the selection algorithms that I haven't fully worked out yet. I was also thinking about this, an advantage of a sorting network is that it can be vectorized to be significantly faster than an insertion sort. Handling NaN's should also be directly possible. The issue is that its probably too much complicated code for only a very small gain. * Using Cilk Plus to parallelize merge sort can make it significantly faster than quicksort at the expense of only a few lines of code, but I haven't checked whether Cilk Plus plays nicely with multiprocessing and other parallelism options (remember the trouble with OpenMP-ified OpenBLAS?). you should also be able to do this with openmp tasks, though it will be a little less efficient as cilk+ has a better scheduler for this type of work. But I assume you will get the same trouble as openmp but that needs testing, also cilk+ in gcc is not really production ready yet, I got lots of crashes when I last tried it (it will be in 5.0 though). This isn't really an answer to your questions, more like a brain dump from someone who's stared at the same code for a while and did some experiments. I'm not saying we should implement all of this, but keep in mind that there are some interesting options besides implementing timsort. [1] https://en.wikipedia.org/wiki/Introsort [2] http://www.sorting-algorithms.com/static/QuicksortIsOptimal.pdf [3] https://github.com/larsmans/numpy/tree/sorting-nets ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Sorting refactor
I agree; an np.setnumthreads to manage a numpy-global threadpool makes
sense to me.
Of course there are a great many cases where just spawning as many threads
as cores is a sensible default, but if this kind of behavior could not be
overridden I could see that greatly reduce performance for some of my more
complex projects
On Fri, Jan 16, 2015 at 4:11 PM, Julian Taylor
jtaylor.deb...@googlemail.com wrote:
On 01/16/2015 03:14 PM, Lars Buitinck wrote:
2015-01-16 13:29 GMT+01:00 numpy-discussion-requ...@scipy.org:
Date: Fri, 16 Jan 2015 12:43:43 +0100
From: Julian Taylor jtaylor.deb...@googlemail.com
Subject: Re: [Numpy-discussion] Sorting refactor
To: Discussion of Numerical Python numpy-discussion@scipy.org
Message-ID: 54b8f96f.7090...@googlemail.com
Content-Type: text/plain; charset=windows-1252
On 16.01.2015 12:33, Lars Buitinck wrote:
* Quicksort has a quadratic time worst case. I think it should be
turned into an introsort [1] for O(n log n) worst case; we have the
heapsort needed to do that.
This probably rarely happens in numeric data, and we do have guaranteed
nlog runtime algorithms available.
It's no more likely or unlikely than in any other type of data
(AFAIK), but it's possible for an adversary to DOS any (web server)
code that uses np.sort.
if you are using numpy where an arbitrary user is allowed to control the
data passed to a non isolated environment you have a problem anyway.
numpy is far from secure software and there are likely hundreds of ways
to produce DOS and dozens of ways to do code execution in any nontrivial
numpy using application.
I was also thinking about this, an advantage of a sorting network is
that it can be vectorized to be significantly faster than an insertion
sort. Handling NaN's should also be directly possible.
Tried that, and it didn't give any speedup, at least not without
explicit vector instructions.
Just moving the NaNs aside didn't cost anything in my preliminary
benchmarks (without sorting nets), the cost of the operation was
almost exactly compensated by simpler comparisons.
an SSE2 implementation a 16 entry bitonic sort is available here:
https://github.com/mischasan/sse2/blob/master/ssesort.c
there is also a benchmark, on my machine its 6 times faster than
insertion sort.
But again this would only gain us 5-10% improvement at best as the
partition part of quicksort is still the major time consuming part.
The issue is that its probably too much complicated code for only a very
small gain.
Maybe. The thing is that the selection algorithms are optimized for
NaNs and seem to depend on the current comparison code. We'd need
distinct TYPE_LT and TYPE_LT_NONAN for each TYPE.
The sorting nets themselves aren't complicated, just lengthy. My
branch has the length-optimal (not optimally parallel) ones for n =
16.
But I assume you will get the same trouble as openmp but that needs
testing, also cilk+ in gcc is not really production ready yet, I got
lots of crashes when I last tried it (it will be in 5.0 though).
The data parallel constructs tend to crash the compiler, but task
spawning seems to be stable in 4.9.2. I've still to see how it handles
multiprocessing/fork.
What do you mean by will be in 5.0, did they do a big push?
gcc 5.0 changelog reports full support for cilk plus.
Also all bugs I have filed have been fixed in 5.0.
Date: Fri, 16 Jan 2015 13:28:56 +0100
From: Da?id davidmen...@gmail.com
Subject: Re: [Numpy-discussion] Sorting refactor
To: Discussion of Numerical Python numpy-discussion@scipy.org
Message-ID:
CAJhcF=1O5Own_5ydzu+To8HHbm3e66k=
iunqreiasdy23dn...@mail.gmail.com
Content-Type: text/plain; charset=UTF-8
On 16 January 2015 at 13:15, Eelco Hoogendoorn
hoogendoorn.ee...@gmail.com wrote:
Perhaps an algorithm can be made faster, but often these multicore
algorithms are also less efficient, and a less data-dependent way of
putting
my cores to good use would have been preferable. Potentially, other
code
could slow down due to cache trashing if too many parallel tasks run in
parallel. Id rather be in charge of such matters myself; but I imagine
adding a keyword arg for these matters would not be much trouble?
As I understand it, that is where the strength of Cilk+ lies. It does
not force parallelisation, just suggests it. The decision to actually
spawn parallel is decided at runtime depending on the load of the
other cores.
cilk+ guarantees that the amount of space used by a pool of P threads
is at most P times the stack space used by the sequential version (+ a
constant). The idea is that you can say
for (i = 0; i 100; i++) {
cilk_spawn f(a[i]);
}
and it will never create more than P work items in memory, rather than
1e6, even if each f() spawns a bunch itself. Of course, it won't
guarantee that OpenMP will not also spawn P threads
Re: [Numpy-discussion] Question about dtype
This is a general problem in trying to use JSON to send arbitrary python
objects. Its not made for that purpose, JSON itself only supports a very
limited grammar (only one sequence type for instance, as you noticed), so
in general you will need to specify your own encoding/decoding for more
complex objects you want to send over JSON.
In the case of an object dtype, dtypestr = str(dtype) gives you a nice
JSONable string representation, which you can convert back into a dtype
using np.dtype(eval(dtypestr))
On Sat, Dec 13, 2014 at 9:25 AM, Sebastian se...@sebix.at wrote:
Hi,
I'll just comment on the creation of your dtype:
dt = [(f8, f8)]
You are creating a dtype with one field called 'f8' and with type 'f8':
dt = [(f8, f8)]
dty = np.dtype(dt)
dty.names
('f8',)
What you may want are two fields with type 'f8' and without fieldname:
dt = [(f8, f8)]
dty = np.dtype(('f8,f8'))
dty.names
('f0', 'f1')
dty.descr
[('f0', 'f8'), ('f1', 'f8')]
I can't help you with the json-module and what it's doing there. As the
output is unequal to the input, I suspect JSON to be misbehaving here.
If you need to store the dtype as strings, not as binary pickle, you can
use pickle.dumps and pickle.loads
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Should ndarray be a context manager?
My impression is that this level of optimization does and should not fall within the scope of numpy.. -Original Message- From: Sturla Molden sturla.mol...@gmail.com Sent: 9-12-2014 16:02 To: numpy-discussion@scipy.org numpy-discussion@scipy.org Subject: [Numpy-discussion] Should ndarray be a context manager? I wonder if ndarray should be a context manager so we can write something like this: with np.zeros(n) as x: [...] The difference should be that __exit__ should free the memory in x (if owned by x) and make x a zero size array. Unlike the current ndarray, which does not have an __exit__ method, this would give precise control over when the memory is freed. The timing of the memory release would not be dependent on the Python implementation, and a reference cycle or reference leak would not accidentally produce a memory leak. It would allow us to deterministically decide when the memory should be freed, which e.g. is useful when we work with large arrays. A problem with this is that the memory in the ndarray would be volatile with respect to other Python threads and view arrays. However, there are dozens of other ways to produce segfaults or buffer overflows with NumPy (cf. stride_tricks or wrapping external buffers). Below is a Cython class that does something similar, but we would need to e.g. write something like with Heapmem(n * np.double().itemsize) as hm: x = hm.doublearray [...] instead of just with np.zeros(n) as x: [...] Sturla # (C) 2014 Sturla Molden from cpython cimport PyMem_Malloc, PyMem_Free from libc.string cimport memset cimport numpy as cnp cnp.init_array() cdef class Heapmem: cdef: void *_pointer cnp.intp_t _size def __cinit__(Heapmem self, Py_ssize_t n): self._pointer = NULL self._size = cnp.intp_t n def __init__(Heapmem self, Py_ssize_t n): self.allocate() def allocate(Heapmem self): if self._pointer != NULL: raise RuntimeError(Memory already allocated) else: self._pointer = PyMem_Malloc(self._size) if (self._pointer == NULL): raise MemoryError() memset(self._pointer, 0, self._size) def __dealloc__(Heapmem self): if self._pointer != NULL: PyMem_Free(self._pointer) self._pointer = NULL property pointer: def __get__(Heapmem self): return cnp.intp_t self._pointer property doublearray: def __get__(Heapmem self): cdef cnp.intp_t n = self._size//sizeof(double) if self._pointer != NULL: return cnp.PyArray_SimpleNewFromData(1, n, cnp.NPY_DOUBLE, self._pointer) else: raise RuntimeError(Memory not allocated) property chararray: def __get__(Heapmem self): if self._pointer != NULL: return cnp.PyArray_SimpleNewFromData(1, self._size, cnp.NPY_CHAR, self._pointer) else: raise RuntimeError(Memory not allocated) def __enter__(self): if self._pointer != NULL: raise RuntimeError(Memory not allocated) def __exit__(Heapmem self, type, value, traceback): if self._pointer != NULL: PyMem_Free(self._pointer) self._pointer = NULL ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] FFTS for numpy's FFTs (was: Re: Choosing between NumPy and SciPy functions)
My point isn't about speed; its about the scope of numpy. typing np.fft.fft isn't more or less convenient than using some other symbol from the scientific python stack. Numerical algorithms should be part of the stack, for sure; but should they be part of numpy? I think its cleaner to have them in a separate package. Id rather have us discuss how to facilitate the integration of as many possible fft libraries with numpy behind a maximally uniform interface, rather than having us debate which fft library is 'best'. On Tue, Oct 28, 2014 at 6:21 PM, Sturla Molden sturla.mol...@gmail.com wrote: Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Perhaps the 'batteries included' philosophy made sense in the early days of numpy; but given that there are several fft libraries with their own pros and cons, and that most numpy projects will use none of them at all, why should numpy bundle any of them? Because sometimes we just need to compute a DFT, just like we sometimes need to compute a sine or an exponential. It does that job perfectly well. It is not always about speed. Just typing np.fft.fft(x) is convinient. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] help using np.einsum for stacked matrix multiplication
You need to specify your input format. Also, if your output matrix misses the NY dimension, that implies you wish to contract (sum) over it, which contradicts your statement that the 2x2 subblocks form the matrices to multiply with. In general, I think it would help if you give a little more background on your problem. On Wed, Oct 29, 2014 at 10:39 AM, Andrew Nelson andyf...@gmail.com wrote: Dear list, I have a 4D array, A, that has the shape (NX, NY, 2, 2). I wish to perform matrix multiplication of the 'NY' 2x2 matrices, resulting in the matrix B. B would have shape (NX, 2, 2). I believe that np.einsum would be up to the task, but I'm not quite sure of the subscripts I would need to achieve this. Can anyone help, please? cheers, Andrew. -- _ Dr. Andrew Nelson _ ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] FFTS for numpy's FFTs (was: Re: Choosing between NumPy and SciPy functions)
If I may 'hyjack' the discussion back to the meta-point: should we be having this discussion on the numpy mailing list at all? Perhaps the 'batteries included' philosophy made sense in the early days of numpy; but given that there are several fft libraries with their own pros and cons, and that most numpy projects will use none of them at all, why should numpy bundle any of them? To have a scipy.linalg and scipy.fft makes sense to me, although import pyfftw or import pyFFTPACK would arguably be better still. Just as in the case of linear algebra, those different libraries represent meaningful differences, and if the user wants to paper over those differences with a named import they are always free to do so themselves, explicitly. To be sure, the maintenance of quality fft libraries should be part of the numpy/scipy-stack in some way or another. But I would argue that the core thing that numpy should do is ndarrays alone. On Tue, Oct 28, 2014 at 11:11 AM, Sturla Molden sturla.mol...@gmail.com wrote: David Cournapeau courn...@gmail.com wrote: The real issue with fftw (besides the license) is the need for plan computation, which are expensive (but are not needed for each transform). This is not a problem if you thell FFTW to guess a plan instead of making measurements. FFTPACK needs to set up a look-up table too. I made some experiments with the Bluestein transform to handle prime transforms on fftpack, but the precision seemed to be an issue. Maybe I should revive this work (if I still have it somewhere). You have it in a branch on Github. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Choosing between NumPy and SciPy functions
The same occurred to me when reading that question. My personal opinion is that such functionality should be deprecated from numpy. I don't know who said this, but it really stuck with me: but the power of numpy is first and foremost in it being a fantastic interface, not in being a library. There is nothing more annoying than every project having its own array type. The fact that the whole scientific python stack can so seamlessly communicate is where all good things begin. In my opinion, that is what numpy should focus on; basic data structures, and tools for manipulating them. Linear algebra is way too high level for numpy imo, and used by only a small subsets of its 'matlab-like' users. When I get serious about linear algebra or ffts or what have you, id rather import an extra module that wraps a specific library. On Mon, Oct 27, 2014 at 2:26 PM, D. Michael McFarland dm...@dmmcf.net wrote: A recent post raised a question about differences in results obtained with numpy.linalg.eigh() and scipy.linalg.eigh(), documented at http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html#numpy.linalg.eigh and http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigh.html#scipy.linalg.eigh , respectively. It is clear that these functions address different mathematical problems (among other things, the SciPy routine can solve the generalized as well as standard eigenproblems); I am not concerned here with numerical differences in the results for problems both should be able to solve (the author of the original post received useful replies in that thread). What I would like to ask about is the situation this illustrates, where both NumPy and SciPy provide similar functionality (sometimes identical, to judge by the documentation). Is there some guidance on which is to be preferred? I could argue that using only NumPy when possible avoids unnecessary dependence on SciPy in some code, or that using SciPy consistently makes for a single interface and so is less error prone. Is there a rule of thumb for cases where SciPy names shadow NumPy names? I've used Python for a long time, but have only recently returned to doing serious numerical work with it. The tools are very much improved, but sometimes, like now, I feel I'm missing the obvious. I would appreciate pointers to any relevant documentation, or just a summary of conventional wisdom on the topic. Regards, Michael ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Memory efficient alternative for np.loadtxt and np.genfromtxt
Im not sure why the memory doubling is necessary. Isnt it possible to preallocate the arrays and write to them? I suppose this might be inefficient though, in case you end up reading only a small subset of rows out of a mostly corrupt file? But that seems to be a rather uncommon corner case. Either way, id say a doubling of memory use is fair game for numpy. Generality is more important than absolute performance. The most important thing is that temporary python datastructures are avoided. That shouldn't be too hard to accomplish, and would realize most of the performance and memory gains, I imagine. On Sun, Oct 26, 2014 at 12:54 PM, Jeff Reback jeffreb...@gmail.com wrote: you should have a read here/ http://wesmckinney.com/blog/?p=543 going below the 2x memory usage on read in is non trivial and costly in terms of performance On Oct 26, 2014, at 4:46 AM, Saullo Castro saullogiov...@gmail.com wrote: I would like to start working on a memory efficient alternative for np.loadtxt and np.genfromtxt that uses arrays instead of lists to store the data while the file iterator is exhausted. The motivation came from this SO question: http://stackoverflow.com/q/26569852/832621 where for huge arrays the current NumPy ASCII readers are really slow and require ~6 times more memory. This case I tested with Pandas' read_csv() and it required 2 times more memory. I would be glad if you could share your experience on this matter. Greetings, Saullo ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Request for enhancement to numpy.random.shuffle
Thanks Warren, I think these are sensible additions. I would argue to treat the None-False condition as an error. Indeed I agree one might argue the correcr behavior is to 'shuffle' the singleton block of data, which does nothing; but its more likely to come up as an unintended error than as a natural outcome of parametrized behavior. On Sun, Oct 12, 2014 at 3:31 AM, John Zwinck jzwi...@gmail.com wrote: On Sun, Oct 12, 2014 at 6:51 AM, Warren Weckesser warren.weckes...@gmail.com wrote: I created an issue on github for an enhancement to numpy.random.shuffle: https://github.com/numpy/numpy/issues/5173 I like this idea. I was a bit surprised there wasn't something like this already. A small wart in this API is the meaning of shuffle(a, independent=False, axis=None) It could be argued that the correct behavior is to leave the array unchanged. (The current behavior can be interpreted as shuffling a 1-d sequence of monolithic blobs; the axis argument specifies which axis of the array corresponds to the sequence index. Then `axis=None` means the argument is a single monolithic blob, so there is nothing to shuffle.) Or an error could be raised. Let's think about it from the other direction: if a user wants to shuffle all the elements as if it were 1-d, as you point out they could do this: shuffle(a, axis=None, independent=True) But that's a lot of typing. Maybe we should just let this do the same thing: shuffle(a, axis=None) That seems to be in keeping with the other APIs taking axis as you mentioned. To me, independent has no relevance when the array is 1-d, it can simply be ignored. John Zwinck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Request for enhancement to numpy.random.shuffle
yeah, a shuffle function that does not shuffle indeed seems like a major source of bugs to me. Indeed one could argue that setting axis=None should suffice to give a clear enough declaration of intent; though I wouldn't mind typing the extra bit to ensure consistent semantics. On Sun, Oct 12, 2014 at 10:56 AM, Stefan van der Walt ste...@sun.ac.za wrote: Hi Warren On 2014-10-12 00:51:56, Warren Weckesser warren.weckes...@gmail.com wrote: A small wart in this API is the meaning of shuffle(a, independent=False, axis=None) It could be argued that the correct behavior is to leave the array unchanged. I like the suggested changes. Since independent loses its meaning when axis is None, I would expect this to have the same effect as `shuffle(a, independent=True, axis=None)`. I think a shuffle function that doesn't shuffle will confuse a lot of people! Stéfan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] 0/0 == 0?
slightly OT; but fwiw, its all ill-thought out nonsense from the start anyway. ALL numbers satisfy the predicate 0*x=0. what the IEEE calls 'not a number' would be more accurately called 'not a specific number', or 'a number'. whats a logical negation among computer scientists? On Fri, Oct 3, 2014 at 6:13 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Thu, Oct 2, 2014 at 10:12 PM, Charles R Harris charlesr.har...@gmail.com wrote: On Thu, Oct 2, 2014 at 9:29 PM, Nathaniel Smith n...@pobox.com wrote: On Fri, Oct 3, 2014 at 3:20 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Thu, Oct 2, 2014 at 7:06 PM, Benjamin Root ben.r...@ou.edu wrote: Out[1] has an integer divided by an integer, and you can't represent nan as an integer. Perhaps something weird was happening with type promotion between versions? Also note that in python3 the '/' operator does float rather than integer division. np.array(0) / np.array(0) __main__:1: RuntimeWarning: invalid value encountered in true_divide nan Floor division still acts the same though: np.array(0) // np.array(0) __main__:1: RuntimeWarning: divide by zero encountered in floor_divide 0 The seterr warning system makes a lot of sense for IEEE754 floats, which are specifically designed so that 0/0 has a unique well-defined answer. For ints though this seems really broken to me. 0 / 0 = 0 is just the wrong answer. It would be nice if we had something reasonable to return, but we don't, and I'd rather raise an error than return the wrong answer. That's an option, although arguable for arrays of numbers. However, the fact that we don't know *which* numbers caused the problem strengthens the argument for an error. Plus the g*dawful warning default to only warn once. That has always bothered me, it just seems useless. Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Proposal: add ndarray.keys() to return dtype.names
Well, the method will have to be present on all ndarrays, since structured arrays do not have a different type from regular arrays, only a different dtype. Thus the attribute has to be present regardless, but some Exception will have to be raised depending on the dtype, to make it quack like the kind of duck it is, so to speak. Indeed it seems like an atypical design pattern; but I don't see a problem with it. On Wed, Oct 1, 2014 at 4:08 PM, John Zwinck jzwi...@gmail.com wrote: On 1 Oct 2014 04:30, Stephan Hoyer sho...@gmail.com wrote: On Tue, Sep 30, 2014 at 1:22 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On more careful reading of your words, I think we agree; indeed, if keys() is present is should return an iterable; but I don't think it should be present for non-structured arrays. Indeed, I think we do agree. The attribute can simply be missing (e.g., accessing it raises AttributeError) for non-structured arrays. I'm generally fine with this, though I would like to know if there is precedent for methods being present on structured arrays only. Even if there is no precedent I am still OK with the idea, I just think we should understand if how novel this will be. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Proposal: add ndarray.keys() to return dtype.names
Sounds fair to me. Indeed the ducktyping argument makes sense, and I have a hard time imagining any namespace conflicts or other confusion. Should this attribute return none for non-structured arrays, or simply be undefined? On Tue, Sep 30, 2014 at 12:49 PM, John Zwinck jzwi...@gmail.com wrote: I first proposed this on GitHub: https://github.com/numpy/numpy/issues/5134 ; jaimefrio requested that I bring it to this list for discussion. My proposal is to add a keys() method to NumPy's array class ndarray. The behavior would be to return self.dtype.names, i.e. the column names for a structured array (and None when dtype.names is None, which it is for pure numeric arrays without named columns). I originally proposed to add a values() method also, but I am tabling that for now so we needn't discuss it in this thread. The motivation is to enhance the ability to use duck typing with NumPy arrays, Python dicts, and other types like Pandas DataFrames, h5py Files, and more. It's a fairly common thing to want to get the keys of a container, where keys is understood to be a sequence of values one can pass to __getitem__(), and this is exactly what I'm aiming at. Thoughts? John Zwinck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Proposal: add ndarray.keys() to return dtype.names
So a non-structured array should return an empty list/iterable as its keys? That doesn't seem right to me, but perhaps you have a compelling example to the contrary. I mean, wouldn't we want the duck-typing to fail if it isn't a structured array? Throwing an attributeError seems like the best thing to do, from a duck-typing perspective. On Tue, Sep 30, 2014 at 8:05 PM, Stephan Hoyer sho...@gmail.com wrote: I like this idea. But I am -1 on returning None if the array is unstructured. I expect .keys(), if present, to always return an iterable. In fact, this would break some of my existing code, which checks for the existence of keys as a way to do duck typed checks for dictionary like objects (e.g., including pandas.DataFrame): https://github.com/xray/xray/blob/v0.3/xray/core/utils.py#L165 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Proposal: add ndarray.keys() to return dtype.names
On more careful reading of your words, I think we agree; indeed, if keys() is present is should return an iterable; but I don't think it should be present for non-structured arrays. On Tue, Sep 30, 2014 at 10:21 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: So a non-structured array should return an empty list/iterable as its keys? That doesn't seem right to me, but perhaps you have a compelling example to the contrary. I mean, wouldn't we want the duck-typing to fail if it isn't a structured array? Throwing an attributeError seems like the best thing to do, from a duck-typing perspective. On Tue, Sep 30, 2014 at 8:05 PM, Stephan Hoyer sho...@gmail.com wrote: I like this idea. But I am -1 on returning None if the array is unstructured. I expect .keys(), if present, to always return an iterable. In fact, this would break some of my existing code, which checks for the existence of keys as a way to do duck typed checks for dictionary like objects (e.g., including pandas.DataFrame): https://github.com/xray/xray/blob/v0.3/xray/core/utils.py#L165 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Tracking and inspecting numpy objects
On Mon, Sep 15, 2014 at 11:55 AM, Sebastian Berg sebast...@sipsolutions.net wrote: On Mo, 2014-09-15 at 10:16 +0200, Mads Ipsen wrote: Hi, I am trying to inspect the reference count of numpy arrays generated by my application. Initially, I thought I could inspect the tracked objects using gc.get_objects(), but, with respect to numpy objects, the returned container is empty. For example: import numpy import gc data = numpy.ones(1024).reshape((32,32)) objs = [o for o in gc.get_objects() if isinstance(o, numpy.ndarray)] print objs# Prints empty list print gc.is_tracked(data) # Print False Why is this? Also, is there some other technique I can use to inspect all numpy generated objects? Two reasons. First of all, unless your array is an object arrays (or a structured one with objects in it), there are no objects to track. The array is a single python object without any referenced objects (except possibly its `arr.base`). Second of all -- and this is an issue -- numpy doesn't actually implement the traverse slot, so it won't even work for object arrays (numpy object arrays do not support circular garbage collection at this time, please feel free to implement it ;)). - Sebastian Does this answer why the ndarray object itself isn't tracked though? I must say I find this puzzling; the only thing I can think of is that the python compiler notices that data isn't used anymore after its creation, and deletes it right after its creation as an optimization, but that conflicts with my own experience of the GC. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Tracking and inspecting numpy objects
I figured the reference to the object through the local scope would also be tracked by the GC somehow, since the entire stack frame can be regarded as a separate object itself, but apparently not. On Mon, Sep 15, 2014 at 1:06 PM, Mads Ipsen mads.ip...@gmail.com wrote: Thanks to everybody for taking time to answer! Best regards, Mads On 15/09/14 12:11, Sebastian Berg wrote: On Mo, 2014-09-15 at 12:05 +0200, Eelco Hoogendoorn wrote: On Mon, Sep 15, 2014 at 11:55 AM, Sebastian Berg sebast...@sipsolutions.net wrote: On Mo, 2014-09-15 at 10:16 +0200, Mads Ipsen wrote: Hi, I am trying to inspect the reference count of numpy arrays generated by my application. Initially, I thought I could inspect the tracked objects using gc.get_objects(), but, with respect to numpy objects, the returned container is empty. For example: import numpy import gc data = numpy.ones(1024).reshape((32,32)) objs = [o for o in gc.get_objects() if isinstance(o, numpy.ndarray)] print objs# Prints empty list print gc.is_tracked(data) # Print False Why is this? Also, is there some other technique I can use to inspect all numpy generated objects? Two reasons. First of all, unless your array is an object arrays (or a structured one with objects in it), there are no objects to track. The array is a single python object without any referenced objects (except possibly its `arr.base`). Second of all -- and this is an issue -- numpy doesn't actually implement the traverse slot, so it won't even work for object arrays (numpy object arrays do not support circular garbage collection at this time, please feel free to implement it ;)). - Sebastian Does this answer why the ndarray object itself isn't tracked though? I must say I find this puzzling; the only thing I can think of is that the python compiler notices that data isn't used anymore after its creation, and deletes it right after its creation as an optimization, but that conflicts with my own experience of the GC. Not sure if it does, but my quick try and error says: In [15]: class T(tuple): : pass : In [16]: t = T() In [17]: objs = [o for o in gc.get_objects() if isinstance(o, T)] In [18]: objs Out[18]: [()] In [19]: a = 123. In [20]: objs = [o for o in gc.get_objects() if isinstance(o, float)] In [21]: objs Out[21]: [] So I guess nothing is tracked, unless it contains things, and numpy arrays don't say they can contain things (i.e. no traverse). - Sebastian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- +-+ | Mads Ipsen | +--+--+ | Gåsebæksvej 7, 4. tv | phone: +45-29716388 | | DK-2500 Valby| email: mads.ip...@gmail.com | | Denmark | map : www.tinyurl.com/ns52fpa | +--+--+ ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] why does u.resize return None?
agreed; I never saw the logic in returning none either. On Thu, Sep 11, 2014 at 4:27 PM, Neal Becker ndbeck...@gmail.com wrote: It would be useful if u.resize returned the new array, so it could be used for chaining operations -- -- Those who don't understand recursion are doomed to repeat it ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Generalize hstack/vstack -- stack; Block matrices like in matlab
Sturla: im not sure if the intention is always unambiguous, for such more flexible arrangements. Also, I doubt such situations arise often in practice; if the arrays arnt a grid, they are probably a nested grid, and the code would most naturally concatenate them with nested calls to a stacking function. However, some form of nd-stack function would be neat in my opinion. On Mon, Sep 8, 2014 at 6:10 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Mon, Sep 8, 2014 at 7:41 AM, Sturla Molden sturla.mol...@gmail.com wrote: Stefan Otte stefan.o...@gmail.com wrote: stack([[a, b], [c, d]]) In my case `stack` replaced `hstack` and `vstack` almost completely. If you're interested in including it in numpy I created a pull request [1]. I'm looking forward to getting some feedback! As far as I can see, it uses hstack and vstack. But that means a and b have to have the same number of rows, c and d must have the same rumber of rows, and hstack((a,b)) and hstack((c,d)) must have the same number of columns. Thus it requires a regularity like this: BB BB CCCDDD CCCDDD CCCDDD CCCDDD What if we just ignore this constraint, and only require the output to be rectangular? Now we have a 'tetris game': BB BB BB BB DD DD or BB BB BB BB BB BB This should be 'stackable', yes? Or perhaps we need another stacking function for this, say numpy.tetris? And while we're at it, what about higher dimensions? should there be an ndstack function too? This is starting to look like the second time in a row Stefan tries to extend numpy with a simple convenience function, and he gets tricked into implementing some sophisticated algorithm... For his next PR I expect nothing less than an NP-complete problem. ;-) Jaime -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Generalize hstack/vstack -- stack; Blockmatrices like in matlab
Blaze aims to do something like that; to make the notion of an array and how it stores it's data far more flexible. But if it isn't a single strided ND array, it isn't numpy. This concept lies at its very heart; and for good reasons I would add. -Original Message- From: Benjamin Root ben.r...@ou.edu Sent: 8-9-2014 19:00 To: Discussion of Numerical Python numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] Generalize hstack/vstack -- stack; Blockmatrices like in matlab Btw, on a somewhat related note, whoever can implement ndarray to be able to use views from other ndarrays stitched together would get a fruit basket from me come the holidays and possibly naming rights for the next kid... Cheers! Ben Root On Mon, Sep 8, 2014 at 12:55 PM, Benjamin Root ben.r...@ou.edu wrote: A use case would be image stitching or even data tiling. I have had to implement something like this at work (so, I can't share it, unfortunately) and it even goes so far as to allow the caller to specify how much the tiles can overlap and such. The specification is ungodly hideous and I doubt I would be willing to share it even if I could lest I release code-thulu upon the world... I think just having this generalize stack feature would be nice start. Tetris could be built on top of that later. (Although, I do vote for at least 3 or 4 dimensional stacking, if possible). Cheers! Ben Root On Mon, Sep 8, 2014 at 12:41 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Sturla: im not sure if the intention is always unambiguous, for such more flexible arrangements. Also, I doubt such situations arise often in practice; if the arrays arnt a grid, they are probably a nested grid, and the code would most naturally concatenate them with nested calls to a stacking function. However, some form of nd-stack function would be neat in my opinion. On Mon, Sep 8, 2014 at 6:10 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Mon, Sep 8, 2014 at 7:41 AM, Sturla Molden sturla.mol...@gmail.com wrote: Stefan Otte stefan.o...@gmail.com wrote: stack([[a, b], [c, d]]) In my case `stack` replaced `hstack` and `vstack` almost completely. If you're interested in including it in numpy I created a pull request [1]. I'm looking forward to getting some feedback! As far as I can see, it uses hstack and vstack. But that means a and b have to have the same number of rows, c and d must have the same rumber of rows, and hstack((a,b)) and hstack((c,d)) must have the same number of columns. Thus it requires a regularity like this: BB BB CCCDDD CCCDDD CCCDDD CCCDDD What if we just ignore this constraint, and only require the output to be rectangular? Now we have a 'tetris game': BB BB BB BB DD DD or BB BB BB BB BB BB This should be 'stackable', yes? Or perhaps we need another stacking function for this, say numpy.tetris? And while we're at it, what about higher dimensions? should there be an ndstack function too? This is starting to look like the second time in a row Stefan tries to extend numpy with a simple convenience function, and he gets tricked into implementing some sophisticated algorithm... For his next PR I expect nothing less than an NP-complete problem. ;-) Jaime -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
On Wed, Sep 3, 2014 at 6:46 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 9:33 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:41 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Not sure about the hashing. Indeed one can also build an index of a set by means of a hash table, but its questionable if this leads to improved performance over performing an argsort. Hashing may have better asymptotic time complexity in theory, but many datasets used in practice are very easy to sort (O(N)-ish), and the time-constant of hashing is higher. But more importantly, using a hash-table guarantees poor cache behavior for many operations using this index. By contrast, sorting may (but need not) make one random access pass to build the index, and may (but need not) perform one random access to reorder values for grouping. But insofar as the keys are better behaved than pure random, this coherence will be exploited. If you want to give it a try, these branch of my numpy fork has hash table based implementations of unique (with no extra indices) and in1d: https://github.com/jaimefrio/numpy/tree/hash-unique A use cases where the hash table is clearly better: In [1]: import numpy as np In [2]: from numpy.lib._compiled_base import _unique, _in1d In [3]: a = np.random.randint(10, size=(1,)) In [4]: %timeit np.unique(a) 1000 loops, best of 3: 258 us per loop In [5]: %timeit _unique(a) 1 loops, best of 3: 143 us per loop In [6]: %timeit np.sort(_unique(a)) 1 loops, best of 3: 149 us per loop It typically performs between 1.5x and 4x faster than sorting. I haven't profiled it properly to know, but there may be quite a bit of performance to dig out: have type specific comparison functions, optimize the starting hash table size based on the size of the array to avoid reinsertions... If getting the elements sorted is a necessity, and the array contains very few or no repeated items, then the hash table approach may even perform worse,: In [8]: a = np.random.randint(1, size=(5000,)) In [9]: %timeit np.unique(a) 1000 loops, best of 3: 277 us per loop In [10]: %timeit np.sort(_unique(a)) 1000 loops, best of 3: 320 us per loop But the hash table still wins in extracting the unique items only: In [11]: %timeit _unique(a) 1 loops, best of 3: 187 us per loop Where the hash table shines is in more elaborate situations. If you keep the first index where it was found, and the number of repeats, in the hash table, you can get return_index and return_counts almost for free, which means you are performing an extra 3x faster than with sorting. return_inverse requires a little more trickery, so I won;t attempt to quantify the improvement. But I wouldn't be surprised if, after fine tuning it, there is close to an order of magnitude overall improvement The spped-up for in1d is also nice: In [16]: a = np.random.randint(1000, size=(1000,)) In [17]: b = np.random.randint(1000, size=(500,)) In [18]: %timeit np.in1d(a, b) 1000 loops, best of 3: 178 us per loop In [19]: %timeit _in1d(a, b) 1 loops, best of 3: 30.1 us per loop Of course, there is no point in Ooops!!! Hit the send button too quick. Not to extend myself too long: if we are going to rethink all of this, we should approach it with an open mind. Still, and this post is not helping with that either, I am afraid that we are discussing implementation details, but are missing a broader vision of what we want to accomplish and why. That vision of what numpy's grouping functionality, if any, should be, and how it complements or conflicts with what pandas is providing, should precede anything else. I now I haven't, but has anyone looked at how pandas implements grouping? Their documentation on the subject is well worth a read: http://pandas.pydata.org/pandas-docs/stable/groupby.html Does numpy need to replicate this? What/why/how can we add to that? Jaime -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion I would certainly not be opposed to having a hashing based indexing mechanism; I think it would make sense design-wise to have a HashIndex class with the same interface as the rest, and use that subclass in those arraysetops where it makes sense. The 'how to' of indexing and its applications are largely orthogonal I think (with some tiny performance compromises which are worth the abstraction imo). For datasets which are not purely random, have many unique items, and which do not fit into cache, I would expect sorting to come out on top, but indeed it depends on the dataset. Yeah, the question how pandas does grouping, and whether we can do better, is a relevant one. From what
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
On Thu, Sep 4, 2014 at 10:31 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:46 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 9:33 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:41 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Not sure about the hashing. Indeed one can also build an index of a set by means of a hash table, but its questionable if this leads to improved performance over performing an argsort. Hashing may have better asymptotic time complexity in theory, but many datasets used in practice are very easy to sort (O(N)-ish), and the time-constant of hashing is higher. But more importantly, using a hash-table guarantees poor cache behavior for many operations using this index. By contrast, sorting may (but need not) make one random access pass to build the index, and may (but need not) perform one random access to reorder values for grouping. But insofar as the keys are better behaved than pure random, this coherence will be exploited. If you want to give it a try, these branch of my numpy fork has hash table based implementations of unique (with no extra indices) and in1d: https://github.com/jaimefrio/numpy/tree/hash-unique A use cases where the hash table is clearly better: In [1]: import numpy as np In [2]: from numpy.lib._compiled_base import _unique, _in1d In [3]: a = np.random.randint(10, size=(1,)) In [4]: %timeit np.unique(a) 1000 loops, best of 3: 258 us per loop In [5]: %timeit _unique(a) 1 loops, best of 3: 143 us per loop In [6]: %timeit np.sort(_unique(a)) 1 loops, best of 3: 149 us per loop It typically performs between 1.5x and 4x faster than sorting. I haven't profiled it properly to know, but there may be quite a bit of performance to dig out: have type specific comparison functions, optimize the starting hash table size based on the size of the array to avoid reinsertions... If getting the elements sorted is a necessity, and the array contains very few or no repeated items, then the hash table approach may even perform worse,: In [8]: a = np.random.randint(1, size=(5000,)) In [9]: %timeit np.unique(a) 1000 loops, best of 3: 277 us per loop In [10]: %timeit np.sort(_unique(a)) 1000 loops, best of 3: 320 us per loop But the hash table still wins in extracting the unique items only: In [11]: %timeit _unique(a) 1 loops, best of 3: 187 us per loop Where the hash table shines is in more elaborate situations. If you keep the first index where it was found, and the number of repeats, in the hash table, you can get return_index and return_counts almost for free, which means you are performing an extra 3x faster than with sorting. return_inverse requires a little more trickery, so I won;t attempt to quantify the improvement. But I wouldn't be surprised if, after fine tuning it, there is close to an order of magnitude overall improvement The spped-up for in1d is also nice: In [16]: a = np.random.randint(1000, size=(1000,)) In [17]: b = np.random.randint(1000, size=(500,)) In [18]: %timeit np.in1d(a, b) 1000 loops, best of 3: 178 us per loop In [19]: %timeit _in1d(a, b) 1 loops, best of 3: 30.1 us per loop Of course, there is no point in Ooops!!! Hit the send button too quick. Not to extend myself too long: if we are going to rethink all of this, we should approach it with an open mind. Still, and this post is not helping with that either, I am afraid that we are discussing implementation details, but are missing a broader vision of what we want to accomplish and why. That vision of what numpy's grouping functionality, if any, should be, and how it complements or conflicts with what pandas is providing, should precede anything else. I now I haven't, but has anyone looked at how pandas implements grouping? Their documentation on the subject is well worth a read: http://pandas.pydata.org/pandas-docs/stable/groupby.html Does numpy need to replicate this? What/why/how can we add to that? Jaime -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion I would certainly not be opposed to having a hashing based indexing mechanism; I think it would make sense design-wise to have a HashIndex class with the same interface as the rest, and use that subclass in those arraysetops where it makes sense. The 'how to' of indexing and its applications are largely orthogonal I think (with some tiny performance compromises which are worth the abstraction imo). For datasets which are not purely random, have many unique items, and which do not fit into cache, I would expect sorting to come out on top, but indeed it depends on the dataset
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
I should clarify: I am speaking about my implementation, I havnt looked at the numpy implementation for a while so im not sure what it is up to. Note that by 'almost free', we are still talking about three passes over the whole array plus temp allocations, but I am assuming a use-case where the various sorts involved are the dominant cost, which I imagine they are, for out-of-cache sorts. Perhaps this isn't too realistic an assumption about the average use case though, I don't know. Though I suppose its a reasonable guideline to assume that either the dataset is big, or performance isn't that big a concern in the first place. On Thu, Sep 4, 2014 at 7:55 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Thu, Sep 4, 2014 at 10:39 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On Thu, Sep 4, 2014 at 10:31 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:46 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 9:33 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:41 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Not sure about the hashing. Indeed one can also build an index of a set by means of a hash table, but its questionable if this leads to improved performance over performing an argsort. Hashing may have better asymptotic time complexity in theory, but many datasets used in practice are very easy to sort (O(N)-ish), and the time-constant of hashing is higher. But more importantly, using a hash-table guarantees poor cache behavior for many operations using this index. By contrast, sorting may (but need not) make one random access pass to build the index, and may (but need not) perform one random access to reorder values for grouping. But insofar as the keys are better behaved than pure random, this coherence will be exploited. If you want to give it a try, these branch of my numpy fork has hash table based implementations of unique (with no extra indices) and in1d: https://github.com/jaimefrio/numpy/tree/hash-unique A use cases where the hash table is clearly better: In [1]: import numpy as np In [2]: from numpy.lib._compiled_base import _unique, _in1d In [3]: a = np.random.randint(10, size=(1,)) In [4]: %timeit np.unique(a) 1000 loops, best of 3: 258 us per loop In [5]: %timeit _unique(a) 1 loops, best of 3: 143 us per loop In [6]: %timeit np.sort(_unique(a)) 1 loops, best of 3: 149 us per loop It typically performs between 1.5x and 4x faster than sorting. I haven't profiled it properly to know, but there may be quite a bit of performance to dig out: have type specific comparison functions, optimize the starting hash table size based on the size of the array to avoid reinsertions... If getting the elements sorted is a necessity, and the array contains very few or no repeated items, then the hash table approach may even perform worse,: In [8]: a = np.random.randint(1, size=(5000,)) In [9]: %timeit np.unique(a) 1000 loops, best of 3: 277 us per loop In [10]: %timeit np.sort(_unique(a)) 1000 loops, best of 3: 320 us per loop But the hash table still wins in extracting the unique items only: In [11]: %timeit _unique(a) 1 loops, best of 3: 187 us per loop Where the hash table shines is in more elaborate situations. If you keep the first index where it was found, and the number of repeats, in the hash table, you can get return_index and return_counts almost for free, which means you are performing an extra 3x faster than with sorting. return_inverse requires a little more trickery, so I won;t attempt to quantify the improvement. But I wouldn't be surprised if, after fine tuning it, there is close to an order of magnitude overall improvement The spped-up for in1d is also nice: In [16]: a = np.random.randint(1000, size=(1000,)) In [17]: b = np.random.randint(1000, size=(500,)) In [18]: %timeit np.in1d(a, b) 1000 loops, best of 3: 178 us per loop In [19]: %timeit _in1d(a, b) 1 loops, best of 3: 30.1 us per loop Of course, there is no point in Ooops!!! Hit the send button too quick. Not to extend myself too long: if we are going to rethink all of this, we should approach it with an open mind. Still, and this post is not helping with that either, I am afraid that we are discussing implementation details, but are missing a broader vision of what we want to accomplish and why. That vision of what numpy's grouping functionality, if any, should be, and how it complements or conflicts with what pandas is providing, should precede anything else. I now I haven't, but has anyone looked at how pandas implements grouping? Their documentation on the subject is well worth a read: http://pandas.pydata.org/pandas-docs/stable/groupby.html Does numpy need to replicate this? What/why/how can we add to that? Jaime -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
On Thu, Sep 4, 2014 at 8:14 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: I should clarify: I am speaking about my implementation, I havnt looked at the numpy implementation for a while so im not sure what it is up to. Note that by 'almost free', we are still talking about three passes over the whole array plus temp allocations, but I am assuming a use-case where the various sorts involved are the dominant cost, which I imagine they are, for out-of-cache sorts. Perhaps this isn't too realistic an assumption about the average use case though, I don't know. Though I suppose its a reasonable guideline to assume that either the dataset is big, or performance isn't that big a concern in the first place. On Thu, Sep 4, 2014 at 7:55 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Thu, Sep 4, 2014 at 10:39 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On Thu, Sep 4, 2014 at 10:31 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:46 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 9:33 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:41 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Not sure about the hashing. Indeed one can also build an index of a set by means of a hash table, but its questionable if this leads to improved performance over performing an argsort. Hashing may have better asymptotic time complexity in theory, but many datasets used in practice are very easy to sort (O(N)-ish), and the time-constant of hashing is higher. But more importantly, using a hash-table guarantees poor cache behavior for many operations using this index. By contrast, sorting may (but need not) make one random access pass to build the index, and may (but need not) perform one random access to reorder values for grouping. But insofar as the keys are better behaved than pure random, this coherence will be exploited. If you want to give it a try, these branch of my numpy fork has hash table based implementations of unique (with no extra indices) and in1d: https://github.com/jaimefrio/numpy/tree/hash-unique A use cases where the hash table is clearly better: In [1]: import numpy as np In [2]: from numpy.lib._compiled_base import _unique, _in1d In [3]: a = np.random.randint(10, size=(1,)) In [4]: %timeit np.unique(a) 1000 loops, best of 3: 258 us per loop In [5]: %timeit _unique(a) 1 loops, best of 3: 143 us per loop In [6]: %timeit np.sort(_unique(a)) 1 loops, best of 3: 149 us per loop It typically performs between 1.5x and 4x faster than sorting. I haven't profiled it properly to know, but there may be quite a bit of performance to dig out: have type specific comparison functions, optimize the starting hash table size based on the size of the array to avoid reinsertions... If getting the elements sorted is a necessity, and the array contains very few or no repeated items, then the hash table approach may even perform worse,: In [8]: a = np.random.randint(1, size=(5000,)) In [9]: %timeit np.unique(a) 1000 loops, best of 3: 277 us per loop In [10]: %timeit np.sort(_unique(a)) 1000 loops, best of 3: 320 us per loop But the hash table still wins in extracting the unique items only: In [11]: %timeit _unique(a) 1 loops, best of 3: 187 us per loop Where the hash table shines is in more elaborate situations. If you keep the first index where it was found, and the number of repeats, in the hash table, you can get return_index and return_counts almost for free, which means you are performing an extra 3x faster than with sorting. return_inverse requires a little more trickery, so I won;t attempt to quantify the improvement. But I wouldn't be surprised if, after fine tuning it, there is close to an order of magnitude overall improvement The spped-up for in1d is also nice: In [16]: a = np.random.randint(1000, size=(1000,)) In [17]: b = np.random.randint(1000, size=(500,)) In [18]: %timeit np.in1d(a, b) 1000 loops, best of 3: 178 us per loop In [19]: %timeit _in1d(a, b) 1 loops, best of 3: 30.1 us per loop Of course, there is no point in Ooops!!! Hit the send button too quick. Not to extend myself too long: if we are going to rethink all of this, we should approach it with an open mind. Still, and this post is not helping with that either, I am afraid that we are discussing implementation details, but are missing a broader vision of what we want to accomplish and why. That vision of what numpy's grouping functionality, if any, should be, and how it complements or conflicts with what pandas is providing, should precede anything else. I now I haven't, but has anyone looked at how pandas implements grouping? Their documentation on the subject is well worth a read: http://pandas.pydata.org/pandas-docs/stable/groupby.html Does numpy need to replicate
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
Naturally, youd want to avoid redoing the indexing where you can, which is another good reason to factor out the indexing mechanisms into separate classes. A factor two performance difference does not get me too excited; again, I think it would be the other way around for an out-of-cache dataset being grouped. But this by itself is ofcourse another argument for factoring out the indexing behind a uniform interface, so we can play around with those implementation details later, and specialize the indexing to serve different scenarios. Also, it really helps with code maintainability; most arraysetops are almost trivial to implement once you have abstracted away the indexing machinery. On Thu, Sep 4, 2014 at 8:36 PM, Jeff Reback jeffreb...@gmail.com wrote: FYI pandas DOES use a very performant hash table impl for unique (and value_counts). Sorted state IS maintained by underlying Index implmentation. https://github.com/pydata/pandas/blob/master/pandas/hashtable.pyx In [8]: a = np.random.randint(10, size=(1,)) In [9]: %timeit np.unique(a) 1000 loops, best of 3: 284 µs per loop In [10]: %timeit Series(a).unique() 1 loops, best of 3: 161 µs per loop In [11]: s = Series(a) # without the creation overhead In [12]: %timeit s.unique() 1 loops, best of 3: 75.3 µs per loop On Thu, Sep 4, 2014 at 2:29 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On Thu, Sep 4, 2014 at 8:14 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: I should clarify: I am speaking about my implementation, I havnt looked at the numpy implementation for a while so im not sure what it is up to. Note that by 'almost free', we are still talking about three passes over the whole array plus temp allocations, but I am assuming a use-case where the various sorts involved are the dominant cost, which I imagine they are, for out-of-cache sorts. Perhaps this isn't too realistic an assumption about the average use case though, I don't know. Though I suppose its a reasonable guideline to assume that either the dataset is big, or performance isn't that big a concern in the first place. On Thu, Sep 4, 2014 at 7:55 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Thu, Sep 4, 2014 at 10:39 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On Thu, Sep 4, 2014 at 10:31 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:46 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 9:33 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Sep 3, 2014 at 6:41 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Not sure about the hashing. Indeed one can also build an index of a set by means of a hash table, but its questionable if this leads to improved performance over performing an argsort. Hashing may have better asymptotic time complexity in theory, but many datasets used in practice are very easy to sort (O(N)-ish), and the time-constant of hashing is higher. But more importantly, using a hash-table guarantees poor cache behavior for many operations using this index. By contrast, sorting may (but need not) make one random access pass to build the index, and may (but need not) perform one random access to reorder values for grouping. But insofar as the keys are better behaved than pure random, this coherence will be exploited. If you want to give it a try, these branch of my numpy fork has hash table based implementations of unique (with no extra indices) and in1d: https://github.com/jaimefrio/numpy/tree/hash-unique A use cases where the hash table is clearly better: In [1]: import numpy as np In [2]: from numpy.lib._compiled_base import _unique, _in1d In [3]: a = np.random.randint(10, size=(1,)) In [4]: %timeit np.unique(a) 1000 loops, best of 3: 258 us per loop In [5]: %timeit _unique(a) 1 loops, best of 3: 143 us per loop In [6]: %timeit np.sort(_unique(a)) 1 loops, best of 3: 149 us per loop It typically performs between 1.5x and 4x faster than sorting. I haven't profiled it properly to know, but there may be quite a bit of performance to dig out: have type specific comparison functions, optimize the starting hash table size based on the size of the array to avoid reinsertions... If getting the elements sorted is a necessity, and the array contains very few or no repeated items, then the hash table approach may even perform worse,: In [8]: a = np.random.randint(1, size=(5000,)) In [9]: %timeit np.unique(a) 1000 loops, best of 3: 277 us per loop In [10]: %timeit np.sort(_unique(a)) 1000 loops, best of 3: 320 us per loop But the hash table still wins in extracting the unique items only: In [11]: %timeit _unique(a) 1 loops, best of 3: 187 us per loop Where the hash table shines is in more elaborate situations. If you keep the first index where it was found, and the number of repeats
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
On Wed, Sep 3, 2014 at 4:07 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Tue, Sep 2, 2014 at 5:40 PM, Charles R Harris charlesr.har...@gmail.com wrote: What do you think about the suggestion of timsort? One would need to concatenate the arrays before sorting, but it should be fairly efficient. Timsort is very cool, and it would definitely be fun to implement in numpy. It is also a lot more work that merging two sorted arrays! I think +1 if someone else does it, but although I would love to be part of it, I am not sure I will be able to find time to look into it seriously in the next couple of months. From a setops point of view, merging two sorted arrays makes it very straightforward to compute, together with (or instead of) the result of the merge, index arrays that let you calculate things like `in1d` faster. Although perhaps an `argtimsort` could provide the same functionality, not sure. I will probably wrap up what I have, put a lace on it, and submit it as a PR. Even if it is not destined to be merged, it may serve as a warning to others. I have also been thinking lately that one of the problems with all these unique-related computations may be a case of having a hammer and seeing everything as nails. Perhaps the algorithm that needs to be ported from Python is not the sorting one, but the hash table... Jaime Not sure about the hashing. Indeed one can also build an index of a set by means of a hash table, but its questionable if this leads to improved performance over performing an argsort. Hashing may have better asymptotic time complexity in theory, but many datasets used in practice are very easy to sort (O(N)-ish), and the time-constant of hashing is higher. But more importantly, using a hash-table guarantees poor cache behavior for many operations using this index. By contrast, sorting may (but need not) make one random access pass to build the index, and may (but need not) perform one random access to reorder values for grouping. But insofar as the keys are better behaved than pure random, this coherence will be exploited. Also, getting the unique values/keys in sorted order is a nice side-benefit for many applications. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Give Jaime Fernandez commit rights.
+1; though I am relatively new to the scene, Jaime's contributions have always stood out to me as thoughtful. On Thu, Sep 4, 2014 at 12:42 AM, Ralf Gommers ralf.gomm...@gmail.com wrote: On Wed, Sep 3, 2014 at 11:48 PM, Robert Kern robert.k...@gmail.com wrote: On Wed, Sep 3, 2014 at 10:47 PM, Charles R Harris charlesr.har...@gmail.com wrote: Hi All, I'd like to give Jaime commit rights. Having at three active developers with commit rights is the goal and Jaime has been pretty consistent with code submissions and discussion participation. +1 +1 excellent idea Ralf ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
Sure, id like to do the hashing things out, but I would also like some preliminary feedback as to whether this is going in a direction anyone else sees the point of, if it conflicts with other plans, and indeed if we can agree that numpy is the right place for it; a point which I would very much like to defend. If there is some obvious no-go that im missing, I can do without the drudgery of writing proper documentation ;). As for whether this belongs in numpy: yes, I would say so. There are the extension of functionality to functions already in numpy, which are a no-brainer (it need not cost anything performance wise, and ive needed unique graph edges many many times), and there is the grouping functionality, which is the main novelty. However, note that the grouping functionality itself is a very small addition, just a few 100 lines of pure python, given that the indexing logic has been factored out of the classic arraysetops. At least from a developers perspective, it very much feels like a logical extension of the same 'thing'. But also from a conceptual numpy perspective, grouping is really more an 'elementary manipulation of an ndarray' than a 'special purpose algorithm'. It is useful for literally all kinds of programming; hence there is similar functionality in the python standard library (itertools.groupby); so why not have an efficient vectorized equivalent in numpy? It belongs there more than the linalg module, arguably. Also, from a community perspective, a significant fraction of all stackoverflow numpy questions are (unknowingly) exactly about 'how to do grouping in numpy'. On Mon, Sep 1, 2014 at 4:36 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Sun, Aug 31, 2014 at 1:48 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Ive organized all code I had relating to this subject in a github repository https://github.com/EelcoHoogendoorn/Numpy_arraysetops_EP. That should facilitate shooting around ideas. Ive also added more documentation and structure to make it easier to see what is going on. Hopefully we can converge on a common vision, and then improve the documentation and testing to make it worthy of including in the numpy master. Note that there is also a complete rewrite of the classic numpy.arraysetops, such that they are also generalized to more complex input, such as finding unique graph edges, and so on. You mentioned getting the numpy core developers involved; are they not subscribed to this mailing list? I wouldn't be surprised; youd hope there is a channel of discussion concerning development with higher signal to noise There are only about 2.5 of us at the moment. Those for whom this is an itch that need scratching should hash things out and make a PR. The main question for me is if it belongs in numpy, scipy, or somewhere else. Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
On Mon, Sep 1, 2014 at 2:05 PM, Charles R Harris charlesr.har...@gmail.com wrote: On Mon, Sep 1, 2014 at 1:49 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Sure, id like to do the hashing things out, but I would also like some preliminary feedback as to whether this is going in a direction anyone else sees the point of, if it conflicts with other plans, and indeed if we can agree that numpy is the right place for it; a point which I would very much like to defend. If there is some obvious no-go that im missing, I can do without the drudgery of writing proper documentation ;). As for whether this belongs in numpy: yes, I would say so. There are the extension of functionality to functions already in numpy, which are a no-brainer (it need not cost anything performance wise, and ive needed unique graph edges many many times), and there is the grouping functionality, which is the main novelty. However, note that the grouping functionality itself is a very small addition, just a few 100 lines of pure python, given that the indexing logic has been factored out of the classic arraysetops. At least from a developers perspective, it very much feels like a logical extension of the same 'thing'. But also from a conceptual numpy perspective, grouping is really more an 'elementary manipulation of an ndarray' than a 'special purpose algorithm'. It is useful for literally all kinds of programming; hence there is similar functionality in the python standard library (itertools.groupby); so why not have an efficient vectorized equivalent in numpy? It belongs there more than the linalg module, arguably. Also, from a community perspective, a significant fraction of all stackoverflow numpy questions are (unknowingly) exactly about 'how to do grouping in numpy'. What I'm trying to say is that numpy is a community project. We don't have a central planning committee, the only difference between developers and everyone else is activity and commit rights. Which is to say if you develop and push this topic it is likely to go in. There certainly seems to be interest in this functionality. The reason that I brought up scipy is that there are some graph algorithms there that went in a couple of years ago. Note that the convention on the list is bottom posting. snip Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion I understand that numpy is a community project, so that the decision isn't up to any one particular person; but some early stage feedback from those active in the community would be welcome. I am generally confident that this addition makes sense, but I have not contributed to numpy before, and you don't know what you don't know and all... given that there are multiple suggestions for changing arraysetops, some coordination would be useful I think. Note that I use graph edges merely as an example; the proposed functionality is much more general than graphing algorithms specifically. The radial reduction https://github.com/EelcoHoogendoorn/Numpy_arraysetops_EP/blob/master/examples.pyexample I included on github is particularly illustrative of the general utility of grouping functionality I think. Operations like radial reductions are rather common, and a custom implementation is quite lengthy, very bug prone, and potentially very slow. Thanks for the heads up on posting convention; ive always let gmail do my thinking for me, which works well enough for me, but I can see how not following this convention is annoying to others. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
Ive organized all code I had relating to this subject in a github repository https://github.com/EelcoHoogendoorn/Numpy_arraysetops_EP. That should facilitate shooting around ideas. Ive also added more documentation and structure to make it easier to see what is going on. Hopefully we can converge on a common vision, and then improve the documentation and testing to make it worthy of including in the numpy master. Note that there is also a complete rewrite of the classic numpy.arraysetops, such that they are also generalized to more complex input, such as finding unique graph edges, and so on. You mentioned getting the numpy core developers involved; are they not subscribed to this mailing list? I wouldn't be surprised; youd hope there is a channel of discussion concerning development with higher signal to noise On Thu, Aug 28, 2014 at 1:49 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: I just checked the docs on ufuncs, and it appears that's a solved problem now, since ufunc.reduceat now comes with an axis argument. Or maybe it already did when I wrote that, but I simply wasn't paying attention. Either way, the code is fully vectorized now, in both grouped and non-grouped axes. Its a lot of code, but all that happens for a grouping other than some O(1) and O(n) stuff is an argsort of the keys, and then the reduction itself, all fully vectorized. Note that I sort the values first, and then use ufunc.reduceat on the groups. It would seem to me that ufunc.at should be more efficient, by avoiding this indirection, but testing very much revealed the opposite, for reasons unclear to me. Perhaps that's changed now as well. On Wed, Aug 27, 2014 at 11:32 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: Yes, I was aware of that. But the point would be to provide true vectorization on those operations. The way I see it, numpy may not have to have a GroupBy implementation, but it should at least enable implementing one that is fast and efficient over any axis. On Wed, Aug 27, 2014 at 12:38 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: i.e, if the grouped axis is small but the other axes are not, you could write this, which avoids the python loop over the long axis that np.vectorize would otherwise perform. import numpy as np from grouping import group_by keys = np.random.randint(0,4,10) values = np.random.rand(10,2000) for k,g in zip(*group_by(keys)(values)): print k, g.mean(0) On Wed, Aug 27, 2014 at 9:29 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: f.i., this works as expected as well (100 keys of 1d int arrays and 100 values of 1d float arrays): group_by(randint(0,4,(100,2))).mean(rand(100,2)) On Wed, Aug 27, 2014 at 9:27 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: If I understand you correctly, the current implementation supports these operations. All reductions over groups (except for median) are performed through the corresponding ufunc (see GroupBy.reduce). This works on multidimensional arrays as well, although this broadcasting over the non-grouping axes is accomplished using np.vectorize. Actual vectorization only happens over the axis being grouped over, but this is usually a long axis. If it isn't, it is more efficient to perform a reduction by means of splitting the array by its groups first, and then map the iterable of groups over some reduction operation (as noted in the docstring of GroupBy.reduce). On Wed, Aug 27, 2014 at 8:29 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: Hi Eelco, I took a deeper look into your code a couple of weeks back. I don't think I have fully grasped what it allows completely, but I agree that some form of what you have there is highly desirable. Along the same lines, for sometime I have been thinking that the right place for a `groupby` in numpy is as a method of ufuncs, so that `np.add.groupby(arr, groups)` would do a multidimensional version of `np.bincount(groups, weights=arr)`. You would then need a more powerful version of `np.unique` to produce the `groups`, but that is something that Joe Kington's old PR was very close to achieving, that should probably be resurrected as well. But yes, there seems to be material for a NEP here, and some guidance from one of the numpy devs would be helpful in getting this somewhere. Jaime On Wed, Aug 27, 2014 at 10:35 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: It wouldn't hurt to have this function, but my intuition is that its use will be minimal. If you are already working with sorted arrays, you already have a flop cost on that order of magnitude, and the optimized merge saves you a factor two at the very most. Using numpy means you are sacrificing factors of two and beyond relative to pure C left right and center anyway, so if this kind of thing matters to you, you probably wont be working in numpy in the first place. That said
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
It wouldn't hurt to have this function, but my intuition is that its use will be minimal. If you are already working with sorted arrays, you already have a flop cost on that order of magnitude, and the optimized merge saves you a factor two at the very most. Using numpy means you are sacrificing factors of two and beyond relative to pure C left right and center anyway, so if this kind of thing matters to you, you probably wont be working in numpy in the first place. That said, I share your interest in overhauling arraysetops. I see many opportunities for expanding its functionality. There is a question that amounts to 'how do I do group-by in numpy' on stackoverflow almost every week. That would have my top priority, but also things like extending np.unique to things like graph edges, or other more complex input, is very often useful to me. Ive written up a draft http://pastebin.com/c5WLWPbpa while ago which accomplishes all of the above and more. It reimplements functions like np.unique around a common Index object. This index object encapsulates the precomputation (sorting) required for efficient set-ops on different datatypes, and provides a common interface to obtain the kind of information you are talking about (which is used extensively internally in the implementation of group_by, for instance). ie, this functionality allows you to write neat things like group_by(randint(0,9,(100,2))).median(rand(100)) But I have the feeling much more could be done in this direction, and I feel this draft could really use a bit of back and forth. If we are going to completely rewrite arraysetops, we might as well do it right. On Wed, Aug 27, 2014 at 7:02 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: A request was open in github to add a `merge` function to numpy that would merge two sorted 1d arrays into a single sorted 1d array. I have been playing around with that idea for a while, and have a branch in my numpy fork that adds a `mergesorted` function to `numpy.lib`: https://github.com/jaimefrio/numpy/commit/ce5d480afecc989a36e5d2bf4ea1d1ba58a83b0a I drew inspiration from C++ STL algorithms, and merged into a single function what merge, set_union, set_intersection, set_difference and set_symmetric_difference do there. My first thought when implementing this was to not make it a public function, but use it under the hood to speed-up some of the functions of `arraysetops.py`, which are now merging two already sorted functions by doing `np.sort(np.concatenate((a, b)))`. I would need to revisit my testing, but the speed-ups weren't that great. One other thing I saw value in for some of the `arraysetops.py` functions, but couldn't fully figure out, was in providing extra output aside from the merged arrays, either in the form of indices, or of boolean masks, indicating which items of the original arrays made it into the merged one, and/or where did they end up in it. Since there is at least one other person out there that likes it, is there any more interest in such a function? If yes, any comments on what the proper interface for extra output should be? Although perhaps the best is to leave that out for starters and see what use people make of it, if any. Jaime -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
If I understand you correctly, the current implementation supports these operations. All reductions over groups (except for median) are performed through the corresponding ufunc (see GroupBy.reduce). This works on multidimensional arrays as well, although this broadcasting over the non-grouping axes is accomplished using np.vectorize. Actual vectorization only happens over the axis being grouped over, but this is usually a long axis. If it isn't, it is more efficient to perform a reduction by means of splitting the array by its groups first, and then map the iterable of groups over some reduction operation (as noted in the docstring of GroupBy.reduce). On Wed, Aug 27, 2014 at 8:29 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: Hi Eelco, I took a deeper look into your code a couple of weeks back. I don't think I have fully grasped what it allows completely, but I agree that some form of what you have there is highly desirable. Along the same lines, for sometime I have been thinking that the right place for a `groupby` in numpy is as a method of ufuncs, so that `np.add.groupby(arr, groups)` would do a multidimensional version of `np.bincount(groups, weights=arr)`. You would then need a more powerful version of `np.unique` to produce the `groups`, but that is something that Joe Kington's old PR was very close to achieving, that should probably be resurrected as well. But yes, there seems to be material for a NEP here, and some guidance from one of the numpy devs would be helpful in getting this somewhere. Jaime On Wed, Aug 27, 2014 at 10:35 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: It wouldn't hurt to have this function, but my intuition is that its use will be minimal. If you are already working with sorted arrays, you already have a flop cost on that order of magnitude, and the optimized merge saves you a factor two at the very most. Using numpy means you are sacrificing factors of two and beyond relative to pure C left right and center anyway, so if this kind of thing matters to you, you probably wont be working in numpy in the first place. That said, I share your interest in overhauling arraysetops. I see many opportunities for expanding its functionality. There is a question that amounts to 'how do I do group-by in numpy' on stackoverflow almost every week. That would have my top priority, but also things like extending np.unique to things like graph edges, or other more complex input, is very often useful to me. Ive written up a draft http://pastebin.com/c5WLWPbpa while ago which accomplishes all of the above and more. It reimplements functions like np.unique around a common Index object. This index object encapsulates the precomputation (sorting) required for efficient set-ops on different datatypes, and provides a common interface to obtain the kind of information you are talking about (which is used extensively internally in the implementation of group_by, for instance). ie, this functionality allows you to write neat things like group_by(randint(0,9,(100,2))).median(rand(100)) But I have the feeling much more could be done in this direction, and I feel this draft could really use a bit of back and forth. If we are going to completely rewrite arraysetops, we might as well do it right. On Wed, Aug 27, 2014 at 7:02 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: A request was open in github to add a `merge` function to numpy that would merge two sorted 1d arrays into a single sorted 1d array. I have been playing around with that idea for a while, and have a branch in my numpy fork that adds a `mergesorted` function to `numpy.lib`: https://github.com/jaimefrio/numpy/commit/ce5d480afecc989a36e5d2bf4ea1d1ba58a83b0a I drew inspiration from C++ STL algorithms, and merged into a single function what merge, set_union, set_intersection, set_difference and set_symmetric_difference do there. My first thought when implementing this was to not make it a public function, but use it under the hood to speed-up some of the functions of `arraysetops.py`, which are now merging two already sorted functions by doing `np.sort(np.concatenate((a, b)))`. I would need to revisit my testing, but the speed-ups weren't that great. One other thing I saw value in for some of the `arraysetops.py` functions, but couldn't fully figure out, was in providing extra output aside from the merged arrays, either in the form of indices, or of boolean masks, indicating which items of the original arrays made it into the merged one, and/or where did they end up in it. Since there is at least one other person out there that likes it, is there any more interest in such a function? If yes, any comments on what the proper interface for extra output should be? Although perhaps the best is to leave that out for starters and see what use people make of it, if any. Jaime -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
i.e, if the grouped axis is small but the other axes are not, you could write this, which avoids the python loop over the long axis that np.vectorize would otherwise perform. import numpy as np from grouping import group_by keys = np.random.randint(0,4,10) values = np.random.rand(10,2000) for k,g in zip(*group_by(keys)(values)): print k, g.mean(0) On Wed, Aug 27, 2014 at 9:29 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: f.i., this works as expected as well (100 keys of 1d int arrays and 100 values of 1d float arrays): group_by(randint(0,4,(100,2))).mean(rand(100,2)) On Wed, Aug 27, 2014 at 9:27 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: If I understand you correctly, the current implementation supports these operations. All reductions over groups (except for median) are performed through the corresponding ufunc (see GroupBy.reduce). This works on multidimensional arrays as well, although this broadcasting over the non-grouping axes is accomplished using np.vectorize. Actual vectorization only happens over the axis being grouped over, but this is usually a long axis. If it isn't, it is more efficient to perform a reduction by means of splitting the array by its groups first, and then map the iterable of groups over some reduction operation (as noted in the docstring of GroupBy.reduce). On Wed, Aug 27, 2014 at 8:29 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: Hi Eelco, I took a deeper look into your code a couple of weeks back. I don't think I have fully grasped what it allows completely, but I agree that some form of what you have there is highly desirable. Along the same lines, for sometime I have been thinking that the right place for a `groupby` in numpy is as a method of ufuncs, so that `np.add.groupby(arr, groups)` would do a multidimensional version of `np.bincount(groups, weights=arr)`. You would then need a more powerful version of `np.unique` to produce the `groups`, but that is something that Joe Kington's old PR was very close to achieving, that should probably be resurrected as well. But yes, there seems to be material for a NEP here, and some guidance from one of the numpy devs would be helpful in getting this somewhere. Jaime On Wed, Aug 27, 2014 at 10:35 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: It wouldn't hurt to have this function, but my intuition is that its use will be minimal. If you are already working with sorted arrays, you already have a flop cost on that order of magnitude, and the optimized merge saves you a factor two at the very most. Using numpy means you are sacrificing factors of two and beyond relative to pure C left right and center anyway, so if this kind of thing matters to you, you probably wont be working in numpy in the first place. That said, I share your interest in overhauling arraysetops. I see many opportunities for expanding its functionality. There is a question that amounts to 'how do I do group-by in numpy' on stackoverflow almost every week. That would have my top priority, but also things like extending np.unique to things like graph edges, or other more complex input, is very often useful to me. Ive written up a draft http://pastebin.com/c5WLWPbpa while ago which accomplishes all of the above and more. It reimplements functions like np.unique around a common Index object. This index object encapsulates the precomputation (sorting) required for efficient set-ops on different datatypes, and provides a common interface to obtain the kind of information you are talking about (which is used extensively internally in the implementation of group_by, for instance). ie, this functionality allows you to write neat things like group_by(randint(0,9,(100,2))).median(rand(100)) But I have the feeling much more could be done in this direction, and I feel this draft could really use a bit of back and forth. If we are going to completely rewrite arraysetops, we might as well do it right. On Wed, Aug 27, 2014 at 7:02 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: A request was open in github to add a `merge` function to numpy that would merge two sorted 1d arrays into a single sorted 1d array. I have been playing around with that idea for a while, and have a branch in my numpy fork that adds a `mergesorted` function to `numpy.lib`: https://github.com/jaimefrio/numpy/commit/ce5d480afecc989a36e5d2bf4ea1d1ba58a83b0a I drew inspiration from C++ STL algorithms, and merged into a single function what merge, set_union, set_intersection, set_difference and set_symmetric_difference do there. My first thought when implementing this was to not make it a public function, but use it under the hood to speed-up some of the functions of `arraysetops.py`, which are now merging two already sorted functions by doing `np.sort(np.concatenate((a, b)))`. I would need to revisit my testing, but the speed-ups weren't that great
Re: [Numpy-discussion] Does a `mergesorted` function make sense?
I just checked the docs on ufuncs, and it appears that's a solved problem now, since ufunc.reduceat now comes with an axis argument. Or maybe it already did when I wrote that, but I simply wasn't paying attention. Either way, the code is fully vectorized now, in both grouped and non-grouped axes. Its a lot of code, but all that happens for a grouping other than some O(1) and O(n) stuff is an argsort of the keys, and then the reduction itself, all fully vectorized. Note that I sort the values first, and then use ufunc.reduceat on the groups. It would seem to me that ufunc.at should be more efficient, by avoiding this indirection, but testing very much revealed the opposite, for reasons unclear to me. Perhaps that's changed now as well. On Wed, Aug 27, 2014 at 11:32 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: Yes, I was aware of that. But the point would be to provide true vectorization on those operations. The way I see it, numpy may not have to have a GroupBy implementation, but it should at least enable implementing one that is fast and efficient over any axis. On Wed, Aug 27, 2014 at 12:38 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: i.e, if the grouped axis is small but the other axes are not, you could write this, which avoids the python loop over the long axis that np.vectorize would otherwise perform. import numpy as np from grouping import group_by keys = np.random.randint(0,4,10) values = np.random.rand(10,2000) for k,g in zip(*group_by(keys)(values)): print k, g.mean(0) On Wed, Aug 27, 2014 at 9:29 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: f.i., this works as expected as well (100 keys of 1d int arrays and 100 values of 1d float arrays): group_by(randint(0,4,(100,2))).mean(rand(100,2)) On Wed, Aug 27, 2014 at 9:27 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: If I understand you correctly, the current implementation supports these operations. All reductions over groups (except for median) are performed through the corresponding ufunc (see GroupBy.reduce). This works on multidimensional arrays as well, although this broadcasting over the non-grouping axes is accomplished using np.vectorize. Actual vectorization only happens over the axis being grouped over, but this is usually a long axis. If it isn't, it is more efficient to perform a reduction by means of splitting the array by its groups first, and then map the iterable of groups over some reduction operation (as noted in the docstring of GroupBy.reduce). On Wed, Aug 27, 2014 at 8:29 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: Hi Eelco, I took a deeper look into your code a couple of weeks back. I don't think I have fully grasped what it allows completely, but I agree that some form of what you have there is highly desirable. Along the same lines, for sometime I have been thinking that the right place for a `groupby` in numpy is as a method of ufuncs, so that `np.add.groupby(arr, groups)` would do a multidimensional version of `np.bincount(groups, weights=arr)`. You would then need a more powerful version of `np.unique` to produce the `groups`, but that is something that Joe Kington's old PR was very close to achieving, that should probably be resurrected as well. But yes, there seems to be material for a NEP here, and some guidance from one of the numpy devs would be helpful in getting this somewhere. Jaime On Wed, Aug 27, 2014 at 10:35 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: It wouldn't hurt to have this function, but my intuition is that its use will be minimal. If you are already working with sorted arrays, you already have a flop cost on that order of magnitude, and the optimized merge saves you a factor two at the very most. Using numpy means you are sacrificing factors of two and beyond relative to pure C left right and center anyway, so if this kind of thing matters to you, you probably wont be working in numpy in the first place. That said, I share your interest in overhauling arraysetops. I see many opportunities for expanding its functionality. There is a question that amounts to 'how do I do group-by in numpy' on stackoverflow almost every week. That would have my top priority, but also things like extending np.unique to things like graph edges, or other more complex input, is very often useful to me. Ive written up a draft http://pastebin.com/c5WLWPbpa while ago which accomplishes all of the above and more. It reimplements functions like np.unique around a common Index object. This index object encapsulates the precomputation (sorting) required for efficient set-ops on different datatypes, and provides a common interface to obtain the kind of information you are talking about (which is used extensively internally in the implementation of group_by, for instance). ie, this functionality allows you to write neat things like group_by
Re: [Numpy-discussion] np.unique with structured arrays
It does not sound like an issue with unique, but rather like a matter of floating point equality and representation. Do the ' identical' elements pass an equality test? -Original Message- From: Nicolas P. Rougier nicolas.roug...@inria.fr Sent: 22-8-2014 15:21 To: Discussion of Numerical Python numpy-discussion@scipy.org Subject: [Numpy-discussion] np.unique with structured arrays Hello, I've found a strange behavior or I'm missing something obvious (or np.unique is not supposed to work with structured arrays). I'm trying to extract unique values from a simple structured array but it does not seem to work as expected. Here is a minimal script showing the problem: import numpy as np V = np.zeros(4, dtype=[(v, np.float32, 3)]) V[v] = [ [0.5,0.0, 1.0], [0.5, -1.e-16, 1.0], # [0.5, +1.e-16, 1.0] works [0.5,0.0, -1.0], [0.5, -1.e-16, -1.0]] # [0.5, +1.e-16, -1.0]] works V_ = np.zeros_like(V) V_[v][:,0] = V[v][:,0].round(decimals=3) V_[v][:,1] = V[v][:,1].round(decimals=3) V_[v][:,2] = V[v][:,2].round(decimals=3) print np.unique(V_) [([0.5, 0.0, 1.0],) ([0.5, 0.0, -1.0],) ([0.5, -0.0, 1.0],) ([0.5, -0.0, -1.0],)] While I would have expected: [([0.5, 0.0, 1.0],) ([0.5, 0.0, -1.0],)] Can anyone confirm ? Nicolas___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] np.unique with structured arrays
Oh yeah this could be. Floating point equality and bitwise equality are not the
same thing.
-Original Message-
From: Jaime Fernández del Río jaime.f...@gmail.com
Sent: 22-8-2014 16:22
To: Discussion of Numerical Python numpy-discussion@scipy.org
Subject: Re: [Numpy-discussion] np.unique with structured arrays
I can confirm, the issue seems to be in sorting:
np.sort(V_)
array([([0.5, 0.0, 1.0],), ([0.5, 0.0, -1.0],), ([0.5, -0.0, 1.0],),
([0.5, -0.0, -1.0],)],
dtype=[('v', 'f4', (3,))])
These I think are handled by the generic sort functions, and it looks like the
comparison function being used is the one for a VOID dtype with no fields, so
it is being done byte-wise, hence the problems with 0.0 and -0.0. Not sure
where exactly the bug is, though...
Jaime
On Fri, Aug 22, 2014 at 6:20 AM, Nicolas P. Rougier nicolas.roug...@inria.fr
wrote:
Hello,
I've found a strange behavior or I'm missing something obvious (or np.unique is
not supposed to work with structured arrays).
I'm trying to extract unique values from a simple structured array but it does
not seem to work as expected.
Here is a minimal script showing the problem:
import numpy as np
V = np.zeros(4, dtype=[(v, np.float32, 3)])
V[v] = [ [0.5,0.0, 1.0],
[0.5, -1.e-16, 1.0], # [0.5, +1.e-16, 1.0] works
[0.5,0.0, -1.0],
[0.5, -1.e-16, -1.0]] # [0.5, +1.e-16, -1.0]] works
V_ = np.zeros_like(V)
V_[v][:,0] = V[v][:,0].round(decimals=3)
V_[v][:,1] = V[v][:,1].round(decimals=3)
V_[v][:,2] = V[v][:,2].round(decimals=3)
print np.unique(V_)
[([0.5, 0.0, 1.0],) ([0.5, 0.0, -1.0],) ([0.5, -0.0, 1.0],) ([0.5, -0.0,
-1.0],)]
While I would have expected:
[([0.5, 0.0, 1.0],) ([0.5, 0.0, -1.0],)]
Can anyone confirm ?
Nicolas
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
--
(\__/)
( O.o)
( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de
dominación mundial. ___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Proposed new feature for numpy.einsum: repeated output subscripts as diagonal
Agreed; this addition occurred to me as well. Note that the implemenatation
should be straightforward: just allocate an enlarged array, use some
striding logic to construct the relevant view, and let einsums internals
act on the view. hopefully, you wont even have to touch the guts of einsum
at the C level, because id say that isn't for the faint of heart...
On Fri, Aug 15, 2014 at 3:53 PM, Sebastian Berg sebast...@sipsolutions.net
wrote:
On Do, 2014-08-14 at 12:42 -0700, Stephan Hoyer wrote:
I think this would be very nice addition.
On Thu, Aug 14, 2014 at 12:21 PM, Benjamin Root ben.r...@ou.edu
wrote:
You had me at Kronecker delta... :-) +1
Sounds good to me. I don't see a reason for not relaxing the
restriction, unless there is some technical issue, but I doubt that.
- Sebastian
On Thu, Aug 14, 2014 at 3:07 PM, Pierre-Andre Noel
noel.pierre.an...@gmail.com wrote:
(I created issue 4965 earlier today on this topic, and
I have been
advised to email to this mailing list to discuss
whether it is a good
idea or not. I include my original post as-is,
followed by additional
comments.)
I think that the following new feature would make
`numpy.einsum` even
more powerful/useful/awesome than it already is.
Moreover, the change
should not interfere with existing code, it would
preserve the
minimalistic spirit of `numpy.einsum`, and the new
functionality would
integrate in a seamless/intuitive manner for the
users.
In short, the new feature would allow for repeated
subscripts to appear
in the output part of the `subscripts` parameter
(i.e., on the
right-hand side of `-`). The corresponding dimensions
in the resulting
`ndarray` would only be filled along their diagonal,
leaving the off
diagonal entries to the default value for this `dtype`
(typically zero).
Note that the current behavior is to raise an
exception when repeated
output subscripts are being used.
This is simplest to describe with an example involving
the dual behavior
of `numpy.diag`.
```python
# Extracting the diagonal of a 2-D array.
A = arange(16).reshape(4,4)
print(diag(A)) # Output: [ 0 5 10 15 ]
print(einsum('ii-i', A)) # Same as previous line
(current behavior).
# Constructing a diagonal 2-D array.
v = arange(4)
print(diag(v)) # Output: [[0 0 0 0] [0 1 0 0] [0 0 2
0] [0 0 0 3]]
print(einsum('i-ii', v)) # New behavior would be same
as previous line.
# The current behavior of the previous line is to
raise an exception.
```
By opposition to `numpy.diag`, the approach
generalizes to higher
dimensions: `einsum('iii-i', A)` extracts the
diagonal of a 3-D array,
and `einsum('i-iii', v)` would build a diagonal 3-D
array.
The proposed behavior really starts to shine in more
intricate cases.
```python
# Dummy values, these should be probabilities to make
sense below.
P_w_ab = arange(24).reshape(3,2,4)
P_y_wxab = arange(144).reshape(3,3,2,2,4)
# With the proposed behavior, the following two lines
should be equivalent.
P_xyz_ab = einsum('wab,xa,ywxab,zy-xyzab', P_w_ab,
eye(2), P_y_wxab,
eye(3))
also_P_xyz_ab = einsum('wab,ywaab-ayyab', P_w_ab,
P_y_wxab)
```
If this is not convincing enough, replace `eye(2)` by
`eye(P_w_ab.shape[1])` and replace `eye(3)` by
`eye(P_y_wxab.shape[0])`,
then imagine more dimensions and repeated indices...
The new notation
would allow for crisper codes and reduce the
opportunities for dumb
mistakes.
For those who wonder, the above computation amounts to
$P(X=x,Y=y,Z=z|A=a,B=b) = \sum_w P(W=w|A=a,B=b) P(X=x|
A=a)
Re: [Numpy-discussion] Proposed new feature for numpy.einsum: repeated output subscripts as diagonal
Well, there is the numpy-API C level, and then there is the arcane macro C
level. The two might as well be a completely different language.
Indeed, it should be doing something similar for the inputs. Actually, I
think I wrote a wrapper around einsum/numexpr once that performed this
generalized indexing once... ill see if I can dig that up.
On Fri, Aug 15, 2014 at 5:01 PM, Sebastian Berg sebast...@sipsolutions.net
wrote:
On Fr, 2014-08-15 at 16:42 +0200, Eelco Hoogendoorn wrote:
Agreed; this addition occurred to me as well. Note that the
implemenatation should be straightforward: just allocate an enlarged
array, use some striding logic to construct the relevant view, and let
einsums internals act on the view. hopefully, you wont even have to
touch the guts of einsum at the C level, because id say that isn't for
the faint of heart...
I am not sure if einsum isn't pure C :). But even if, it should be doing
something identical already for duplicate indices on the inputs...
- Sebastian
On Fri, Aug 15, 2014 at 3:53 PM, Sebastian Berg
sebast...@sipsolutions.net wrote:
On Do, 2014-08-14 at 12:42 -0700, Stephan Hoyer wrote:
I think this would be very nice addition.
On Thu, Aug 14, 2014 at 12:21 PM, Benjamin Root
ben.r...@ou.edu
wrote:
You had me at Kronecker delta... :-) +1
Sounds good to me. I don't see a reason for not relaxing the
restriction, unless there is some technical issue, but I doubt
that.
- Sebastian
On Thu, Aug 14, 2014 at 3:07 PM, Pierre-Andre Noel
noel.pierre.an...@gmail.com wrote:
(I created issue 4965 earlier today on this
topic, and
I have been
advised to email to this mailing list to
discuss
whether it is a good
idea or not. I include my original post
as-is,
followed by additional
comments.)
I think that the following new feature would
make
`numpy.einsum` even
more powerful/useful/awesome than it already
is.
Moreover, the change
should not interfere with existing code, it
would
preserve the
minimalistic spirit of `numpy.einsum`, and
the new
functionality would
integrate in a seamless/intuitive manner for
the
users.
In short, the new feature would allow for
repeated
subscripts to appear
in the output part of the `subscripts`
parameter
(i.e., on the
right-hand side of `-`). The corresponding
dimensions
in the resulting
`ndarray` would only be filled along their
diagonal,
leaving the off
diagonal entries to the default value for
this `dtype`
(typically zero).
Note that the current behavior is to raise
an
exception when repeated
output subscripts are being used.
This is simplest to describe with an example
involving
the dual behavior
of `numpy.diag`.
```python
# Extracting the diagonal of a 2-D array.
A = arange(16).reshape(4,4)
print(diag(A)) # Output: [ 0 5 10 15 ]
print(einsum('ii-i', A)) # Same as previous
line
(current behavior).
# Constructing a diagonal 2-D array.
v = arange(4)
print(diag(v)) # Output: [[0 0 0 0] [0 1 0
0] [0 0 2
0] [0 0 0 3]]
print(einsum('i-ii', v)) # New behavior
would be same
as previous line.
# The current behavior of the previous line
is to
raise an exception.
```
By opposition to `numpy.diag`, the approach
generalizes to higher
Re: [Numpy-discussion] Proposed new feature for numpy.einsum: repeated output subscripts as diagonal
here is a snippet I extracted from a project with similar aims (integrating
the functionality of einsum and numexpr, actually)
Not much to it, but in case someone needs a reminder on how to use striding
tricks:
http://pastebin.com/kQNySjcj
On Fri, Aug 15, 2014 at 5:20 PM, Eelco Hoogendoorn
hoogendoorn.ee...@gmail.com wrote:
Well, there is the numpy-API C level, and then there is the arcane macro C
level. The two might as well be a completely different language.
Indeed, it should be doing something similar for the inputs. Actually, I
think I wrote a wrapper around einsum/numexpr once that performed this
generalized indexing once... ill see if I can dig that up.
On Fri, Aug 15, 2014 at 5:01 PM, Sebastian Berg
sebast...@sipsolutions.net wrote:
On Fr, 2014-08-15 at 16:42 +0200, Eelco Hoogendoorn wrote:
Agreed; this addition occurred to me as well. Note that the
implemenatation should be straightforward: just allocate an enlarged
array, use some striding logic to construct the relevant view, and let
einsums internals act on the view. hopefully, you wont even have to
touch the guts of einsum at the C level, because id say that isn't for
the faint of heart...
I am not sure if einsum isn't pure C :). But even if, it should be doing
something identical already for duplicate indices on the inputs...
- Sebastian
On Fri, Aug 15, 2014 at 3:53 PM, Sebastian Berg
sebast...@sipsolutions.net wrote:
On Do, 2014-08-14 at 12:42 -0700, Stephan Hoyer wrote:
I think this would be very nice addition.
On Thu, Aug 14, 2014 at 12:21 PM, Benjamin Root
ben.r...@ou.edu
wrote:
You had me at Kronecker delta... :-) +1
Sounds good to me. I don't see a reason for not relaxing the
restriction, unless there is some technical issue, but I doubt
that.
- Sebastian
On Thu, Aug 14, 2014 at 3:07 PM, Pierre-Andre Noel
noel.pierre.an...@gmail.com wrote:
(I created issue 4965 earlier today on this
topic, and
I have been
advised to email to this mailing list to
discuss
whether it is a good
idea or not. I include my original post
as-is,
followed by additional
comments.)
I think that the following new feature would
make
`numpy.einsum` even
more powerful/useful/awesome than it already
is.
Moreover, the change
should not interfere with existing code, it
would
preserve the
minimalistic spirit of `numpy.einsum`, and
the new
functionality would
integrate in a seamless/intuitive manner for
the
users.
In short, the new feature would allow for
repeated
subscripts to appear
in the output part of the `subscripts`
parameter
(i.e., on the
right-hand side of `-`). The corresponding
dimensions
in the resulting
`ndarray` would only be filled along their
diagonal,
leaving the off
diagonal entries to the default value for
this `dtype`
(typically zero).
Note that the current behavior is to raise
an
exception when repeated
output subscripts are being used.
This is simplest to describe with an example
involving
the dual behavior
of `numpy.diag`.
```python
# Extracting the diagonal of a 2-D array.
A = arange(16).reshape(4,4)
print(diag(A)) # Output: [ 0 5 10 15 ]
print(einsum('ii-i', A)) # Same as previous
line
(current behavior).
# Constructing a diagonal 2-D array.
v = arange(4)
print(diag(v)) # Output: [[0 0 0 0] [0 1 0
0] [0 0 2
0] [0 0 0 3]]
print(einsum('i-ii', v)) # New behavior
would be same
Re: [Numpy-discussion] New function `count_unique` to generate contingency tables.
Its pretty easy to implement this table functionality and more on top of the code I linked above. I still think such a comprehensive overhaul of arraysetops is worth discussing. import numpy as np import grouping x = [1, 1, 1, 1, 2, 2, 2, 2, 2] y = [3, 4, 3, 3, 3, 4, 5, 5, 5] z = np.random.randint(0,2,(9,2)) def table(*keys): desired table implementation, building on the index object cleaner, and more functionality performance should be the same indices = [grouping.as_index(k, axis=0) for k in keys] uniques = [i.unique for i in indices] inverses = [i.inverse for i in indices] shape= [i.groups for i in indices] t = np.zeros(shape, np.int) np.add.at(t, inverses, 1) return tuple(uniques), t #here is how to use print table(x,y) #but we can use fancy keys as well; here a composite key and a row-key print table((x,y), z) #this effectively creates a sparse matrix equivalent of your desired table print grouping.count((x,y)) On Wed, Aug 13, 2014 at 11:25 PM, Warren Weckesser warren.weckes...@gmail.com wrote: On Wed, Aug 13, 2014 at 5:15 PM, Benjamin Root ben.r...@ou.edu wrote: The ever-wonderful pylab mode in matplotlib has a table function for plotting a table of text in a plot. If I remember correctly, what would happen is that matplotlib's table() function will simply obliterate the numpy's table function. This isn't a show-stopper, I just wanted to point that out. Personally, while I wasn't a particular fan of count_unique because I wouldn't necessarially think of it when needing a contingency table, I do like that it is verb-ish. table(), in this sense, is not a verb. That said, I am perfectly fine with it if you are fine with the name collision in pylab mode. Thanks for pointing that out. I only changed it to have something that sounded more table-ish, like the Pandas, R and Matlab functions. I won't update it right now, but if there is interest in putting it into numpy, I'll rename it to avoid the pylab conflict. Anything along the lines of `crosstab`, `xtable`, etc., would be fine with me. Warren On Wed, Aug 13, 2014 at 4:57 PM, Warren Weckesser warren.weckes...@gmail.com wrote: On Tue, Aug 12, 2014 at 12:51 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: ah yes, that's also an issue I was trying to deal with. the semantics I prefer in these type of operators, is (as a default), to have every array be treated as a sequence of keys, so if calling unique(arr_2d), youd get unique rows, unless you pass axis=None, in which case the array is flattened. I also agree that the extension you propose here is useful; but ideally, with a little more discussion on these subjects we can converge on an even more comprehensive overhaul On Tue, Aug 12, 2014 at 6:33 PM, Joe Kington joferking...@gmail.com wrote: On Tue, Aug 12, 2014 at 11:17 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Thanks. Prompted by that stackoverflow question, and similar problems I had to deal with myself, I started working on a much more general extension to numpy's functionality in this space. Like you noted, things get a little panda-y, but I think there is a lot of panda's functionality that could or should be part of the numpy core, a robust set of grouping operations in particular. see pastebin here: http://pastebin.com/c5WLWPbp On a side note, this is related to a pull request of mine from awhile back: https://github.com/numpy/numpy/pull/3584 There was a lot of disagreement on the mailing list about what to call a unique slices along a given axis function, so I wound up closing the pull request pending more discussion. At any rate, I think it's a useful thing to have in base numpy. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion Update: I renamed the function to `table` in the pull request: https://github.com/numpy/numpy/pull/4958 Warren ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] New function `count_unique` to generate contingency tables.
Thanks. Prompted by that stackoverflow question, and similar problems I had to deal with myself, I started working on a much more general extension to numpy's functionality in this space. Like you noted, things get a little panda-y, but I think there is a lot of panda's functionality that could or should be part of the numpy core, a robust set of grouping operations in particular. see pastebin here: http://pastebin.com/c5WLWPbp Ive posted about it on this list before, but without apparent interest; and I havnt gotten around to getting this up to professional standards yet either. But there is a lot more that could be done in this direction. Note that the count functionality in the stackoverflow answer is relatively indirect and inefficient, using the inverse_index and such. A much more efficient method is obtained by the code used here. On Tue, Aug 12, 2014 at 5:57 PM, Warren Weckesser warren.weckes...@gmail.com wrote: On Tue, Aug 12, 2014 at 11:35 AM, Warren Weckesser warren.weckes...@gmail.com wrote: I created a pull request (https://github.com/numpy/numpy/pull/4958) that defines the function `count_unique`. `count_unique` generates a contingency table from a collection of sequences. For example, In [7]: x = [1, 1, 1, 1, 2, 2, 2, 2, 2] In [8]: y = [3, 4, 3, 3, 3, 4, 5, 5, 5] In [9]: (xvals, yvals), counts = count_unique(x, y) In [10]: xvals Out[10]: array([1, 2]) In [11]: yvals Out[11]: array([3, 4, 5]) In [12]: counts Out[12]: array([[3, 1, 0], [1, 1, 3]]) It can be interpreted as a multi-argument generalization of `np.unique(x, return_counts=True)`. It overlaps with Pandas' `crosstab`, but I think this is a pretty fundamental counting operation that fits in numpy. Matlab's `crosstab` (http://www.mathworks.com/help/stats/crosstab.html) and R's `table` perform the same calculation (with a few more bells and whistles). For comparison, here's Pandas' `crosstab` (same `x` and `y` as above): In [28]: import pandas as pd In [29]: xs = pd.Series(x) In [30]: ys = pd.Series(y) In [31]: pd.crosstab(xs, ys) Out[31]: col_0 3 4 5 row_0 1 3 1 0 2 1 1 3 And here is R's `table`: x - c(1,1,1,1,2,2,2,2,2) y - c(3,4,3,3,3,4,5,5,5) table(x, y) y x 3 4 5 1 3 1 0 2 1 1 3 Is there any interest in adding this (or some variation of it) to numpy? Warren While searching StackOverflow in the numpy tag for count unique, I just discovered that I basically reinvented Eelco Hoogendoorn's code in his answer to http://stackoverflow.com/questions/10741346/numpy-frequency-counts-for-unique-values-in-an-array. Nice one, Eelco! Warren ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] New function `count_unique` to generate contingency tables.
ah yes, that's also an issue I was trying to deal with. the semantics I prefer in these type of operators, is (as a default), to have every array be treated as a sequence of keys, so if calling unique(arr_2d), youd get unique rows, unless you pass axis=None, in which case the array is flattened. I also agree that the extension you propose here is useful; but ideally, with a little more discussion on these subjects we can converge on an even more comprehensive overhaul On Tue, Aug 12, 2014 at 6:33 PM, Joe Kington joferking...@gmail.com wrote: On Tue, Aug 12, 2014 at 11:17 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Thanks. Prompted by that stackoverflow question, and similar problems I had to deal with myself, I started working on a much more general extension to numpy's functionality in this space. Like you noted, things get a little panda-y, but I think there is a lot of panda's functionality that could or should be part of the numpy core, a robust set of grouping operations in particular. see pastebin here: http://pastebin.com/c5WLWPbp On a side note, this is related to a pull request of mine from awhile back: https://github.com/numpy/numpy/pull/3584 There was a lot of disagreement on the mailing list about what to call a unique slices along a given axis function, so I wound up closing the pull request pending more discussion. At any rate, I think it's a useful thing to have in base numpy. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Calculation of a hessian
Do it in pure numpy? How about copying the source of numdifftools? What exactly is the obstacle to using numdifftools? There seem to be no licensing issues. In my experience, its a crafty piece of work; and calculating a hessian correctly, accounting for all kinds of nasty floating point issues, is no walk in the park. Even if an analytical derivative isn't too big a pain in the ass to implement, there is a good chance that what numdifftools does is more numerically stable (though in all likelihood much slower). The only good reason for a specialized solution I can think of is speed; but be aware what you are trading it in for. If speed is your major concern though, you really cant go wrong with Theano. http://deeplearning.net/software/theano/library/gradient.html#theano.gradient.hessian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Preliminary thoughts on implementing __matmul__
I don't expect stacked matrices/vectors to be used often, although there are some areas that might make heavy use of them, so I think we could live with the simple implementation, it's just a bit of a wart when there is broadcasting of arrays. Just to be clear, the '@' broadcasting differs from the dot broadcasting, agreed? This lack of elegance and unity combined with frankly; a lack of utility, made me plead @ is a bad idea in the first place; but I guess I lost that debate... On Thu, Aug 7, 2014 at 2:00 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Aug 6, 2014 at 5:51 PM, Nathaniel Smith n...@pobox.com wrote: On 7 Aug 2014 00:41, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Aug 6, 2014 at 5:33 PM, Nathaniel Smith n...@pobox.com wrote: On Thu, Aug 7, 2014 at 12:24 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Aug 6, 2014 at 4:57 PM, Nathaniel Smith n...@pobox.com wrote: On Wed, Aug 6, 2014 at 4:32 PM, Charles R Harris charlesr.har...@gmail.com wrote: Should also mention that we don't have the ability to operate on stacked vectors because they can't be identified by dimension info. One workaround is to add dummy dimensions where needed, another is to add two flags, row and col, and set them appropriately. Two flags are needed for backward compatibility, i.e., both false is a traditional array. It's possible I could be convinced to like this, but it would take a substantial amount of convincing :-). It seems like a pretty big violation of orthogonality/one obvious way/etc. to have two totally different ways of representing row/column vectors. The '@' operator supports matrix stacks, so it would seem we also need to support vector stacks. The new addition would only be effective with the '@' operator. The main problem I see with flags is that adding them would require an extensive audit of the C code to make sure they were preserved. Another option, already supported to a large extent, is to have row and col classes inheriting from ndarray that add nothing, except for a possible new transpose type function/method. I did mock up such a class just for fun, and also added a 'dyad' function. If we really don't care to support stacked vectors we can get by without adding anything. It's possible you could convince me that this is a good idea, but I'm starting at like -0.95 :-). Wouldn't it be vastly simpler to just have np.linalg.matvec, matmat, vecvec or something (each of which are single-liners in terms of @), rather than deal with two different ways of representing row/column vectors everywhere? Sure, but matvec and vecvec would not be supported by '@' except when vec was 1d because there is no way to distinguish a stack of vectors from a matrix or a stack of matrices. Yes. But @ can never be magic - either people will have to write something extra to flip these flags on their array objects, or they'll have to write something extra to describe which operation they want. @ was never intended to cover every case, just the simple-but-super-common ones that dot covers, plus a few more (simple broadcasting). We have np.add even though + exists too... I don't expect stacked matrices/vectors to be used often, although there are some areas that might make heavy use of them, so I think we could live with the simple implementation, it's just a bit of a wart when there is broadcasting of arrays. Just to be clear, the '@' broadcasting differs from the dot broadcasting, agreed? Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Array2 subset of array1
np.all(np.in1d(array1,array2)) On Tue, Aug 5, 2014 at 2:58 PM, Jurgens de Bruin debrui...@gmail.com wrote: Hi, I am new to numpy so any help would be greatly appreciated. I have two arrays: array1 = np.arange(1,100+1) array2 = np.arange(1,50+1) How can I calculate/determine if array2 is a subset of array1 (falls within array 1) Something like : array2 in array1 = TRUE for the case above. Thank ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Array2 subset of array1
ah yes, that may indeed be what you want. depending on your datatype, you could access the underlying raw data as a string. b.tostring() in a.tostring() sort of works; but isn't entirely safe, as you may have false positive matches which arnt aligned to your datatype using str.find in combination with dtype.itemsize could solve that problem; though it isn't the most elegant solution id say. also note that you need to check for identical datatypes and memory layout for this to guarantee correct results. On Tue, Aug 5, 2014 at 6:33 PM, Sebastian Berg sebast...@sipsolutions.net wrote: On Di, 2014-08-05 at 14:58 +0200, Jurgens de Bruin wrote: Hi, I am new to numpy so any help would be greatly appreciated. I have two arrays: array1 = np.arange(1,100+1) array2 = np.arange(1,50+1) How can I calculate/determine if array2 is a subset of array1 (falls within array 1) Something like : array2 in array1 = TRUE for the case above. Just to be clear. You are looking for the whole of array1 (as a block/subarray) as far as I understand. And there is no obvious numpy way to do this. Depending on your array sizes, you could blow up the first array from (N,) to (N-M+1,M) and then check if any row matches completely. There may be better tricks available though, especially if array1 is large. - Sebastian Thank ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
To rephrase my most pressing question: may np.ones((N,2)).mean(0) and np.ones((2,N)).mean(1) produce different results with the implementation in the current master? If so, I think that would be very much regrettable; and if this is a minority opinion, I do hope that at least this gets documented in a most explicit manner. On Sun, Jul 27, 2014 at 8:26 PM, Sturla Molden sturla.mol...@gmail.com wrote: Nathaniel Smith n...@pobox.com wrote: The problem here is that when summing up the values, the sum gets large enough that after rounding, x + 1 = x and the sum stops increasing. Interesting. That explains why the divide-and-conquer reduction is much more robust. Thanks :) Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
Sebastian: Those are good points. Indeed iteration order may already produce different results, even though the semantics of numpy suggest identical operations. Still, I feel this different behavior without any semantical clues is something to be minimized. Indeed copying might have large speed implications. But on second thought, does it? Either the data is already aligned and no copy is required, or it isn't aligned, and we need one pass of cache inefficient access to the data anyway. Infact, if we had one low level function which does cache-intelligent transposition of numpy arrays (using some block strategy), this might be faster even than the current simple reduction operations when forced to work on awkwardly aligned data. Ideally, this intelligent access and intelligent reduction would be part of a single pass of course; but that wouldn't really fit within the numpy design, and merely such an intelligent transpose would provide most of the benefit I think. Or is the mechanism behind ascontiguousarray already intelligent in this sense? On Mon, Jul 28, 2014 at 4:06 PM, Sebastian Berg sebast...@sipsolutions.net wrote: On Mo, 2014-07-28 at 15:35 +0200, Sturla Molden wrote: On 28/07/14 15:21, alex wrote: Are you sure they always give different results? Notice that np.ones((N,2)).mean(0) np.ones((2,N)).mean(1) compute means of different axes on transposed arrays so these differences 'cancel out'. They will be if different algorithms are used. np.ones((N,2)).mean(0) will have larger accumulated rounding error than np.ones((2,N)).mean(1), if only the latter uses the divide-and-conquer summation. What I wanted to point out is that to some extend the algorithm does not matter. You will not necessarily get identical results already if you use a different iteration order, and we have been doing that for years for speed reasons. All libs like BLAS do the same. Yes, the new changes make this much more dramatic, but they only make some paths much better, never worse. It might be dangerous, but only in the sense that you test it with the good path and it works good enough, but later (also) use the other one in some lib. I am not even sure if I I would suggest that in the first case we try to copy the array to a temporary contiguous buffer and use the same divide-and-conquer algorithm, unless some heuristics on memory usage fails. Sure, but you have to make major changes to the buffered iterator to do that without larger speed implications. It might be a good idea, but it requires someone who knows this stuff to spend a lot of time and care in the depths of numpy. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
I see, thanks for the clarification. Just for the sake of argument, since unfortunately I don't have the time to go dig in the guts of numpy myself: a design which always produces results of the same (high) accuracy, but only optimizes the common access patterns in a hacky way, and may be inefficient in case it needs to fall back on dumb iteration or array copying, is the best compromise between features and the ever limiting amount of time available, I would argue, no? Its preferable if your code works, but may be hacked to work more efficiently, than that it works efficiently, but may need hacking to work correctly under all circumstances. But fun as it is to think about what ought to be, i suppose the people who do actually pour in the effort have thought about these things already. A numpy 2.0 could probably borrow/integrate a lot from numexpr, I suppose. By the way, the hierarchical summation would make it fairly easy to erase (and in any case would minimize) summation differences due to differences between logical and actual ordering in memory of the data, no? On Mon, Jul 28, 2014 at 5:22 PM, Sebastian Berg sebast...@sipsolutions.net wrote: On Mo, 2014-07-28 at 16:31 +0200, Eelco Hoogendoorn wrote: Sebastian: Those are good points. Indeed iteration order may already produce different results, even though the semantics of numpy suggest identical operations. Still, I feel this different behavior without any semantical clues is something to be minimized. Indeed copying might have large speed implications. But on second thought, does it? Either the data is already aligned and no copy is required, or it isn't aligned, and we need one pass of cache inefficient access to the data anyway. Infact, if we had one low level function which does cache-intelligent transposition of numpy arrays (using some block strategy), this might be faster even than the current simple reduction operations when forced to work on awkwardly aligned data. Ideally, this intelligent access and intelligent reduction would be part of a single pass of course; but that wouldn't really fit within the numpy design, and merely such an intelligent transpose would provide most of the benefit I think. Or is the mechanism behind ascontiguousarray already intelligent in this sense? The iterator is currently smart in the sense that it will (obviously low level), do something like [1]. Most things in numpy use that iterator, ascontiguousarray does so as well. Such a blocked cache aware iterator is what I mean by, someone who knows it would have to spend a lot of time on it :). [1] Appendix: arr = np.ones((100, 100)) arr.sum(1) # being equivalent (iteration order wise) to: res = np.zeros(100) for i in range(100): res += arr[i, :] # while arr.sum(0) would be: for i in range(100): res[i] = arr[i, :].sum() On Mon, Jul 28, 2014 at 4:06 PM, Sebastian Berg sebast...@sipsolutions.net wrote: On Mo, 2014-07-28 at 15:35 +0200, Sturla Molden wrote: On 28/07/14 15:21, alex wrote: Are you sure they always give different results? Notice that np.ones((N,2)).mean(0) np.ones((2,N)).mean(1) compute means of different axes on transposed arrays so these differences 'cancel out'. They will be if different algorithms are used. np.ones((N,2)).mean(0) will have larger accumulated rounding error than np.ones((2,N)).mean(1), if only the latter uses the divide-and-conquer summation. What I wanted to point out is that to some extend the algorithm does not matter. You will not necessarily get identical results already if you use a different iteration order, and we have been doing that for years for speed reasons. All libs like BLAS do the same. Yes, the new changes make this much more dramatic, but they only make some paths much better, never worse. It might be dangerous, but only in the sense that you test it with the good path and it works good enough, but later (also) use the other one in some lib. I am not even sure if I I would suggest that in the first case we try to copy the array to a temporary contiguous buffer and use the same divide-and-conquer algorithm, unless some heuristics on memory usage fails. Sure, but you have to make major changes to the buffered iterator to do that without larger speed implications. It might be a good idea, but it requires someone who knows this stuff to spend a lot of time and care in the depths of numpy. Sturla
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
Cool, sounds like great improvements. I can imagine that after some loop unrolling one becomes memory bound pretty soon. Is the summation guaranteed to traverse the data in its natural order? And do you happen to know what the rules for choosing accumulator dtypes are? -Original Message- From: Julian Taylor jtaylor.deb...@googlemail.com Sent: 26-7-2014 00:58 To: Discussion of Numerical Python numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays On 25.07.2014 23:51, Eelco Hoogendoorn wrote: Ray: I'm not working with Hubble data, but yeah these are all issues I've run into with my terrabytes of microscopy data as well. Given that such raw data comes as uint16, its best to do your calculations as much as possible in good old ints. What you compute is what you get, no obscure shenanigans. integers are dangerous too, they overflow quickly and signed overflow is even undefined in C the standard. It just occurred to me that pairwise summation will lead to highly branchy code, and you can forget about any vector extensions. Tradeoffs indeed. Any such hierarchical summation is probably best done by aggregating naively summed blocks. pairwise summation is usually implemented with a naive sum cutoff large enough so the recursion does not matter much. In numpy 1.9 this cutoff is 128 elements, but the inner loop is unrolled 8 times which makes it effectively 16 elements. the unrolling factor of 8 was intentionally chosen to allow using AVX in the inner loop without changing the summation ordering, but last I tested actually using AVX here only gave mediocre speedups (10%-20% on an i5). From: RayS mailto:r...@blue-cove.com Sent: 25-7-2014 23:26 To: Discussion of Numerical Python mailto:numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays At 11:29 AM 7/25/2014, you wrote: On Fri, Jul 25, 2014 at 5:56 PM, RayS r...@blue-cove.com wrote: The important point was that it would be best if all of the methods affected by summing 32 bit floats with 32 bit accumulators had the same Notes as numpy.mean(). We went through a lot of code yesterday, assuming that any numpy or Scipy.stats functions that use accumulators suffer the same issue, whether noted or not, and found it true. Do you have a list of the functions that are affected? We only tested a few we used, but scipy.stats.nanmean, or any .*mean() numpy.sum, mean, average, std, var,... via something like: import numpy import scipy.stats print numpy.__version__ print scipy.__version__ onez = numpy.ones((2**25, 1), numpy.float32) step = 2**10 func = scipy.stats.nanmean for s in range(2**24-step, 2**25, step): if func(onez[:s+step])!=1.: print '\nbroke', s, func(onez[:s+step]) break else: print '\r',s, That said, it does seem that np.mean could be implemented better than it is, even given float32's inherent limitations. If anyone wants to implement better algorithms for computing the mean, variance, sums, etc., then we would love to add them to numpy. Others have pointed out the possible tradeoffs in summation algos, perhaps a method arg would be appropriate, better depending on your desire for speed vs. accuracy. It just occurred to me that if the STSI folks (who count photons) took the mean() or other such func of an image array from Hubble sensors to find background value, they'd better always be using float64. - Ray ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
I was wondering the same thing. Are there any known tradeoffs to this method of reduction? On Sat, Jul 26, 2014 at 12:39 PM, Sturla Molden sturla.mol...@gmail.com wrote: Sebastian Berg sebast...@sipsolutions.net wrote: chose more stable algorithms for such statistical functions. The pairwise summation that is in master now is very awesome, but it is not secure enough in the sense that a new user will have difficulty understanding when he can be sure it is used. Why is it not always used? ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
A context manager makes sense. I very much appreciate the time constraints and the effort put in this far, but if we can not make something work uniformly, I wonder if we should include it in the master at all. I don't have a problem with customizing algorithms where fp accuracy demands it; I have more of a problem with hard to predict behavior. If np.ones(bigN).sum() gives different results than np.ones((bigN,2)).sum(0), aside from the obvious differences, that would be one hard to catch source of bugs. Wouldn't per-axis reduction, as a limited form of nested reduction, provide most of the benefits, without any of the drawbacks? On Sat, Jul 26, 2014 at 3:53 PM, Julian Taylor jtaylor.deb...@googlemail.com wrote: On 26.07.2014 15:38, Eelco Hoogendoorn wrote: Why is it not always used? for 1d reduction the iterator blocks by 8192 elements even when no buffering is required. There is a TODO in the source to fix that by adding additional checks. Unfortunately nobody knows hat these additional tests would need to be and Mark Wiebe who wrote it did not reply to a ping yet. Also along the non-fast axes the iterator optimizes the reduction to remove the strided access, see: https://github.com/numpy/numpy/pull/4697#issuecomment-42752599 Instead of having a keyword argument to mean I would prefer a context manager that changes algorithms for different requirements. This would easily allow changing the accuracy and performance of third party functions using numpy without changing the third party library as long as they are using numpy as the base. E.g. with np.precisionstate(sum=kahan): scipy.stats.nanmean(d) We also have case where numpy uses algorithms that are far more precise than most people needs them. E.g. np.hypot and the related complex absolute value and division. These are very slow with glibc as it provides 1ulp accuracy, this is hardly ever needed. Another case that could use dynamic changing is flushing subnormals to zero. But this api is like Nathaniels parameterizable dtypes just an idea floating in my head which needs proper design and implementation written down. The issue is as usual ENOTIME. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
Perhaps I in turn am missing something; but I would suppose that any algorithm that requires multiple passes over the data is off the table? Perhaps I am being a little old fashioned and performance oriented here, but to make the ultra-majority of use cases suffer a factor two performance penalty for an odd use case which already has a plethora of fine and dandy solutions? Id vote against, fwiw... On Sat, Jul 26, 2014 at 6:34 PM, Sturla Molden sturla.mol...@gmail.com wrote: Sturla Molden sturla.mol...@gmail.com wrote: Sebastian Berg sebast...@sipsolutions.net wrote: Yes, it is much more complicated and incompatible with naive ufuncs if you want your memory access to be optimized. And optimizing that is very much worth it speed wise... Why? Couldn't we just copy the data chunk-wise to a temporary buffer of say 2**13 numbers and then reduce that? I don't see why we need another iterator for that. I am sorry if this is a stupid suggestion. My knowledge of how NumPy ufuncs works could have been better. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
To elaborate on that point; knowing that numpy accumulates in a simple first-to-last sweep, and does not implicitly upcast, the original problem can be solved in several ways; specifying a higher precision to sum with, or by a nested summation, like X.mean(0).mean(0)==1.0. I personally like this explicitness, and am wary of numpy doing overly clever things behind the scenes, as I can think of other code that might become broken if things change too radically. For instance, I often sort large arrays with a large spread in magnitudes before summation, relying on the fact that summing the smallest values first gives best precision. Any changes made to reduction behavior should try and be backwards compatible with such properties of straightforward reductions, or else a lot of code is going to be broken without warning. I suppose using maximum precision internally, and nesting all reductions over multiple axes of an ndarray, are both easy to implement improvements that do not come with any drawbacks that I can think of. Actually the maximum precision I am not so sure of, as I personally prefer to make an informed decision about precision used, and get an error on a platform that does not support the specified precision, rather than obtain subtly or horribly broken results without warning when moving your code to a different platform/compiler whatever. On Fri, Jul 25, 2014 at 5:37 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Perhaps it is a slightly semantical discussion; but all fp calculations have errors, and there are always strategies for making them smaller. We just don't happen to like the error for this case; but rest assured it won't be hard to find new cases of 'blatantly wrong' results, no matter what accumulator is implemented. That's no reason to not try and be clever about it, but there isn't going to be an algorithm that is best for all possible inputs, and in the end the most important thing is that the algorithm used is specified in the docs. -- From: Alan G Isaac alan.is...@gmail.com Sent: 25-7-2014 00:10 To: Discussion of Numerical Python numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays On 7/24/2014 4:42 PM, Eelco Hoogendoorn wrote: This isn't a bug report, but rather a feature request. I'm not sure statement this is correct. The mean of a float32 array can certainly be computed as a float32. Currently this is not necessarily what happens, not even approximately. That feels a lot like a bug, even if we can readily understand how the algorithm currently used produces it. To say whether it is a bug or not, don't we have to ask about the intent of `mean`? If the intent is to sum and divide, then it is not a bug. If the intent is to produce the mean, then it is a bug. Alan Isaac ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
Arguably, the whole of floating point numbers and their related shenanigans is not very pythonic in the first place. The accuracy of the output WILL depend on the input, to some degree or another. At the risk of repeating myself: explicit is better than implicit -Original Message- From: RayS r...@blue-cove.com Sent: 25-7-2014 19:56 To: Discussion of Numerical Python numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays At 07:22 AM 7/25/2014, you wrote: We were talking on this in the office, as we realized it does affect a couple of lines dealing with large arrays, including complex64. As I expect Python modules to work uniformly cross platform unless documented otherwise, to me that includes 32 vs 64 bit platforms, implying that the modules should automatically use large enough accumulators for the data type input. The 32/64-bitness of your platform has nothing to do with floating point. As a naive end user, I can, and do, download different binaries for different CPUs/Windows versions and will get different results http://mail.scipy.org/pipermail/numpy-discussion/2014-July/070747.html Nothing discussed in this thread is platform-specific (modulo some minor details about the hardware FPU, but that should be taken as read). And compilers, apparently. The important point was that it would be best if all of the methods affected by summing 32 bit floats with 32 bit accumulators had the same Notes as numpy.mean(). We went through a lot of code yesterday, assuming that any numpy or Scipy.stats functions that use accumulators suffer the same issue, whether noted or not, and found it true. Depending on the input data, this can cause the results to be inaccurate, especially for float32 (see example below). Specifying a higher-precision accumulator using the dtype keyword can alleviate this issue. seems rather un-Pythonic. - Ray___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
It need not be exactly representable as such; take the mean of [1, 1+eps] for instance. Granted, there are at most two number in the range of the original dtype which are closest to the true mean; but im not sure that computing them exactly is a tractable problem for arbitrary input. Im not sure what is considered best practice for these problems; or if there is one, considering the hetrogenity of the problem. As noted earlier, summing a list of floating point values is a remarkably multifaceted problem, once you get down into the details. I think it should be understood that all floating point algorithms are subject to floating point errors. As long as the algorithm used is specified, one can make an informed decision if the given algorithm will do what you expect of it. That's the best we can hope for. If we are going to advocate doing 'clever' things behind the scenes, we have to take backwards compatibility (not creating a possibility of producing worse results on the same input) and platform independence in mind. Funny summation orders could violate the former depending on the implementation details, and 'using the highest machine precision available' violates the latter (and is horrible practice in general, imo. Either you don't need the extra accuracy, or you do, and the absence on a given platform should be an error) Perhaps pairwise summation in the original order of the data is the best option: q = np.ones((2,)*26, np.float32) print q.mean() while q.ndim 0: q = q.mean(axis=-1, dtype=np.float32) print q This only requires log(N) space on the stack if properly implemented, and is not platform dependent, nor should have any backward compatibility issues that I can think of. But im not sure how easy it would be to implement, given the current framework. The ability to specify different algorithms per kwarg wouldn't be a bad idea either, imo; or the ability to explicitly specify a separate output and accumulator dtype. On Fri, Jul 25, 2014 at 8:00 PM, Alan G Isaac alan.is...@gmail.com wrote: On 7/25/2014 1:40 PM, Eelco Hoogendoorn wrote: At the risk of repeating myself: explicit is better than implicit This sounds like an argument for renaming the `mean` function `naivemean` rather than `mean`. Whatever numpy names `mean`, shouldn't it implement an algorithm that produces the mean? And obviously, for any float data type, the mean value of the values in the array is representable as a value of the same type. Alan Isaac ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
Ray: I'm not working with Hubble data, but yeah these are all issues I've run into with my terrabytes of microscopy data as well. Given that such raw data comes as uint16, its best to do your calculations as much as possible in good old ints. What you compute is what you get, no obscure shenanigans. It just occurred to me that pairwise summation will lead to highly branchy code, and you can forget about any vector extensions. Tradeoffs indeed. Any such hierarchical summation is probably best done by aggregating naively summed blocks. -Original Message- From: RayS r...@blue-cove.com Sent: 25-7-2014 23:26 To: Discussion of Numerical Python numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays At 11:29 AM 7/25/2014, you wrote: On Fri, Jul 25, 2014 at 5:56 PM, RayS r...@blue-cove.com wrote: The important point was that it would be best if all of the methods affected by summing 32 bit floats with 32 bit accumulators had the same Notes as numpy.mean(). We went through a lot of code yesterday, assuming that any numpy or Scipy.stats functions that use accumulators suffer the same issue, whether noted or not, and found it true. Do you have a list of the functions that are affected? We only tested a few we used, but scipy.stats.nanmean, or any .*mean() numpy.sum, mean, average, std, var,... via something like: import numpy import scipy.stats print numpy.__version__ print scipy.__version__ onez = numpy.ones((2**25, 1), numpy.float32) step = 2**10 func = scipy.stats.nanmean for s in range(2**24-step, 2**25, step): if func(onez[:s+step])!=1.: print '\nbroke', s, func(onez[:s+step]) break else: print '\r',s, That said, it does seem that np.mean could be implemented better than it is, even given float32's inherent limitations. If anyone wants to implement better algorithms for computing the mean, variance, sums, etc., then we would love to add them to numpy. Others have pointed out the possible tradeoffs in summation algos, perhaps a method arg would be appropriate, better depending on your desire for speed vs. accuracy. It just occurred to me that if the STSI folks (who count photons) took the mean() or other such func of an image array from Hubble sensors to find background value, they'd better always be using float64. - Ray ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for large float32 arrays
Arguably, this isn't a problem of numpy, but of programmers being trained to think of floating point numbers as 'real' numbers, rather than just a finite number of states with a funny distribution over the number line. np.mean isn't broken; your understanding of floating point number is. What you appear to wish for is a silent upcasting of the accumulated result. This is often performed in reducing operations, but I can imagine it runs into trouble for nd-arrays. After all, if I have a huge array that I want to reduce over a very short axis, upcasting might be very undesirable; it wouldn't buy me any extra precision, but it would increase memory use from 'huge' to 'even more huge'. np.mean has a kwarg that allows you to explicitly choose the dtype of the accumulant. X.mean(dtype=np.float64)==1.0. Personally, I have a distaste for implicit behavior, unless the rule is simple and there really can be no negative downsides; which doesn't apply here I would argue. Perhaps when reducing an array completely to a single value, there is no harm in upcasting to the maximum machine precision; but that becomes a rather complex rule which would work out differently for different machines. Its better to be confronted with the limitations of floating point numbers earlier, rather than later when you want to distribute your work and run into subtle bugs on other peoples computers. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for large float32arrays
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.is...@gmail.com Sent: 24-7-2014 18:09 To: Discussion of Numerical Python numpy-discussion@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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for large float32arrays
Inaccurate and utterly wrong are subjective. If You want To Be sufficiently strict, floating point calculations are almost always 'utterly wrong'. Granted, It would Be Nice if the docs specified the algorithm used. But numpy does not produce anything different than what a standard c loop or c++ std lib func would. This isn't a bug report, but rather a feature request. That said, support for fancy reduction algorithms would certainly be nice, if implementing it in numpy in a coherent manner is feasible. -Original Message- From: Joseph Martinot-Lagarde joseph.martinot-laga...@m4x.org Sent: 24-7-2014 20:04 To: numpy-discussion@scipy.org numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] numpy.mean still broken for large float32arrays Le 24/07/2014 12:55, Thomas Unterthiner a écrit : I don't agree. The problem is that I expect `mean` to do something reasonable. The documentation mentions that the results can be inaccurate, which is a huge understatement: the results can be utterly wrong. That is not reasonable. At the very least, a warning should be issued in cases where the dtype might not be appropriate. Maybe the problem is the documentation, then. If this is a common error, it could be explicitly documented in the function documentation. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays
Perhaps it is a slightly semantical discussion; but all fp calculations have errors, and there are always strategies for making them smaller. We just don't happen to like the error for this case; but rest assured it won't be hard to find new cases of 'blatantly wrong' results, no matter what accumulator is implemented. That's no reason to not try and be clever about it, but there isn't going to be an algorithm that is best for all possible inputs, and in the end the most important thing is that the algorithm used is specified in the docs. -Original Message- From: Alan G Isaac alan.is...@gmail.com Sent: 25-7-2014 00:10 To: Discussion of Numerical Python numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] numpy.mean still broken for largefloat32arrays On 7/24/2014 4:42 PM, Eelco Hoogendoorn wrote: This isn't a bug report, but rather a feature request. I'm not sure statement this is correct. The mean of a float32 array can certainly be computed as a float32. Currently this is not necessarily what happens, not even approximately. That feels a lot like a bug, even if we can readily understand how the algorithm currently used produces it. To say whether it is a bug or not, don't we have to ask about the intent of `mean`? If the intent is to sum and divide, then it is not a bug. If the intent is to produce the mean, then it is a bug. Alan Isaac ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Find n closest values
Well, if the spacing is truly uniform, then of course you don't really need the search, and you can do away with the extra log-n, and there is a purely linear solution: def find_closest_direct(start, end, count, A): Q = (A-start)/(end-start)*count mid = ((Q[1:]+Q[:-1]+1)/2).astype(np.int) boundary = np.zeros(count, np.int) boundary[mid] = 1 return np.add.accumulate(boundary) I expect this to be a bit faster, but nothing dramatic, unless your datasets are huge. It isn't really more or less elegant either, id say. Note that the output isn't 100% identical; youd need to do a little tinkering to figure out the correct/desired rounding behavior. On Sun, Jun 22, 2014 at 5:16 PM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Thanks for the answer. I was secretly hoping for some kind of hardly-known numpy function that would make things faster auto-magically... Nicolas On 22 Jun 2014, at 10:30, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Perhaps you could simplify some statements, but at least the algorithmic complexity is fine, and everything is vectorized, so I doubt you will get huge gains. You could take a look at the functions in scipy.spatial, and see how they perform for your problem parameters. On Sun, Jun 22, 2014 at 10:22 AM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Hi, I have an array L with regular spaced values between 0 and width. I have a (sorted) array I with irregular spaced values between 0 and width. I would like to find the closest value in I for any value in L. Currently, I'm using the following script but I wonder if I missed an obvious (and faster) solution: import numpy as np def find_closest(A, target): idx = A.searchsorted(target) idx = np.clip(idx, 1, len(A) - 1) left = A[idx - 1] right = A[idx] idx -= target - left right - target return idx n, width = 256, 100.0 # 10 random sorted values in [0,width] I = np.sort(np.random.randint(0,width,10)) # n regular spaced values in [0,width] L = np.linspace(0, width, n) print I[find_closest(I,L)] Nicolas ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Find n closest values
Also, if you use scipy.spatial.KDTree, make sure to use cKDTree; the native python kdtree is sure to be slow as hell. On Sun, Jun 22, 2014 at 7:05 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Well, if the spacing is truly uniform, then of course you don't really need the search, and you can do away with the extra log-n, and there is a purely linear solution: def find_closest_direct(start, end, count, A): Q = (A-start)/(end-start)*count mid = ((Q[1:]+Q[:-1]+1)/2).astype(np.int) boundary = np.zeros(count, np.int) boundary[mid] = 1 return np.add.accumulate(boundary) I expect this to be a bit faster, but nothing dramatic, unless your datasets are huge. It isn't really more or less elegant either, id say. Note that the output isn't 100% identical; youd need to do a little tinkering to figure out the correct/desired rounding behavior. On Sun, Jun 22, 2014 at 5:16 PM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Thanks for the answer. I was secretly hoping for some kind of hardly-known numpy function that would make things faster auto-magically... Nicolas On 22 Jun 2014, at 10:30, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Perhaps you could simplify some statements, but at least the algorithmic complexity is fine, and everything is vectorized, so I doubt you will get huge gains. You could take a look at the functions in scipy.spatial, and see how they perform for your problem parameters. On Sun, Jun 22, 2014 at 10:22 AM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Hi, I have an array L with regular spaced values between 0 and width. I have a (sorted) array I with irregular spaced values between 0 and width. I would like to find the closest value in I for any value in L. Currently, I'm using the following script but I wonder if I missed an obvious (and faster) solution: import numpy as np def find_closest(A, target): idx = A.searchsorted(target) idx = np.clip(idx, 1, len(A) - 1) left = A[idx - 1] right = A[idx] idx -= target - left right - target return idx n, width = 256, 100.0 # 10 random sorted values in [0,width] I = np.sort(np.random.randint(0,width,10)) # n regular spaced values in [0,width] L = np.linspace(0, width, n) print I[find_closest(I,L)] Nicolas ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Find n closest values
Protip: if you are writing your own rasterization code in python, be prepared to forget about performance altogether. Something like numba or other c-like extension will be necessary unless you are willing to leave big gobs of performance on the table; and even with pure C you will get nowhere close to the performance of super-duper optimized library code you are used to. But before you go down that rabbit hole, its probably worth thinking about whether you can get an existing rendering framework to do what you want to do. On Sun, Jun 22, 2014 at 8:30 PM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Thanks, I'll try your solution. Data (L) is not so big actually, it represents pixels on screen and (I) represents line position (for grids). I need to compute this quantity everytime the user zoom in or out. Nicolas On 22 Jun 2014, at 19:05, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Well, if the spacing is truly uniform, then of course you don't really need the search, and you can do away with the extra log-n, and there is a purely linear solution: def find_closest_direct(start, end, count, A): Q = (A-start)/(end-start)*count mid = ((Q[1:]+Q[:-1]+1)/2).astype(np.int) boundary = np.zeros(count, np.int) boundary[mid] = 1 return np.add.accumulate(boundary) I expect this to be a bit faster, but nothing dramatic, unless your datasets are huge. It isn't really more or less elegant either, id say. Note that the output isn't 100% identical; youd need to do a little tinkering to figure out the correct/desired rounding behavior. On Sun, Jun 22, 2014 at 5:16 PM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Thanks for the answer. I was secretly hoping for some kind of hardly-known numpy function that would make things faster auto-magically... Nicolas On 22 Jun 2014, at 10:30, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Perhaps you could simplify some statements, but at least the algorithmic complexity is fine, and everything is vectorized, so I doubt you will get huge gains. You could take a look at the functions in scipy.spatial, and see how they perform for your problem parameters. On Sun, Jun 22, 2014 at 10:22 AM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Hi, I have an array L with regular spaced values between 0 and width. I have a (sorted) array I with irregular spaced values between 0 and width. I would like to find the closest value in I for any value in L. Currently, I'm using the following script but I wonder if I missed an obvious (and faster) solution: import numpy as np def find_closest(A, target): idx = A.searchsorted(target) idx = np.clip(idx, 1, len(A) - 1) left = A[idx - 1] right = A[idx] idx -= target - left right - target return idx n, width = 256, 100.0 # 10 random sorted values in [0,width] I = np.sort(np.random.randint(0,width,10)) # n regular spaced values in [0,width] L = np.linspace(0, width, n) print I[find_closest(I,L)] Nicolas ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Find n closest values
That's pretty cool; and it makes sense that way. Still, couldn't you fold this kind of computation into a shader? Have you looked at vispy btw? I think its a really nice initiative; having a high quality vector graphics module in there would make it even better. Would be nice if those projects could be merged. On Sun, Jun 22, 2014 at 9:51 PM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Actually, it's already working pretty well but it slows down when you're doing a lot of zoom in/out. The trick is that rendering is done using shader (OpenGL) and this computation is used to give information to the shader to where to draw antialiased lines. In the end, this shader is able to draw any amiunt of grids/ticks (as in matplotlib). Some old example are available from here: https://github.com/rougier/gl-agg I tested your solution and it is faster by only a tiny amount but the way you wrote it might open the door for other improvements. Thanks. Nicolas On 22 Jun 2014, at 21:14, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Protip: if you are writing your own rasterization code in python, be prepared to forget about performance altogether. Something like numba or other c-like extension will be necessary unless you are willing to leave big gobs of performance on the table; and even with pure C you will get nowhere close to the performance of super-duper optimized library code you are used to. But before you go down that rabbit hole, its probably worth thinking about whether you can get an existing rendering framework to do what you want to do. On Sun, Jun 22, 2014 at 8:30 PM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Thanks, I'll try your solution. Data (L) is not so big actually, it represents pixels on screen and (I) represents line position (for grids). I need to compute this quantity everytime the user zoom in or out. Nicolas On 22 Jun 2014, at 19:05, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Well, if the spacing is truly uniform, then of course you don't really need the search, and you can do away with the extra log-n, and there is a purely linear solution: def find_closest_direct(start, end, count, A): Q = (A-start)/(end-start)*count mid = ((Q[1:]+Q[:-1]+1)/2).astype(np.int) boundary = np.zeros(count, np.int) boundary[mid] = 1 return np.add.accumulate(boundary) I expect this to be a bit faster, but nothing dramatic, unless your datasets are huge. It isn't really more or less elegant either, id say. Note that the output isn't 100% identical; youd need to do a little tinkering to figure out the correct/desired rounding behavior. On Sun, Jun 22, 2014 at 5:16 PM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Thanks for the answer. I was secretly hoping for some kind of hardly-known numpy function that would make things faster auto-magically... Nicolas On 22 Jun 2014, at 10:30, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Perhaps you could simplify some statements, but at least the algorithmic complexity is fine, and everything is vectorized, so I doubt you will get huge gains. You could take a look at the functions in scipy.spatial, and see how they perform for your problem parameters. On Sun, Jun 22, 2014 at 10:22 AM, Nicolas P. Rougier nicolas.roug...@inria.fr wrote: Hi, I have an array L with regular spaced values between 0 and width. I have a (sorted) array I with irregular spaced values between 0 and width. I would like to find the closest value in I for any value in L. Currently, I'm using the following script but I wonder if I missed an obvious (and faster) solution: import numpy as np def find_closest(A, target): idx = A.searchsorted(target) idx = np.clip(idx, 1, len(A) - 1) left = A[idx - 1] right = A[idx] idx -= target - left right - target return idx n, width = 256, 100.0 # 10 random sorted values in [0,width] I = np.sort(np.random.randint(0,width,10)) # n regular spaced values in [0,width] L = np.linspace(0, width, n) print I[find_closest(I,L)] Nicolas ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org
Re: [Numpy-discussion] Easter Egg or what I am missing here?
I agree; this 'wart' has also messed with my code a few times. I didn't find it to be the case two years ago, but perhaps I should reevaluate if the scientific python stack has sufficiently migrated to python 3. On Thu, May 22, 2014 at 7:35 AM, Siegfried Gonzi siegfried.go...@ed.ac.ukwrote: On 22/05/2014 00:37, numpy-discussion-requ...@scipy.org wrote: Message: 4 Date: Wed, 21 May 2014 18:32:30 -0400 From: Warren Weckesser warren.weckes...@gmail.com Subject: Re: [Numpy-discussion] Easter Egg or what I am missing here? To: Discussion of Numerical Python numpy-discussion@scipy.org Message-ID: cagzf1udkdap+yd2sqy9cca6rm4zzfyjjzfv0tieesv_driv...@mail.gmail.com Content-Type: text/plain; charset=UTF-8 On 5/21/14, Siegfried Gonzi siegfried.go...@ed.ac.uk wrote: Please would anyone tell me the following is an undocumented bug otherwise I will lose faith in everything: == import numpy as np years = [2004,2005,2006,2007] dates = [20040501,20050601,20060801,20071001] for x in years: print 'year ',x xy = np.array([x*1.0e-4 for x in dates]).astype(np.int) print 'year ',x == Or is this a recipe to blow up a power plant? This is a wart of Python 2.x. The dummy variable used in a list comprehension remains defined with its final value in the enclosing scope. For example, this is Python 2.7: x = 100 w = [x*x for x in range(4)] x 3 This behavior has been changed in Python 3. Here's the same sequence in Python 3.4: x = 100 w = [x*x for x in range(4)] x 100 Guido van Rossum gives a summary of this issue near the end of this blog: http://python-history.blogspot.com/2010/06/from-list-comprehensions-to-generator.html Warren [I still do not know how to properly use the reply function here. I apologise.] Hi all and thanks to all the respondes. I think I would have expected my code to be behaving like you said version 3.4 will do. I would never have thought 'x' is being changed during execution. I took me nearly 2 hours in my code to figure out what was going on (it was a lenghty piece of code an not so easy to spot). Siegfried -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] repeat an array without allocation
If b is indeed big I don't see a problem with the python loop, elegance aside; but Cython will not beat it on that front. On Mon, May 5, 2014 at 9:34 AM, srean srean.l...@gmail.com wrote: Great ! thanks. I should have seen that. Is there any way array multiplication (as opposed to matrix multiplication) can be sped up without forming A and (A * b) explicitly. A = np.repeat(x, [4, 2, 1, 3], axis = 0)# A.shape == 10,10 c = sum(b * A, axis = 1)# b.shape == 10,10 In my actual setting b is pretty big, so I would like to avoid creating another array the same size. I would also like to avoid a Python loop. st = 0 for (i,rep) in enumerate([4, 2, 1, 3]): end = st + rep c[st : end] = np.dot(b[st : end, :], a[i,:]) st = end Is Cython the only way ? On Mon, May 5, 2014 at 1:20 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Sun, May 4, 2014 at 9:34 PM, srean srean.l...@gmail.com wrote: Hi all, is there an efficient way to do the following without allocating A where A = np.repeat(x, [4, 2, 1, 3], axis=0) c = A.dot(b)# b.shape If x is a 2D array you can call repeat **after** dot, not before, which will save you some memory and a few operations: a = np.random.rand(4, 5) b = np.random.rand(5, 6) np.allclose(np.repeat(a, [4, 2, 1, 3], axis=0).dot(b), ... np.repeat(a.dot(b), [4, 2, 1, 3], axis=0)) True Similarly, if x is a 1D array, you can sum the corresponding items of b before calling dot: a = np.random.rand(4) b = np.random.rand(10) idx = np.concatenate(([0], np.cumsum([4,2,1,3])[:-1])) np.allclose(np.dot(np.repeat(a, [4,2,1,3] ,axis=0), b), ... np.dot(a, np.add.reduceat(b, idx))) ... ) True Jaime -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] repeat an array without allocation
nope; its impossible to express A as a strided view on x, for the repeats you have. even if you had uniform repeats, it still would not work. that would make it easy to add an extra axis to x without a new allocation; but reshaping/merging that axis with axis=0 would again trigger a copy, as it would require a non-integer stride. On Mon, May 5, 2014 at 6:34 AM, srean srean.l...@gmail.com wrote: Hi all, is there an efficient way to do the following without allocating A where A = np.repeat(x, [4, 2, 1, 3], axis=0) c = A.dot(b)# b.shape thanks -- srean ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] arrays and : behaviour
You problem isn't with colon indexing, but with the interpretation of the arguments to plot. multiple calls to plot with scalar arguments do not have the same result as a single call with array arguments. For this to work as intended, you would need plt.hold(True), for starters, and maybe there are other subtleties. On Thu, May 1, 2014 at 1:31 PM, did did 21di...@gmx.com wrote: Hello all and sorry for my bad english, i am a beginner with python and i try to save a lot of data in several folders in a 4D matrix and then to plot two columns of this 4D matrix. Bellow, i have the code to fill my 4D matrix, it works very well : [CODE]matrix4D=[] for i in Numbers: readInFolder=folderPrefixe+i+/ matrix3D=[] for j in listeOfdata: nameOfFile=filePrefixe+i+-+j+extensionTXT nameOfFile=readInFolder+nameOfFile matrix2D=np.loadtxt(nameOfFile,delimiter=,,skiprows=1) matrix3D.append(matrix2D) matrix4D.append(matrix3D) array4D = np.asarray(matrix4D)[/CODE] But now, i want to plot the third column as function of the third too (just for trying) and i use this stupid manner that works well : [CODE]plt.figure(1) temp=plt.plot(array4D[0][0][0][0],array4D[0][0][0][0],'bo') temp=plt.plot(array4D[0][0][1][0],array4D[0][0][1][0],'bo') temp=plt.plot(array4D[0][0][2][0],array4D[0][0][2][0],'bo') temp=plt.plot(array4D[0][0][3][0],array4D[0][0][3][0],'bo') plt.show()[/CODE] Now, i want to use a more smart manner and i use : like this [CODE]plt.figure(1) temp=plt.plot(array4D[0][0][0:3][0],array4D[0][0][0:3][0],'bo') plt.show()[/CODE] The result should be the same but i don't got the same results!!! In attachement you have the two corresponding plots, can you explain to me with i don't have the same plots ?? thanks for all ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numerical gradient, Jacobian, and Hessian
I was going to suggest numdifftools; its a very capable package in my experience. Indeed it would be nice to have it integrated into scipy. Also, in case trying to calculate a numerical gradient is a case of 'the math getting too bothersome' rather than no closed form gradient actually existing: Theano may be your best bet; I have very good experiences with it as well. As far as I can tell, it is actually the only tensor/ndarray aware differentiator out there (maple and mathematica don't appear to support this) On Sun, Apr 20, 2014 at 4:55 PM, Alan G Isaac alan.is...@gmail.com wrote: Awhile back there were good signs that SciPy would end up with a `diff` module: https://github.com/scipy/scipy/issues/2035 Is this still moving forward? It would certainly be nice for SciPy to have intuitive numerical gradients, Jacobians, and Hessians. The last two are I think missing altogether. The first exists as scipy.optimize.approx_fprime. `approx_fprime` seems to work fine, but I suggest it has the following drawbacks: - it is hard to find (e.g., try doing a Google search on scipy gradient or scipy numerical gradient - related, it is in the wrong location (scipy.optimize) - the signature is odd: (x,f,dx) instead of (f,x,dx) (This matters for ease of recall and for teaching.) In any case, as I understand it, the author's of numdifftools http://code.google.com/p/numdifftools/ expressed willingness to have their code moved into SciPy. This seems like an excellent way forward. There was talk of making this a summer of code project, but that seems to have sputtered. Alan Isaac ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] string replace
Indeed this isn't numpy, and I don't see how your collegues opinions have
bearing on that issue; but anyway..
There isn't a 'python' way to do this, the best method involves some form
of parsing library. Undoubtly there is a one-line regex to do this kind of
thing, but regexes are themselves the antithesis of python.
Here is how id do it, using pyparsing:
from pyparsing import *
line_left = './erfo/restart.ST010.EN0001-EN0090.MMDDhh'
n1 = Word(nums,exact=3).setParseAction(replaceWith('000'))
n2 = Word(nums,exact=4).setParseAction(replaceWith('0001'))
n3 = Word(nums,exact=4).setParseAction(replaceWith('0092'))
pattern = Literal('ST')+ n1 +'.EN'+n2+'-EN'+n3
print pattern.transformString(line_left)
On Sun, Apr 20, 2014 at 5:25 PM, Siegfried Gonzi
siegfried.go...@ed.ac.ukwrote:
Hi all
I know this is not numpy related but a colleague insists the following
is supposed to work. But it doesn't:
==
line_left = './erfo/restart.ST010.EN0001-EN0090.MMDDhh'
enafix = 'ST000.EN0001-EN0092'
line_left = line_left.replace('STYYY.EN-EN', enafix)
print 'line_left',line_left
==
[right answer would be: ./erfo/restart.ST000.EN0001-EN0092.MMDDhh' ]
I think it works in Fortran but not in Python. What would be the easiest
way to replacing this kind of pattern in a variable lenght string?
Thanks,
Siegfried
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] min depth to nonzero in 3d array
I agree; argmax would the best option here; though I would hardly call it abuse. It seems perfectly readable and idiomatic to me. Though the != comparison requires an extra pass over the array, that's the kind of tradeoff you make in using numpy. On Thu, Apr 17, 2014 at 7:45 PM, Stephan Hoyer sho...@gmail.com wrote: Hi Alan, You can abuse np.argmax to calculate the first nonzero element in a vectorized manner: import numpy as np A = (2 * np.random.rand(100, 50, 50)).astype(int) Compare: np.argmax(A != 0, axis=0) np.array([[np.flatnonzero(A[:,i,j])[0] for j in range(50)] for i in range(50)]) You'll also want to check for all zero arrays with np.all: np.all(A == 0, axis=0) Cheers, Stephan On Thu, Apr 17, 2014 at 9:32 AM, Alan G Isaac alan.is...@gmail.comwrote: Given an array A of shape m x n x n (i.e., a stack of square matrices), I want an n x n array that gives the minimum depth to a nonzero element. E.g., the 0,0 element of the result is np.flatnonzero(A[:,0,0])[0] Can this be vectorized? (Assuming a nonzero element exists is ok, but dealing nicely with its absence is even better.) Thanks, Alan Isaac ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Wiki page for building numerical stuff on Windows
I wonder: how hard would it be to create a more 21th-century oriented BLAS, relying more on code generation tools, and perhaps LLVM/JITting? Wouldn't we get ten times the portability with one-tenth the lines of code? Or is there too much dark magic going on in BLAS for such an approach to come close enough to hand-tuned performance? On Sat, Apr 12, 2014 at 12:15 AM, Sturla Molden sturla.mol...@gmail.comwrote: On 11/04/14 04:44, Matthew Brett wrote: I've been working on a general wiki page on building numerical stuff on Windows: https://github.com/numpy/numpy/wiki/Numerical-software-on-Windows I am worried that the conclusion will be that there is no viable BLAS alternative on Windows... Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Wiki page for building numerical stuff onWindows
BLIS seems like a nice project as well. I like the arbitrary striding; BLAS lacking this has always annoyed me. -Original Message- From: Sturla Molden sturla.mol...@gmail.com Sent: 12-4-2014 13:12 To: numpy-discussion@scipy.org numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] Wiki page for building numerical stuff onWindows Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: I wonder: how hard would it be to create a more 21th-century oriented BLAS, relying more on code generation tools, and perhaps LLVM/JITting? Wouldn't we get ten times the portability with one-tenth the lines of code? Or is there too much dark magic going on in BLAS for such an approach to come close enough to hand-tuned performance? The dark magic in OpenBLAS is mostly to place prefetch instructions strategically, to make sure hierarchical memory is used optimally. This is very hard for the compiler to get correctly, because it doesn't know matrix algebra like we do. The reason prefetching is needed, is because when two matrices are multiplied, one of them will have strided memory access. On the other hand, putting in other SIMD instructions than _mm_prefetch is something a compiler might be able to vectorize without a lot of help today. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Standard Deviation (std): Suggested change for ddof default value
I agree; breaking code over this would be ridiculous. Also, I prefer the zero default, despite the mean/std combo probably being more common. On Tue, Apr 1, 2014 at 10:02 PM, Sturla Molden sturla.mol...@gmail.comwrote: Haslwanter Thomas thomas.haslwan...@fh-linz.at wrote: Personally I cannot think of many applications where it would be desired to calculate the standard deviation with ddof=0. In addition, I feel that there should be consistency between standard modules such as numpy, scipy, and pandas. ddof=0 is the maxiumum likelihood estimate. It is also needed in Bayesian estimation. If you are not eatimating from a sample, but rather calculating for the whole population, you always want ddof=0. What does Matlab do by default? (Yes, it is a retorical question.) I am wondering if there is a good reason to stick to ddof=0 as the default for std, or if others would agree with my suggestion to change the default to ddof=1? It is a bad idea to suddenly break everyone's code. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Is there a pure numpy recipe for this?
Id recommend taking a look at pytables as well. It has support for out-of-core array computations on large arrays. On Thu, Mar 27, 2014 at 9:00 PM, RayS r...@blue-cove.com wrote: Thanks for all of the suggestions; we are migrating to 64bit Python soon as well. The environments are Win7 and Mac Maverics. carray sounds like what you said Chris - more I just found at http://kmike.ru/python-data-structures/ - Ray Schumacher At 12:31 PM 3/27/2014, you wrote: On Thu, Mar 27, 2014 at 7:42 AM, RayS r...@blue-cove.com wrote: I find this interesting, since I work with medical data sets of 100s of MB, and regularly run into memory allocation problems when doing a lot of Fourrier analysis, waterfalls etc. The per-process limit seems to be about 1.3GB on this 6GB quad-i7 with Win7. This sounds like 32 bit -- have you tried a 64 bit Python_numpy? Nt that you wont have issues anyway, but you should be abel to do better than 1.3GB...   memmaps are also limited to RAM, I don't think so, no -- but are limited to 2GB (I think)  if you're using a 32 bit process There is also a compressed array package out there -- I can't remember what it's called -- but if you have large compressible arrays -- that might help.  -CHB -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/ORR       (206) 526-6959  voice 7600 Sand Point Way NE   (206) 526-6329  fax Seattle, WA  98115     (206) 526-6317  main reception chris.bar...@noaa.gov ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Is there a pure numpy recipe for this?
Without looking ahead, here is what I came up with; but I see more elegant solutions have been found already. import numpy as np def as_dense(f, length): i = np.zeros(length+1, np.int) i[f[0]] = 1 i[f[1]] = -1 return np.cumsum(i)[:-1] def as_sparse(d): diff = np.diff(np.concatenate(([0], d))) on, = np.nonzero(diff) on = on if on.size%2==0 else np.append(on, len(d)) return on.reshape(-1,2).T def join(f, g): on = np.sort(np.concatenate((f[0], g[0]))) off = np.sort(np.concatenate((f[1], g[1]))) I = np.argsort( np.concatenate((on, off)) ).argsort().reshape(2,-1) Q = -np.ones((2,I.size), np.int) Q[0,I[0]] = on Q[1,I[1]] = off idx_on = np.logical_and( Q[0,1:]*Q[0,:-1] 0, Q[0,:-1]!=-1) idx_off = np.logical_and( Q[1,1:]*Q[1,:-1] 0, Q[1,1:]!=-1) idx_on = np.concatenate( (idx_on, [False])) idx_off = np.concatenate( ([False], idx_off)) return np.array(( Q[0,idx_on], Q[1,idx_off])) length = 150 f_2_changes_at = np.array( [ 2 , 3 , 39, 41 , 58 , 59, 65 , 66 , 93 ,102, 145, length]) g_2_changes_at = np.array( [ 2 , 94 ,101, 146, 149, length]) f = f_2_changes_at.reshape(-1,2).T g = g_2_changes_at.reshape(-1,2).T dense_result = as_sparse( np.logical_and( as_dense(f, length), as_dense(g,length))) sparse_result = join(f,g) print np.allclose(dense_result, sparse_result) On Wed, Mar 26, 2014 at 10:50 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: On Wed, Mar 26, 2014 at 2:23 PM, Slaunger slaun...@gmail.com wrote: Jaime Fernández del Río wrote You saved my evening! Actually, my head has been spinning about this problem the last three evenings without having been able to nail it down. I had to quit Project Euler about 5 years ago because it was taking a huge toll on my mental health. I did learn/remember a ton of math, but was staying up all night banging my head against the problems much too often. Every now and then I do peek back and sometimes attempt a problem or two, but try to stay away for my own good. If you want to be projecteuler friends, I'm jfrio over there... Jaime ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Implementing elementary matrices
Sounds (marginally) useful; although elementary row/column operations are
in practice usually better implemented directly by indexing rather than in
an operator form. Though I can see a use for the latter.
My suggestion: its not a common enough operation to deserve a 4 letter
acronym (assuming those are good things in any context). A full
'elementary' would be much preferable I think.
On Mon, Mar 24, 2014 at 4:06 AM, Matt Pagan m...@pagan.io wrote:
Greetings!
I made a patch for NumPy that adds a function for easily creating
elementary matrices. Sorry for not knowing the process for submitting
patches.
Is this function something the NumPy community could see adding to the
codebase? Are there ways I can improve on this?
diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py
index 12c0f9b..10073af 100644
--- a/numpy/lib/twodim_base.py
+++ b/numpy/lib/twodim_base.py
@@ -967,3 +967,85 @@ def triu_indices_from(arr, k=0):
if not (arr.ndim == 2 and arr.shape[0] == arr.shape[1]):
raise ValueError(input array must be 2-d and square)
return triu_indices(arr.shape[0], k)
+
+def elem(N, i, j=None, t=None, dtype=float):
+
+Return an elementary matrix.
+
+Parameters
+--
+N : int
+The size of the NxN array to be returned. Elementary matrices
+should be square.
+i : int
+The index of the first row on which operations are to be
+performed.
+j : int
+If set, the index of the second row of which operations are to
+be performed.
+t : scalar
+If set, the factor by which a given row will be multiplied.
+
+Returns
+---
+m: ndarray of shape (NxN)
+ The identity matrix after a single row operation has been
+ performed on it.
+
+See also
+
+eye, identity
+
+Examples
+---
+To swap the the first and third rows of a 4x4 identity matirx:
+
+ L = elem(4, 0, 2)
+ L
+array([[ 0., 0., 1., 0.],
+ [ 0., 1., 0., 0.],
+ [ 1., 0., 0., 0.],
+ [ 0., 0., 0., 1.]])
+
+This array then becomes quite useful for matrix multiplication.
+
+ H = np.matrix([[ 2, 3, 5, 7],
+ [11, 13, 17, 19],
+ [23, 29, 31, 37],
+ [41, 43, 47, 53]])
+ L*H
+matrix([[ 23., 29., 31., 37.],
+[ 11., 13., 17., 19.],
+[ 2., 3., 5., 7.],
+[ 41., 43., 47., 53.]])
+
+When the elemntary matrix is multiplied by the given matrix, the
+result is the given matrix with it's first and third rows swapped.
+
+If the given matrix is multiplied by the elementary matrix (i.e.,
+the multiplication takes place in reverse order, the result is
+ the given matrix with its first and third columns swapped.
+
+ H*L
+matrix([[ 5., 3., 2., 7.],
+[ 17., 13., 11., 19.],
+[ 31., 29., 23., 37.],
+[ 47., 43., 41., 53.]])
+
+
+m=eye(N, dtype=dtype)
+if j==None and t==None:
+raise ValueError(One or more of %s and %s must be set. % \
+('j', 't'))
+return None
+elif t==None:
+swap = np.array(m[i])
+m[i] = m[j]
+m[j] = swap
+return m
+elif j==None:
+m[i] *= t
+return m
+else:
+m[j] += (t * m[i])
+return m
--
Matt Pagan
m...@pagan.io
PGP: 0xE9284418E360583C
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [help needed] associativity and precedence of '@'
Perhaps this a bit of a thread hyjack; but this discussion got me thinking about how to arrive at a more vectorized/tensorified way of specifying linear algebra operations, in an elegant manner. I probably got a little carried away, but what about this syntax? - indexing/calling an ndarray with a string returns a TensorExpression object - these TensorExpression objects can be combined into a graph using operator overloads - and these graphs are translated to calls to BLAS or einsum, as is appropriate #declare some symbols i,j,ij,k = 'i','j','ij','k' #we may force evaluation of a (sub) TensorExpression by calling it #this is trivial to translate to call to einsum #but such special cases could be dispatched to BLAS as well b = (A(ij) * x(j)) (i) #alternatively, we can predeclare a LHS which is automatically sized later #note that this translates into the same call as the above; just some syntactic sugar b = np.empty(()) b[i] = A(ij) * x(j) #more complex TensorExpression graphs of this form are trivial to translate to a call to einsum as well a(i)*b(j)*c(k) #conceptually, there is no need to limit this scheme to multiplications only! #although such generalizations would require a more complex execution engine #however, the revamped nditer should make this quite managable to implement a(i)*b(j) + c(k) #if axes strings are omitted, standard numpy broadcasting rules are applied to the expressiongraph created #this is identical to a*b+c; except that we have the opportunity to eliminate temporaries a()*b()+c() Note that such an approach kills quite some birds with one stone it allows for the elimination of temporaries along the lines of numexpr But if i could write: b[i] = A[ij] * x[j] I would much prefer that over b = A @ x even though the latter is shorter Now if i had n input and output vectors, it would be easy what to do with them: b[ni] = A[ij] * x[nj] As i argued earlier, I much prefer this form of explicitness over conventions about what constitutes a row or column vector. And vectorization of linear algebra is a trivial extension in this manner, which in itself is just a subset of even more general multilinear products, which themselves are a subset of more general expression involving things other than products Its a somewhat ambitious idea, and there are probably reasons why it isnt a good idea as well, but it does not require python language modifications, and it does not clash with any other functionality or syntax of numpy, as far as i can tell. Calling of arrays is not yet defined, and alternatively array indexing could be overloaded on string type. Either way, something to chew on when deciding on the best way to go forward. On Tue, Mar 18, 2014 at 4:28 AM, Mark Daoust daoust...@gmail.com wrote: On Mon, Mar 17, 2014 at 8:54 PM, Nathaniel Smith n...@pobox.com wrote: But, this is actually a feature! Because obviously what *should* be returned in this case is *not* (Mat @ vec) @ Mat, *or* Mat @ (vec @ Mat). Both of those answers are terrible; it's just, if you have an ordinary left-/right-associative operator, those are your only options. What *should* be returned is an error. And in this scheme we get to see the whole @ expression at once, so we actually can raise an error for such things. Sorry if this is a little off topic. But there's still something about the vector examples that bugs me, matrix@vector and vector@@2, keep popping up (this also applies to the matrix@matrix examples to a lesser extent). I'm a little unconformable looking at the shape to to decide what's a matrix and what's a vector. (Matlab has some problems like this) If it only has one or two dimensions it's easy, but I always find that if I've written code that works for 1 matrix or vector, 5 minutes later I want it to work for fields of matrices or vectors. If we're just going by shape there's no way to distinguish between a 2d field of matrices and a 3d field of vectors. I guess this is a repeat of part of what Eelco Hoogendoorn saying a few posts back I was just wondering if anyone sees a place, to get @ a little closer to Einsum, for some sort of array class that understands the difference between a 4D array of scalars, a 3D array of vectors, and a 2D array of matrices... The difference between the axes that broad-cast and the axes that can sum when you hit them with an @ ... or something like that. Just a thought. Einsum is fantastic by the way, totally worth learning and using. Mark Daoust ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [help needed] associativity and precedence of '@'
To elaborate a little on such a more general and explicit method of specifying linear operations (perhaps 'expressions with named axes' is a good nomer to cover this topic). I think indexing rather than calling is preferable. I worried at first about the performance overhead of checking for strings at every indexing op, but get ndarray__getitem__ is already quite a complex beast anyway, and adding (yet another) type test on its args isn't a significant difference. For those who disagree; we could also approach strings with a 'forgiveness is better then permission' attitude. The general rules could be: if no string args, everything works as normal. In case of string args, we may think of the effect of __getitem__ as indexing with strings replaced by colons first, and then creating a NamedAxisIndexExpression (NAIE), associating the given string label with each corresponding axis. Thus, we can write things like A[0:3,'i'] As some additional rules; string arguments can be 'expanded', the string is split on commas if present, and otherwise split into characters, which are then the axis labels. In expressions, all non-labeled axes are treated in sequential order, similar to the ... construct, and have standard numpy broadcasting semantics. The only problem with [] notation is field name lookup; though I have always felt that tables with named columns should be an ndarray subtype, given their fundamentally different indexing semantics. Realizing the full potential of such an approach would be a complex undertaking, but to start with, a more elegant interface to np.einsum would be rather easy to implement. On Tue, Mar 18, 2014 at 9:46 AM, Sebastian Haase seb.ha...@gmail.comwrote: Just add one vote: I am for * right association * because 1) I'm thinking of matrix multiplication more like operators, which I also learned to work from right to left and because 2) I would put a vector to the right, which would result in better performance. I don't have an opinion on tight/same/ or weak (maybe that means then 'same' because it's easier to remember !?) My two cents, Sebastian Haase On Tue, Mar 18, 2014 at 7:13 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Perhaps this a bit of a thread hyjack; but this discussion got me thinking about how to arrive at a more vectorized/tensorified way of specifying linear algebra operations, in an elegant manner. I probably got a little carried away, but what about this syntax? indexing/calling an ndarray with a string returns a TensorExpression object these TensorExpression objects can be combined into a graph using operator overloads and these graphs are translated to calls to BLAS or einsum, as is appropriate #declare some symbols i,j,ij,k = 'i','j','ij','k' #we may force evaluation of a (sub) TensorExpression by calling it #this is trivial to translate to call to einsum #but such special cases could be dispatched to BLAS as well b = (A(ij) * x(j)) (i) #alternatively, we can predeclare a LHS which is automatically sized later #note that this translates into the same call as the above; just some syntactic sugar b = np.empty(()) b[i] = A(ij) * x(j) #more complex TensorExpression graphs of this form are trivial to translate to a call to einsum as well a(i)*b(j)*c(k) #conceptually, there is no need to limit this scheme to multiplications only! #although such generalizations would require a more complex execution engine #however, the revamped nditer should make this quite managable to implement a(i)*b(j) + c(k) #if axes strings are omitted, standard numpy broadcasting rules are applied to the expressiongraph created #this is identical to a*b+c; except that we have the opportunity to eliminate temporaries a()*b()+c() Note that such an approach kills quite some birds with one stone it allows for the elimination of temporaries along the lines of numexpr But if i could write: b[i] = A[ij] * x[j] I would much prefer that over b = A @ x even though the latter is shorter Now if i had n input and output vectors, it would be easy what to do with them: b[ni] = A[ij] * x[nj] As i argued earlier, I much prefer this form of explicitness over conventions about what constitutes a row or column vector. And vectorization of linear algebra is a trivial extension in this manner, which in itself is just a subset of even more general multilinear products, which themselves are a subset of more general expression involving things other than products Its a somewhat ambitious idea, and there are probably reasons why it isnt a good idea as well, but it does not require python language modifications, and it does not clash with any other functionality or syntax of numpy, as far as i can tell. Calling of arrays is not yet defined, and alternatively array indexing could be overloaded on string type. Either way, something
Re: [Numpy-discussion] It looks like Py 3.5 will include a dedicated infix matrix multiply operator
Note that I am not opposed to extra operators in python, and only mildly opposed to a matrix multiplication operator in numpy; but let me lay out the case against, for your consideration. First of all, the use of matrix semantics relative to arrays semantics is extremely rare; even in linear algebra heavy code, arrays semantics often dominate. As such, the default of array semantics for numpy has been a great choice. Ive never looked back at MATLAB semantics. Secondly, I feel the urge to conform to a historical mathematical notation is misguided, especially for the problem domain of linear algebra. Perhaps in the world of mathematics your operation is associative or commutes, but on your computer, the order of operations will influence both outcomes and performance. Even for products, we usually care not only about the outcome, but also how that outcome is arrived at. And along the same lines, I don't suppose I need to explain how I feel about A@@-1 and the likes. Sure, it isn't to hard to learn or infer this implies a matrix inverse, but why on earth would I want to pretend the rich complexity of numerical matrix inversion can be mangled into one symbol? Id much rather write inv or pinv, or whatever particular algorithm happens to be called for given the situation. Considering this isn't the num-lisp discussion group, I suppose I am hardly the only one who feels so. On the whole, I feel the @ operator is mostly superfluous. I prefer to be explicit about where I place my brackets. I prefer to be explicit about the data layout and axes that go into a (multi)linear product, rather than rely on obtuse row/column conventions which are not transparent across function calls. When I do linear algebra, it is almost always vectorized over additional axes; how does a special operator which is only well defined for a few special cases of 2d and 1d tensors help me with that? On the whole, the linear algebra conventions inspired by the particular constraints of people working with blackboards, are a rather ugly and hacky beast in my opinion, which I feel no inclination to emulate. As a sidenote to the contrary; I love using broadcasting semantics when writing papers. Sure, your reviewers will balk at it, but it wouldn't do to give the dinosaurs the last word on what any given formal language ought to be like. We get to define the future, and im not sure the set of conventions that goes under the name of 'matrix multiplication' is one of particular importance to the future of numerical linear algebra. Note that I don't think there is much harm in an @ operator; but I don't see myself using it either. Aside from making textbook examples like a gram-schmidt orthogonalization more compact to write, I don't see it having much of an impact in the real world. On Sat, Mar 15, 2014 at 3:52 PM, Charles R Harris charlesr.har...@gmail.com wrote: On Fri, Mar 14, 2014 at 6:51 PM, Nathaniel Smith n...@pobox.com wrote: Well, that was fast. Guido says he'll accept the addition of '@' as an infix operator for matrix multiplication, once some details are ironed out: https://mail.python.org/pipermail/python-ideas/2014-March/027109.html http://legacy.python.org/dev/peps/pep-0465/ Specifically, we need to figure out whether we want to make an argument for a matrix power operator (@@), and what precedence/associativity we want '@' to have. I'll post two separate threads to get feedback on those in an organized way -- this is just a heads-up. Surprisingly little discussion on python-ideas, or so it seemed to me. Guido came out in favor less than halfway through. Congratulations on putting together a successful proposal, many of us had given up on ever seeing a matrix multiplication operator. Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] It looks like Py 3.5 will include a dedicated infix matrix multiply operator
Different people work on different code and have different experiences here -- yours may or may be typical yours. Pauli did some quick checks on scikit-learn nipy scipy, and found that in their test suites, uses of np.dot and uses of elementwise-multiplication are ~equally common: https://github.com/numpy/numpy/pull/4351#issuecomment-37717330h Yeah; these are examples of linalg-heavy packages. Even there, dot does not dominate. My impression from the other thread is that @@ probably won't end up existing, so you're safe here ;-). I know; my point is that the same objections apply to @, albeit in weaker form. Einstein notation is coming up on its 100th birthday and is just as blackboard-friendly as matrix product notation. Yet there's still a huge number of domains where the matrix notation dominates. It's cool if you aren't one of the people who find it useful, but I don't think it's going anywhere soon. Einstein notation is just as blackboard friendly; but also much more computer-future proof. I am not saying matrix multiplication is going anywhere soon; but as far as I can tell that is all inertia; historical circumstance has not accidentially prepared it well for numerical needs, as far as I can tell. The analysis in the PEP found ~780 calls to np.dot, just in the two projects I happened to look at. @ will get tons of use in the real world. Maybe all those people who will be using it would be happier if they were using einsum instead, I dunno, but it's an argument you'll have to convince them of, not me :-). 780 calls is not tons of use, and these projects are outliers id argue. I just read for the first time two journal articles in econometrics that use einsum notation. I have no idea what their formulas are supposed to mean, no sum signs and no matrix algebra. If they could have been expressed more clearly otherwise, of course this is what they should have done; but could they? b_i = A_ij x_j isnt exactly hard to read, but if it was some form of complicated product, its probably tensor notation was their best bet.___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] It looks like Py 3.5 will include a dedicated infix matrix multiply operator
An important distinction between calling dot or @ is that matrix multiplication is a domain where enormous effort has already been spent on algorithms and building fast, scalable libraries. Yes einsum can call these for some subset of calls but it's also trivial to set up a case where it can't. This is a huge pitfall because it hides this complexity. Einsum, despite the brevity that it can provide, is too general to make a basic building block. There isn't a good way to reason about its runtime. I am not arguing in favor of einsum; I am arguing in favor of being explicit, rather than hiding semantically meaningful information from the code. Whether using @ or dot or einsum, you are not explicitly specifying the type of algorithm used, so on that front, its a wash, really. But at least dot and einsum have room for keyword arguments. '@' is in my perception simply too narrow an interface to cram in all meaningful information that you might want to specify concerning a linear product. Matrix-matrix and matrix-vector products are the fundamental operations, generalized multilinear products etc are not. Perhaps from a library perspective, but from a conceptual perspective, it is very much the other way around. If we keep going in the direction that numba/theano/loopy take, such library functionality will soon be moot. Id argue that the priority of the default semantics should be in providing a unified conceptual scheme, rather than maximum performance considerations. Ideally, the standard operator would pick a sensible default which can be inferred from the arguments, while allowing for explicit specification of the kind of algorithm used where this verbosity is worth the hassle. On Sun, Mar 16, 2014 at 5:33 PM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Different people work on different code and have different experiences here -- yours may or may be typical yours. Pauli did some quick checks on scikit-learn nipy scipy, and found that in their test suites, uses of np.dot and uses of elementwise-multiplication are ~equally common: https://github.com/numpy/numpy/pull/4351#issuecomment-37717330h Yeah; these are examples of linalg-heavy packages. Even there, dot does not dominate. My impression from the other thread is that @@ probably won't end up existing, so you're safe here ;-). I know; my point is that the same objections apply to @, albeit in weaker form. Einstein notation is coming up on its 100th birthday and is just as blackboard-friendly as matrix product notation. Yet there's still a huge number of domains where the matrix notation dominates. It's cool if you aren't one of the people who find it useful, but I don't think it's going anywhere soon. Einstein notation is just as blackboard friendly; but also much more computer-future proof. I am not saying matrix multiplication is going anywhere soon; but as far as I can tell that is all inertia; historical circumstance has not accidentially prepared it well for numerical needs, as far as I can tell. The analysis in the PEP found ~780 calls to np.dot, just in the two projects I happened to look at. @ will get tons of use in the real world. Maybe all those people who will be using it would be happier if they were using einsum instead, I dunno, but it's an argument you'll have to convince them of, not me :-). 780 calls is not tons of use, and these projects are outliers id argue. I just read for the first time two journal articles in econometrics that use einsum notation. I have no idea what their formulas are supposed to mean, no sum signs and no matrix algebra. If they could have been expressed more clearly otherwise, of course this is what they should have done; but could they? b_i = A_ij x_j isnt exactly hard to read, but if it was some form of complicated product, its probably tensor notation was their best bet. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Pickling of memory aliasing patterns
I have been working on a general function caching mechanism, and in doing so I stumbled upon the following quirck: @cached def foo(a,b): b[0] = 1 return a[0] a = np.zeros(1) b = a[:] print foo(a, b)#computes and returns 1 print foo(a, b)#gets 1 from cache, as it should a = np.zeros(1) #no aliasing between inputs b = np.zeros(1) print foo(a, b)#should compute and return 0 but instead gets 1 from cache Fundamentaly, this is because it turns out that the memory aliasing patterns that arrays may have are lost during pickling. This leads me to two questions: 1: Is this desirable behavior 2: is this preventable behavior? It seems to me the answer to the first question is no, and the answer to the second question is yes. Here is what I am using at the moment to generate a correct hash under such circumstances; but unpickling along these lines should be possible too, methinks. Or am I missing some subtlety as to why something along these lines couldn't be the default pickling behavior for numpy arrays? class ndarray_own(object): def __init__(self, arr): self.buffer = np.getbuffer(arr) self.dtype = arr.dtype self.shape = arr.shape self.strides= arr.strides class ndarray_view(object): def __init__(self, arr): self.base = arr.base self.offset = self.base.ctypes.data - arr.ctypes.data #so we have a view; but where is it? self.dtype = arr.dtype self.shape = arr.shape self.strides= arr.strides class NumpyDeterministicPickler(DeterministicPickler): Special case for numpy. in general, external C objects may include internal state which does not serialize in a way we want it to ndarray memory aliasing is one of those things def save(self, obj): remap a numpy array to a representation which conserves all semantically relevant information concerning memory aliasing note that this mapping is 'destructive'; we will not get our original numpy arrays back after unpickling; not without custom deserialization code at least but we dont care, since this is only meant to be used to obtain correct keying behavior keys dont need to be deserialized if isinstance(obj, np.ndarray): if obj.flags.owndata: obj = ndarray_own(obj) else: obj = ndarray_view(obj) DeterministicPickler.save(self, obj) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] dtype promotion
The tuple gets cast to an ndarray; which invokes a different codepath than the scalar addition. Somehow, numpy has gotten more aggressive at upcasting to float64 as of 1.8, but I havnt been able to discover the logic behind it either. On Mon, Mar 3, 2014 at 10:06 PM, Nicolas Rougier nicolas.roug...@inria.frwrote: Hi all, I'm using numpy 1.8.0 (osx 10.9, python 2.7.6) and I can't understand dtype promotion in the following case: Z = np.zeros((2,2),dtype=np.float32) + 1 print Z.dtype float32 Z = np.zeros((2,2),dtype=np.float32) + (1,1) print Z.dtype float64 Is this the expected behavior ? What it the difference between the two lines ? Nicolas ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ANN: XDress v0.4
I have; but if I recall correctly, it does not solve the problem of distributing code that uses it, or does it? On Thu, Feb 27, 2014 at 10:51 AM, Toby St Clere Smithe pyvienn...@tsmithe.net wrote: Hi, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com writes: Thanks for the heads up, I wasn't aware of this project. While boost.python is a very nice package, its distributability is nothing short of nonexistent, so its great to have a pure python binding generator. One thing which I have often found frustrating is natural ndarray interop between python and C++. Is there a (planned) mechanism for mapping arbitrary strided python ndarrays to boost arrays? Have you tried boost.numpy? https://github.com/ndarray/Boost.NumPy I have a fork which builds against Python 3, as well -- though it's mainly used for PyViennaCL, and might need a bit of cleaning. Cheers, Toby On Thu, Feb 27, 2014 at 1:24 AM, Anthony Scopatz scop...@gmail.com wrote: Hello All, I am *extremely *pleased to be able to announce the version 0.4 release of xdress. This version contains much anticipated full support for Clang as a parser! This is almost entirely due to the efforts of Geoffrey Irving. Please thank him the next time you get a chance :) This release also contains a lot of other goodies that you can read about in the release notes below. Happy Generating! Anthony XDress 0.4 Release Notes http://xdress.org/previous/0.4_release_notes.html#xdress-0-4-release-notes XDress is a numpy-aware 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: - Clang support! All kudos to Geoffrey Irving! - NumPy dtypes may be created independently of C++ STL vectors - A complete test suite refactor - Arbitrary source code locations - Global run control files - A plethora of useful bug fixes This version of xdress is *not* 100% backwards compatible with previous versions of xdress. We apologize in the name of progress. It represents ans impressive 245 files changed, 44917 aggregate line insertions (+), and 7893 deletions (-). Please visit the website for more information: http://xdress.org/ Ask questions on the mailing list: https://groups.google.com/forum/#!forum/xdress Download the code from GitHub: http://github.com/xdress/xdress XDress is free open source (BSD 2-clause license) and requires Python 2.7+, NumPy 1.5+, Cython 0.19+, and optionally Clang, GCC-XML, pycparser, dOxygen, or lxml. New Features http://xdress.org/previous/0.4_release_notes.html#new-features Clang Support http://xdress.org/previous/0.4_release_notes.html#clang-support Through the herculean efforts of Geoffrey Irving xdress finally has full, first-class Clang/LLVM support! This is major advancement as it allows xdress to wrap more modern versions of C++ than GCC-XML can handle. Because of deficiencies in the existing libclang and Python bindings it was necessary for us to fork libclang for xdress in the short term. We hope to integrate these changes upstream. Clang versions 3.2 - 3.4 are supported. Independent NumPy Dtypes http://xdress.org/previous/0.4_release_notes.html#independent-numpy-dtypes In previous versions of xdress, to create a dtype of type T the user needed to declare the desire for a wrapper of an STL vector of type T. These two desires have now been separated. It is now possible to create a dtype via the dtypes run control parameter. STL vectors are still wrapped via dtypes. See the dtypes module for more information. Shiny New Test Suite http://xdress.org/previous/0.4_release_notes.html#shiny-new-test-suite The xdress test suite has been completely revamped to include both unit and integration tests which are run for all available parsers. The integration tests are accomplished though two fake projects - cproj and cppproj - on which the xdress CLI is run. These tests are now fully platform independent, unlike the previous BASH-based test suite. Source Paths http://xdress.org/previous/0.4_release_notes.html#source-paths Source file paths are now given by either their absolute or relative path. This allows source code to be located anywhere on the user's file system and enable the wrapping of dependencies or externally supplied libraries as needed. The run control parametersourcedir has been deprecated. Global Run Control Files http://xdress.org/previous/0.4_release_notes.html#global-run-control-files
Re: [Numpy-discussion] ANN: XDress v0.4
I have a file numpy_boost_python.hpp in one of my projects by Michael Droettboom (can seem to find an online source anymore!), which adds mappings between numpy.ndarray and boost.ndarray, which is very neat and seemless. But like boost.python, it tightly couples with the clusterfuck that is bjam. However, something conceptually like that but integrated with XDress would be great. Indeed it does not sound too complicated; though I don't think I will get around to it anytime soon, unfortunately... On Thu, Feb 27, 2014 at 1:36 PM, Anthony Scopatz scop...@gmail.com wrote: On Thu, Feb 27, 2014 at 1:51 AM, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: Thanks for the heads up, I wasn't aware of this project. While boost.python is a very nice package, its distributability is nothing short of nonexistent, so its great to have a pure python binding generator. Thanks! One thing which I have often found frustrating is natural ndarray interop between python and C++. Is there a (planned) mechanism for mapping arbitrary strided python ndarrays to boost arrays? Not yet! The architecture is very modular (it is just a series of plugins) so I would welcome anyone who wants to tackle this to take a look into it. I don't think that it would be *that* hard. You'd just need to write the Py-to-C++ and C++-to-Py converters for the boost array type. This shouldn't be too hard since std::vector goes through pretty much exactly the same mechanism for exposing to numpy. So there are already a couple of examples of this workflow. Please feel free to jump on the xdress mailing list if you want to discuss this in more depth! Also I didn't know about ndarray/Boost.NumPy. This seems like it could be useful! Be Well Anthony On Thu, Feb 27, 2014 at 1:24 AM, Anthony Scopatz scop...@gmail.comwrote: Hello All, I am *extremely *pleased to be able to announce the version 0.4 release of xdress. This version contains much anticipated full support for Clang as a parser! This is almost entirely due to the efforts of Geoffrey Irving. Please thank him the next time you get a chance :) This release also contains a lot of other goodies that you can read about in the release notes below. Happy Generating! Anthony XDress 0.4 Release Noteshttp://xdress.org/previous/0.4_release_notes.html#xdress-0-4-release-notes XDress is a numpy-aware 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: - Clang support! All kudos to Geoffrey Irving! - NumPy dtypes may be created independently of C++ STL vectors - A complete test suite refactor - Arbitrary source code locations - Global run control files - A plethora of useful bug fixes This version of xdress is *not* 100% backwards compatible with previous versions of xdress. We apologize in the name of progress. It represents ans impressive 245 files changed, 44917 aggregate line insertions (+), and 7893 deletions (-). Please visit the website for more information: http://xdress.org/ Ask questions on the mailing list: https://groups.google.com/forum/#!forum/xdress Download the code from GitHub: http://github.com/xdress/xdress XDress is free open source (BSD 2-clause license) and requires Python 2.7+, NumPy 1.5+, Cython 0.19+, and optionally Clang, GCC-XML, pycparser, dOxygen, or lxml. New Featureshttp://xdress.org/previous/0.4_release_notes.html#new-features Clang Supporthttp://xdress.org/previous/0.4_release_notes.html#clang-support Through the herculean efforts of Geoffrey Irving xdress finally has full, first-class Clang/LLVM support! This is major advancement as it allows xdress to wrap more modern versions of C++ than GCC-XML can handle. Because of deficiencies in the existing libclang and Python bindings it was necessary for us to fork libclang for xdress in the short term. We hope to integrate these changes upstream. Clang versions 3.2 - 3.4 are supported. Independent NumPy Dtypeshttp://xdress.org/previous/0.4_release_notes.html#independent-numpy-dtypes In previous versions of xdress, to create a dtype of type T the user needed to declare the desire for a wrapper of an STL vector of type T. These two desires have now been separated. It is now possible to create a dtype via the dtypes run control parameter. STL vectors are still wrapped via dtypes. See the dtypes module for more information. Shiny New Test Suitehttp://xdress.org/previous/0.4_release_notes.html#shiny-new-test-suite The xdress test suite has been completely revamped to include both
Re: [Numpy-discussion] ANN: XDress v0.4
That is good to know. The boost documentation makes it appear as if bjam is the only way to build boost.python, but good to see examples to the contrary! On Thu, Feb 27, 2014 at 2:19 PM, Toby St Clere Smithe pyvienn...@tsmithe.net wrote: Eelco Hoogendoorn hoogendoorn.ee...@gmail.com writes: I have a file numpy_boost_python.hpp in one of my projects by Michael Droettboom (can seem to find an online source anymore!), which adds mappings between numpy.ndarray and boost.ndarray, which is very neat and seemless. But like boost.python, it tightly couples with the clusterfuck that is bjam. However, something conceptually like that but integrated with XDress would be great. Indeed it does not sound too complicated; though I don't think I will get around to it anytime soon, unfortunately... You don't have to use bjam! I have built my projects with distutils and CMake, and never once touched bjam; CMake provides find_package scripts for Boost, Python and NumPy, and for distutils, I just include the relevant files and flags in my project. See [1] for a CMake example, and [2] for a distutils example. [1] https://github.com/tsmithe/viennacl-dev/blob/pyviennacl/CMakeLists.txt [2] https://github.com/viennacl/pyviennacl-dev/blob/master/setup.py Cheers, Toby ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ANN: XDress v0.4
Thanks for the heads up, I wasn't aware of this project. While boost.python is a very nice package, its distributability is nothing short of nonexistent, so its great to have a pure python binding generator. One thing which I have often found frustrating is natural ndarray interop between python and C++. Is there a (planned) mechanism for mapping arbitrary strided python ndarrays to boost arrays? On Thu, Feb 27, 2014 at 1:24 AM, Anthony Scopatz scop...@gmail.com wrote: Hello All, I am *extremely *pleased to be able to announce the version 0.4 release of xdress. This version contains much anticipated full support for Clang as a parser! This is almost entirely due to the efforts of Geoffrey Irving. Please thank him the next time you get a chance :) This release also contains a lot of other goodies that you can read about in the release notes below. Happy Generating! Anthony XDress 0.4 Release Noteshttp://xdress.org/previous/0.4_release_notes.html#xdress-0-4-release-notes XDress is a numpy-aware 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: - Clang support! All kudos to Geoffrey Irving! - NumPy dtypes may be created independently of C++ STL vectors - A complete test suite refactor - Arbitrary source code locations - Global run control files - A plethora of useful bug fixes This version of xdress is *not* 100% backwards compatible with previous versions of xdress. We apologize in the name of progress. It represents ans impressive 245 files changed, 44917 aggregate line insertions (+), and 7893 deletions (-). Please visit the website for more information: http://xdress.org/ Ask questions on the mailing list: https://groups.google.com/forum/#!forum/xdress Download the code from GitHub: http://github.com/xdress/xdress XDress is free open source (BSD 2-clause license) and requires Python 2.7+, NumPy 1.5+, Cython 0.19+, and optionally Clang, GCC-XML, pycparser, dOxygen, or lxml. New Featureshttp://xdress.org/previous/0.4_release_notes.html#new-features Clang Supporthttp://xdress.org/previous/0.4_release_notes.html#clang-support Through the herculean efforts of Geoffrey Irving xdress finally has full, first-class Clang/LLVM support! This is major advancement as it allows xdress to wrap more modern versions of C++ than GCC-XML can handle. Because of deficiencies in the existing libclang and Python bindings it was necessary for us to fork libclang for xdress in the short term. We hope to integrate these changes upstream. Clang versions 3.2 - 3.4 are supported. Independent NumPy Dtypeshttp://xdress.org/previous/0.4_release_notes.html#independent-numpy-dtypes In previous versions of xdress, to create a dtype of type T the user needed to declare the desire for a wrapper of an STL vector of type T. These two desires have now been separated. It is now possible to create a dtype via the dtypes run control parameter. STL vectors are still wrapped via dtypes. See the dtypes module for more information. Shiny New Test Suitehttp://xdress.org/previous/0.4_release_notes.html#shiny-new-test-suite The xdress test suite has been completely revamped to include both unit and integration tests which are run for all available parsers. The integration tests are accomplished though two fake projects - cproj and cppproj - on which the xdress CLI is run. These tests are now fully platform independent, unlike the previous BASH-based test suite. Source Pathshttp://xdress.org/previous/0.4_release_notes.html#source-paths Source file paths are now given by either their absolute or relative path. This allows source code to be located anywhere on the user's file system and enable the wrapping of dependencies or externally supplied libraries as needed. The run control parametersourcedir has been deprecated. Global Run Control Fileshttp://xdress.org/previous/0.4_release_notes.html#global-run-control-files It is sometimes useful to be able to set system-wide run control parameters. XDress will now search the following files in order of increasing precedence. - $HOME/.xdressrc - $HOME/.xdressrc.py - $HOME/.config/xdressrc - $HOME/.config/xdressrc.py $HOME is the user's home directory. Settings in the project run control file take precedence over the values here. Major Bug Fixeshttp://xdress.org/previous/0.4_release_notes.html#major-bug-fixes - Debug file now always written when in debug mode. - STL sets of custom types now allowed. - Template parameters now allowed to be enum
Re: [Numpy-discussion] Help Understanding Indexing Behavior
To elaborate on what Julian wrote: it is indeed simply a convention; slices/ranges in python are from the start to one-past-the-end. The reason for the emergence of this convention is that C code using iterators looks most natural this way. This manifests in a simple for (i = 0; i 5; i++), but also when specifying a slice of a linked list, for instance. We don't want to terminate the loop when we are just arriving at the last item; we want to terminate a loop when we have gone past the last item. Also, the length of a range is simply end-start under this convention; no breaking your head over -1 or +1. Such little nudges of elegance pop up all over C code; and that's where the convention comes from. Same as zero-based indexing; just a convention, and if you are going to pick a convention you might as well pick one that minimizes the number of required operations. Anything but zero-based indexing will require additional integer math to find an array element, given its base pointer. On Wed, Feb 26, 2014 at 12:15 AM, Julian Taylor jtaylor.deb...@googlemail.com wrote: On 26.02.2014 00:04, JB wrote: At the risk of igniting a flame war...can someone please help me understand the indexing behavior of NumPy? I will readily I admit I come from a Matlab background, but I appreciate the power of Python and am trying to learn more. From a Matlab user's perspective, the behavior of indexing in NumPy seems very bizarre. For example, if I define an array: x = np.array([1,2,3,4,5,6,7,8,9,10]) If I want the first 5 elements, what do I do? Well, I say to myself, Python is zero-based, whereas Matlab is one-based, so if I want the values 1 - 5, then I want to index 0 - 4. So I type: x[0:4] And get in return: array([1, 2, 3, 4]). So I got the first value of my array, but I did not get the 5th value of the array. So the start index needs to be zero-based, but the end index needs to be one-based. Or to put it another way, if I type x[4] and x[0:4], the 4 means different things depending on which set of brackets you're looking at! It's hard for me to see this as anything by extremely confusing. Can someone explain this more clearly. Feel free to post links if you'd like. I know this has been discussed ad nauseam online; I just haven't found any of the explanations satisfactory (or sufficiently clear, at any rate). numpy indexing is like conventional C indexing beginning from inclusive 0 to exclusive upper bound: [0, 5[. So the selection length is upper bound - lower bound. as a for loop: for (i = 0; i 5; i++) select(i); This is the same way Python treats slices. in comparison one based indexing is usually inclusive 1 to inclusive upper bound: [1, 4]. So the selection length is upper bound - lower bound + 1. for (i = 1; i = 4; i++) select(i); ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Proposal: Chaining np.dot with mdot helper function
If the standard semantics are not affected, and the most common two-argument scenario does not take more than a single if-statement overhead, I don't see why it couldn't be a replacement for the existing np.dot; but others mileage may vary. On Thu, Feb 20, 2014 at 11:34 AM, Stefan Otte stefan.o...@gmail.com wrote: Hey, so I propose the following. I'll implement a new function `mdot`. Incorporating the changes in `dot` are unlikely. Later, one can still include the features in `dot` if desired. `mdot` will have a default parameter `optimize`. If `optimize==True` the reordering of the multiplication is done. Otherwise it simply chains the multiplications. I'll test and benchmark my implementation and create a pull request. Cheers, Stefan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Proposal: Chaining np.dot with mdot helper function
Erik; take a look at np.einsum The only reason against such dot semantics is that there isn't much to be gained in elegance that np.einsum already provides, For a plain chaining, multiple arguments to dot would be an improvement; but if you want to go for more complex products, the elegance of np.einsum will be hard to beat On Thu, Feb 20, 2014 at 3:27 PM, Eric Moore e...@redtetrahedron.org wrote: On Thursday, February 20, 2014, Eelco Hoogendoorn hoogendoorn.ee...@gmail.com wrote: If the standard semantics are not affected, and the most common two-argument scenario does not take more than a single if-statement overhead, I don't see why it couldn't be a replacement for the existing np.dot; but others mileage may vary. On Thu, Feb 20, 2014 at 11:34 AM, Stefan Otte stefan.o...@gmail.comwrote: Hey, so I propose the following. I'll implement a new function `mdot`. Incorporating the changes in `dot` are unlikely. Later, one can still include the features in `dot` if desired. `mdot` will have a default parameter `optimize`. If `optimize==True` the reordering of the multiplication is done. Otherwise it simply chains the multiplications. I'll test and benchmark my implementation and create a pull request. Cheers, Stefan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion Another consideration here is that we need a better way to work with stacked matrices such as np.linalg handles now. Ie I want to compute the matrix product of two (k, n, n) arrays producing a (k,n,n) result. Near as I can tell there isn't a way to do this right now that doesn't involve an explicit loop. Since dot will return a (k, n, k, n) result. Yes this output contains what I want but it also computes a lot of things that I don't want too. It would also be nice to be able to do a matrix product reduction, (k, n, n) - (n, n) in a single line too. Eric ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Proposal: Chaining np.dot with mdot helper function
considering np.dot takes only its binary positional args and a single defaulted kwarg, passing in a variable number of positional args as a list makes sense. Then just call the builtin reduce on the list, and there you go. I also generally approve of such semantics for binary associative operations. On Mon, Feb 17, 2014 at 11:27 PM, Jaime Fernández del Río jaime.f...@gmail.com wrote: Perhaps you could reuse np.dot, by giving its second argument a default None value, and passing a tuple as first argument, i.e. np.dot((a, b, c)) would compute a.dot(b).dot(c), possibly not in that order. As is suggested in the matlab thread linked by Josef, if you do implement an optimal ordering algorithm, then precalculating the ordering and passing it in as an argument should be an option. If I get a vote, I am definitely +1 on this, especially the more sophisticated version. On Feb 17, 2014 1:40 PM, Stefan Otte stefan.o...@gmail.com wrote: ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Requesting Code Review of nanmedian ENH
hi david, I havnt run the code; but the _replace_nan(0) call worries me; especially considering that the unit tests seem to deal with positive numbers exclusively. Have you tested with mixed positive/negative inputs? On Sun, Feb 16, 2014 at 6:13 PM, David Freese dfre...@stanford.edu wrote: Hi everyone, I put together a np.nanmedian function to extend np.median to handle nans. Could someone review this code and give me some feedback on it before I submit a pull request for it? https://github.com/dfreese/numpy/compare/master...feature;nanmedian Thanks, David ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] argsort speed
My guess; First of all, you are actually manipulating twice as much data as opposed to an inplace sort. Moreover, an inplace sort gains locality as it is being sorted, whereas the argsort is continuously making completely random memory accesses. -Original Message- From: josef.p...@gmail.com Sent: Sunday, February 16, 2014 11:43 PM To: Discussion of Numerical Python Subject: [Numpy-discussion] argsort speed currently using numpy 1.6.1 What's the fastest argsort for a 1d array with around 28 Million elements, roughly uniformly distributed, random order? Is there a reason that np.argsort is almost 3 times slower than np.sort? I'm doing semi-systematic timing for a stats(models) algorithm. Josef ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] libflatarray
As usual, 'it depends', but a struct of arrays layout (which is a virtual necessity on GPU's), can also be advantageous on the CPU. One rarely acts on only a single object at a time; but quite often, you only work on a subset of the objects attributes at a time. In an array of structs layout, you are always pulling the whole objects from main memory into the cache, even if you only use a single attribute. In a struct of arrays layout, we can do efficient prefetching on a single attribute when looping over all objects. On Thu, Feb 13, 2014 at 3:37 PM, Sturla Molden sturla.mol...@gmail.comwrote: Neal Becker ndbeck...@gmail.com wrote: I thought this was interesting: http://www.libgeodecomp.org/libflatarray.html This is mostly flawed thinking. Nowadays, CPUs are much faster than memory access, and the gap is just increasing. In addition, CPUs have hierarchical memory (several layers of cache). Most algorithms therefore benefit from doing as much computation on the data as possible, before reading more data from RAM. That means that an interleaved memory layout is usually the more effective. This of course deepends on the algorithm, but an array of structs is usually better than a struct of arrays. Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] deprecate numpy.matrix
My 2pc: I personally hardly ever use matrix, even in linear algebra dense code. It can be nice though, to use matrix semantics within a restricted scope. When I first came to numpy, the ability to choose linear algebra versus array semantics seemed like a really neat thing to me; though in practice the array semantics are so much more useful you really don't mind having to write the occasional np.dot. There seems to be some resistance form the developer side in having to maintain this architecture, which I cannot comment on, but from a user perspective, I am perfectly fine with the way things are. On Tue, Feb 11, 2014 at 11:16 AM, Todd toddr...@gmail.com wrote: On Feb 11, 2014 3:23 AM, Alan G Isaac alan.is...@gmail.com wrote: On 2/10/2014 7:39 PM, Pauli Virtanen wrote: The issue here is semantics for basic linear algebra operations, such as matrix multiplication, that work for different matrix objects, including ndarrays. I'll see if I can restate my suggestion in another way, because I do not feel you are responding to it. (I might be wrong.) What is a duck? If you ask it to quack, it quacks. OK, but what is it to quack? Here, quacking is behaving like an ndarray (in your view, as I understand it) when asked. But how do we ask? Your view (if I understand) is we ask via the operations supported by ndarrays. But maybe that is the wrong way for the library to ask this question. If so, then scipy libraries could ask an object to behave like an an ndarray by calling, e.g., __asarray__ on it. It becomes the responsibility of the object to return something appropriate when __asarray__ is called. Objects that know how to do this will provide __asarray__ and respond appropriately. Other types can be coerced if that is the documented behavior (e.g., lists). The libraries will then be written for a single type of behavior. What it means to quack is pretty easily documented, and a matrix object already knows how (e.g., m.A). Presumably in this scenario __asarray__ would return an object that behaves like an ndarray and a converter for turning the final result into the desired object type (e.g., into a `matrix` if necessary). Hope that clearer, even if it proves a terrible idea. Alan Isaac I don't currently use the matrix class, but having taken many linear algebra classes I can see the appeal, and if I end up teaching the subject I think I would appreciate having it available. On the other hand, I certainly can see the possibility for confusion, and I don't think it is something that should be used unless someone has a really good reason. So I come out somewhere in the middle here. So, although this may end up being a terrible idea, I would like to purpose what I think is a compromise: instead of just removing matrix, split it out into a scikit. That way, it it's still available for those who need it, but will be less likely to be used accidentally, and won't be interfering with the rest of numpy and scipy development. Specifically, I would split matrix into a scikit, while in the same release deprecate np.matrix. They can then exist in parallel for a few releases to allow code to be ported away from it. However, I would suggest that before the split, all linear algebra routines should be available as functions or methods in numpy proper. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
