[Yt-svn] yt-commit r397 - in trunk: tests tests/regression_scripts yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri Apr 11 00:20:32 PDT 2008
Author: mturk
Date: Fri Apr 11 00:20:30 2008
New Revision: 397
URL: http://yt.spacepope.org/changeset/397
Log:
Moved the contour finding into proper C.
There is a bit of a regression, as noted in the regression test file I am
committing. It's just a single cell. I am going to look into it tomorrow.
Added:
trunk/tests/regression_scripts/gal.py
Modified:
trunk/tests/test_lagos.py
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/ContourFinder.py
trunk/yt/lagos/PointCombine.c
Added: trunk/tests/regression_scripts/gal.py
==============================================================================
--- (empty file)
+++ trunk/tests/regression_scripts/gal.py Fri Apr 11 00:20:30 2008
@@ -0,0 +1,117 @@
+""" Original with Weave:
+Doing contour 9 / 10 (4.56276e-14 9.91625e-14)
+yt.lagos INFO 2008-04-10 17:04:45,085 Getting field Density from 197
+yt.lagos INFO 2008-04-10 17:04:46,284 Contouring over 1061337 cells with 30115 candidates
+yt.lagos INFO 2008-04-10 17:05:17,880 Getting field tempContours from 197
+yt.lagos INFO 2008-04-10 17:05:18,305 Identified 5 contours between 4.56276e-14 and 2.15510e-13
+yt.lagos INFO 2008-04-10 17:05:18,312 Getting field Contours from 197
+yt.lagos INFO 2008-04-10 17:05:19,330 Getting field GridIndices from 197
+ Contour id 0.0 has: 4.56297e-14 9.88512e-14 (2351) (4 grids, 11081.0 11087.0)
+ Contour id 1.0 has: 4.58381e-14 5.41707e-14 (13) (1 grids, 11079.0 11079.0)
+ Contour id 2.0 has: 4.60523e-14 4.79728e-14 (19) (2 grids, 11079.0 11082.0)
+ Contour id 3.0 has: 4.56287e-14 2.15439e-13 (26063) (32 grids, 11040.0 11103.0)
+ Contour id 4.0 has: 4.56339e-14 7.18032e-14 (1669) (3 grids, 11039.0 11051.0)
+"""
+
+""" New with C extension:
+yt.lagos INFO 2008-04-11 00:18:03,943 Getting field Density from 197
+yt.lagos INFO 2008-04-11 00:18:05,109 Contouring over 1061337 cells with 30115 candidates
+yt.lagos INFO 2008-04-11 00:18:36,243 Getting field tempContours from 197
+yt.lagos INFO 2008-04-11 00:18:36,626 Identified 5 contours between 4.56276e-14 and 2.15510e-13
+yt.lagos INFO 2008-04-11 00:18:36,633 Getting field Contours from 197
+yt.lagos INFO 2008-04-11 00:18:37,601 Getting field GridIndices from 197
+ Contour id 0.0 has: 3.43716e-14 9.88512e-14 (2352) (4 grids, 11081.0 11087.0)
+ Contour id 1.0 has: 4.58381e-14 5.41707e-14 (13) (1 grids, 11079.0 11079.0)
+ Contour id 2.0 has: 4.60523e-14 4.79728e-14 (19) (2 grids, 11079.0 11082.0)
+ Contour id 3.0 has: 4.56287e-14 2.15439e-13 (26063) (32 grids, 11040.0 11103.0)
+ Contour id 4.0 has: 4.56339e-14 7.18032e-14 (1669) (3 grids, 11039.0 11051.0)
+"""
+
+
+import sys
+sys.path.insert(0,"/Users/matthewturk/Development/yt/trunk")
+from yt.config import ytcfg
+
+ytcfg["yt","LogLevel"] = '20'
+ytcfg["yt","logFile"] = "False"
+ytcfg["yt","suppressStreamLogging"] = "True"
+
+con_field = "Temperature"
+con_field = "Density"
+
+import yt.lagos as lagos
+import yt.raven as raven
+import numpy as na
+import scipy.weave as weave
+import scipy.weave.converters as converters
+
+print lagos, raven
+
+#fn = "/Users/matthewturk/Research/data/yt_test_data/galaxy1360.dir/galaxy1360"
+fn = "/Users/matthewturk/Research/data/DataDir0036/DataDump0036"
+a = lagos.EnzoStaticOutput(fn)
+
+v,c = a.h.find_max("Density")
+
+sp = a.h.sphere(c,20000.0/a['au'])
+sp2 = a.h.sphere(c,200.0/a['au'])
+
+#lagos.identify_contours(sp, "Density", 1e-27, 1e-24)
+#lagos.identify_contours(sp, "Density", 1e-33, 1e-22)
+nc = 10
+cons = na.logspace(na.log10(sp2[con_field].min()*0.9),
+ na.log10(sp2[con_field].max()),nc+1)
+do_plot = 1
+k = []
+k_bad = []
+for i in [8]:#range(nc):
+ [grid.clear_data() for grid in sp._grids]
+ mi, ma = cons[i], cons[i+1]
+ print "Doing contour %s / %s (%0.5e %0.5e)" % (i+1,nc,mi,ma)
+ mq=lagos.identify_contours(sp, con_field, mi, sp[con_field].max())
+ for cid in mq:
+ sp["Contours"][mq[cid]] = cid
+ k.append(na.unique(sp["Contours"]).size)
+ my_con = sp["Contours"]
+ sp._flush_data_to_grids("Contours",-1)
+ for cid in na.unique(sp["Contours"])[1:]:
+ cid_ind = na.where(sp["Contours"] == cid)
+ print "\tContour id %s has: %0.5e %0.5e (%s) (%s grids, %s %s)" % \
+ (cid, sp[con_field][cid_ind].min(), sp[con_field][cid_ind].max(),
+ sp["Density"][cid_ind].size,
+ na.unique(sp["GridIndices"][cid_ind]).size,
+ na.unique(sp["GridIndices"][cid_ind]).min(),
+ na.unique(sp["GridIndices"][cid_ind]).max(),
+ )
+ #sp = a.h.sphere(c,20000.0/a['au'])
+ sp.get_data("Contours",in_grids=True)
+ if not na.all(sp["Contours"] == my_con):
+ khdoihiueh0
+ g_bad = 0
+ for grid in sp._grids:
+ grid['dx'] = na.ones(grid.ActiveDimensions)*grid.dx
+ grid['dy'] = na.ones(grid.ActiveDimensions)*grid.dy
+ grid['dz'] = na.ones(grid.ActiveDimensions)*grid.dz
+ g_bad+=lagos.check_neighbors(grid)
+ print "GRID:",g_bad
+ k_bad.append(lagos.check_neighbors(sp))
+ print "SOURCE",k_bad[-1]
+ if do_plot:
+ #c=na.array([sp["x"][ind][-1],sp["y"][ind][-1],sp["z"][ind][-1]])
+ pc = raven.PlotCollection(a, center=c)
+ pc.add_slice("Density",0)
+ pc.add_slice("Density",1)
+ pc.add_slice("Density",2)
+ for field in ["Density","Contours"]:
+ pc.switch_field(field)
+ pc.set_width(10000,'au')
+ pc.save("test_bds_%02i_10000au" % i)
+ pc.set_width(2000,'au')
+ pc.save("test_bds_%02i_02000au" % i)
+ pc.set_width(1000,'au')
+ pc.save("test_bds_%02i_01000au" % i)
+ pc.set_width(100,'au')
+ pc.save("test_bds_%02i_00100au" % i)
+ del sp.data["Contours"]
+print k
+print k_bad
Modified: trunk/tests/test_lagos.py
==============================================================================
--- trunk/tests/test_lagos.py (original)
+++ trunk/tests/test_lagos.py Fri Apr 11 00:20:30 2008
@@ -108,6 +108,27 @@
p = self.hierarchy.proj(axis, "Ones") # Derived field
self.assertTrue(na.all(p["Ones"] == 1.0))
+ def testContours(self):
+ # As a note, unfortunately this dataset only has one sphere.
+ # Frownie face.
+ dx = self.hierarchy.grids[0].dx
+ reg = self.hierarchy.region([0.5]*3, [2*dx]*3,[1.0-2*dx]*3)
+ cid = yt.lagos.identify_contours(reg, "Density",
+ reg["Density"].min()*0.99, reg["Density"].max()*1.01)
+ self.assertEqual(len(cid), 1)
+ v1 = reg["Density"].max()*0.99
+ v2 = reg["Density"].max()*1.01
+ cid = yt.lagos.identify_contours(reg, "Density", v1, v2)
+ self.assertTrue(na.all(v1 < reg["Density"][cid[0]])
+ and na.all(v2 > reg["Density"][cid[0]]))
+ self.assertEqual(len(cid), 1)
+ v1 = reg["Density"].min()*0.99
+ v2 = reg["Density"].min()*1.01
+ cid = yt.lagos.identify_contours(reg, "Density", v1, v2)
+ self.assertTrue(na.all(v1 < reg["Density"][cid[0]])
+ and na.all(v2 > reg["Density"][cid[0]]))
+ self.assertEqual(len(cid), 1)
+
# Now we test each datatype in turn
def _returnFieldFunction(field):
@@ -190,9 +211,9 @@
for g in self.hierarchy.grids:
cube = g.retrieve_ghost_zones(0, "Density")
self.assertTrue(na.all(cube["Density"] == g["Density"]))
- cube["Density"] = na.ones(cube["Density"].shape)
+ cube["Density"] = na.arange(cube["Density"].size).reshape(cube["Density"].shape)
cube.flush_data(field="Density")
- self.assertTrue(na.all(g["Density"] == 1.0))
+ self.assertTrue(na.all(g["Density"] == cube["Density"]))
def testTwoGhost(self):
for g in self.hierarchy.grids:
@@ -204,7 +225,9 @@
cg.flush_data(field="Ones")
for g in self.hierarchy.grids:
self.assertTrue(g["Ones"].max() == 2.0)
-
+ cg2 = self.hierarchy.covering_grid(3, [0.0]*3, [1.0]*3, [64,64,64])
+ self.assertTrue(na.all(cg["Ones"] == cg2["Ones"]))
+
def testAllCover(self):
cg = self.hierarchy.covering_grid(0, [0.0]*3, [1.0]*3, [32,32,32])
self.assertTrue(cg["Density"].max() \
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Fri Apr 11 00:20:30 2008
@@ -1457,7 +1457,7 @@
continue
mylog.debug("Getting field %s from %s possible grids",
field, len(self._grids))
- self[field] = na.zeros(self.ActiveDimensions, dtype='float64') - 999
+ self[field] = na.zeros(self.ActiveDimensions, dtype='float64') -999
if self._use_pbar: pbar = \
get_pbar('Searching grids for values ', len(self._grids))
for i,grid in enumerate(self._grids):
Modified: trunk/yt/lagos/ContourFinder.py
==============================================================================
--- trunk/yt/lagos/ContourFinder.py (original)
+++ trunk/yt/lagos/ContourFinder.py Fri Apr 11 00:20:30 2008
@@ -91,7 +91,10 @@
my_queue.add(data_source._grids)
for i,grid in enumerate(my_queue):
max_before = grid["tempContours"].max()
- cg = grid.retrieve_ghost_zones(1,["tempContours","GridIndices"], all_levels=False)
+ if na.all(grid.LeftEdge == 0.0) and na.all(grid.RightEdge == 1.0):
+ cg = grid.retrieve_ghost_zones(0,["tempContours","GridIndices"])
+ else:
+ cg = grid.retrieve_ghost_zones(1,["tempContours","GridIndices"])
local_ind = na.where( (cg[field] > min_val)
& (cg[field] < max_val)
& (cg["tempContours"] == -1) )
@@ -106,13 +109,7 @@
xi = xi_u[cor_order]
yi = yi_u[cor_order]
zi = zi_u[cor_order]
- np = xi.size
- mi, mj, mk = fd.shape
- weave.inline(iterate_over_contours,
- ['xi','yi','zi','np','fd', 'mi','mj','mk'],
- compiler='gcc', type_converters=converters.blitz,
- auto_downcast=0, verbose=2, force=0,
- support_code=recursively_find_contours)
+ PointCombine.FindContours(fd, xi, yi, zi)
cg["tempContours"] = fd.copy()
cg.flush_data("tempContours")
my_queue.add(cg._grids)
Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c (original)
+++ trunk/yt/lagos/PointCombine.c Fri Apr 11 00:20:30 2008
@@ -583,7 +583,111 @@
return to_return;
}
+static PyObject *_findContoursError;
+static PyObject *
+Py_FindContours(PyObject *obj, PyObject *args)
+{
+ PyObject *ocon_ids, *oxi, *oyi, *ozi;
+ PyArrayObject *con_ids, *xi, *yi, *zi;
+ int i, j, k, n;
+
+ i = 0;
+ if (!PyArg_ParseTuple(args, "OOOO",
+ &ocon_ids, &oxi, &oyi, &ozi))
+ return PyErr_Format(_findContoursError,
+ "FindContours: Invalid parameters.");
+
+ con_ids = (PyArrayObject *) PyArray_FromAny(ocon_ids,
+ PyArray_DescrFromType(NPY_FLOAT64), 3, 3,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((con_ids==NULL) || (con_ids->nd != 3)) {
+ PyErr_Format(_findContoursError,
+ "FindContours: Three dimensions required for con_ids.");
+ goto _fail;
+ }
+
+ xi = (PyArrayObject *) PyArray_FromAny(oxi,
+ PyArray_DescrFromType(NPY_INT32), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if(xi==NULL) {
+ PyErr_Format(_findContoursError,
+ "Bin2DProfile: One dimension required for xi.");
+ goto _fail;
+ }
+
+ yi = (PyArrayObject *) PyArray_FromAny(oyi,
+ PyArray_DescrFromType(NPY_INT32), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((yi==NULL) || (PyArray_SIZE(xi) != PyArray_SIZE(yi))) {
+ PyErr_Format(_findContoursError,
+ "Bin2DProfile: One dimension required for yi, same size as xi.");
+ goto _fail;
+ }
+
+ zi = (PyArrayObject *) PyArray_FromAny(ozi,
+ PyArray_DescrFromType(NPY_INT32), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((zi==NULL) || (PyArray_SIZE(xi) != PyArray_SIZE(zi))) {
+ PyErr_Format(_findContoursError,
+ "Bin2DProfile: One dimension required for zi, same size as xi.");
+ goto _fail;
+ }
+
+ for(n=0;n<PyArray_SIZE(xi);n++) {
+ i=xi->data[n];j=yi->data[n];k=zi->data[n];
+ process_neighbors(con_ids, i, j, k);
+ }
+
+ Py_DECREF(con_ids);
+ Py_DECREF(xi);
+ Py_DECREF(yi);
+ Py_DECREF(zi);
+
+ PyObject *status = PyInt_FromLong(1);
+ return status;
+
+ _fail:
+ Py_XDECREF(con_ids);
+ Py_XDECREF(xi);
+ Py_XDECREF(yi);
+ Py_XDECREF(zi);
+ return NULL;
+}
+
+int process_neighbors(PyArrayObject *con_ids, int i, int j, int k)
+{
+ int off_i, off_j, off_k, spawn_check;
+ int mi, mj, mk;
+ npy_float64 val1, val2, new_cid;
+ mi = con_ids->dimensions[0];
+ mj = con_ids->dimensions[1];
+ mk = con_ids->dimensions[2];
+ do {
+ spawn_check = 0;
+ for (off_i=max(i-1,0);off_i<=min(i+1,mi-1);off_i++)
+ for (off_j=max(j-1,0);off_j<=min(j+1,mj-1);off_j++)
+ for (off_k=max(k-1,0);off_k<=min(k+1,mk-1);off_k++) {
+ if((off_i==i)&&(off_j==j)&&(off_k==k)) continue;
+ val1 = *((npy_float64*)PyArray_GETPTR3(con_ids, off_i, off_j, off_k));
+ val2 = *((npy_float64*)PyArray_GETPTR3(con_ids, i, j, k));
+ if(val1!=-1) {
+ if(val1 > val2){
+ val2 = *((npy_float64*)PyArray_GETPTR3(con_ids, i,j,k)) = val1;
+ spawn_check += 1;
+ }
+ if(val1 < val2){
+ *((npy_float64*)PyArray_GETPTR3(con_ids, off_i,off_j,off_k)) = val2;
+ new_cid = process_neighbors(con_ids, off_i, off_j, off_k);
+ if (new_cid != val2) spawn_check += 1;
+ *((npy_float64*)PyArray_GETPTR3(con_ids, i,j,k)) = new_cid;
+ }
+ }
+ }
+ } while (spawn_check > 0);
+ val2 = *((npy_float64*)PyArray_GETPTR3(con_ids, i, j, k));
+ return val2;
+}
static PyObject *_interpolateError;
@@ -689,6 +793,7 @@
{"DataCubeRefine", Py_DataCubeRefine, METH_VARARGS},
{"DataCubeReplace", Py_DataCubeReplace, METH_VARARGS},
{"Bin2DProfile", Py_Bin2DProfile, METH_VARARGS},
+ {"FindContours", Py_FindContours, METH_VARARGS},
{NULL, NULL} /* Sentinel */
};
@@ -710,6 +815,8 @@
PyDict_SetItemString(d, "error", _dataCubeError);
_profile2DError = PyErr_NewException("PointCombine.Profile2DError", NULL, NULL);
PyDict_SetItemString(d, "error", _profile2DError);
+ _findContoursError = PyErr_NewException("PointCombine.FindContoursError", NULL, NULL);
+ PyDict_SetItemString(d, "error", _findContoursError);
import_array();
}
More information about the yt-svn
mailing list