Hello, I am relatively new to numpy and am having trouble with the speed of a specific array based calculation that I'm trying to do. What I'm trying to do is to calculate the total total potential energy and coordination number of each atom within a relatively large simulation. Each atom is at a position (x,y,z) given by a row in a large array (approximately 1e6 by 3) and presently I have no information about its nearest neighbours so each its position must be checked against all others before cutting the list down prior to calculating the energy. My current numpy code is below and in it I have tried to push as much of the work for this computation into compiled C (ideally) of numpy. However it is still taking a very long time at approx 1 second per row. At this speed even some simple speed ups like weave won't really help. Are there any other numpy ways that it would be possible to calculate this, or should I start looking at going to C/C++? I am also left wondering if there is something wrong with my installation of numpy (i.e. not making use of lapack/ATLAS). Other than that if it is relevant - I am running 32 bit x86 Centos 5.1 linux on a dual Xeon 3.0 GHz Dell tower with 4 GB of memory. Any suggestions of things to try would be great. Dan ======================================================================== ======================= def epairANDcoord(xyz,cutoff=4.0): box_length = 160.0 cut2 = cutoff**2 ljTotal = numpy.zeros(xyz.shape[0]) coord = numpy.zeros(xyz.shape[0]) i = 0 for r0 in xyz: r2 = r0-xyz # application of minimum image convention so each atom interacts with nearest periodic image w/i box_length/2 r2 -= box_length*numpy.rint(r2/box_length) r2 = numpy.power(r2,2).sum(axis=1) r2 = numpy.extract(r2<cut2, r2) coord[i] = r2.shape[0]-1 r2 = numpy.extract(r2 != 0, r2) # avoid problems of inf when considering own atom. r2 = numpy.divide(1.0,r2) ljTotal[i] = numpy.multiply(2,(numpy.power(r2,6)- numpy.power(r2,3))).sum(axis=0) # ljTotal[i] = ljPerAtom.sum(axis=0) i += 1 del r2 # if i>25: break # used to cut down the length of calculation as needed return ljTotal,coord
On Thu, Jul 10, 2008 at 10:38 AM, Dan Lussier <dtlussier@gmail.com> wrote:
Hello,
I am relatively new to numpy and am having trouble with the speed of a specific array based calculation that I'm trying to do.
What I'm trying to do is to calculate the total total potential energy and coordination number of each atom within a relatively large simulation. Each atom is at a position (x,y,z) given by a row in a large array (approximately 1e6 by 3) and presently I have no information about its nearest neighbours so each its position must be checked against all others before cutting the list down prior to calculating the energy.
This looks to be O(n^2) and might well be the bottle neck. There are various ways to speed up such things but more information would help determine the method, i.e., is this operation within a loop so that the values change a lot. However, one quick thing to try is a sort on one of the coordinates so you only need to check a subset of the vectors. Searchsorted could be useful here also. Chuck
Charles R Harris wrote:
On Thu, Jul 10, 2008 at 10:38 AM, Dan Lussier <dtlussier@gmail.com <mailto:dtlussier@gmail.com>> wrote:
Hello,
I am relatively new to numpy and am having trouble with the speed of a specific array based calculation that I'm trying to do.
What I'm trying to do is to calculate the total total potential energy and coordination number of each atom within a relatively large simulation. Each atom is at a position (x,y,z) given by a row in a large array (approximately 1e6 by 3) and presently I have no information about its nearest neighbours so each its position must be checked against all others before cutting the list down prior to calculating the energy.
This looks to be O(n^2) and might well be the bottle neck. There are various ways to speed up such things but more information would help determine the method, i.e., is this operation within a loop so that the values change a lot. However, one quick thing to try is a sort on one of the coordinates so you only need to check a subset of the vectors. Searchsorted could be useful here also.
Chuck
------------------------------------------------------------------------
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
If I understand correctly, I notice that you are doing many scalar multiplications where some are the same or constants. For example, you compute r2*r2 a few times: numpy.power(r2,2)= r2*r2 numpy.power(r2,3)= r2*r2*r2 numpy.power(r2,6)= r2*r2*r2*r2*r2*r2 But you don't need this numpy.power(r2,6) as it can be factored. Also the division is unnecessary. So the liTotal is really: ljTotal[i] = 2.0*((numpy.power(r2,-3)*(numpy.power(r2,-3)-1)).sum(axis=0)) However, the real problem is getting the coordinates that I do not follow and limits how you can remove the loop over xyz. Like why the need for the repeated (r0-xyz)? This is a huge cost that you need to avoid perhaps by just extracting those elements within your criteria. Bruce
2008/7/10 Dan Lussier <dtlussier@gmail.com>:
I am relatively new to numpy and am having trouble with the speed of a specific array based calculation that I'm trying to do.
What I'm trying to do is to calculate the total total potential energy and coordination number of each atom within a relatively large simulation. Each atom is at a position (x,y,z) given by a row in a large array (approximately 1e6 by 3) and presently I have no information about its nearest neighbours so each its position must be checked against all others before cutting the list down prior to calculating the energy.
The way you've implemented this goes as the square of the number of atoms. This is going to absolutely kill performance, and you can spend weeks trying to optimize the code for a factor of two speedup. I would look hard for algorithms that do this in less than O(n^2). This problem of finding the neighbours of an object has seen substantial research, and there are a variety of well-established techniques covering many possible situations. My knowledge of it is far from up-to-date, but since you have only a three-dimensional problem, you could try a three-dimensional grid (if your atoms aren't too clumpy) or octrees (if they are somewhat clumpy); kd-trees are probably overkill (since they're designed for higher-dimensional problems).
My current numpy code is below and in it I have tried to push as much of the work for this computation into compiled C (ideally) of numpy. However it is still taking a very long time at approx 1 second per row. At this speed even some simple speed ups like weave won't really help.
Are there any other numpy ways that it would be possible to calculate this, or should I start looking at going to C/C++?
I am also left wondering if there is something wrong with my installation of numpy (i.e. not making use of lapack/ATLAS). Other than that if it is relevant - I am running 32 bit x86 Centos 5.1 linux on a dual Xeon 3.0 GHz Dell tower with 4 GB of memory.
Unfortunately, implementing most of the algorithms in the literature within numpy is going to be fairly cumbersome. But I think there are better algorithms you could try: * Put all your atoms in a grid. Ideally the cells would be about the size of your cutoff radius, so that for each atom you need to pull out all atoms in at most eight cells for checking against the cutoff. I think this can be done fairly efficiently in numpy. * Sort the atoms along a coordinate and use searchsorted to pull out those for which that coordinate is within the cutoff. This should get you down to reasonably modest lists fairly quickly. There is of course a tradeoff between preprocessing time and lookup time. We seem to get quite a few posts from people wanting some kind of spatial data structure (whether they know it or not). Would it make sense to come up with some kind of compiled spatial data structure to include in scipy? Anne
On Thu, Jul 10, 2008 at 2:55 PM, Anne Archibald <peridot.faceted@gmail.com> wrote:
2008/7/10 Dan Lussier <dtlussier@gmail.com>:
I am relatively new to numpy and am having trouble with the speed of a specific array based calculation that I'm trying to do.
What I'm trying to do is to calculate the total total potential energy and coordination number of each atom within a relatively large simulation. Each atom is at a position (x,y,z) given by a row in a large array (approximately 1e6 by 3) and presently I have no information about its nearest neighbours so each its position must be checked against all others before cutting the list down prior to calculating the energy.
The way you've implemented this goes as the square of the number of atoms. This is going to absolutely kill performance, and you can spend weeks trying to optimize the code for a factor of two speedup. I would look hard for algorithms that do this in less than O(n^2).
This problem of finding the neighbours of an object has seen substantial research, and there are a variety of well-established techniques covering many possible situations. My knowledge of it is far from up-to-date, but since you have only a three-dimensional problem, you could try a three-dimensional grid (if your atoms aren't too clumpy) or octrees (if they are somewhat clumpy); kd-trees are probably overkill (since they're designed for higher-dimensional problems).
My current numpy code is below and in it I have tried to push as much of the work for this computation into compiled C (ideally) of numpy. However it is still taking a very long time at approx 1 second per row. At this speed even some simple speed ups like weave won't really help.
Are there any other numpy ways that it would be possible to calculate this, or should I start looking at going to C/C++?
I am also left wondering if there is something wrong with my installation of numpy (i.e. not making use of lapack/ATLAS). Other than that if it is relevant - I am running 32 bit x86 Centos 5.1 linux on a dual Xeon 3.0 GHz Dell tower with 4 GB of memory.
Unfortunately, implementing most of the algorithms in the literature within numpy is going to be fairly cumbersome. But I think there are better algorithms you could try:
* Put all your atoms in a grid. Ideally the cells would be about the size of your cutoff radius, so that for each atom you need to pull out all atoms in at most eight cells for checking against the cutoff. I think this can be done fairly efficiently in numpy.
* Sort the atoms along a coordinate and use searchsorted to pull out those for which that coordinate is within the cutoff. This should get you down to reasonably modest lists fairly quickly.
There is of course a tradeoff between preprocessing time and lookup time.
We seem to get quite a few posts from people wanting some kind of spatial data structure (whether they know it or not). Would it make sense to come up with some kind of compiled spatial data structure to include in scipy?
I think there should be a "computer science" module containing such things as r-b trees, k-d trees, equivalence relations, and so forth. For this problem one could probably make a low resolution cube of list objects and index atoms into the appropriate list by severe rounding. I have done similar things for indexing stars in a (fairly) wide field image and it worked well. But I think the sorting approach would be a good first try here. Argsort on the proper column followed by take would be the way to go for that. Chuck
On Fri, Jul 11, 2008 at 5:55 AM, Anne Archibald <peridot.faceted@gmail.com> wrote:
2008/7/10 Dan Lussier <dtlussier@gmail.com>:
We seem to get quite a few posts from people wanting some kind of spatial data structure (whether they know it or not). Would it make sense to come up with some kind of compiled spatial data structure to include in scipy?
Yes! http://en.wikipedia.org/wiki/Nearest_neighbor_search I've been using cover trees recently. http://hunch.net/~jl/projects/cover_tree/cover_tree.html Though I haven't hammered on it very hard, it seems to be working well so far. --bb
2008/7/10 Anne Archibald <peridot.faceted@gmail.com>:
Unfortunately, implementing most of the algorithms in the literature within numpy is going to be fairly cumbersome. But I think there are better algorithms you could try:
There's also http://pypi.python.org/pypi/scikits.ann Regards Stéfan
Anne Archibald wrote:
Unfortunately, implementing most of the algorithms in the literature within numpy is going to be fairly cumbersome.
Maybe this rtree implementation would help: http://pypi.python.org/pypi/Rtree -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov
On Jul 10, 2008, at 6:38 PM, Dan Lussier wrote:
What I'm trying to do is to calculate the total total potential energy and coordination number of each atom within a relatively large simulation.
Anne Archibald already responded:
you could try a three-dimensional grid (if your atoms aren't too clumpy) or octrees (if they are somewhat clumpy); kd-trees are probably overkill (since they're designed for higher-dimensional problems).
I implemented something like what you did in VMD about 14(!) years ago. (VMD is a molecular structure visualization program.) I needed to find which atoms might be bonded to each other. I assume you have a cutoff value, which means you don't need to worry about the general nearest-neighbor problem. Molecules are nice because the distributions are not clumpy. Atoms don't get that close to nor that far from other atoms. I implemented a grid. It goes something like this: import collections # Search within 3 A d = 3.0 coordinates = [ ("C1", 34.287, 50.970, 115.006), ("O1", 34.972, 51.144, 113.870), ("C2", 33.929, 52.255, 115.739), ("N2", 34.753, 52.387, 116.954), ("C3", 32.448, 52.219, 116.121), ("O3", 32.033, 50.877, 116.336), ("C4", 31.528, 52.817, 115.054), ("C5", 30.095, 53.013, 115.558), ("C6", 29.226, 53.835, 114.609), ("C7", 29.807, 55.217, 114.304), ("C8", 29.092, 55.920, 113.147), ("C9", 29.525, 57.375, 112.971), ("C10", 28.409, 58.267, 112.422), ("C11", 28.828, 59.734, 112.294), ("C12", 27.902, 60.542, 111.385), ("C13", 26.617, 60.996, 112.085), ("C14", 26.182, 62.401, 111.667), ("C15", 24.739, 62.731, 112.054), ("C16", 24.251, 64.046, 111.441), ("C17", 23.026, 64.624, 112.150), ("C18", 22.631, 66.007, 111.623),] def dict_of_list(): return collections.defaultdict(list) def dict_of_dict(): return collections.defaultdict(dict_of_list) grid = collections.defaultdict(dict_of_dict) for atom in coordinates: name,x,y,z = atom i = int(x/d) j = int(y/d) k = int(z/d) grid[i][j][k].append(atom) query_name, query_x, query_y, query_z = coordinates[8] query_i = int(query_x/d) query_j = int(query_y/d) query_k = int(query_z/d) # Given the search distance 'd', only need to check cells up to 1 unit away within_names = set() for i in range(query_i-1, query_i+2): for j in range(query_j-1, query_j+2): for k in range(query_k-1, query_k+2): for atom in grid[i][j][k]: name, x, y, z = atom print "Check", atom, query_d2 = (x-query_x)**2+(y-query_y)**2+(z-query_z)**2 if query_d2 < d*d: print "Within!", query_d2 within_names.add(name) else: print "Too far", query_d2 print len(within_names), "were close enough" # Linear search to verify count = 0 for name, x, y, z in coordinates: query_d2 = (x-query_x)**2+(y-query_y)**2+(z-query_z)**2 if query_d2 < d*d: print "Within", name if name not in within_names: raise AssertionError(name) count += 1 if count != len(within_names): raise AssertionError("count problem") You can also grab the KDTree from Biopython, which is implemented in C. http://www.biopython.org/DIST/docs/api/Bio.KDTree.KDTree'-module.html It was designed for just this task. Andrew dalke@dalkescientific.com
On Thu, Jul 10, 2008 at 5:48 PM, Andrew Dalke <dalke@dalkescientific.com> wrote: <snip>
You can also grab the KDTree from Biopython, which is implemented in C.
It was designed for just this task.
Looks nice, has a BSD type license, but uses Numeric :( Oh well, a little fixing up should do the trick. Chuck
Hi, scikits.learn.machine/manifold_learning.regression.neighbor contains a kd-tree for neighbors search as well (implemented in C++). Matthieu 2008/7/11 Charles R Harris <charlesr.harris@gmail.com>:
On Thu, Jul 10, 2008 at 5:48 PM, Andrew Dalke <dalke@dalkescientific.com> wrote:
<snip>
You can also grab the KDTree from Biopython, which is implemented in C.
http://www.biopython.org/DIST/docs/api/Bio.KDTree.KDTree'-module.html
It was designed for just this task.
Looks nice, has a BSD type license, but uses Numeric :( Oh well, a little fixing up should do the trick.
Chuck
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
-- French PhD student Website : http://matthieu-brucher.developpez.com/ Blogs : http://matt.eifelle.com and http://blog.developpez.com/?blog=92 LinkedIn : http://www.linkedin.com/in/matthieubrucher
If your positions are static (I'm not clear on that from your message), then you might want to check the technique of "slice searching". It only requires one sort of the data for each dimension initially, then uses a simple, but clever look up to find neighbors within some epsilon of a chosen point. Speeds appear to be about equal to k-d trees. Programming is vastly simpler than k-d trees, however. See, [1] "A Simple Algorithm for Nearest Neighbor Search in High Dimensions," Sameer A. Nene and Shree K. Nayar, IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE 19 (9), 989 (1997). -- Lou Pecora, my views are my own. --- On Thu, 7/10/08, Dan Lussier <dtlussier@gmail.com> wrote:
From: Dan Lussier <dtlussier@gmail.com> Subject: [Numpy-discussion] huge array calculation speed To: numpy-discussion@scipy.org Date: Thursday, July 10, 2008, 12:38 PM Hello,
I am relatively new to numpy and am having trouble with the speed of a specific array based calculation that I'm trying to do.
What I'm trying to do is to calculate the total total potential energy and coordination number of each atom within a relatively large simulation. Each atom is at a position (x,y,z) given by a row in a large array (approximately 1e6 by 3) and presently I have no information about its nearest neighbours so each its position must be checked against all others before cutting the list down prior to calculating the energy.
On Fri, Jul 11, 2008 at 11:04 AM, Lou Pecora <lou_boog2000@yahoo.com> wrote:
If your positions are static (I'm not clear on that from your message), then you might want to check the technique of "slice searching". It only requires one sort of the data for each dimension initially, then uses a simple, but clever look up to find neighbors within some epsilon of a chosen point. Speeds appear to be about equal to k-d trees. Programming is vastly simpler than k-d trees, however.
This one is actually easy to implement in numpy using argsort. I'm not sure how much speed the integer comparisons buy as opposed to straight floating comparisons; they probably did it that way for the hardware implementation. It might be interesting to make a comparison. Chuck
[ Sorry I never sent this, I just found it in my drafts folder. Just in case it's useful to the OP, here it is. ] On Thu, Jul 10, 2008 at 9:38 AM, Dan Lussier <dtlussier@gmail.com> wrote:
r2 = numpy.power(r2,2).sum(axis=1)
r2 = numpy.extract(r2<cut2, r2) coord[i] = r2.shape[0]-1 r2 = numpy.extract(r2 != 0, r2) # avoid problems of inf when considering own atom. r2 = numpy.divide(1.0,r2) ljTotal[i] = numpy.multiply(2,(numpy.power(r2,6)- numpy.power(r2,3))).sum(axis=0)
That's a Lennard-Jones potential, and there seems to be a fair amount of work on FMM (fast multipole method) for that type of potential. Depending on what the final set of quantities you're after is, the FMM may be worth looking into. Cheers, f
participants (11)
-
Andrew Dalke -
Anne Archibald -
Bill Baxter -
Bruce Southey -
Charles R Harris -
Christopher Barker -
Dan Lussier -
Fernando Perez -
Lou Pecora -
Matthieu Brucher -
Stéfan van der Walt