[Yt-svn] yt-commit r1263 - in trunk/yt: lagos lagos/hop raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Mon Apr 20 10:49:11 PDT 2009
Author: mturk
Date: Mon Apr 20 10:49:10 2009
New Revision: 1263
URL: http://yt.spacepope.org/changeset/1263
Log:
A bunch of changes:
* find_max now works in parallel, but not super efficiently
* IOError from nonexistent filenames now returns the name
* Slices now work in parallel but do not accept source= anymore.
* Cutting planes now work in parallel
* Got rid of cruft code for old HaloFinder
* Improved error handling in _MPL.c for CPixelize
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/DerivedQuantities.py
trunk/yt/lagos/HierarchyType.py
trunk/yt/lagos/OutputTypes.py
trunk/yt/lagos/hop/HopOutput.py
trunk/yt/raven/_MPL.c
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Mon Apr 20 10:49:10 2009
@@ -206,7 +206,7 @@
def __delitem__(self, key):
"""
- Sets a field to be some other value.
+ Deletes a field
"""
try:
del self.fields[self.fields.index(key)]
@@ -629,7 +629,7 @@
_con_args = ('axis', 'coord')
#@time_execution
def __init__(self, axis, coord, fields = None, center=None, pf=None,
- node_name = False, source = None, **kwargs):
+ node_name = False, **kwargs):
"""
Slice along *axis*:ref:`axis-specification`, at the coordinate *coord*.
Optionally supply fields.
@@ -637,21 +637,12 @@
AMR2DData.__init__(self, axis, fields, pf, **kwargs)
self.center = center
self.coord = coord
- self._initialize_source(source)
if node_name is False:
self._refresh_data()
else:
if node_name is True: self._deserialize()
else: self._deserialize(node_name)
- def _initialize_source(self, source = None):
- if source is None:
- check, source = self._partition_hierarchy_2d(self.axis)
- self._check_region = check
- else:
- self._check_region = True
- self.source = source
-
def reslice(self, coord):
"""
Change the entire dataset, clearing out the current data and slicing at
@@ -699,17 +690,12 @@
self.ActiveDimensions = (t.shape[0], 1, 1)
def _get_list_of_grids(self):
- goodI = ((self.source.gridRightEdge[:,self.axis] > self.coord)
- & (self.source.gridLeftEdge[:,self.axis] <= self.coord ))
- self._grids = self.source._grids[goodI] # Using sources not hierarchy
+ goodI = ((self.pf.h.gridRightEdge[:,self.axis] > self.coord)
+ & (self.pf.h.gridLeftEdge[:,self.axis] <= self.coord ))
+ self._grids = self.pf.h.grids[goodI] # Using sources not hierarchy
def __cut_mask_child_mask(self, grid):
mask = grid.child_mask.copy()
- if self._check_region:
- cut_mask = self.source._get_cut_mask(grid)
- if mask is False: mask *= False
- elif mask is True: pass
- else: mask &= cut_mask
return mask
def _generate_grid_coords(self, grid):
@@ -1704,8 +1690,6 @@
self.right_edge)
def _is_fully_enclosed(self, grid):
- offsets = na.array([-1,0,1])
-
for off_x, off_y, off_z in self.offsets:
region_left = [self.left_edge[0]+off_x,
self.left_edge[1]+off_y,self.left_edge[2]+off_z]
Modified: trunk/yt/lagos/DerivedQuantities.py
==============================================================================
--- trunk/yt/lagos/DerivedQuantities.py (original)
+++ trunk/yt/lagos/DerivedQuantities.py Mon Apr 20 10:49:10 2009
@@ -382,7 +382,26 @@
return [(na.min(mins[:,i]), na.max(maxs[:,i])) for i in range(n_fields)]
add_quantity("Extrema", function=_Extrema, combine_function=_combExtrema,
n_ret=3)
-
+
+def _MaxLocation(data, field):
+ """
+ This function returns the location of the maximum of a set
+ of fields.
+ """
+ ma, maxi, mx, my, mz, mg = -1e90, -1, -1, -1, -1, -1
+ if data[field].size > 0:
+ maxi = na.argmax(data[field])
+ ma = data[field][maxi]
+ mx, my, mz = [data[ax][maxi] for ax in 'xyz']
+ mg = data["GridIndices"][maxi]
+ return (ma, maxi, mx, my, mz, mg)
+def _combMaxLocation(data, *args):
+ args = [na.atleast_1d(arg) for arg in args]
+ i = na.argmax(args[0]) # ma is arg[0]
+ return [arg[i] for arg in args]
+add_quantity("MaxLocation", function=_MaxLocation,
+ combine_function=_combMaxLocation, n_ret = 6)
+
def _TotalQuantity(data, fields):
"""
This function sums up a given field over the entire region
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Mon Apr 20 10:49:10 2009
@@ -305,25 +305,22 @@
def find_max_cell_location(self, field, finestLevels = True):
if finestLevels:
- gI = na.where(self.gridLevels >= self.maxLevel - NUMTOCHECK)
+ gi = (self.gridLevels >= self.max_level - NUMTOCHECK).ravel()
+ source = self.grid_collection([0.0]*3,
+ self.grids[gi])
else:
- gI = na.where(self.gridLevels >= 0) # Slow but pedantic
- maxVal = -1e100
- for grid in self.grids[gI[0]]:
- mylog.debug("Checking %s (level %s)", grid.id, grid.Level)
- val, coord = grid.find_max(field)
- if val > maxVal:
- maxCoord = coord
- maxVal = val
- maxGrid = grid
- mc = na.array(maxCoord)
- pos=maxGrid.get_position(mc)
+ source = self.all_data()
+ mylog.debug("Searching %s grids for maximum value of %s",
+ len(source._grids), field)
+ max_val, maxi, mx, my, mz, mg = source.quantities["MaxLocation"](
+ field, lazy_reader=True)
+ max_grid = self.grids[mg]
+ mc = na.unravel_index(maxi, max_grid.ActiveDimensions)
mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s %s", \
- maxVal, pos[0], pos[1], pos[2], maxGrid, maxGrid.Level, mc)
- self.center = pos
- self.parameters["Max%sValue" % (field)] = maxVal
- self.parameters["Max%sPos" % (field)] = "%s" % (pos)
- return maxGrid, maxCoord, maxVal, pos
+ max_val, mx, my, mz, max_grid, max_grid.Level, mc)
+ self.parameters["Max%sValue" % (field)] = max_val
+ self.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
+ return max_grid, mc, max_val, (mx,my,mz)
@time_execution
def find_min(self, field):
Modified: trunk/yt/lagos/OutputTypes.py
==============================================================================
--- trunk/yt/lagos/OutputTypes.py (original)
+++ trunk/yt/lagos/OutputTypes.py Mon Apr 20 10:49:10 2009
@@ -45,7 +45,7 @@
def __new__(cls, filename, *args, **kwargs):
apath = os.path.abspath(filename)
- if not os.path.exists(apath): raise IOError
+ if not os.path.exists(apath): raise IOError(filename)
if apath not in _cached_pfs:
obj = object.__new__(cls)
obj.__init__(filename, *args, **kwargs)
Modified: trunk/yt/lagos/hop/HopOutput.py
==============================================================================
--- trunk/yt/lagos/hop/HopOutput.py (original)
+++ trunk/yt/lagos/hop/HopOutput.py Mon Apr 20 10:49:10 2009
@@ -231,44 +231,3 @@
def get_size(self):
return self.indices.size
-
-class HaloFinder(ParallelAnalysisInterface):
- def __init__(self, pf, threshold=160.0, dm_only=True):
- self.pf = pf
- self.hierarchy = pf.hierarchy
- self.padding = 0.2 * pf["unitary"]
- LE, RE, self.source = self._partition_hierarchy_3d(padding=self.padding)
- self.bounds = (LE, RE)
- self._reposition_particles((LE, RE))
- # For scaling the threshold, note that it's a passthrough
- total_mass = self._mpi_allsum(self.source["ParticleMassMsun"].sum())
- hop_list = HopList(self.source, threshold, dm_only)
- self._join_hoplists(hop_list)
-
- @parallel_passthrough
- def _join_hoplists(self, hop_list):
- # First we get the total number of halos the entire collection
- # has identified
- nhalos = self._mpi_allsum(len(hop_list))
- # Now we identify our padding-region particles
- LE, RE = self.bounds
- ind = na.zeros(hop_list.particle_fields["particle_position_x"].size, dtype='bool')
- for i, ax in enumerate('xyz'):
- arr = hop_list.particle_fields["particle_position_%s" % ax]
- ind |= (arr < LE[i]-self.padding)
- ind |= (arr > RE[i]+self.padding)
- # This is a one-d array of all buffer particle indices
- indices = self._mpi_catarray(hop_list.particle_fields["particle_index"][ind])
- # This is a one-d array of halo IDs
- halos = self._mpi_catarray(hop_list.tags[ind])
-
- @parallel_passthrough
- def _reposition_particles(self, bounds):
- # This only does periodicity. We do NOT want to deal with anything
- # else. The only reason we even do periodicity is the
- LE, RE = bounds
- dw = self.pf["DomainRightEdge"] - self.pf["DomainLeftEdge"]
- for i, ax in enumerate('xyz'):
- arr = self.source["particle_position_%s" % ax]
- arr[arr < LE[i]-self.padding] += dw[i]
- arr[arr > RE[i]+self.padding] -= dw[i]
Modified: trunk/yt/raven/_MPL.c
==============================================================================
--- trunk/yt/raven/_MPL.c (original)
+++ trunk/yt/raven/_MPL.c Mon Apr 20 10:49:10 2009
@@ -193,6 +193,16 @@
PyObject *xp, *yp, *zp, *pxp, *pyp,
*dxp, *dyp, *dzp, *dp,
*centerp, *inv_matp, *indicesp;
+
+ xp = yp = zp = pxp = pyp = dxp = dyp = dzp = dp = NULL;
+ centerp = inv_matp = indicesp = NULL;
+
+ PyArrayObject *x, *y, *z, *px, *py, *d,
+ *dx, *dy, *dz, *center, *inv_mat, *indices;
+
+ x = y = z = px = py = dx = dy = dz = d = NULL;
+ center = inv_mat = indices = NULL;
+
unsigned int rows, cols;
double px_min, px_max, py_min, py_max;
@@ -211,79 +221,79 @@
PyErr_Format( _pixelizeError, "Cannot scale to zero size.");
// Get numeric arrays
- PyArrayObject *x = (PyArrayObject *) PyArray_FromAny(xp,
+ x = (PyArrayObject *) PyArray_FromAny(xp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if (x == NULL) {
PyErr_Format( _pixelizeError, "x is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *y = (PyArrayObject *) PyArray_FromAny(yp,
+ y = (PyArrayObject *) PyArray_FromAny(yp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((y == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "y is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *z = (PyArrayObject *) PyArray_FromAny(zp,
+ z = (PyArrayObject *) PyArray_FromAny(zp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((z == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "z is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *px = (PyArrayObject *) PyArray_FromAny(pxp,
+ px = (PyArrayObject *) PyArray_FromAny(pxp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((px == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "px is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *py = (PyArrayObject *) PyArray_FromAny(pyp,
+ py = (PyArrayObject *) PyArray_FromAny(pyp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((py == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "py is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *d = (PyArrayObject *) PyArray_FromAny(dp,
+ d = (PyArrayObject *) PyArray_FromAny(dp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((d == NULL) || (PyArray_SIZE(d) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "data is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *dx = (PyArrayObject *) PyArray_FromAny(dxp,
+ dx = (PyArrayObject *) PyArray_FromAny(dxp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((dx == NULL) || (PyArray_SIZE(dx) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "dx is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *dy = (PyArrayObject *) PyArray_FromAny(dyp,
+ dy = (PyArrayObject *) PyArray_FromAny(dyp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((dy == NULL) || (PyArray_SIZE(dy) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "dy is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *dz = (PyArrayObject *) PyArray_FromAny(dzp,
+ dz = (PyArrayObject *) PyArray_FromAny(dzp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((dz == NULL) || (PyArray_SIZE(dz) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "dz is of incorrect type (wanted 1D float)");
goto _fail;
}
- PyArrayObject *center = (PyArrayObject *) PyArray_FromAny(centerp,
+ center = (PyArrayObject *) PyArray_FromAny(centerp,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
if ((dz == NULL) || (PyArray_SIZE(center) != 3)) {
PyErr_Format( _pixelizeError, "Center must have three points");
goto _fail;
}
- PyArrayObject *inv_mat = (PyArrayObject *) PyArray_FromAny(inv_matp,
+ inv_mat = (PyArrayObject *) PyArray_FromAny(inv_matp,
PyArray_DescrFromType(NPY_FLOAT64), 2, 2, 0, NULL);
if ((inv_mat == NULL) || (PyArray_SIZE(inv_mat) != 9)) {
PyErr_Format( _pixelizeError, "inv_mat must be three by three");
goto _fail;
}
- PyArrayObject *indices = (PyArrayObject *) PyArray_FromAny(indicesp,
+ indices = (PyArrayObject *) PyArray_FromAny(indicesp,
PyArray_DescrFromType(NPY_INT64), 1, 1, 0, NULL);
if ((indices == NULL) || (PyArray_SIZE(indices) != PyArray_SIZE(dx))) {
PyErr_Format( _pixelizeError, "indices must be same length as dx");
More information about the yt-svn
mailing list