[yt-svn] commit/yt-3.0: 6 new changesets

Bitbucket commits-reply at bitbucket.org
Fri Feb 8 14:00:24 PST 2013


6 new commits in yt-3.0:

https://bitbucket.org/yt_analysis/yt-3.0/commits/bc401eff6878/
changeset:   bc401eff6878
branch:      yt-3.0
user:        MatthewTurk
date:        2013-02-04 23:55:43
summary:     Initial import of a smoothing_kernel routine.
affected #:  3 files

diff -r 4a75438eafcb349e21f9fe0df651b9f3c1d86179 -r bc401eff6878310015520dec8fc03329ca0f2d05 yt/frontends/setup.py
--- a/yt/frontends/setup.py
+++ b/yt/frontends/setup.py
@@ -7,17 +7,18 @@
     config = Configuration('frontends', parent_package, top_path)
     config.make_config_py()  # installs __config__.py
     #config.make_svn_version_py()
+    config.add_subpackage("art")
     config.add_subpackage("athena")
-    config.add_subpackage("gdf")
+    config.add_subpackage("castro")
     config.add_subpackage("chombo")
     config.add_subpackage("enzo")
     config.add_subpackage("flash")
+    config.add_subpackage("gdf")
+    config.add_subpackage("maestro")
     config.add_subpackage("nyx")
     config.add_subpackage("orion")
     config.add_subpackage("ramses")
+    config.add_subpackage("sph")
+    config.add_subpackage("stream")
     config.add_subpackage("tiger")
-    config.add_subpackage("art")
-    config.add_subpackage("maestro")
-    config.add_subpackage("castro")
-    config.add_subpackage("stream")
     return config

diff -r 4a75438eafcb349e21f9fe0df651b9f3c1d86179 -r bc401eff6878310015520dec8fc03329ca0f2d05 yt/frontends/sph/setup.py
--- a/yt/frontends/sph/setup.py
+++ b/yt/frontends/sph/setup.py
@@ -3,11 +3,18 @@
 import os
 import sys
 import os.path
-
+import glob
 
 def configuration(parent_package='', top_path=None):
     from numpy.distutils.misc_util import Configuration
-    config = Configuration('gadget', parent_package, top_path)
+    config = Configuration('sph', parent_package, top_path)
+    config.add_extension("smoothing_kernel",
+        ["yt/frontends/sph/smoothing_kernel.pyx"],
+        include_dirs=["yt/frontends/sph/",
+                      "yt/geometry/",
+                      "yt/utilities/lib/"],
+        depends=glob.glob("yt/geometry/*.pxd"),
+        )
     config.make_config_py()  # installs __config__.py
     #config.make_svn_version_py()
     return config

diff -r 4a75438eafcb349e21f9fe0df651b9f3c1d86179 -r bc401eff6878310015520dec8fc03329ca0f2d05 yt/frontends/sph/smoothing_kernel.pyx
--- /dev/null
+++ b/yt/frontends/sph/smoothing_kernel.pyx
@@ -0,0 +1,99 @@
+# This set of routines was originally written by Robert Thompson.
+cimport numpy as np
+import numpy as np
+from libc.stdlib cimport free, malloc
+from libc.math cimport sqrt
+cimport cython
+
+#@cython.boundscheck(False)
+#@cython.wraparound(False)
+ at cython.cdivision(True)
+def grid(input_fields, output_grids,
+         np.ndarray[np.float64_t, ndim=1] left_edge,
+         np.ndarray[np.float64_t, ndim=1] right_edge,
+         np.ndarray[np.float64_t, ndim=2] ppos,
+         np.ndarray[np.float64_t, ndim=1] hsml):
+
+    cdef np.int64_t ngas
+    cdef np.float64_t dds[3], pos[3], idist[3], kern
+    cdef int p, i, j, k, d, ind[3], ib0[3], ib1[3], dims[3]
+    cdef int nf, half_len, skip
+    cdef np.float64_t dist
+    cdef np.ndarray[np.float64_t, ndim=1] gas_arr
+    cdef np.ndarray[np.float64_t, ndim=3] grid_arr
+
+    ngas = input_fields[0].size
+    for i in range(3):
+        dims[i] = output_grids[0].shape[i]
+        dds[i] = (right_edge[i] - left_edge[i])/dims[i]
+
+    cdef np.float64_t *kernel_sum = \
+        <np.float64_t *>malloc(ngas * sizeof(np.float64_t))
+
+    nf = len(input_fields)
+    assert(nf == len(output_grids))
+
+    cdef np.float64_t **pdata = <np.float64_t**> \
+            malloc(sizeof(np.float64_t *) * nf)
+    cdef np.float64_t **gdata = <np.float64_t**> \
+            malloc(sizeof(np.float64_t *) * nf)
+
+    for i in range(nf):
+        gas_arr = input_fields[i]
+        grid_arr = output_grids[i]
+        pdata[i] = <np.float64_t*> gas_arr.data
+        gdata[i] = <np.float64_t*> grid_arr.data
+
+    for p in range(ngas):
+        kernel_sum[p] = 0.0
+        skip = 0
+        for i in range(3):
+            pos[i] = ppos[p, i]
+            ind[i] = <int>((pos[i] - left_edge[i]) / dds[i])
+            half_len = <int>(hsml[p]/dds[i])
+            ib0[i] = ind[i] - half_len
+            ib1[i] = ind[i] + half_len
+            if ib0[i] >= dims[i] or ib1[i] < 0:
+                skip = 1
+        if skip == 1: continue
+
+        for i from ib0[0] <= i <= ib1[0]:
+            idist[0] = (ind[0] - i) * (ind[0] - i) * dds[0]
+            for j from ib0[1] <= j <= ib1[1]:
+                idist[1] = (ind[1] - j) * (ind[1] - j) * dds[1] 
+                idist[1] += idist[0]
+                for k from ib0[2] <= k <= ib1[2]:
+                    idist[2] = (ind[2] - k) * (ind[2] - k) * dds[2]
+                    idist[2] += idist[1]
+                    dist = sqrt(idist[2]) / hsml[p]
+                    kernel_sum[p] += sph_kernel(dist)
+
+        for i from ib0[0] <= i <= ib1[0]:
+            idist[0] = (ind[0] - i) * (ind[0] - i) * dds[0]
+            for j from ib0[1] <= j <= ib1[1]:
+                idist[1] = (ind[1] - j) * (ind[1] - j) * dds[1] 
+                idist[1] += idist[0]
+                for k from ib0[2] <= k <= ib1[2]:
+                    idist[2] = (ind[2] - k) * (ind[2] - k) * dds[2]
+                    idist[2] += idist[1]
+                    dist = sqrt(idist[2]) / hsml[p]
+                    kern = sph_kernel(dist)
+                    gi = ((i * dims[1] + j) * dims[2]) + k
+                    for d in range(nf):
+                        gdata[d][gi] += pdata[d][p] * kern / kernel_sum[p]
+
+    free(kernel_sum)
+    free(gdata)
+    free(pdata)
+
+##############################################
+#Standard SPH kernel for use with the Grid method
+cdef float sph_kernel(float x) nogil:
+    cdef float kernel
+    if x <= 0.5:
+        kernel = 1.-6.*x*x*(1.-x)
+    elif x>0.5 and x<=1.0:
+        kernel = 2.*(1.-x)*(1.-x)*(1.-x)
+    else:
+        kernel = 0.
+    return kernel


https://bitbucket.org/yt_analysis/yt-3.0/commits/9a8cdca24326/
changeset:   9a8cdca24326
branch:      yt-3.0
user:        MatthewTurk
date:        2013-02-05 02:41:37
summary:     This fixes a missing factor of dds[i] in calculating smoothing kernel radius.
affected #:  1 file

diff -r bc401eff6878310015520dec8fc03329ca0f2d05 -r 9a8cdca24326a31954ca7828cc2ff337ece64ef6 yt/frontends/sph/smoothing_kernel.pyx
--- a/yt/frontends/sph/smoothing_kernel.pyx
+++ b/yt/frontends/sph/smoothing_kernel.pyx
@@ -3,19 +3,21 @@
 import numpy as np
 from libc.stdlib cimport free, malloc
 from libc.math cimport sqrt
+from fp_utils cimport iclip
 cimport cython
 
 #@cython.boundscheck(False)
 #@cython.wraparound(False)
 @cython.cdivision(True)
-def grid(input_fields, output_grids,
-         np.ndarray[np.float64_t, ndim=1] left_edge,
-         np.ndarray[np.float64_t, ndim=1] right_edge,
-         np.ndarray[np.float64_t, ndim=2] ppos,
-         np.ndarray[np.float64_t, ndim=1] hsml):
+def smooth_particles(
+        input_fields, output_grids,
+        np.ndarray[np.float64_t, ndim=1] left_edge,
+        np.ndarray[np.float64_t, ndim=1] right_edge,
+        np.ndarray[np.float64_t, ndim=2] ppos,
+        np.ndarray[np.float64_t, ndim=1] hsml):
 
     cdef np.int64_t ngas
-    cdef np.float64_t dds[3], pos[3], idist[3], kern
+    cdef np.float64_t dds[3], sdds[3], pos[3], idist[3], kern
     cdef int p, i, j, k, d, ind[3], ib0[3], ib1[3], dims[3]
     cdef int nf, half_len, skip
     cdef np.float64_t dist
@@ -26,6 +28,7 @@
     for i in range(3):
         dims[i] = output_grids[0].shape[i]
         dds[i] = (right_edge[i] - left_edge[i])/dims[i]
+        sdds[i] = dds[i] * dds[i]
 
     cdef np.float64_t *kernel_sum = \
         <np.float64_t *>malloc(ngas * sizeof(np.float64_t))
@@ -50,38 +53,35 @@
         for i in range(3):
             pos[i] = ppos[p, i]
             ind[i] = <int>((pos[i] - left_edge[i]) / dds[i])
-            half_len = <int>(hsml[p]/dds[i])
+            half_len = <int>(hsml[p]/dds[i]) + 1
             ib0[i] = ind[i] - half_len
             ib1[i] = ind[i] + half_len
             if ib0[i] >= dims[i] or ib1[i] < 0:
                 skip = 1
+            ib0[i] = iclip(ib0[i], 0, dims[i] - 1)
+            ib1[i] = iclip(ib1[i], 0, dims[i] - 1)
         if skip == 1: continue
-
         for i from ib0[0] <= i <= ib1[0]:
-            idist[0] = (ind[0] - i) * (ind[0] - i) * dds[0]
+            idist[0] = (ind[0] - i) * (ind[0] - i) * sdds[0]
             for j from ib0[1] <= j <= ib1[1]:
-                idist[1] = (ind[1] - j) * (ind[1] - j) * dds[1] 
-                idist[1] += idist[0]
+                idist[1] = (ind[1] - j) * (ind[1] - j) * sdds[1] 
                 for k from ib0[2] <= k <= ib1[2]:
-                    idist[2] = (ind[2] - k) * (ind[2] - k) * dds[2]
-                    idist[2] += idist[1]
-                    dist = sqrt(idist[2]) / hsml[p]
+                    idist[2] = (ind[2] - k) * (ind[2] - k) * sdds[2]
+                    dist = idist[0] + idist[1] + idist[2]
+                    dist = sqrt(dist) / hsml[p]
                     kernel_sum[p] += sph_kernel(dist)
-
         for i from ib0[0] <= i <= ib1[0]:
-            idist[0] = (ind[0] - i) * (ind[0] - i) * dds[0]
+            idist[0] = (ind[0] - i) * (ind[0] - i) * sdds[0]
             for j from ib0[1] <= j <= ib1[1]:
-                idist[1] = (ind[1] - j) * (ind[1] - j) * dds[1] 
-                idist[1] += idist[0]
+                idist[1] = (ind[1] - j) * (ind[1] - j) * sdds[1] 
                 for k from ib0[2] <= k <= ib1[2]:
-                    idist[2] = (ind[2] - k) * (ind[2] - k) * dds[2]
-                    idist[2] += idist[1]
-                    dist = sqrt(idist[2]) / hsml[p]
+                    idist[2] = (ind[2] - k) * (ind[2] - k) * sdds[2]
+                    dist = idist[0] + idist[1] + idist[2]
+                    dist = sqrt(dist) / hsml[p]
                     kern = sph_kernel(dist)
                     gi = ((i * dims[1] + j) * dims[2]) + k
                     for d in range(nf):
                         gdata[d][gi] += pdata[d][p] * kern / kernel_sum[p]
-
     free(kernel_sum)
     free(gdata)
     free(pdata)


https://bitbucket.org/yt_analysis/yt-3.0/commits/83c602e62523/
changeset:   83c602e62523
branch:      yt-3.0
user:        MatthewTurk
date:        2013-02-05 06:52:36
summary:     Speed up the kernel computation considerably by storing distance.

This also includes a different calculation of the half_len, but commented out.
I need to think carefully what the correct half_len is.
affected #:  1 file

diff -r 9a8cdca24326a31954ca7828cc2ff337ece64ef6 -r 83c602e62523446d613e7f55d3087400b652f8e3 yt/frontends/sph/smoothing_kernel.pyx
--- a/yt/frontends/sph/smoothing_kernel.pyx
+++ b/yt/frontends/sph/smoothing_kernel.pyx
@@ -6,8 +6,8 @@
 from fp_utils cimport iclip
 cimport cython
 
-#@cython.boundscheck(False)
-#@cython.wraparound(False)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
 @cython.cdivision(True)
 def smooth_particles(
         input_fields, output_grids,
@@ -32,6 +32,9 @@
 
     cdef np.float64_t *kernel_sum = \
         <np.float64_t *>malloc(ngas * sizeof(np.float64_t))
+    cdef np.float64_t *pdist = \
+        <np.float64_t *>malloc(dims[0]*dims[1]*dims[2]*
+                               sizeof(np.float64_t))
 
     nf = len(input_fields)
     assert(nf == len(output_grids))
@@ -56,6 +59,10 @@
             half_len = <int>(hsml[p]/dds[i]) + 1
             ib0[i] = ind[i] - half_len
             ib1[i] = ind[i] + half_len
+            #pos[i] = ppos[p, i] - left_edge[i]
+            #ind[i] = <int>(pos[i] / dds[i])
+            #ib0[i] = <int>((pos[i] - hsml[i]) / dds[i]) - 1
+            #ib1[i] = <int>((pos[i] + hsml[i]) / dds[i]) + 1
             if ib0[i] >= dims[i] or ib1[i] < 0:
                 skip = 1
             ib0[i] = iclip(ib0[i], 0, dims[i] - 1)
@@ -69,20 +76,18 @@
                     idist[2] = (ind[2] - k) * (ind[2] - k) * sdds[2]
                     dist = idist[0] + idist[1] + idist[2]
                     dist = sqrt(dist) / hsml[p]
-                    kernel_sum[p] += sph_kernel(dist)
+                    gi = ((i * dims[1] + j) * dims[2]) + k
+                    pdist[gi] = sph_kernel(dist)
+                    kernel_sum[p] += pdist[gi]
         for i from ib0[0] <= i <= ib1[0]:
-            idist[0] = (ind[0] - i) * (ind[0] - i) * sdds[0]
             for j from ib0[1] <= j <= ib1[1]:
-                idist[1] = (ind[1] - j) * (ind[1] - j) * sdds[1] 
                 for k from ib0[2] <= k <= ib1[2]:
-                    idist[2] = (ind[2] - k) * (ind[2] - k) * sdds[2]
-                    dist = idist[0] + idist[1] + idist[2]
-                    dist = sqrt(dist) / hsml[p]
-                    kern = sph_kernel(dist)
                     gi = ((i * dims[1] + j) * dims[2]) + k
+                    dist = pdist[gi] / kernel_sum[p]
                     for d in range(nf):
-                        gdata[d][gi] += pdata[d][p] * kern / kernel_sum[p]
+                        gdata[d][gi] += pdata[d][p] * dist
     free(kernel_sum)
+    free(pdist)
     free(gdata)
     free(pdata)
 


https://bitbucket.org/yt_analysis/yt-3.0/commits/32db5f147394/
changeset:   32db5f147394
branch:      yt-3.0
user:        MatthewTurk
date:        2013-02-05 15:46:20
summary:     Considerable speedup in OWLS IO.  Slight Smoothing Kernel speedup.

By reading the data from disk *before* getting the mask, we can perform the
subselection in memory, rather than creating the HDF5 H5S objects.
Additionally, we can do all three vector field components at once inside Cython
for even further speedup.
affected #:  3 files

diff -r 83c602e62523446d613e7f55d3087400b652f8e3 -r 32db5f147394385d6b001462c18734b82ae9a55f yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -27,6 +27,7 @@
 import numpy as np
 from yt.funcs import *
 from yt.utilities.exceptions import *
+from yt.geometry.selection_routines import particle_mask_fill
 
 from yt.utilities.io_handler import \
     BaseIOHandler
@@ -80,12 +81,19 @@
                     del coords
                     if mask is None: continue
                     for field in field_list:
-                        data = g[field][mask,...]
                         my_ind = ind[ptype, field]
-                        mylog.debug("Filling from %s to %s with %s",
-                            my_ind, my_ind+data.shape[0], field)
-                        rv[ptype, field][my_ind:my_ind + data.shape[0],...] = data
-                        ind[ptype, field] += data.shape[0]
+                        if field in _vector_fields:
+                            mylog.debug("Filling (vectors) from %s with %s",
+                                my_ind, field)
+                            data = g[field][:].astype("float64")
+                            ind[ptype, field] = particle_mask_fill(
+                                mask, data, rv[ptype, field], my_ind)
+                        else:
+                            data = g[field][:][mask]
+                            mylog.debug("Filling (scalars) from %s to %s with %s",
+                                my_ind, my_ind+data.shape[0], field)
+                            rv[ptype, field][my_ind:my_ind + data.shape[0]] = data
+                            ind[ptype, field] += data.shape[0]
                 f.close()
         return rv
 

diff -r 83c602e62523446d613e7f55d3087400b652f8e3 -r 32db5f147394385d6b001462c18734b82ae9a55f yt/frontends/sph/smoothing_kernel.pyx
--- a/yt/frontends/sph/smoothing_kernel.pyx
+++ b/yt/frontends/sph/smoothing_kernel.pyx
@@ -19,7 +19,7 @@
     cdef np.int64_t ngas
     cdef np.float64_t dds[3], sdds[3], pos[3], idist[3], kern
     cdef int p, i, j, k, d, ind[3], ib0[3], ib1[3], dims[3]
-    cdef int nf, half_len, skip
+    cdef int nf, half_len, skip, gi
     cdef np.float64_t dist
     cdef np.ndarray[np.float64_t, ndim=1] gas_arr
     cdef np.ndarray[np.float64_t, ndim=3] grid_arr
@@ -93,8 +93,8 @@
 
 ##############################################
 #Standard SPH kernel for use with the Grid method
-cdef float sph_kernel(float x) nogil:
-    cdef float kernel
+cdef np.float64_t sph_kernel(np.float64_t x) nogil:
+    cdef np.float64_t kernel
     if x <= 0.5:
         kernel = 1.-6.*x*x*(1.-x)
     elif x>0.5 and x<=1.0:

diff -r 83c602e62523446d613e7f55d3087400b652f8e3 -r 32db5f147394385d6b001462c18734b82ae9a55f yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -114,6 +114,21 @@
     else:
         raise RuntimeError
 
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def particle_mask_fill(np.ndarray[np.uint8_t, ndim=1, cast=True] mask,
+                       np.ndarray[np.float64_t, ndim=2] input,
+                       np.ndarray[np.float64_t, ndim=2] output,
+                       np.int64_t offset):
+    cdef np.int64_t i, j
+    for i in range(mask.shape[0]):
+        if mask[i] == 0: continue
+        for j in range(3):
+            output[offset, j] = input[i, j]
+        offset += 1
+    return offset
+
 # Inclined Box
 
 cdef class SelectorObject:


https://bitbucket.org/yt_analysis/yt-3.0/commits/51519d4c147e/
changeset:   51519d4c147e
branch:      yt-3.0
user:        MatthewTurk
date:        2013-02-05 15:48:03
summary:     Considerable speedup for OWLS data by removing implicit H5S creation.
affected #:  2 files
Diff not available.

https://bitbucket.org/yt_analysis/yt-3.0/commits/82ad2207eeed/
changeset:   82ad2207eeed
branch:      yt-3.0
user:        MatthewTurk
date:        2013-02-05 15:49:11
summary:     Merging (no-op)
affected #:  2 files
Diff not available.

Repository URL: https://bitbucket.org/yt_analysis/yt-3.0/

--

This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.



More information about the yt-svn mailing list