Re: [Numpy-discussion] Bug fix or behavior change?
On Sat, Sep 27, 2008 at 12:51 AM, Travis E. Oliphant <[EMAIL PROTECTED] > wrote: > Charles R Harris wrote: > > Hi All, > > > > Currently subtract for boolean arrays is defined in > > > > /**begin repeat > > * Arithmetic operators > > * > > * # OP = ||, ^, & > > * #kind = add, subtract, multiply# > > */ > > static void > > [EMAIL PROTECTED]@(char **args, intp *dimensions, intp *steps, void *func) > > { > > register intp i; > > intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; > > char *i1=args[0], *i2=args[1], *op=args[2]; > > for(i=0; i > *((Bool *)op)=*((Bool *)i1) @OP@ *((Bool *)i2); > > } > > } > > /**end repeat**/ > > > > > > Note that this might yield unexpected results if the boolean value is > > not 0 or 1, which is possible using views. Note also that bitwise_xor > > converts the boolean to 0 or 1 before using ^. I think subtract should > > work the same way, but that would change current behavior. So... do I > > make the change and call it a bug fix or leave it as it is? > > I think it's a bugfix and so am +1 on the change. > Done. I've made all the boolean ufuncs convert args to 0 or 1 before applying the operators. Next question, what about logical_and applied to complex numbers. Currently it gives In [16]: a = complex(1,0) In [17]: b = complex(0,1) In [18]: logical_and(a,b) Out[18]: False i.e., it && all the components together. I suspect what it should be doing is anding first and second arguments being non-zero. Also, there is this bit for absolute of signed numbers static void @[EMAIL PROTECTED](char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1 = *((@type@ *)ip1); *((@type@ *)op) = (in1 < 0) ? -in1 : in1; *((@type@ *)op) += 0; /* clear sign-bit */ } } I don't understand the clearing the sign-bit line, it looks unnecessary. Was there some problem that this once fixed? Also, the signature for the ones_like ufuncs has 1 input argument, but only the output argument is used. I suppose changing this would be a pain, but it's a bit odd. Chuck ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Bug fix or behavior change?
On Sat, Sep 27, 2008 at 1:18 AM, Charles R Harris <[EMAIL PROTECTED] > wrote: > > > On Sat, Sep 27, 2008 at 12:51 AM, Travis E. Oliphant < > [EMAIL PROTECTED]> wrote: > >> Charles R Harris wrote: >> > Hi All, >> > >> > Currently subtract for boolean arrays is defined in >> > >> > /**begin repeat >> > * Arithmetic operators >> > * >> > * # OP = ||, ^, & >> > * #kind = add, subtract, multiply# >> > */ >> > static void >> > [EMAIL PROTECTED]@(char **args, intp *dimensions, intp *steps, void *func) >> > { >> > register intp i; >> > intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; >> > char *i1=args[0], *i2=args[1], *op=args[2]; >> > for(i=0; i> > *((Bool *)op)=*((Bool *)i1) @OP@ *((Bool *)i2); >> > } >> > } >> > /**end repeat**/ >> > >> > >> > Note that this might yield unexpected results if the boolean value is >> > not 0 or 1, which is possible using views. Note also that bitwise_xor >> > converts the boolean to 0 or 1 before using ^. I think subtract should >> > work the same way, but that would change current behavior. So... do I >> > make the change and call it a bug fix or leave it as it is? >> >> I think it's a bugfix and so am +1 on the change. >> > > Done. I've made all the boolean ufuncs convert args to 0 or 1 before > applying the operators. Next question, what about logical_and applied to > complex numbers. Currently it gives > I'm also thinking that bitwise operators should not be allowed for booleans. But that might break some code... Chuck ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ANN: NumPy 1.2.0
Thomas Heller wrote: > > Well, the first question is: What does happen when I install the SSE3 version > (or how it's called) on my machine, use py2exe to build an app, and this > app runs on a SSE2 machine - degraded performance, or hard crashes? Hard crash. That's the whole point of the installer, actually: install the right version. The problem is two fold: - it is practically impossible to force ATLAS to configure itself to use an instruction set different than the one used to build ATLAS. - we don't have a system to dynamically load the right ATLAS when importing numpy. Note that this is not new: you had the problem before, because before, numpy was *only* built with SSE2 support, and any machine wo SSE2 would crash when using numpy. > > So, maybe the gui could allow to select whether to install the > high-performance > version specialized for the current cpu, or a more portable but a little bit > slower version (I assume there is one included) that can be safely used for > py2exe. Is it really complicated to decompress the .exe to get the installers and choose the one you want ? I am reluctant to add a GUI option because nsis is primitive, and adding gui is no fun (I use it only because it is open source and has a plug-in system; the scripting language to build the binary is awful). Also, the nsis installer itself has no GUI on purpose, to avoid confusing people with a two stages installer. cheers, David ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] profiling line by line
On Sat, Sep 27, 2008 at 00:50, Fernando Perez <[EMAIL PROTECTED]> wrote: > Hey, > > On Fri, Sep 26, 2008 at 5:36 PM, Robert Kern <[EMAIL PROTECTED]> wrote: > >> There is now! Add >> >>import line_profiler >>ip.expose_magic('lprun', line_profiler.magic_lprun) >> >> to your ipy_user_conf.py . > > Mmh, should we ship this out of the box in ipython? The C dependency > is the only thing that gives me a bit of pause. There goes our pure > C, easy to build binaries for... On the other hand, me wants :) > > What's your take on that? I'd prefer not to. The exclusion of C modules is an eminently reasonable one. You would inevitably make it optional, and my opinion is that optional functionality is best handled by separate packages. Feel free to add it as "highly recommended" in the IPython installation instructions, though. :-) -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Minimum distance between 2 paths in 3D
Hi All, I was wondering if someone had any suggestions/references/snippets of code on how to find the minimum distance between 2 paths in 3D. Basically, for every path, I have I series of points (x, y, z) and I would like to know if there is a better way, other than comparing point by point the 2 paths, to find the minimum distance between them. Any suggestion is welcome :-) Thank you. Andrea. "Imagination Is The Only Weapon In The War Against Reality." http://xoomer.alice.it/infinity77/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Minimum distance between 2 paths in 3D
Hi Andrea, >I was wondering if someone had any suggestions/references/snippets > of code on how to find the minimum distance between 2 paths in 3D. > Basically, for every path, I have I series of points (x, y, z) and I > would like to know if there is a better way, other than comparing > point by point the 2 paths, to find the minimum distance between them. In 2D there would be a few tricks you could use, but in higher dimensions anything smart that you could attempt might cost you more computation time than just comparing the points (unless N is huge). At least make sure to put the "looping" over points into a vectorized form to avoid python for loops. e.g. two curves given by 3xN arrays c and d: from numpy import concatenate, argmin from numpy.linalg import norm distvec = concatenate([c[:,i]-d.T for i in range(N)]) # all N**2 distance vectors ns = [norm(a) for a in distvec] # all N**2 norms of the distance vectors cix, dix = divmod(argmin(ns), N) # find the index of the minimum norm from [0 .. N**2] and decode which point this corresponds to The minimum is between the points c[:,cix] and d[:,dix]. A numpy wonk might be able to squeeze a bit more optimization out of this, but I think this code works OK. Unfortunately, unless you know more mathematical properties of your curves in advance (info about their maximum curvature, for instance) you'll end up needing to check every pair of points. If N is really big, like over 10**4 maybe, it might be worth trying to break the curves up into pieces contained in bounding boxes which you can eliminate from a full search if they don't intersect. -Rob ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] numpy.random.logseries - incorrect convergence for k=1, k=2
random numbers generated by numpy.random.logseries do not converge to theoretical distribution: for probability paramater pr = 0.8, the random number generator converges to a frequency for k=1 at 39.8 %, while the theoretical probability mass is 49.71 k=2 is oversampled, other k's look ok check frequency of k=1 and k=2 at N = 100 0.398406 0.296465 pmf at k = 1 and k=2 with formula [ 0.4971 0.1988] for probability paramater pr = 0.3, the results are not as bad, but still off: frequency for k=1 at 82.6 %, while the theoretical probability mass is 84.11 check frequency of k=1 and k=2 at N = 100 0.826006 0.141244 pmf at k = 1 and k=2 with formula [ 0.8411 0.1262] below is a quick script for checking this Josef {{{ import numpy as np from scipy import stats pr = 0.8 np.set_printoptions(precision=2, suppress=True) # calculation for N=1million takes some time for N in [1000, 1, 1, 100]: rvsn=np.random.logseries(pr,size=N) fr=stats.itemfreq(rvsn) pmfs=stats.logser.pmf(fr[:,0],pr)*100 print 'log series sample frequency and pmf (in %) with N = ', N print np.column_stack((fr[:,0],fr[:,1]*100.0/N,pmfs)) np.set_printoptions(precision=4, suppress=True) print 'check frequency of k=1 and k=2 at N = ', N print np.sum(rvsn==1)/float(N), print np.sum(rvsn==2)/float(N) k = np.array([1,2]) print 'pmf at k = 1 and k=2 with formula' print -pr**k * 1.0 / k / np.log(1-pr) }}} ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] profiling line by line
On Sat, Sep 27, 2008 at 1:31 AM, Robert Kern <[EMAIL PROTECTED]> wrote: >> Mmh, should we ship this out of the box in ipython? The C dependency >> is the only thing that gives me a bit of pause. There goes our pure >> C, easy to build binaries for... On the other hand, me wants :) s/pure C/pure Python/ grr... >> What's your take on that? > > I'd prefer not to. The exclusion of C modules is an eminently > reasonable one. You would inevitably make it optional, and my opinion > is that optional functionality is best handled by separate packages. Agreed. Not having C dependencies simplifies lots of things for us. If we ever gain an unavoidable one, we can then revisit this. > Feel free to add it as "highly recommended" in the IPython > installation instructions, though. :-) Most certainly will. It would be better to have a more stable release/download link than your hg repo. Do you have a page or a Pypi link for it yet? A little readme? Cheers, f ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Interesting bug with numpy.ma
All, I've just been pointed to a strange bug in numpy.ma import numpy as np import numpy.ma as ma assert(ma.masked*float(1) is ma.masked) assert(float(1)*ma.masked is ma.masked) assert(ma.masked*np.float(1) is ma.masked) assert(np.float(1)*ma.masked is ma.masked) assert(ma.masked*np.float_(1) is ma.masked) assert(np.float_(1)*ma.masked is ma.masked) The last assertion fails, and I have no idea why. What's the difference between assertions 3 and 5 ? FYI, ma.masked=ma.MaskedArray(0.,dtype=np.float_,mask=True), and MaskedArray has a __array_priority__ of 15. Implementing a __rmul__ doesn't work either. Any idea welcome. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Cython docs including numpy info
Hi all, I just found out something that might interest others on this list. Gabriel Gellner (forgive me if I'm not duly crediting others, I don't know the detailed history of this) has been leading the effort to produce self-contained Cython documentation, in sphix, that includes the new numpy support added recently to Cython: http://www.mudskipper.ca/cython-doc/ I'm looking at these docs right now, but I'm sure Gabriel and the rest of the Cython team would benefit from more feedback. Cheers, f ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Minimum distance between 2 paths in 3D
2008/9/27 Andrea Gavana <[EMAIL PROTECTED]>: >I was wondering if someone had any suggestions/references/snippets > of code on how to find the minimum distance between 2 paths in 3D. > Basically, for every path, I have I series of points (x, y, z) and I > would like to know if there is a better way, other than comparing > point by point the 2 paths, to find the minimum distance between them. If you treat the points as simply two clouds of points (that is, ignore the fact that they represent points on paths), spatial data structures can make this kind of question faster. For example a kd-tree can be constructed in something like O(n log n) time, and you can answer questions like "what is the closest point in this set to (a,b,c)?" in something like O(log n) time. So you could do this by putting one set of points in a kd-tree and just running queries against each point in the other. There exist other spatial data structures - octrees, voxel grids, bounding-volume hierarchies, other binary space partitioning trees, as well as various more specialized ones. As the number of dimensions becomes large they become substantially more difficult to work with, and you do need to balance construction time against lookup time, but they are definitely worth thinking about in your problem. Sadly, no such data structure exists in either numpy or scipy, though biopython apparently has an implementation of kd-trees. Good luck, Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ANN: NumPy 1.2.0
David Cournapeau schrieb: > Thomas Heller wrote: >> >> Well, the first question is: What does happen when I install the SSE3 >> version >> (or how it's called) on my machine, use py2exe to build an app, and this >> app runs on a SSE2 machine - degraded performance, or hard crashes? > > Hard crash. That's the whole point of the installer, actually: install > the right version. The problem is two fold: > - it is practically impossible to force ATLAS to configure itself to > use an instruction set different than the one used to build ATLAS. > - we don't have a system to dynamically load the right ATLAS when > importing numpy. > > Note that this is not new: you had the problem before, because before, > numpy was *only* built with SSE2 support, and any machine wo SSE2 would > crash when using numpy. I see. >> So, maybe the gui could allow to select whether to install the >> high-performance >> version specialized for the current cpu, or a more portable but a little bit >> slower version (I assume there is one included) that can be safely used for >> py2exe. > > Is it really complicated to decompress the .exe to get the installers > and choose the one you want ? I am reluctant to add a GUI option because > nsis is primitive, and adding gui is no fun (I use it only because it is > open source and has a plug-in system; the scripting language to build > the binary is awful). Also, the nsis installer itself has no GUI on > purpose, to avoid confusing people with a two stages installer. No, it isn't complicated. I searched a little about command line options, and didn't find any. But 7-zip did unpack the installers. I can live with that. -- Thanks, Thomas ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Minimum distance between 2 paths in 3D
On Sep 27, 2008, at 10:23 PM, Anne Archibald wrote: > Sadly, no such data structure exists > in either numpy or scipy, though biopython apparently has an > implementation of kd-trees. It's not part of scipy but I saw there was a scikit which supported ANN: http://scipy.org/scipy/scikits/wiki/AnnWrapper scikits.ann exposes a single class, kdtree that wraps the Approximate Nearest Neighbor library's kd-tree implementation. kdtree has a single (non-constructor) method, knn that finds the indecies (of the points used to construct the kdtree) of the k-nearest neighbors and the squared distances to those points Does anyone have any experience with it? Though perhaps I should ask on the scipy list instead of numpy Andrew [EMAIL PROTECTED] ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Minimum distance between 2 paths in 3D
On Sat, Sep 27, 2008 at 15:23, Anne Archibald <[EMAIL PROTECTED]> wrote: > 2008/9/27 Andrea Gavana <[EMAIL PROTECTED]>: > >>I was wondering if someone had any suggestions/references/snippets >> of code on how to find the minimum distance between 2 paths in 3D. >> Basically, for every path, I have I series of points (x, y, z) and I >> would like to know if there is a better way, other than comparing >> point by point the 2 paths, to find the minimum distance between them. > > If you treat the points as simply two clouds of points (that is, > ignore the fact that they represent points on paths), spatial data > structures can make this kind of question faster. For example a > kd-tree can be constructed in something like O(n log n) time, and you > can answer questions like "what is the closest point in this set to > (a,b,c)?" in something like O(log n) time. So you could do this by > putting one set of points in a kd-tree and just running queries > against each point in the other. There exist other spatial data > structures - octrees, voxel grids, bounding-volume hierarchies, other > binary space partitioning trees, as well as various more specialized > ones. As the number of dimensions becomes large they become > substantially more difficult to work with, and you do need to balance > construction time against lookup time, but they are definitely worth > thinking about in your problem. Sadly, no such data structure exists > in either numpy or scipy, though biopython apparently has an > implementation of kd-trees. Sturla Molden has just contributed a pure Python+numpy implementation on the Scipy Cookbook. http://www.scipy.org/Cookbook/KDTree -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Minimum distance between 2 paths in 3D
2008/9/27 Robert Kern <[EMAIL PROTECTED]>: > On Sat, Sep 27, 2008 at 15:23, Anne Archibald <[EMAIL PROTECTED]> wrote: >> 2008/9/27 Andrea Gavana <[EMAIL PROTECTED]>: >> >>>I was wondering if someone had any suggestions/references/snippets >>> of code on how to find the minimum distance between 2 paths in 3D. >>> Basically, for every path, I have I series of points (x, y, z) and I >>> would like to know if there is a better way, other than comparing >>> point by point the 2 paths, to find the minimum distance between them. >> >> If you treat the points as simply two clouds of points (that is, >> ignore the fact that they represent points on paths), spatial data >> structures can make this kind of question faster. For example a >> kd-tree can be constructed in something like O(n log n) time, and you >> can answer questions like "what is the closest point in this set to >> (a,b,c)?" in something like O(log n) time. So you could do this by >> putting one set of points in a kd-tree and just running queries >> against each point in the other. There exist other spatial data >> structures - octrees, voxel grids, bounding-volume hierarchies, other >> binary space partitioning trees, as well as various more specialized >> ones. As the number of dimensions becomes large they become >> substantially more difficult to work with, and you do need to balance >> construction time against lookup time, but they are definitely worth >> thinking about in your problem. Sadly, no such data structure exists >> in either numpy or scipy, though biopython apparently has an >> implementation of kd-trees. > > Sturla Molden has just contributed a pure Python+numpy implementation > on the Scipy Cookbook. > > http://www.scipy.org/Cookbook/KDTree I think a kd-tree implementation would be a valuable addition to scipy, perhaps in a submodule scipy.spatial that might eventually contain other spatial data structures and algorithms. What do you think? Should we have one? Should it be based on Sturla Molden's code, if the license permits? I am willing to contribute one, if not. Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Minimum distance between 2 paths in 3D
On Sat, Sep 27, 2008 at 11:18 PM, Anne Archibald <[EMAIL PROTECTED]> wrote: > > I think a kd-tree implementation would be a valuable addition to > scipy, perhaps in a submodule scipy.spatial that might eventually > contain other spatial data structures and algorithms. What do you > think? Should we have one? Should it be based on Sturla Molden's code, > if the license permits? I am willing to contribute one, if not. +1 If you're implementing one, I would highly recommend the "left-balanced" kd-tree. http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/2535/pdf/imm2535.pdf Here's one implementation (undetermined license): http://www2.imm.dtu.dk/~jab/software.html#kdtrees I have a C++ implementation based on the one above that does a few more operations (like nearest-n neighbors) that you're welcome to use. The use of n_th() in the STL makes left-balancing fairly straightforward. FWIW I also have a pure python implementation here: http://code.google.com/p/pydec/source/browse/trunk/pydec/math/kd_tree.py -- Nathan Bell [EMAIL PROTECTED] http://graphics.cs.uiuc.edu/~wnbell/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] nonuniform scatter operations
Hello, Is there an efficient way to implement a nonuniform gather operation in numpy? Specifically, I want to do something like n,m = 100,1000 X = random.uniform(size=n) K = random.randint(n, size=m) Y = random.uniform(size=m) for k,y in zip(K,Y): X[k] += y but I want it to be fast. The naive attempt "X[K] += Y" does not work, since the slice assumes the indices don't repeat. Thanks, Geoffrey ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] nonuniform scatter operations
On Sun, Sep 28, 2008 at 12:34 AM, Geoffrey Irving <[EMAIL PROTECTED]> wrote: > > Is there an efficient way to implement a nonuniform gather operation > in numpy? Specifically, I want to do something like > > n,m = 100,1000 > X = random.uniform(size=n) > K = random.randint(n, size=m) > Y = random.uniform(size=m) > > for k,y in zip(K,Y): >X[k] += y > > but I want it to be fast. The naive attempt "X[K] += Y" does not > work, since the slice assumes the indices don't repeat. > I don't know of numpy solution, but in scipy you could use a sparse matrix to perform the operation. I think the following does what you want. from scipy.sparse import coo_matrix X += coo_matrix( (Y, (K,zeros(m,dtype=int)), shape=(n,1)).sum(axis=1) This reduces to a simple C++ loop, so speed should be good: http://projects.scipy.org/scipy/scipy/browser/trunk/scipy/sparse/sparsetools/coo.h#L139 -- Nathan Bell [EMAIL PROTECTED] http://graphics.cs.uiuc.edu/~wnbell/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ANN: NumPy 1.2.0
Thomas Heller wrote: > > No, it isn't complicated. I searched a little about command line options, > and didn't find any. But 7-zip did unpack the installers. I can live with > that. I will look into adding a command line switch. I don't want to add a GUI, but adding some command line switches for advanced uses and scripting makes a lot of sense, without risks of disrupting the whole install "experience". cheers, David ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] nonuniform scatter operations
2008/9/28 Geoffrey Irving <[EMAIL PROTECTED]>: > Is there an efficient way to implement a nonuniform gather operation > in numpy? Specifically, I want to do something like > > n,m = 100,1000 > X = random.uniform(size=n) > K = random.randint(n, size=m) > Y = random.uniform(size=m) > > for k,y in zip(K,Y): >X[k] += y > > but I want it to be fast. The naive attempt "X[K] += Y" does not > work, since the slice assumes the indices don't repeat. I believe histogram can be persuaded to do this. Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion