[Yt-svn] yt: Moving the joining phase for contour finding into Cython

[hg at spacepope.org]

hg Repository: yt
details:   yt/rev/1b81f9d8847c
changeset: 3746:1b81f9d8847c
user:      Matthew Turk <matthewturk at gmail.com>
date:
Wed Feb 16 04:38:44 2011 -0800
description:
Moving the joining phase for contour finding into Cython

diffstat:

 yt/analysis_modules/level_sets/contour_finder.py |  28 +++++++++++++++++------
 yt/utilities/_amr_utils/ContourFinding.pyx       |  17 ++++++++++++++
 2 files changed, 37 insertions(+), 8 deletions(-)

diffs (78 lines):

diff -r 73bbf5a5c797 -r 1b81f9d8847c yt/analysis_modules/level_sets/contour_finder.py
--- a/yt/analysis_modules/level_sets/contour_finder.py	Tue Feb 15 16:49:27 2011 -0700
+++ b/yt/analysis_modules/level_sets/contour_finder.py	Wed Feb 16 04:38:44 2011 -0800
@@ -256,7 +256,11 @@
                     s1.update(joins.pop(k2))
                     s1.update([k2])
                     updated += 1
-    return joins
+    tr = []
+    for k in joins.keys():
+        v = joins.pop(k)
+        tr.append((k, na.array(list(v), dtype="int64")))
+    return tr
 
 def identify_contours(data_source, field, min_val, max_val,
                           cached_fields=None):
@@ -300,15 +304,23 @@
     sort_new = na.array(list(set(tree)), dtype='int64')
     mylog.info("Coalescing %s joins", sort_new.shape[0])
     joins = coalesce_join_tree(sort_new)
+    #joins = [(i, na.array(list(j), dtype="int64")) for i, j in sorted(joins.items())]
     pbar = get_pbar("Joining ", len(joins))
     # This process could and should be done faster
-    for i, new in enumerate(sorted(joins.keys())):
-        pbar.update(i)
-        old_set = joins[new]
-        for old in old_set:
-            if old == new: continue
-            i1 = (data_source["tempContours"] == old)
-            data_source["tempContours"][i1] = new
+    print "Joining..."
+    t1 = time.time()
+    ff = data_source["tempContours"].astype("int64")
+    amr_utils.update_joins(joins, ff)
+    data_source["tempContours"] = ff.astype("float64")
+    #for i, new in enumerate(sorted(joins.keys())):
+    #    pbar.update(i)
+    #    old_set = joins[new]
+    #    for old in old_set:
+    #        if old == new: continue
+    #        i1 = (data_source["tempContours"] == old)
+    #        data_source["tempContours"][i1] = new
+    t2 = time.time()
+    print "Finished joining in %0.2e seconds" % (t2-t1)
     pbar.finish()
     data_source._flush_data_to_grids("tempContours", -1, dtype='int64')
     del data_source.data["tempContours"] # Force a reload from the grids
diff -r 73bbf5a5c797 -r 1b81f9d8847c yt/utilities/_amr_utils/ContourFinding.pyx
--- a/yt/utilities/_amr_utils/ContourFinding.pyx	Tue Feb 15 16:49:27 2011 -0700
+++ b/yt/utilities/_amr_utils/ContourFinding.pyx	Wed Feb 16 04:38:44 2011 -0800
@@ -39,6 +39,7 @@
     return i1
 
 @cython.boundscheck(False)
+ at cython.wraparound(False)
 def construct_boundary_relationships(
         np.ndarray[dtype=np.int64_t, ndim=3] contour_ids):
     # We only look at the boundary and one cell in
@@ -174,3 +175,19 @@
         for j in proto_contour:
             contours[j] = proto_contour
     return contours
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+def update_joins(joins, np.ndarray[np.int64_t, ndim=1] contour_ids):
+    cdef np.int64_t new, old, i, oi
+    cdef int n, on
+    cdef np.ndarray[np.int64_t, ndim=1] old_set
+    #print contour_ids.shape[0]
+    n = contour_ids.shape[0]
+    for new, old_set in joins:
+        #print new
+        on = old_set.shape[0]
+        for i in range(n):
+            for oi in range(on):
+                old = old_set[oi]
+                if contour_ids[i] == old: contour_ids[i] = new