[Yt-svn] yt-commit r1572 - trunk/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Mon Jan 11 21:58:04 PST 2010
Author: mturk
Date: Mon Jan 11 21:58:03 2010
New Revision: 1572
URL: http://yt.enzotools.org/changeset/1572
Log:
Updating the new ParticleIO method to be the default for data objects.
This new method utilizes a C function that counts the constituent particles,
allocates a single array, and returns that. Avoiding multiple concatenations
speeds things up considerably and reduces memory fragmentation.
However, a validation method has to be written for every data type. Right now
only regions and spheres work. (Spheres are currently non-periodic.)
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/DataReadingFuncs.py
trunk/yt/lagos/HDF5LightReader.c
trunk/yt/lagos/ParticleIO.py
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Mon Jan 11 21:58:03 2010
@@ -1609,7 +1609,7 @@
dx, grid["GridIndices"][pointI].ravel()], 'float64').swapaxes(0,1)
return tr
- def get_data(self, fields=None, in_grids=False):
+ def get_data(self, fields=None, in_grids=False, force_particle_read = False):
if self._grids == None:
self._get_list_of_grids()
points = []
@@ -1625,6 +1625,14 @@
if field not in self.hierarchy.field_list and not in_grids:
if self._generate_field(field):
continue # True means we already assigned it
+ # There are a lot of 'ands' here, but I think they are all
+ # necessary.
+ if force_particle_read == False and \
+ self.pf.field_info.has_key(field) and \
+ self.pf.field_info[field].particle_type and \
+ self.pf.h.io._particle_reader:
+ self[field] = self.particles[field]
+ continue
self[field] = na.concatenate(
[self._get_data_from_grid(grid, field)
for grid in self._grids])
Modified: trunk/yt/lagos/DataReadingFuncs.py
==============================================================================
--- trunk/yt/lagos/DataReadingFuncs.py (original)
+++ trunk/yt/lagos/DataReadingFuncs.py Mon Jan 11 21:58:03 2010
@@ -32,7 +32,7 @@
class BaseIOHandler(object):
_data_style = None
- _particle_reader = True
+ _particle_reader = False
class __metaclass__(type):
def __init__(cls, name, b, d):
@@ -56,13 +56,6 @@
# We only read the one set and do not store it if it isn't pre-loaded
return self._read_data_set(grid, field)
- def _read_particles(self, fields, rtype, args, grid_list, enclosed,
- conv_factors):
- filenames = [g.filename for g in grid_list]
- ids = [g.id for g in grid_list]
- return HDF5LightReader.ReadParticles(
- rtype, fields, filenames, ids, conv_factors, args, 0)
-
def peek(self, grid, field):
return self.queue[grid.id].get(field, None)
@@ -147,6 +140,7 @@
class IOHandlerHDF5(BaseIOHandler):
_data_style = "enzo_hdf5"
+ _particle_reader = True
def _read_field_names(self, grid):
"""
@@ -181,6 +175,13 @@
def _read_exception(self):
return (exceptions.KeyError, HDF5LightReader.ReadingError)
+ def _read_particles(self, fields, rtype, args, grid_list, enclosed,
+ conv_factors):
+ filenames = [g.filename for g in grid_list]
+ ids = [g.id for g in grid_list]
+ return HDF5LightReader.ReadParticles(
+ rtype, fields, filenames, ids, conv_factors, args, 0)
+
class IOHandlerExtracted(BaseIOHandler):
_data_style = 'extracted'
@@ -390,6 +391,7 @@
class IOHandlerPacked2D(IOHandlerPackedHDF5):
_data_style = "enzo_packed_2d"
+ _particle_reader = False
def _read_data_set(self, grid, field):
return HDF5LightReader.ReadData(grid.filename,
@@ -407,6 +409,7 @@
class IOHandlerPacked1D(IOHandlerPackedHDF5):
_data_style = "enzo_packed_1d"
+ _particle_reader = False
def _read_data_set(self, grid, field):
return HDF5LightReader.ReadData(grid.filename,
Modified: trunk/yt/lagos/HDF5LightReader.c
==============================================================================
--- trunk/yt/lagos/HDF5LightReader.c (original)
+++ trunk/yt/lagos/HDF5LightReader.c Mon Jan 11 21:58:03 2010
@@ -51,6 +51,7 @@
char **field_names;
PyArrayObject *conv_factors;
PyArrayObject **return_values;
+ int *npy_types;
int (*count_func)(struct particle_validation_ *data);
int (*count_func_float)(struct particle_validation_ *data);
int (*count_func_double)(struct particle_validation_ *data);
@@ -67,8 +68,15 @@
int periodic;
} region_validation;
+typedef struct sphere_validation_ {
+ /* These cannot contain any pointers */
+ npy_float64 center[3];
+ npy_float64 radius;
+} sphere_validation;
+
/* Forward declarations */
int setup_validator_region(particle_validation *data, PyObject *InputData);
+int setup_validator_sphere(particle_validation *data, PyObject *InputData);
int run_validators(particle_validation *pv, char *filename,
int grid_id, const int read, const int packed,
int grid_index);
@@ -77,6 +85,11 @@
int count_particles_region_FLOAT(particle_validation *data);
int count_particles_region_DOUBLE(particle_validation *data);
+int count_particles_region_LONGDOUBLE(particle_validation *data);
+
+int count_particles_sphere_FLOAT(particle_validation *data);
+int count_particles_sphere_DOUBLE(particle_validation *data);
+int count_particles_sphere_LONGDOUBLE(particle_validation *data);
int get_my_desc_type(hid_t native_type_id){
@@ -574,9 +587,6 @@
/* Similar to the way Enzo does it, we're going to set the file access
property to store bigger bits in RAM. */
- int memory_increment = 1024*1024; /* in bytes */
- int dump_flag = 0;
-
file_id = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0) {
@@ -735,6 +745,7 @@
pv.validation_reqs = NULL;
pv.particle_position[0] = pv.particle_position[1] = pv.particle_position[2] = NULL;
pv.return_values = NULL;
+ pv.npy_types = NULL;
/* Set initial values for pv */
pv.stride_size = stride_size;
@@ -798,6 +809,10 @@
/* Region type */
setup_validator_region(&pv, vargs);
break;
+ case 1:
+ /* Sphere type */
+ setup_validator_sphere(&pv, vargs);
+ break;
default:
PyErr_Format(_hdf5ReadError,
"Unrecognized data source.\n");
@@ -827,10 +842,12 @@
pv.return_values = (PyArrayObject**) malloc(
sizeof(PyArrayObject*) * nfields);
+ pv.npy_types = (int *) malloc(sizeof(int) * nfields);
pv.field_names = (char **) malloc(
sizeof(char *) * nfields);
for (ifield = 0; ifield < nfields; ifield++) {
pv.return_values[ifield] = NULL;
+ pv.npy_types[ifield] = -999;
pv.field_names[ifield] = PyString_AsString(PyList_GetItem(field_list, ifield));
}
@@ -862,6 +879,7 @@
free(pv.mask);
free(pv.field_names);
free(pv.return_values); /* Has to happen after packing our return value */
+ free(pv.npy_types);
for (i = 0; i<3; i++) {
free(pv.particle_position[i]);
}
@@ -887,6 +905,7 @@
}
free(pv.return_values);
}
+ if(pv.npy_types != NULL) free(pv.npy_types);
for(i = 0; i < 3; i++) {
if(pv.particle_position[i] != NULL) free(pv.particle_position[i]);
}
@@ -929,12 +948,39 @@
data->count_func = NULL;
data->count_func_float = count_particles_region_FLOAT;
data->count_func_double = count_particles_region_DOUBLE;
+ data->count_func_longdouble = count_particles_region_LONGDOUBLE;
/* We need to insert more periodic logic here */
return 1;
}
+int setup_validator_sphere(particle_validation *data, PyObject *InputData)
+{
+ int i;
+ /* These are borrowed references */
+ PyArrayObject *center = (PyArrayObject *) PyTuple_GetItem(InputData, 0);
+ PyObject *radius = (PyObject *) PyTuple_GetItem(InputData, 1);
+
+ /* This will get freed in the finalization of particle validation */
+ sphere_validation *sv = (sphere_validation *)
+ malloc(sizeof(sphere_validation));
+ data->validation_reqs = (void *) sv;
+
+ for (i = 0; i < 3; i++){
+ sv->center[i] = *(npy_float64*) PyArray_GETPTR1(center, i);
+ }
+
+ sv->radius = (npy_float64) PyFloat_AsDouble(radius);
+
+ data->count_func = NULL;
+ data->count_func_float = count_particles_sphere_FLOAT;
+ data->count_func_double = count_particles_sphere_DOUBLE;
+ data->count_func_longdouble = count_particles_sphere_LONGDOUBLE;
+
+ return 1;
+}
+
int run_validators(particle_validation *pv, char *filename,
int grid_id, const int read, const int packed,
int grid_index)
@@ -965,8 +1011,6 @@
dataspace = memspace = datatype_id = native_type_id = -1;
rdatatype_id = rnative_type_id = -1;
- int npy_type;
-
if (packed == 1) {
snprintf(name_x, 254, "/Grid%08d/particle_position_x", grid_id);
snprintf(name_y, 254, "/Grid%08d/particle_position_y", grid_id);
@@ -1045,11 +1089,11 @@
npy_intp dims = pv->total_valid_particles;
rdatatype_id = H5Dget_type(dataset_read[i]);
rnative_type_id = H5Tget_native_type(rdatatype_id, H5T_DIR_ASCEND);
- npy_type = get_my_desc_type(rnative_type_id);
+ pv->npy_types[i] = get_my_desc_type(rnative_type_id);
//fprintf(stderr, "Allocating array of size %d\n", (int) dims);
pv->return_values[i] = (PyArrayObject *)
PyArray_SimpleNewFromDescr(
- 1, &dims, PyArray_DescrFromType(npy_type));
+ 1, &dims, PyArray_DescrFromType(pv->npy_types[i]));
H5Tclose(rnative_type_id);
H5Tclose(rdatatype_id);
}
@@ -1103,14 +1147,26 @@
/* Now we multiply our fields by the appropriate conversion factor */
if (cfactors[ifield] != 1.0) {
for(p_ind = 0; p_ind < read_here; p_ind++)
- if (npy_type == 11) { // floats
+ if (pv->npy_types[ifield] == NPY_FLOAT) { // floats
*(npy_float32 *) PyArray_GETPTR1(
pv->return_values[ifield], p_ind + pv->nread)
*= cfactors[ifield];
- } else { // doubles
+ } else if (pv->npy_types[ifield] == NPY_DOUBLE) { // doubles
*(npy_float64 *) PyArray_GETPTR1(
pv->return_values[ifield], p_ind + pv->nread)
*= cfactors[ifield];
+ } else if (pv->npy_types[ifield] == NPY_LONGDOUBLE) {
+ *(npy_float128 *) PyArray_GETPTR1 (
+ pv->return_values[ifield], p_ind + pv->nread)
+ *= cfactors[ifield];
+ } else if (pv->npy_types[ifield] == NPY_INT) {
+ *(npy_int *) PyArray_GETPTR1 (
+ pv->return_values[ifield], p_ind + pv->nread)
+ *= cfactors[ifield];
+ } else if (pv->npy_types[ifield] == NPY_LONG) {
+ *(npy_long *) PyArray_GETPTR1 (
+ pv->return_values[ifield], p_ind + pv->nread)
+ *= cfactors[ifield];
}
}
}
@@ -1170,7 +1226,7 @@
/* First is our validation requirements, which are a set of three items: */
- int ind, n;
+ int ind, n=0;
region_validation *vdata;
vdata = (region_validation*) data->validation_reqs;
@@ -1304,6 +1360,188 @@
return n;
}
+int count_particles_region_LONGDOUBLE(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n=0;
+ region_validation *vdata;
+
+ vdata = (region_validation*) data->validation_reqs;
+
+ long double **particle_data = (long double **) data->particle_position;
+
+ long double *particle_position_x = particle_data[0];
+ long double *particle_position_y = particle_data[1];
+ long double *particle_position_z = particle_data[2];
+ long double tempx, tempy, tempz;
+
+ if (vdata->periodic == 0) {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ if ( (particle_position_x[ind] >= vdata->left_edge[0])
+ && (particle_position_x[ind] <= vdata->right_edge[0])
+ && (particle_position_y[ind] >= vdata->left_edge[1])
+ && (particle_position_y[ind] <= vdata->right_edge[1])
+ && (particle_position_z[ind] >= vdata->left_edge[2])
+ && (particle_position_z[ind] <= vdata->right_edge[2])) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ } else {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ tempx = particle_position_x[ind];
+ tempy = particle_position_y[ind];
+ tempz = particle_position_z[ind];
+ if ( (tempx < vdata->left_edge[0]) && (tempx < vdata->right_edge[0]) ) {
+ tempx += vdata->period[0];
+ } else if ( (tempx > vdata->left_edge[0]) && (tempx > vdata->right_edge[0]) ) {
+ tempx -= vdata->period[0];
+ }
+ if ( (tempy < vdata->left_edge[1]) && (tempx < vdata->right_edge[1]) ) {
+ tempy += vdata->period[1];
+ } else if ( (tempy > vdata->left_edge[1]) && (tempx > vdata->right_edge[1]) ) {
+ tempy -= vdata->period[1];
+ }
+ if ( (tempz < vdata->left_edge[2]) && (tempx < vdata->right_edge[2]) ) {
+ tempz += vdata->period[2];
+ } else if ( (tempz > vdata->left_edge[2]) && (tempx > vdata->right_edge[2]) ) {
+ tempz -= vdata->period[2];
+ }
+ if ( (tempx >= vdata->left_edge[0])
+ && (tempx <= vdata->right_edge[0])
+ && (tempy >= vdata->left_edge[1])
+ && (tempy <= vdata->right_edge[1])
+ && (tempz >= vdata->left_edge[2])
+ && (tempz <= vdata->right_edge[2])) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ }
+ return n;
+}
+
+int count_particles_sphere_FLOAT(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n=0;
+ sphere_validation *vdata;
+
+ vdata = (sphere_validation*) data->validation_reqs;
+
+ float **particle_data = (float **) data->particle_position;
+
+ float *particle_position_x = particle_data[0];
+ float *particle_position_y = particle_data[1];
+ float *particle_position_z = particle_data[2];
+ float tempr;
+
+ double pradius;
+
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0;
+ tempr = (particle_position_x[ind] - vdata->center[0]); pradius += tempr*tempr;
+ tempr = (particle_position_y[ind] - vdata->center[1]); pradius += tempr*tempr;
+ tempr = (particle_position_z[ind] - vdata->center[2]); pradius += tempr*tempr;
+ pradius = pow(pradius, 0.5);
+ if (pradius <= vdata->radius) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ return n;
+}
+
+int count_particles_sphere_DOUBLE(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n=0;
+ sphere_validation *vdata;
+
+ vdata = (sphere_validation*) data->validation_reqs;
+
+ double **particle_data = (double **) data->particle_position;
+
+ double *particle_position_x = particle_data[0];
+ double *particle_position_y = particle_data[1];
+ double *particle_position_z = particle_data[2];
+ double tempr;
+
+ double pradius;
+
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0;
+ tempr = (particle_position_x[ind] - vdata->center[0]); pradius += tempr*tempr;
+ tempr = (particle_position_y[ind] - vdata->center[1]); pradius += tempr*tempr;
+ tempr = (particle_position_z[ind] - vdata->center[2]); pradius += tempr*tempr;
+ pradius = pow(pradius, 0.5);
+ if (pradius <= vdata->radius) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ return n;
+}
+
+int count_particles_sphere_LONGDOUBLE(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n=0;
+ sphere_validation *vdata;
+
+ vdata = (sphere_validation*) data->validation_reqs;
+
+ long double **particle_data = (long double **) data->particle_position;
+
+ long double *particle_position_x = particle_data[0];
+ long double *particle_position_y = particle_data[1];
+ long double *particle_position_z = particle_data[2];
+ long double tempr;
+
+ long double pradius;
+
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0;
+ tempr = (particle_position_x[ind] - vdata->center[0]); pradius += tempr*tempr;
+ tempr = (particle_position_y[ind] - vdata->center[1]); pradius += tempr*tempr;
+ tempr = (particle_position_z[ind] - vdata->center[2]); pradius += tempr*tempr;
+ pradius = pow(pradius, 0.5);
+ if (pradius <= vdata->radius) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ return n;
+}
+
+
static PyMethodDef _hdf5LightReaderMethods[] = {
{"ReadData", Py_ReadHDF5DataSet, METH_VARARGS},
{"ReadDataSlice", Py_ReadHDF5DataSetSlice, METH_VARARGS},
Modified: trunk/yt/lagos/ParticleIO.py
==============================================================================
--- trunk/yt/lagos/ParticleIO.py (original)
+++ trunk/yt/lagos/ParticleIO.py Mon Jan 11 21:58:03 2010
@@ -36,6 +36,7 @@
particle_handler_registry[cls._source_type] = cls
_source_type = None
+
def __init__(self, pf, source):
self.pf = pf
self.data = {}
@@ -57,31 +58,20 @@
def get_data(self, fields):
fields = ensure_list(fields)
- self.source.get_data(fields)
+ self.source.get_data(fields, force_particle_read=True)
for field in fields:
self[field] = self.source[field]
particle_handler_registry.default_factory = lambda: ParticleIOHandler
-class ParticleIOHandlerRegion(ParticleIOHandler):
- periodic = False
- _source_type = "region"
- def __init__(self, pf, source):
- self.left_edge = source.left_edge
- self.right_edge = source.right_edge
- ParticleIOHandler.__init__(self, pf, source)
-
+class ParticleIOHandlerImplemented(ParticleIOHandler):
def get_data(self, fields):
mylog.info("Getting %s using ParticleIO" % str(fields))
fields = ensure_list(fields)
if not self.pf.h.io._particle_reader:
mylog.info("not self.pf.h.io._particle_reader")
return self.source.get_data(fields)
- rtype = 0
- DLE = na.array(self.pf["DomainLeftEdge"], dtype='float64')
- DRE = na.array(self.pf["DomainRightEdge"], dtype='float64')
- args = (na.array(self.left_edge), na.array(self.right_edge),
- int(self.periodic), DLE, DRE)
+ rtype, args = self._get_args()
count_list, grid_list = [], []
for grid in self.source._grids:
if grid.NumberOfParticles == 0: continue
@@ -113,6 +103,22 @@
conv_factors)
for field, v in zip(fields, rv): self[field] = v
+class ParticleIOHandlerRegion(ParticleIOHandlerImplemented):
+ periodic = False
+ _source_type = "region"
+
+ def __init__(self, pf, source):
+ self.left_edge = source.left_edge
+ self.right_edge = source.right_edge
+ ParticleIOHandler.__init__(self, pf, source)
+
+ def _get_args(self):
+ DLE = na.array(self.pf["DomainLeftEdge"], dtype='float64')
+ DRE = na.array(self.pf["DomainRightEdge"], dtype='float64')
+ args = (na.array(self.left_edge), na.array(self.right_edge),
+ int(self.periodic), DLE, DRE)
+ return (0, args)
+
class ParticleIOHandlerRegionStrict(ParticleIOHandlerRegion):
_source_type = "region_strict"
@@ -122,3 +128,14 @@
class ParticleIOHandlerPeriodicRegionStrict(ParticleIOHandlerPeriodicRegion):
_source_type = "periodic_region_strict"
+
+class ParticleIOHandlerSphere(ParticleIOHandlerImplemented):
+ _source_type = "sphere"
+
+ def __init__(self, pf, source):
+ self.center = source.center
+ self.radius = source.radius
+ ParticleIOHandler.__init__(self, pf, source)
+
+ def _get_args(self):
+ return (1, (na.array(self.center, dtype='float64'), self.radius))
More information about the yt-svn
mailing list