[Yt-svn] yt-commit r396 - in trunk: tests yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Thu Apr 10 15:24:44 PDT 2008
Author: mturk
Date: Thu Apr 10 15:24:43 2008
New Revision: 396
URL: http://yt.spacepope.org/changeset/396
Log:
Removed unused strings from Weave.
Added 2D profiling to PointCombine, which should probably be renamed sometime
soon.
Have not rigorously checked for memory leaks, but from my unit test
observations it should be clean.
Removed weave cruft from the datacube object.
Added a test for making sure all profiles return the sum total of the mass if
you added an accumulating mass field. I would like to include a way to check
for correctness later on.
Modified:
trunk/tests/test_lagos.py
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/PointCombine.c
trunk/yt/lagos/Profiles.py
trunk/yt/lagos/WeaveStrings.py
Modified: trunk/tests/test_lagos.py
==============================================================================
--- trunk/tests/test_lagos.py (original)
+++ trunk/tests/test_lagos.py Thu Apr 10 15:24:43 2008
@@ -147,7 +147,15 @@
LagosTestingBase.setUp(self)
class Data3DBase:
- pass
+ def testProfileAccumulateMass(self):
+ self.data.set_field_parameter("center",[0.5]*3)
+ profile = yt.lagos.BinnedProfile1D(self.data, 8, "RadiusCode", 0, 1.0,
+ False, True)
+ profile.add_fields("CellMassMsun", weight=None, accumulation=True)
+ v1 = profile["CellMassMsun"].max()
+ v2 = self.data["CellMassMsun"].sum()
+ v2 = na.abs(1.0 - v2/v1)
+ self.assertAlmostEqual(v2, 0.0, 7)
for field in yt.lagos.fieldInfo.values():
setattr(DataTypeTestingBase, "test%s" % field.name, _returnFieldFunction(field))
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Thu Apr 10 15:24:43 2008
@@ -1439,27 +1439,6 @@
def extract_region(self, indices):
mylog.error("Sorry, dude, do it yourself, it's already in 3-D.")
- def __setup_weave_dict(self, grid):
- return {
- 'nx_g': int(grid.ActiveDimensions[0]),
- 'ny_g': int(grid.ActiveDimensions[1]),
- 'nz_g': int(grid.ActiveDimensions[2]),
- 'leftEdgeGrid': na.array(grid.LeftEdge),
- 'rf': int(grid.dx / self.dx),
- 'dx_g': float(grid.dx),
- 'dy_g': float(grid.dy),
- 'dz_g': float(grid.dz),
- 'dx_m': float(self.dx),
- 'dy_m': float(self.dy),
- 'dz_m': float(self.dz),
- 'leftEdgeCube': na.array(self.left_edge),
- 'cubeRightEdge': na.array(self.right_edge),
- 'nx_m': int(self.ActiveDimensions[0]),
- 'ny_m': int(self.ActiveDimensions[1]),
- 'nz_m': int(self.ActiveDimensions[2]),
- 'childMask' : grid.child_mask
- }
-
def _refresh_data(self):
Enzo3DData._refresh_data(self)
self['dx'] = self.dx * na.ones(self.ActiveDimensions, dtype='float64')
@@ -1529,29 +1508,3 @@
grid.LeftEdge, g_dx, grid[field], grid.child_mask,
self.left_edge, self.right_edge, c_dx, self[field],
ll)
-
- @restore_grid_state
- def old_get_data_from_grid(self, grid, fields):
- for field in ensure_list(fields):
- locals_dict = self.__setup_weave_dict(grid)
- locals_dict['fieldData'] = grid[field]
- locals_dict['cubeData'] = self[field]
- locals_dict['lastLevel'] = int(grid.Level == self.level)
- weave.inline(DataCubeRefineCoarseData,
- locals_dict.keys(), local_dict=locals_dict,
- compiler='gcc',
- type_converters=converters.blitz,
- auto_downcast=0, verbose=2)
-
- def old_flush_data_to_grid(self, grid, fields):
- for field in ensure_list(fields):
- locals_dict = self.__setup_weave_dict(grid)
- locals_dict['fieldData'] = grid[field].copy()
- locals_dict['cubeData'] = self[field]
- locals_dict['lastLevel'] = int(grid.Level == self.level)
- weave.inline(DataCubeReplaceData,
- locals_dict.keys(), local_dict=locals_dict,
- compiler='gcc',
- type_converters=converters.blitz,
- auto_downcast=0, verbose=2)
- grid[field] = locals_dict['fieldData'].copy()
Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c (original)
+++ trunk/yt/lagos/PointCombine.c Thu Apr 10 15:24:43 2008
@@ -257,8 +257,7 @@
}
/* These functions are both called with
- func(cubedata, griddata)
-*/
+ func(cubedata, griddata) */
static void dcRefine(npy_float64 *val1, npy_float64 *val2) {
*val1 = *val2;
}
@@ -267,13 +266,127 @@
*val2 = *val1;
}
+static PyObject *_profile2DError;
+
+static PyObject *Py_Bin2DProfile(PyObject *obj, PyObject *args)
+{
+ int i, j;
+ PyObject *obins_x, *obins_y, *owsource, *obsource, *owresult, *obresult, *oused;
+ PyArrayObject *bins_x, *bins_y, *wsource, *bsource, *wresult, *bresult, *used;
+
+ if (!PyArg_ParseTuple(args, "OOOOOOO",
+ &obins_x, &obins_y, &owsource, &obsource,
+ &owresult, &obresult, &oused))
+ return PyErr_Format(_profile2DError,
+ "Bin2DProfile: Invalid parameters.");
+ i = 0;
+
+ bins_x = (PyArrayObject *) PyArray_FromAny(obins_x,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if(bins_x==NULL) {
+ PyErr_Format(_profile2DError,
+ "Bin2DProfile: One dimension required for bins_x.");
+ goto _fail;
+ }
+
+ bins_y = (PyArrayObject *) PyArray_FromAny(obins_y,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((bins_y==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(bins_y))) {
+ PyErr_Format(_profile2DError,
+ "Bin2DProfile: One dimension required for bins_y, same size as bins_x.");
+ goto _fail;
+ }
+
+ wsource = (PyArrayObject *) PyArray_FromAny(owsource,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((wsource==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(wsource))) {
+ PyErr_Format(_profile2DError,
+ "Bin2DProfile: One dimension required for wsource, same size as bins_x.");
+ goto _fail;
+ }
+
+ bsource = (PyArrayObject *) PyArray_FromAny(obsource,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((bsource==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(bsource))) {
+ PyErr_Format(_profile2DError,
+ "Bin2DProfile: One dimension required for bsource, same size as bins_x.");
+ goto _fail;
+ }
+
+ wresult = (PyArrayObject *) PyArray_FromAny(owresult,
+ PyArray_DescrFromType(NPY_FLOAT64), 2,2,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if(wresult==NULL){
+ PyErr_Format(_profile2DError,
+ "Bin2DProfile: Two dimensions required for wresult.");
+ goto _fail;
+ }
+
+ bresult = (PyArrayObject *) PyArray_FromAny(obresult,
+ PyArray_DescrFromType(NPY_FLOAT64), 2,2,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((bresult==NULL) ||(PyArray_SIZE(wresult) != PyArray_SIZE(bresult))
+ || (PyArray_DIM(bresult,0) != PyArray_DIM(wresult,0))){
+ PyErr_Format(_profile2DError,
+ "Bin2DProfile: Two dimensions required for bresult, same shape as wresult.");
+ goto _fail;
+ }
+
+ used = (PyArrayObject *) PyArray_FromAny(oused,
+ PyArray_DescrFromType(NPY_FLOAT64), 2,2,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((used==NULL) ||(PyArray_SIZE(used) != PyArray_SIZE(wresult))
+ || (PyArray_DIM(used,0) != PyArray_DIM(wresult,0))){
+ PyErr_Format(_profile2DError,
+ "Bin2DProfile: Two dimensions required for used, same shape as wresult.");
+ goto _fail;
+ }
+
+ npy_float64 wv;
+ int n;
+
+ for(n=0; n<bins_x->dimensions[0]; n++) {
+ i = bins_x->data[n];
+ j = bins_y->data[n];
+ *(npy_float64*)PyArray_GETPTR2(wresult, i, j) +=
+ wsource->data[n];
+ *(npy_float64*)PyArray_GETPTR2(bresult, i, j) +=
+ wsource->data[n] * bsource->data[n];
+ *(npy_float64*)PyArray_GETPTR2(used, i, j) = 1.0;
+ }
+
+ Py_DECREF(bins_x);
+ Py_DECREF(bins_y);
+ Py_DECREF(wsource);
+ Py_DECREF(bsource);
+ Py_DECREF(wresult);
+ Py_DECREF(bresult);
+ Py_DECREF(used);
+
+ PyObject *onum_found = PyInt_FromLong((long)1);
+ return onum_found;
+
+ _fail:
+ Py_XDECREF(bins_x);
+ Py_XDECREF(bins_y);
+ Py_XDECREF(wsource);
+ Py_XDECREF(bsource);
+ Py_XDECREF(wresult);
+ Py_XDECREF(bresult);
+ Py_XDECREF(used);
+ return NULL;
+
+}
+
static PyObject *_dataCubeError;
static PyObject *DataCubeGeneric(PyObject *obj, PyObject *args,
void (*to_call)(npy_float64*,npy_float64*))
{
-#define DUMB
-#ifdef DUMB
/* Standard boilerplate unpacking... */
/*
@@ -454,7 +567,6 @@
Py_XDECREF(c_data);
return NULL;
-#endif
}
static PyObject *
@@ -471,6 +583,8 @@
return to_return;
}
+
+
static PyObject *_interpolateError;
static void
@@ -574,6 +688,7 @@
{"Interpolate", Py_Interpolate, METH_VARARGS},
{"DataCubeRefine", Py_DataCubeRefine, METH_VARARGS},
{"DataCubeReplace", Py_DataCubeReplace, METH_VARARGS},
+ {"Bin2DProfile", Py_Bin2DProfile, METH_VARARGS},
{NULL, NULL} /* Sentinel */
};
@@ -593,6 +708,8 @@
PyDict_SetItemString(d, "error", _interpolateError);
_dataCubeError = PyErr_NewException("PointCombine.DataCubeError", NULL, NULL);
PyDict_SetItemString(d, "error", _dataCubeError);
+ _profile2DError = PyErr_NewException("PointCombine.Profile2DError", NULL, NULL);
+ PyDict_SetItemString(d, "error", _profile2DError);
import_array();
}
Modified: trunk/yt/lagos/Profiles.py
==============================================================================
--- trunk/yt/lagos/Profiles.py (original)
+++ trunk/yt/lagos/Profiles.py Thu Apr 10 15:24:43 2008
@@ -207,19 +207,8 @@
self.total_cells += bin_indices_x.size
nx = bin_indices_x.size
mylog.debug("Binning %s / %s times", source_data.size, nx)
- try:
- weave.inline(_2d_profile_code, ['nx','bin_indices_x','bin_indices_y',
- 'weight_field','weight_data',
- 'binned_field','source_data', 'used_field'],
- compiler='gcc', type_converters=converters.blitz,
- auto_downcast = 0)
- except:
- mylog.debug("SciPy weaving failed; non-fatal, but slow.")
- for k in range(nx):
- i,j = bin_indices_x[k]-1, bin_indices_y[k]-1
- weight_field[i,j] += weight_data[k]
- binned_field[i,j] += source_data[k]*weight_data[k]
- used_field[i,j] = 1.0
+ PointCombine.Bin2DProfile(bin_indices_x, bin_indices_y, weight_data, source_data,
+ weight_field, binned_field, used_field)
if accumulation: # Fix for laziness
if not iterable(accumulation):
raise SyntaxError("Accumulation needs to have length 2")
Modified: trunk/yt/lagos/WeaveStrings.py
==============================================================================
--- trunk/yt/lagos/WeaveStrings.py (original)
+++ trunk/yt/lagos/WeaveStrings.py Thu Apr 10 15:24:43 2008
@@ -23,143 +23,6 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
-ProfileBinningWeighted = \
-"""
-int bin;
-int element;
-int mybin;
-npy_float64 myweight; // should be higher precision
-npy_float64 weights[num_bins]; // should be higher precision
-
-for (bin = 0; bin < num_bins; bin++) {
- weights[bin] = 0;
- profilevalues(bin) = 0;
-}
-
-for (element = 0; element < num_elements ; element++) {
- mybin = binindices(element);
- profilevalues(mybin) += weightvalues(mybin)*fieldvalues(element);
- weights[mybin] += weightvalues(mybin);
-}
-
-for (bin = 0 ; bin < num_bins ; bin++ ) {
- profilevalues(bin) = profilevalues(bin) / weights[bin];
-}
-"""
-
-
-ProfileBinningAccumulation = \
-"""
-int bin;
-int element;
-int mybin;
-npy_float64 myweight; // should be higher precision
-npy_float64 weights[num_bins]; // should be higher precision
-
-for (bin = 0; bin < num_bins; bin++) {
- weights[bin] = 0;
- profilevalues(bin) = 0;
-}
-
-for (element = 0; element < num_elements ; element++) {
- mybin = binindices(element);
- profilevalues(mybin) += fieldvalues(element);
-}
-
-for (bin = 1 ; bin < num_bins ; bin++ ) {
- profilevalues(bin) += profilevalues(bin-1);
-}
-"""
-
-ProjectionCombineSameLevel = \
-r"""
- int g1i, g2i;
-
- for (g1i = 0; g1i < g1points; g1i++) {
- if (g1data_x(g1i) < 0) continue;
- for (g2i = 0; g2i < g2points; g2i++) {
- if (g2data_x(g2i) < 0) continue;
- if ((g1data_x(g1i) == g2data_x(g2i)) &&
- (g1data_y(g1i) == g2data_y(g2i))) {
- g1data_vals(g1i) += g2data_vals(g2i);
- g1data_wgt(g1i) += g2data_wgt(g2i);
- g1data_mask(g1i) = ((g1data_mask(g1i)) && (g2data_mask(g2i)));
- g2data_x(g2i) = -1;
- break; // We map to one g1data point AT MOST
- }
- }
- }
-"""
-
-ProjectionRefineCoarseData = \
-r"""
- int ci, fi, x_off, y_off;
- long int fine_x, fine_y;
-
- for (ci = 0; ci < cpoints; ci++) {
- for (x_off = 0; x_off < rf; x_off++) {
- for (y_off = 0; y_off < rf; y_off++) {
- fine_x = 2*coarsedata_x(ci) + x_off;
- fine_y = 2*coarsedata_y(ci) + y_off;
- for (fi = 0; fi < fpoints; fi++) {
- if ((fine_x == finedata_x(fi)) &&
- (fine_y == finedata_y(fi))) {
- finedata_vals(fi) += coarsedata_vals(ci);
- finedata_wgt(fi) += coarsedata_wgt(ci);
- flagged(ci) += 1;
- }
- }
- }
- }
- }
-
-"""
-
-_baseDataCubeRefineCoarseData = \
-r"""
-
-#define MIN(A,B) ((A) < (B) ? (A) : (B))
-#define MAX(A,B) ((A) < (B) ? (B) : (A))
-
-int bx_m, by_m, bz_m, xg, yg, zg, xm, ym, zm;
-int ixm, iym, izm;
-int max_bx_m, max_by_m, max_bz_m;
-using namespace std;
-
-for (xg = 0; xg < nx_g; xg++) {
- if (leftEdgeGrid(0)+dx_g*xg > cubeRightEdge(0)) continue;
- if (leftEdgeGrid(0)+dx_g*(xg+1) < leftEdgeCube(0)) continue;
- bx_m = MAX(floorl((leftEdgeGrid(0)+dx_g*xg - leftEdgeCube(0))/dx_m),0);
- max_bx_m = MIN(ceill((leftEdgeGrid(0)+dx_g*(xg+1) - leftEdgeCube(0))/dx_m),nx_m);
- for (yg = 0; yg < ny_g; yg++) {
- if (leftEdgeGrid(1)+dy_g*yg > cubeRightEdge(1)) continue;
- if (leftEdgeGrid(1)+dy_g*(yg+1) < leftEdgeCube(1)) continue;
- by_m = MAX(floorl((leftEdgeGrid(1)+dy_g*yg - leftEdgeCube(1))/dy_m),0);
- max_by_m = MIN(ceill((leftEdgeGrid(1)+dy_g*(yg+1) - leftEdgeCube(1))/dy_m),ny_m);
- for (zg = 0; zg < nz_g; zg++) {
- if (!lastLevel)
- if (childMask(xg, yg, zg) == 0) continue;
- if (leftEdgeGrid(2)+dz_g*zg > cubeRightEdge(2)) continue;
- if (leftEdgeGrid(2)+dz_g*(zg+1) < leftEdgeCube(2)) continue;
- bz_m = MAX(floorl((leftEdgeGrid(2)+dz_g*zg - leftEdgeCube(2))/dz_m),0);
- max_bz_m = MIN(ceill((leftEdgeGrid(2)+dz_g*(zg+1) - leftEdgeCube(2))/dz_m),nz_m);
- for (xm = bx_m; xm < max_bx_m ; xm++) {
- for (ym = by_m; ym < max_by_m ; ym++) {
- for (zm = bz_m; zm < max_bz_m ; zm++) {
- %s
- }
- }
- }
- }
- }
-}
-"""
-
-DataCubeRefineCoarseData = _baseDataCubeRefineCoarseData % \
- ("cubeData(xm,ym,zm) = fieldData(xg,yg,zg);")
-DataCubeReplaceData = _baseDataCubeRefineCoarseData % \
- ("fieldData(xg,yg,zg) = cubeData(xm,ym,zm);")
-
iterate_over_contours = r"""
int i,j,k,n;
for(n=0; n<np; n++){
More information about the yt-svn
mailing list