[Yt-svn] yt-commit r487 - trunk/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Sat May 17 21:22:05 PDT 2008
Author: mturk
Date: Sat May 17 21:22:04 2008
New Revision: 487
URL: http://yt.spacepope.org/changeset/487
Log:
Okay, so after much tracking-down, I have determined the issue.
Now, weight_field returns a proper, averaged value, based on weights.
This can be tested pretty easily by looking at a projection of the field 'Ones'
weighted arbitrarily -- it always returns 1.0.
Additionally, non-weighted projections retain correct results.
For an example, give a shot at a Temperature projection weighted by Density.
Thing is, I'm not sure this works for all types of fields, or all types of
weights. I'm going to try to figure out which types of fields it does work for
and which it doesn't, and I will utilize some kind of property in the fieldInfo
object to discriminate between fields that need varying powers of dx to give
correct results.
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/PointCombine.c
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Sat May 17 21:22:04 2008
@@ -758,21 +758,30 @@
self[self.fields[0]] = array[4,:]
return True
+ def __get_dls(self, grid, fields):
+ dls = []
+ convs = []
+ for fi, field in enumerate(fields):
+ if field in fieldInfo and not fieldInfo[field].line_integral:
+ dls.append(1.0)
+ convs.append(1.0)
+ else:
+ dls.append(just_one(grid['d%s' % axis_names[self.axis]]))
+ convs.append(self.pf.units[fieldInfo[field].projection_conversion])
+ return na.array(dls), na.array(convs)
+
def __project_level(self, level, fields):
grids_to_project = self.source.select_grids(level)
+ dls, convs = self.__get_dls(grids_to_project[0], fields)
zero_out = (level != self._max_level)
pbar = get_pbar('Projecting level % 2i / % 2i ' \
% (level, self._max_level), len(grids_to_project))
for pi, grid in enumerate(grids_to_project):
g_coords, g_fields = self._project_grid(grid, fields, zero_out)
- for fi, field in enumerate(fields):
- dl = 1.0
- if field in fieldInfo and fieldInfo[field].line_integral:
- dl = just_one(grid['d%s' % axis_names[self.axis]])
- dl *= self.pf.units[fieldInfo[field].projection_conversion]
- g_fields[fi] *= dl
self.__retval_coords[grid.id] = g_coords
self.__retval_fields[grid.id] = g_fields
+ if self._weight is None:
+ for fi in range(len(fields)): g_fields[fi] *= dls[fi]
pbar.update(pi)
pbar.finish()
self.__combine_grids_on_level(level) # In-place
@@ -789,8 +798,10 @@
self.__retval_fields[grid.id] = [pi[coarse] for pi in self.__retval_fields[grid.id]]
coord_data = na.concatenate(coord_data, axis=1)
field_data = na.concatenate(field_data, axis=1)
- if self._weight != None:
+ if self._weight is not None:
field_data = field_data / coord_data[3,:].reshape((1,coord_data.shape[1]))
+ else:
+ field_data *= convs
mylog.info("Level %s done: %s final", \
level, coord_data.shape[1])
dx = grids_to_project[0].dx# * na.ones(coord_data.shape[0], dtype='float64')
@@ -896,20 +907,20 @@
self.data['pdy'] = self.data['pdx'].copy()
for fi, field in enumerate(fields):
self[field] = field_data[fi,:]
+ self.data['weight_field'] = coord_data[3,:]
def add_fields(self, fields, weight = "CellMassMsun"):
pass
def _project_grid(self, grid, fields, zero_out):
if self._weight is None:
- weight_data = na.ones(grid.ActiveDimensions)
+ weight_data = na.ones(grid.ActiveDimensions, dtype='float64')
else:
- weight_data = self._get_data_from_grid(grid, self._weight)
+ weight_data = self._get_data_from_grid(grid, self._weight).astype('float64')
if zero_out: weight_data[grid.child_indices] = 0
# if we zero it out here, then we only have to zero out the weight!
masked_data = [self._get_data_from_grid(grid, field) * weight_data
for field in fields]
- dl = 1.0
full_proj = [self.func(field,axis=self.axis) for field in masked_data]
weight_proj = self.func(weight_data,axis=self.axis)
if self._check_region and not self.source._is_fully_enclosed(grid):
Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c (original)
+++ trunk/yt/lagos/PointCombine.c Sat May 17 21:22:04 2008
@@ -66,12 +66,12 @@
/* First the regular source arrays */
grid_src_x = (PyArrayObject *) PyArray_FromAny(ogrid_src_x,
- PyArray_DescrFromType(NPY_INT64), 1, 0,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
src_len = PyArray_SIZE(grid_src_x);
grid_src_y = (PyArrayObject *) PyArray_FromAny(ogrid_src_y,
- PyArray_DescrFromType(NPY_INT64), 1, 0,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
if(PyArray_SIZE(grid_src_y) != src_len) {
PyErr_Format(_combineGridsError,
@@ -80,7 +80,7 @@
}
grid_src_mask = (PyArrayObject *) PyArray_FromAny(ogrid_src_mask,
- PyArray_DescrFromType(NPY_INT64), 1, 0,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
if(PyArray_SIZE(grid_src_mask) != src_len) {
PyErr_Format(_combineGridsError,
@@ -89,21 +89,21 @@
}
grid_src_wgt = (PyArrayObject *) PyArray_FromAny(ogrid_src_wgt,
- PyArray_DescrFromType(NPY_FLOAT64), 1, 0,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
- if(PyArray_SIZE(grid_src_wgt) != src_len) {
+ if((grid_src_wgt == NULL) || (PyArray_SIZE(grid_src_wgt) != src_len)) {
PyErr_Format(_combineGridsError,
"CombineGrids: src_x and src_wgt must be the same shape.");
goto _fail;
}
grid_dst_x = (PyArrayObject *) PyArray_FromAny(ogrid_dst_x,
- PyArray_DescrFromType(NPY_INT64), 1, 0,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
dst_len = PyArray_SIZE(grid_dst_x);
grid_dst_y = (PyArrayObject *) PyArray_FromAny(ogrid_dst_y,
- PyArray_DescrFromType(NPY_INT64), 1, 0,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
if(PyArray_SIZE(grid_dst_y) != dst_len) {
PyErr_Format(_combineGridsError,
@@ -112,7 +112,7 @@
}
grid_dst_mask = (PyArrayObject *) PyArray_FromAny(ogrid_dst_mask,
- PyArray_DescrFromType(NPY_INT64), 1, 0,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
if(PyArray_SIZE(grid_dst_mask) != dst_len) {
PyErr_Format(_combineGridsError,
@@ -121,9 +121,9 @@
}
grid_dst_wgt = (PyArrayObject *) PyArray_FromAny(ogrid_dst_wgt,
- PyArray_DescrFromType(NPY_FLOAT64), 1, 0,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
- if(PyArray_SIZE(grid_dst_wgt) != dst_len) {
+ if((grid_dst_wgt == NULL) || (PyArray_SIZE(grid_dst_wgt) != dst_len)) {
PyErr_Format(_combineGridsError,
"CombineGrids: dst_x and dst_wgt must be the same shape.");
goto _fail;
@@ -154,7 +154,7 @@
temp_object,
PyArray_DescrFromType(NPY_FLOAT64), 1, 0,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
- src_vals[i] = (npy_float64 *) grid_src_vals[i]->data;
+ src_vals[i] = (npy_float64 *) PyArray_GETPTR1(grid_src_vals[i],0);
Py_DECREF(temp_object);
temp_object = PySequence_GetItem(ogrid_dst_vals, i);
@@ -162,21 +162,21 @@
temp_object,
PyArray_DescrFromType(NPY_FLOAT64), 1, 0,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
- dst_vals[i] = (npy_float64 *) grid_dst_vals[i]->data;
+ dst_vals[i] = (npy_float64 *) PyArray_GETPTR1(grid_dst_vals[i],0);
Py_DECREF(temp_object);
}
/* Now we're all set to call our sub-function. */
- npy_int64 *src_x = (npy_int64 *) grid_src_x->data;
- npy_int64 *src_y = (npy_int64 *) grid_src_y->data;
- npy_float64 *src_wgt = (npy_float64 *) grid_src_wgt->data;
- npy_int64 *src_mask = (npy_int64 *) grid_src_mask->data;
-
- npy_int64 *dst_x = (npy_int64 *) grid_dst_x->data;
- npy_int64 *dst_y = (npy_int64 *) grid_dst_y->data;
- npy_float64 *dst_wgt = (npy_float64 *) grid_dst_wgt->data;
- npy_int64 *dst_mask = (npy_int64 *) grid_dst_mask->data;
+ npy_int64 *src_x = (npy_int64 *) PyArray_GETPTR1(grid_src_x,0);
+ npy_int64 *src_y = (npy_int64 *) PyArray_GETPTR1(grid_src_y,0);
+ npy_float64 *src_wgt = (npy_float64 *) PyArray_GETPTR1(grid_src_wgt,0);
+ npy_int64 *src_mask = (npy_int64 *) PyArray_GETPTR1(grid_src_mask,0);
+
+ npy_int64 *dst_x = (npy_int64 *) PyArray_GETPTR1(grid_dst_x,0);
+ npy_int64 *dst_y = (npy_int64 *) PyArray_GETPTR1(grid_dst_y,0);
+ npy_float64 *dst_wgt = (npy_float64 *) PyArray_GETPTR1(grid_dst_wgt,0);
+ npy_int64 *dst_mask = (npy_int64 *) PyArray_GETPTR1(grid_dst_mask,0);
int si, di, x_off, y_off;
npy_int64 fine_x, fine_y, init_x, init_y;
@@ -195,12 +195,12 @@
if ((fine_x == dst_x[di]) &&
(fine_y == dst_y[di])) {
num_found++;
- dst_wgt[di] += src_wgt[di];
+ dst_wgt[di] += src_wgt[si];
dst_mask[di] = ((src_mask[si] && dst_mask[di]) ||
((refinement_factor != 1) && (dst_mask[di])));
// So if they are on the same level, then take the logical and
// otherwise, set it to the destination mask
- src_x[si] = -1;
+ src_x[si] = -2;
for (i = 0; i < NumArrays; i++) {
dst_vals[i][di] += src_vals[i][si];
}
More information about the yt-svn
mailing list