Re: [Numpy-discussion] Iterative Matrix Multiplication
2010/3/5 Ian Mallett geometr...@gmail.com: Cool--this works perfectly now :-) :-) Unfortunately, it's actually slower :P Most of the slowest part is in the removing doubles section. Hmm. Let's see ... Can you tell me how I can test the time calls in a script take? I have no idea. #takes 0.04 seconds inner = np.inner(ns, v1s - some_point) I think I can do nothing about that at the moment. #0.0840001106262 sum_1 = sum.reshape((len(sum), 1)).repeat(len(sum), axis = 1) #0.032923706 sum_2 = sum.reshape((1, len(sum))).repeat(len(sum), axis = 0) #0.026504089 comparison_sum = (sum_1 == sum_2) We can leave out the repeat() calls and leave only the reshape() calls there. Numpy will substitute dimi == 1 dimensions with stride == 0, i.e., it will effectively repeat those dimension, just as we did it explicitly. #0.0909998416901 diff_1 = diff.reshape((len(diff), 1)).repeat(len(diff), axis = 1) #0.0340001583099 diff_2 = diff.reshape((1, len(diff))).repeat(len(diff), axis = 0) #0.026504089 comparison_diff = (diff_1 == diff_2) Same here. Delete the repeat() calls, but not the reshape() calls. #0.023019073 same_edges = comparison_sum * comparison_diff Hmm, maybe use numpy.logical_and(comparison_sum, comparison_diff)? I don't know, but I guess it is in some way optimised for such things. #0.12848502 doublet_count = same_edges.sum(axis = 0) Maybe try axis = 1 instead. I wonder why this is so slow. Or maybe it's because he does the conversion to ints on-the-fly, so maybe try same_edges.astype(numpy.int8).sum(axis = 0). Hope this gives some improvement. I attach the modified version. Ah, one thing to mention, have you not accidentally timed also the printout functions? They should be pretty slow. Friedrich ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
2010/3/5 Ian Mallett geometr...@gmail.com: #takes 0.04 seconds inner = np.inner(ns, v1s - some_point) Ok, I don't know why I was able to overlook this: dotprod = (ns * (v1s - some_point)).sum(axis = 1) The things with the inner product have been deleted. Now I will really *attach* it ... Hope it's faster, Friedrich shading.py Description: Binary data ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Building Numpy Windows Superpack
On Sat, Mar 6, 2010 at 1:29 AM, David Cournapeau courn...@gmail.com wrote: On Fri, Mar 5, 2010 at 1:22 PM, Patrick Marsh patrickmars...@gmail.com wrote: I've run the Numpy superpack installer for Python 2.6 built with MinGW through the dependency walker. Unfortunately, outside of checking for some extremely obviously things, I'm in way over my head in interpreting the output (although, I'd like to learn). I've put the output from the program here: http://www.patricktmarsh.com/numpy/20100303.py26.superpack.dependencies.txt . I can also put the binary up somewhere, too if someone wants to check that. I have just attempted to build the super pack installer on windows 7 Ultimate (32 bits), and did not encounter any issue, the testsuite passing everything but a few things unrelated to our problem here. Could you put your binary somewhere so that I can look at it ? You can find the binary here: http://www.patricktmarsh.com/numpy/20100306.py26-1.4.0-win32.superpack.exe I built this one today. Patrick -- Patrick Marsh Ph.D. Student / NSSL Liaison to the HWT School of Meteorology / University of Oklahoma Cooperative Institute for Mesoscale Meteorological Studies National Severe Storms Laboratory http://www.patricktmarsh.com ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] multiprocessing shared arrays and numpy
I did some optimization, and the results are very instructive, although not surprising: javascript:SetCmd(cmdSend); As I wrote before, I processed stereoscopic movie recordings, by making each a memory mapped file and processing it in several steps. By this way I produced extra GB of transient data. Running as one process took 45 seconds, and in dual parallel process ~40 seconds. After rewriting the application to process the recording frame by frame. The code became shorter and the new scores are: One process --- 16 seconds, and dual process --- 9 seconds. What I learned: * Design for multi-procssing from the start, not as afterthought * Shared memory works, but on the expense of code elegance (much like common blocks in fortran) * Memory mapped files can be used much as shared memory. The strange thing is that I got an ignored AttributeError on every frame access to the memory mapped file from the child process. Nadav -Original Message- From: numpy-discussion-boun...@scipy.org on behalf of Brian Granger Sent: Fri 05-Mar-10 21:29 To: Discussion of Numerical Python Subject: Re: [Numpy-discussion] multiprocessing shared arrays and numpy Francesc, Yeah, 10% of improvement by using multi-cores is an expected figure for memory bound problems. This is something people must know: if their computations are memory bound (and this is much more common that one may initially think), then they should not expect significant speed-ups on their parallel codes. +1 Thanks for emphasizing this. This is definitely a big issue with multicore. Cheers, Brian Thanks for sharing your experience anyway, Francesc A Thursday 04 March 2010 18:54:09 Nadav Horesh escrigué: I can not give a reliable answer yet, since I have some more improvement to make. The application is an analysis of a stereoscopic-movie raw-data recording (both channels are recorded in the same file). I treat the data as a huge memory mapped file. The idea was to process each channel (left and right) on a different core. Right now the application is IO bounded since I do classical numpy operation, so each channel (which is handled as one array) is scanned several time. The improvement now over a single process is 10%, but I hope to achieve 10% ore after trivial optimizations. I used this application as an excuse to dive into multi-processing. I hope that the code I posted here would help someone. Nadav. -Original Message- From: numpy-discussion-boun...@scipy.org on behalf of Francesc Alted Sent: Thu 04-Mar-10 15:12 To: Discussion of Numerical Python Subject: Re: [Numpy-discussion] multiprocessing shared arrays and numpy What kind of calculations are you doing with this module? Can you please send some examples and the speed-ups you are getting? Thanks, Francesc A Thursday 04 March 2010 14:06:34 Nadav Horesh escrigué: Extended module that I used for some useful work. Comments: 1. Sturla's module is better designed, but did not work with very large (although sub GB) arrays 2. Tested on 64 bit linux (amd64) + python-2.6.4 + numpy-1.4.0 Nadav. -Original Message- From: numpy-discussion-boun...@scipy.org on behalf of Nadav Horesh Sent: Thu 04-Mar-10 11:55 To: Discussion of Numerical Python Subject: RE: [Numpy-discussion] multiprocessing shared arrays and numpy Maybe the attached file can help. Adpted and tested on amd64 linux Nadav -Original Message- From: numpy-discussion-boun...@scipy.org on behalf of Nadav Horesh Sent: Thu 04-Mar-10 10:54 To: Discussion of Numerical Python Subject: Re: [Numpy-discussion] multiprocessing shared arrays and numpy There is a work by Sturla Molden: look for multiprocessing-tutorial.pdf and sharedmem-feb13-2009.zip. The tutorial includes what is dropped in the cookbook page. I am into the same issue and going to test it today. Nadav On Wed, 2010-03-03 at 15:31 +0100, Jesper Larsen wrote: Hi people, I was wondering about the status of using the standard library multiprocessing module with numpy. I found a cookbook example last updated one year ago which states that: This page was obsolete as multiprocessing's internals have changed. More information will come shortly; a link to this page will then be added back to the Cookbook. http://www.scipy.org/Cookbook/multiprocessing I also found the code that used to be on this page in the cookbook but it does not work any more. So my question is: Is it possible to use numpy arrays as shared arrays in an application using multiprocessing and how do you do it? Best regards, Jesper ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
Hi, On Sat, Mar 6, 2010 at 1:20 AM, Friedrich Romstedt friedrichromst...@gmail.com wrote: Hmm. Let's see ... Can you tell me how I can test the time calls in a script take? I have no idea. I've been doing: t1 = time.time() #code t2 = time.time() print [code]:,t2-t1 #0.0840001106262 sum_1 = sum.reshape((len(sum), 1)).repeat(len(sum), axis = 1) #0.032923706 sum_2 = sum.reshape((1, len(sum))).repeat(len(sum), axis = 0) #0.026504089 comparison_sum = (sum_1 == sum_2) We can leave out the repeat() calls and leave only the reshape() calls there. Numpy will substitute dimi == 1 dimensions with stride == 0, i.e., it will effectively repeat those dimension, just as we did it explicitly. Wow! Drops to immeasurably quick :D #0.0909998416901 diff_1 = diff.reshape((len(diff), 1)).repeat(len(diff), axis = 1) #0.0340001583099 diff_2 = diff.reshape((1, len(diff))).repeat(len(diff), axis = 0) #0.026504089 comparison_diff = (diff_1 == diff_2) Same here. Delete the repeat() calls, but not the reshape() calls. Once again, drops to immeasurably quick :D #0.023019073 same_edges = comparison_sum * comparison_diff Hmm, maybe use numpy.logical_and(comparison_sum, comparison_diff)? I don't know, but I guess it is in some way optimised for such things. It's marginally faster. #0.12848502 doublet_count = same_edges.sum(axis = 0) Maybe try axis = 1 instead. I wonder why this is so slow. Or maybe it's because he does the conversion to ints on-the-fly, so maybe try same_edges.astype(numpy.int8).sum(axis = 0). Actually, it's marginally slower :S Hope this gives some improvement. I attach the modified version. Ah, one thing to mention, have you not accidentally timed also the printout functions? They should be pretty slow. Nope--I've been timing as above. Friedrich On Sat, Mar 6, 2010 at 1:42 AM, Friedrich Romstedt friedrichromst...@gmail.com wrote: 2010/3/5 Ian Mallett geometr...@gmail.com: #takes 0.04 seconds inner = np.inner(ns, v1s - some_point) Ok, I don't know why I was able to overlook this: dotprod = (ns * (v1s - some_point)).sum(axis = 1) Much faster :D So, the most costly lines: comparison_sum = (sum_1 == sum_2) #0.024 sec comparison_diff = (diff_1 == diff_2) #0.024 sec same_edges = np.logical_and(comparison_sum, comparison_diff) #0.029 sec doublet_count = same_edges.sum(axis = 0) #0.147 Thanks again, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
2010/3/6 Ian Mallett geometr...@gmail.com: On Sat, Mar 6, 2010 at 1:20 AM, Friedrich Romstedt friedrichromst...@gmail.com wrote: Hmm. Let's see ... Can you tell me how I can test the time calls in a script take? I have no idea. I've been doing: t1 = time.time() #code t2 = time.time() print [code]:,t2-t1 Ok. I just wonder how you can measure those times. On my five-years old machine, they are 10ms resolution. Wow! Drops to immeasurably quick :D Yeaah! #0.12848502 doublet_count = same_edges.sum(axis = 0) Maybe try axis = 1 instead. I wonder why this is so slow. Or maybe it's because he does the conversion to ints on-the-fly, so maybe try same_edges.astype(numpy.int8).sum(axis = 0). Actually, it's marginally slower :S Hmm. I tried the axis = 1 thing, and it also gave no improvement (maybe you can try too, I'm guessing I'm actually measuring the time Python spends in exeuting my loop to get significant times ...) 2010/3/5 Ian Mallett geometr...@gmail.com: #takes 0.04 seconds inner = np.inner(ns, v1s - some_point) Ok, I don't know why I was able to overlook this: dotprod = (ns * (v1s - some_point)).sum(axis = 1) Much faster :D :-) I'm glad to be able to help! So, the most costly lines: comparison_sum = (sum_1 == sum_2) #0.024 sec comparison_diff = (diff_1 == diff_2) #0.024 sec same_edges = np.logical_and(comparison_sum, comparison_diff) #0.029 sec doublet_count = same_edges.sum(axis = 0) #0.147 At the moment, I can do nothing about that. Seems that we have reached the limit. Anyhow, is it now faster than your Python list implementation, and if yes, how much? How large was your gain by using numpy means at all? I'm just curious. Friedrich ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
On Sat, Mar 6, 2010 at 12:03 PM, Friedrich Romstedt friedrichromst...@gmail.com wrote: At the moment, I can do nothing about that. Seems that we have reached the limit. Anyhow, is it now faster than your Python list implementation, and if yes, how much? How large was your gain by using numpy means at all? I'm just curious. Unfortunately, the pure Python implementation is actually an order of magnitude faster. The fastest solution right now is to use numpy for the transformations, then convert it back into a list (.tolist()) and use Python for the rest. Here's the actual Python code. def glLibInternal_edges(object,lightpos): edge_set = set([]) edges = {} for sublist in xrange(object.number_of_lists): #There's only one sublist here face_data = object.light_volume_face_data[sublist] for indices in face_data: #v1,v2,v3,n normal = object.transformed_normals[sublist][indices[3]] v1,v2,v3 = [ object.transformed_vertices[sublist][indices[i]] for i in xrange(3) ] if abs_angle_between_rad(normal,vec_subt(v1,lightpos))pi_over_two: for p1,p2 in [[indices[0],indices[1]], [indices[1],indices[2]], [indices[2],indices[0]]]: edge = [p1,p2] index = 0 edge2 = list(edge) edge2.sort() edge2 = tuple(edge2) if edge2 in edges: edges[edge2][1] += 1 else: edges[edge2] = [edge,1] edges2 = [] for edge_data in edges.values(): if edge_data[1] == 1: p1 = object.transformed_vertices[sublist][edge_data[0][0]] p2 = object.transformed_vertices[sublist][edge_data[0][1]] edges2.append([p1,p2]) return edges2 Friedrich Thanks, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
I'm a bit unhappy with your code, because it's so hard to read tbh. You don't like objects? 2010/3/6 Ian Mallett geometr...@gmail.com: Unfortunately, the pure Python implementation is actually an order of magnitude faster. The fastest solution right now is to use numpy for the transformations, then convert it back into a list (.tolist()) and use Python for the rest. :-( But well, if it's faster, than do it that way, right? I can only guess, that for large datasets, the comparison to compare each and every vector in parallel makes it slow down. So an iterative approach might be fine. In fact, no one would implement such things in C++ using not lists with pushback(), I guess. A bit of commentary about your code: Here's the actual Python code. def glLibInternal_edges(object,lightpos): edge_set = set([]) Where do you use edge_set? Btw, set() would do it. edges = {} for sublist in xrange(object.number_of_lists): #There's only one sublist here face_data = object.light_volume_face_data[sublist] for indices in face_data: #v1,v2,v3,n Here objects would fit in nicely. for indices in face_data: normal = object.transformed_normals[sublist][indices.nindex] (v1, v2, v3) = [object.transformed_vertices[sublist][vidx] for vidx in indices.vindices] normal = object.transformed_normals[sublist][indices[3]] v1,v2,v3 = [ object.transformed_vertices[sublist][indices[i]] for i in xrange(3) ] if abs_angle_between_rad(normal,vec_subt(v1,lightpos))pi_over_two: v1, lightpos can be stored as numpy ndarrays, no? If you have a list of vectors in an ndarray arr, you can convert to a list of ndarrays by using list(arr). This will only iterate over the first dimension, and not recurse into the vectors in arr. Provided this, you can simply write: if (normal * (v1 - lightpos) / numpy.sqrt(((v1 - lightpos) ** 2).sum())) 0: (...) for p1,p2 in [[indices[0],indices[1]], [indices[1],indices[2]], [indices[2],indices[0]]]: edge = [p1,p2] Why not writing: for edge in numpy.asarray([[indices[0], indices[1]], (...)]): (...) index = 0 Where do you use index? It's a lonely remnant? edge2 = list(edge) Why do you convert the unmodified edge list into a list again? Btw, I found your numbering quite a bit annoying. No one can tell from an appended 2 what purpose that xxx2 has. Furthermore, I think a slight speedup could be reached by: unique = (egde.sum(), abs((egde * [1, -1]).sum())) edge2.sort() edge2 = tuple(edge2) if edge2 in edges: edges[edge2][1] += 1 else: edges[edge2] = [edge,1] edges2 = [] for edge_data in edges.values(): if edge_data[1] == 1: p1 = object.transformed_vertices[sublist][edge_data[0][0]] p2 = object.transformed_vertices[sublist][edge_data[0][1]] edges2.append([p1,p2]) return edges2 My 2 cents: class Edge: def __init__(self, indices, count): self.indices = indices def __hash__(self): return hash(self.unique) edges = {} (...) edges.setdefault() edges[egde2] += 1 for (indices, count) in edges.items(): if count == 1: edges_clean.append(object.transformed_vertices[sublist][indices]) provided that transformed_vertices[sublist] is an ndarray. You can iterate over ndarrays using standard Python iteration, it will provide an iterator object. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
Sent prematurely, typed tab and then space :-(. Real message work in progress. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
I'm a bit unhappy with your code, because it's so hard to read tbh. You don't like objects? I would phrase this differently now: I think you can improve your code by using objects when they are appropriate. :-) 2010/3/6 Ian Mallett geometr...@gmail.com: Unfortunately, the pure Python implementation is actually an order of magnitude faster. The fastest solution right now is to use numpy for the transformations, then convert it back into a list (.tolist()) and use Python for the rest. :-( But well, if it's faster, than do it that way, right? I can only guess, that for large datasets, the comparison to compare each and every vector in parallel makes it slow down. So an iterative approach might be fine. In fact, no one would implement such things in C++ using not lists with pushback(), I guess. A bit of commentary about your code: Here's the actual Python code. def glLibInternal_edges(object,lightpos): edge_set = set([]) Where do you use edge_set? Btw, set() would do it. edges = {} for sublist in xrange(object.number_of_lists): #There's only one sublist here face_data = object.light_volume_face_data[sublist] for indices in face_data: #v1,v2,v3,n Here objects would fit in nicely. for indices in face_data: normal = object.transformed_normals[sublist][indices.nindex] (v1, v2, v3) = [object.transformed_vertices[sublist][vidx] for vidx in indices.vindices] normal = object.transformed_normals[sublist][indices[3]] v1,v2,v3 = [ object.transformed_vertices[sublist][indices[i]] for i in xrange(3) ] if abs_angle_between_rad(normal,vec_subt(v1,lightpos))pi_over_two: v1, lightpos can be stored as numpy ndarrays, no? If you have a list of vectors in an ndarray arr, you can convert to a list of ndarrays by using list(arr). This will only iterate over the first dimension, and not recurse into the vectors in arr. Provided this, you can simply write: if (normal * (v1 - lightpos) / numpy.sqrt(((v1 - lightpos) ** 2).sum())) 0: (...) for p1,p2 in [[indices[0],indices[1]], [indices[1],indices[2]], [indices[2],indices[0]]]: edge = [p1,p2] Why not writing: for edge in numpy.asarray([[indices[0], indices[1]], (...)]): (...) index = 0 Where do you use index? It's a lonely remnant? edge2 = list(edge) Why do you convert the unmodified edge list into a list again? Btw, I found your numbering quite a bit annoying. No one can tell from an appended 2 what purpose that xxx2 has. Furthermore, I think a slight speedup could be reached by: unique = (egde.sum(), abs((egde * [1, -1]).sum())) edge2.sort() edge2 = tuple(edge2) if edge2 in edges: edges[edge2][1] += 1 else: edges[edge2] = [edge,1] edges2 = [] for edge_data in edges.values(): if edge_data[1] == 1: p1 = object.transformed_vertices[sublist][edge_data[0][0]] p2 = object.transformed_vertices[sublist][edge_data[0][1]] edges2.append([p1,p2]) return edges2 My 2 cents: class CountingEdge: def __init__(self, indices): self.indices = indices self.count = 0 edges = {} (...) edges.setdefault(unique, CountingEdge(edge)) edges[unique].count += 1 for counting_edge in edges.values(): if counting_edge.count == 1: edges_clean.append(object.transformed_vertices[sublist][counting_edge.indices]) provided that transformed_vertices[sublist] is an ndarray. You can iterate over ndarrays using standard Python iteration, it will provide an iterator object. I'm always good for fresh ideas, although they nearly never get accepted, but here is another cent: I want to save you calculating the mutual annihilation for several light sources again and again. BEFORE SHADING 1. For the mesh, calculate all edges. Put their unique tuples in a list. Store together with the unique tuple from what triangle they originate and what the original indices were. When you use the sorted indices as unique value, the original indices will be virtually identical with the unique tuple. class EdgeOfTriangle: def __init__(self, unique, vindices, triangleindex): self.unique = unique self.vindices = vindices self.triangleindex = triangleindex def __hash__(self): # For the in statement return hash(self.unique) def __eq__(self, other): return other == self.triangle edges = [] for triangle_descriptor in triangle_descriptors: (...) egde1 = EdgeOfTriangle(unique, vindices1, triangle_descriptor.nindex) edge2 = EdgeOfTriangle(unique, vindices2, traingle_descriptor.nindex) (and so on) for new_egde in (edge1, egde2, egde3): if new_edge not in edges: edges.append(new_edge) edges =
[Numpy-discussion] Why is the shape of a singleton array the empty tuple?
x = numpy.array(3) x array(3) x.shape () My question is: why? DG ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Why is the shape of a singleton array the empty tuple?
On Sat, Mar 6, 2010 at 9:37 PM, Ian Mallett geometr...@gmail.com wrote: x = numpy.array(3) x array(3) x.shape () y = numpy.array([3]) y array([3]) y.shape (1,) Ian Thanks, Ian. I already figured out how to make it not so, but I still want to understand the design reasoning behind it being so in the first place (thus the use of the question why (is it so), not how (to make it different)). DG ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Why is the shape of a singleton array the empty tuple?
On Sat, Mar 6, 2010 at 9:46 PM, David Goldsmith d.l.goldsm...@gmail.comwrote: Thanks, Ian. I already figured out how to make it not so, but I still want to understand the design reasoning behind it being so in the first place (thus the use of the question why (is it so), not how (to make it different)). Well, I can't help you with that. I would also ask why this design even exists? Equating an array with a single number doesn't make sense to me. Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Why is the shape of a singleton array the empty tuple?
On Sat, Mar 6, 2010 at 10:26 PM, Ian Mallett geometr...@gmail.com wrote: On Sat, Mar 6, 2010 at 9:46 PM, David Goldsmith d.l.goldsm...@gmail.comwrote: Thanks, Ian. I already figured out how to make it not so, but I still want to understand the design reasoning behind it being so in the first place (thus the use of the question why (is it so), not how (to make it different)). Well, I can't help you with that. I would also ask why this design even exists? Equating an array with a single number doesn't make sense to me. Ian Here's an (unintended) use case: I wanted to convert anything in an array that's close to zero to be zero (and leave the rest as is), but I want it to be robust so that if it receives a scalar, it can work w/ that, too. Here's my (working) code (I'm sure that once Robert sees it, he'll be able to replace it w/ a one-liner): def convert_close(arg): arg = N.array(arg) if not arg.shape: arg = N.array((arg,)) if arg.size: t = N.array([0 if N.allclose(temp, 0) else temp for temp in arg]) if len(t.shape) - 1: return N.squeeze(t) else: return t else: return N.array() At first I wasn't casting arg to be an array upon entry, but I found that if arg is a scalar, my list comprehension failed, so I had choose _some_ sort of sequence to cast scalars to; since arg will typically be an array and that's what I wanted to return as well, it seemed most appropriate to cast incoming scalars to arrays. So I added the arg = N.array(arg) at the beginning (N.array(array) = array, and N.array(non-array seq) does the right thing as well), but the list comprehension still wouldn't work if arg was a scalar upon entry; after many print statements and much interactive experimenting, I finally figured out that this is because the shape of N.array(scalar) is () (and I thence immediately guessed, correctly of course, that N.array((scalar,)) has shape (1,)). So I added the if not arg.shape: to detect and correct for those zero size arrays, and now it works fine, but I'd still like to know _why_ N.array(scalar).shape == () but N.array((scalar,)).shape == (1,). No biggy, just curious. DG ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion