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
Re: [Numpy-discussion] Minimum distance between 2 paths in 3D
Hi Rob All, On Sat, Sep 27, 2008 at 4:05 PM, Rob Clewley wrote: 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. Thank you very much for your kind suggestions and the code snippet. I think it will do just fine, as my well trajectories may have from 20 to 500 points maximum. I'll try it tomorrow at work and see how it goes. Thank you again. 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
[Numpy-discussion] Vectorization of the product of several matrices ?
Hello, I have two lists of numpy matrices : LM = [M_i, i=1..N] and LN = [N_i, i =1..N] and I would like to compute the list of the products : LP = [M_i * N_i, i=1..N]. I can do : P=[] for i in range(N) : P.append(LM[i]*LN[i]) But this is not vectorized. Is there a faster solution ? Can I vectorize this operation ? How ? Will it be more efficient than the code above ? Thank you for your help, -- O.C. Créez votre adresse électronique [EMAIL PROTECTED] 1 Go d'espace de stockage, anti-spam et anti-virus intégrés. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Help needed with the priority of 0d-arrays and np.scalars
All, In a recent post, I was pointing to a bug in numpy.ma. I think I can understand what happens, but I'd need a couple of clarifications. When multiplying a 0d ndarrray or one of its subclass by a numpy scalar, we end up with a numpy scalar cast to the highest type, right ? So, np.float64(1)*np.array(0.,dtype='g') is a np.float128 object, while np.float64(1)*np.array(0,dtype=np.float32) is a np.float64 object. What surprises me is that numpy scalars have a low priority (listed as -100 or -100 depending on the page of the numpy book), so if I have a 0d subclass of ndarray with a higher priority, I should end up with a 0d subclass of ndarray, but it doesn't seem to be the case. The problem comes and bites me with ma.masked, defined as a 0d MaskedArray with __array_priority__ of 15 and a np.float64 dtype. Doing ma.masked*np.float64(1.) yields ma.masked, as it should. Doing np.float64(1.)*ma.masked yields np.float64(0.) but shouldn't: I should access instead the __rmul__ method of ma.masked. Any comment/suggestion would be welcome. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] nonuniform scatter operations
On Sat, Sep 27, 2008 at 10:01 PM, Nathan Bell [EMAIL PROTECTED] wrote: 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 Thanks. That works great. A slightly cleaner version is X += coo_matrix((Y, (K, zeros_like(K.sum(axis=1) The next question is: is there a similar way that generalizes to the case where X is n by 3 and Y is m by 3 (besides the obvious loop over range(3), that is)? Geoffrey ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
