[Yt-svn] commit/yt: 9 new changesets
Bitbucket
commits-noreply at bitbucket.org
Sun Mar 6 08:25:29 PST 2011
9 new changesets in yt:
http://bitbucket.org/yt_analysis/yt/changeset/19e21a583e6c/
changeset: r3794:19e21a583e6c
branch: yt
user: sskory
date: 2011-03-05 01:07:11
summary: Beginning work on reading in halos off disk. Not functional yet.
affected #: 1 file (7.7 KB)
--- a/yt/analysis_modules/halo_finding/halo_objects.py Fri Mar 04 11:58:40 2011 -0700
+++ b/yt/analysis_modules/halo_finding/halo_objects.py Fri Mar 04 16:07:11 2011 -0800
@@ -32,6 +32,7 @@
import numpy as na
import random
import sys
+from collections import defaultdict
from yt.funcs import *
@@ -786,6 +787,158 @@
r"""Not implemented."""
return self.center_of_mass()
+class LoadedHalo(Halo):
+ def __init__(self, id, size=None, CoM=None,
+ max_dens_point=None, group_total_mass=None, max_radius=None, bulk_vel=None,
+ rms_vel=None, fnames=None):
+ self._max_dens = halo_list._max_dens
+ self.id = id
+ self.size = size
+ self.CoM = CoM
+ self.max_dens_point = max_dens_point
+ self.group_total_mass = group_total_mass
+ self.max_radius = max_radius
+ self.bulk_vel = bulk_vel
+ self.rms_vel = rms_vel
+ # locs=the names of the h5 files that have particle data for this halo
+ self.fnames = fnames
+ self.bin_count = None
+ self.overdensity = None
+ self.saved_fields = defaultdict(None)
+
+ def __getitem__(self, key):
+ if self.saved_fields[key] is not None:
+ # We've already got it.
+ return self.saved_fields[key]
+ else:
+ # Gotta go get it.
+ field_data = self._get_particle_data(self.id, self.locs,
+ self.size)
+ if field_data is not None:
+ self.saved_fields[key] = field_data
+ return self.saved_fields[key]
+ else:
+ # Dynamically create the masking array for particles, and get
+ # the data using standard yt methods.
+
+ def _get_particle_data(self, halo, fnames, size, field):
+ # Given a list of file names, a halo, its size, and the desired field,
+ # this returns the particle data for that halo.
+ # First get the list of fields from the first file. Not all fields
+ # are saved all the time (e.g. creation_time, particle_type).
+ f = h5py.File(fnames[0])
+ fields = f["Halo%08d" % halo].keys()
+ # If we dont have this field, we can give up right now.
+ if field not in fields: return None
+ if field == 'particle_index' or field == 'particle_type':
+ # the only integer field
+ field_data = na.empty(size, dtype='int64')
+ else:
+ field_data = na.empty(size, dtype='float64')
+ f.close()
+ offset = 0
+ for fname in fnames:
+ f = h5py.File(fname)
+ this = f["Halo%08d" % halo][field][:]
+ s = this.size
+ field_data[offset:offset+s] = this
+ offset += s
+ f.close()
+ return field_data
+
+ def center_of_mass(self):
+ r"""Calculate and return the center of mass.
+
+ The center of mass of the halo is directly calculated and returned.
+
+ Examples
+ --------
+ >>> com = halos[0].center_of_mass()
+ """
+ return self.CoM
+
+ def maximum_density_location(self):
+ r"""Return the location HOP identified as maximally dense.
+
+ Return the location HOP identified as maximally dense.
+
+ Examples
+ --------
+ >>> max_dens_loc = halos[0].maximum_density_location()
+ """
+ return self.max_dens_point[1:]
+
+ def maximum_density(self):
+ r"""Return the HOP-identified maximum density.
+
+ Return the HOP-identified maximum density.
+
+ Examples
+ --------
+ >>> max_dens = halos[0].maximum_density()
+ """
+ return self.max_dens_point[0]
+
+ def total_mass(self):
+ r"""Returns the total mass in solar masses of the halo.
+
+ Returns the total mass in solar masses of just the particles in the
+ halo.
+
+ Examples
+ --------
+ >>> halos[0].total_mass()
+ """
+ return self.group_total_mass
+
+ def bulk_velocity(self):
+ r"""Returns the mass-weighted average velocity in cm/s.
+
+ This calculates and returns the mass-weighted average velocity of just
+ the particles in the halo in cm/s.
+
+ Examples
+ --------
+ >>> bv = halos[0].bulk_velocity()
+ """
+ return self.bulk_vel
+
+ def rms_velocity(self):
+ r"""Returns the mass-weighted RMS velocity for the halo
+ particles in cgs units.
+
+ Calculate and return the mass-weighted RMS velocity for just the
+ particles in the halo. The bulk velocity of the halo is subtracted
+ before computation.
+
+ Examples
+ --------
+ >>> rms_vel = halos[0].rms_velocity()
+ """
+ return self.rms_vel
+
+ def maximum_radius(self, center_of_mass=True):
+ r"""Returns the maximum radius in the halo for all particles,
+ either from the point of maximum density or from the
+ center of mass.
+
+ The maximum radius from the most dense point is calculated. This
+ accounts for periodicity.
+
+ Parameters
+ ----------
+ center_of_mass : bool
+ True chooses the center of mass when calculating the maximum radius.
+ False chooses from the maximum density location for HOP halos
+ (it has no effect for FOF halos).
+ Default = True.
+
+ Examples
+ --------
+ >>> radius = halos[0].maximum_radius()
+ """
+ return self.max_radius
+
class HaloList(object):
_fields = ["particle_position_%s" % ax for ax in 'xyz']
@@ -1107,6 +1260,36 @@
"""
HaloList.write_out(self, filename)
+class LoadedHaloList(HaloList):
+ _name = "Loaded"
+
+ def __init__(self, basename):
+ self._groups = []
+ self.basename = basename
+
+ def _get_halo_data(self):
+ # First get the halo particulars.
+ lines = file("%s.out" % self.basename)
+ # The location of particle data for each halo.
+ locations = self._collect_halo_data_locations()
+ for halo, line in enumerate(lines[1:]): # Skip the comment line a top.
+ line = line.split()
+ # get the particle data
+ size = int(line[2])
+ locs = self._collect_halo_data_locations()
+ field_data = self._get_particle_data(halo, locs, size)
+
+
+ def _collect_halo_data_locations(self):
+ # The halos are listed in order in the file.
+ lines = file("%s.txt" % self.basename)
+ locations = []
+ for line in lines:
+ line = line.split()
+ locations.append(line[1:])
+ lines.close()
+ return locations
+
class parallelHOPHaloList(HaloList,ParallelAnalysisInterface):
_name = "parallelHOP"
_halo_class = parallelHOPHalo
@@ -1482,6 +1665,32 @@
if not self._is_mine(halo): continue
halo.write_particle_list(f)
+ def dump(self, basename="HopAnalysis"):
+ r"""Save the full halo data to disk.
+
+ This function will save the halo data in such a manner that it can be
+ easily re-loaded later using `GenericHaloFinder.load`.
+ This is similar in concept to
+ pickling the data, but outputs the data in the already-established
+ data formats. The simple halo data is written to a text file
+ (e.g. "HopAnalysis.out") using
+ write_out(), and the particle data to hdf5 files (e.g. "HopAnalysis.h5")
+ using write_particle_lists().
+
+ Parameters
+ ----------
+ basename : String
+ The base name for the files the data will be written to. Default =
+ "HopAnalysis".
+
+ Examples
+ --------
+ >>> halos.dump("MyHalos")
+ """
+ self.write_out("%s.out" % basename)
+ self.write_particle_lists(basename)
+ self.write_particle_lists_txt(basename)
+
class parallelHF(GenericHaloFinder, parallelHOPHaloList):
def __init__(self, pf, subvolume=None,threshold=160, dm_only=True, \
resize=True, rearrange=True,\
@@ -1928,3 +2137,26 @@
self._join_halolists()
HaloFinder = HOPHaloFinder
+
+class LoadHaloes(GenericHaloFinder, LoadedHaloList):
+ def __init__(self, basename):
+ r"""Load the full halo data into memory.
+
+ This function takes the output of `GenericHaloFinder.dump` and
+ re-establishes the list of halos in memory. This enables the full set
+ of halo analysis features without running the halo finder again.
+
+ Parameters
+ ----------
+ basename : String
+ The base name of the files that will be read in. This should match
+ what was used when `GenericHaloFinder.dump` was called. Default =
+ "HopAnalysis".
+
+ Examples
+ --------
+ >>> halos = LoadHaloes("HopAnalysis")
+ """
+ self.basename = basename
+
+
\ No newline at end of file
http://bitbucket.org/yt_analysis/yt/changeset/adc8939e6bbb/
changeset: r3795:adc8939e6bbb
branch: yt
user: sskory
date: 2011-03-05 18:25:58
summary: Making progress on importing halos. Partially functional!
affected #: 1 file (579 bytes)
--- a/yt/analysis_modules/halo_finding/halo_objects.py Fri Mar 04 16:07:11 2011 -0800
+++ b/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 09:25:58 2011 -0800
@@ -791,7 +791,6 @@
def __init__(self, id, size=None, CoM=None,
max_dens_point=None, group_total_mass=None, max_radius=None, bulk_vel=None,
rms_vel=None, fnames=None):
- self._max_dens = halo_list._max_dens
self.id = id
self.size = size
self.CoM = CoM
@@ -804,22 +803,23 @@
self.fnames = fnames
self.bin_count = None
self.overdensity = None
- self.saved_fields = defaultdict(None)
+ self.saved_fields = {}
def __getitem__(self, key):
- if self.saved_fields[key] is not None:
+ try:
# We've already got it.
return self.saved_fields[key]
- else:
+ except KeyError:
# Gotta go get it.
- field_data = self._get_particle_data(self.id, self.locs,
- self.size)
+ field_data = self._get_particle_data(self.id, self.fnames,
+ self.size, key)
if field_data is not None:
self.saved_fields[key] = field_data
return self.saved_fields[key]
else:
# Dynamically create the masking array for particles, and get
# the data using standard yt methods.
+ pass
def _get_particle_data(self, halo, fnames, size, field):
# Given a list of file names, a halo, its size, and the desired field,
@@ -917,7 +917,7 @@
"""
return self.rms_vel
- def maximum_radius(self, center_of_mass=True):
+ def maximum_radius(self):
r"""Returns the maximum radius in the halo for all particles,
either from the point of maximum density or from the
center of mass.
@@ -1266,19 +1266,33 @@
def __init__(self, basename):
self._groups = []
self.basename = basename
+ self._retrieve_halos()
- def _get_halo_data(self):
+ def _retrieve_halos(self):
# First get the halo particulars.
lines = file("%s.out" % self.basename)
# The location of particle data for each halo.
locations = self._collect_halo_data_locations()
- for halo, line in enumerate(lines[1:]): # Skip the comment line a top.
+ halo = 0
+ for line in lines:
+ # Skip the comment lines at top.
+ if line[0] == "#": continue
line = line.split()
# get the particle data
size = int(line[2])
- locs = self._collect_halo_data_locations()
- field_data = self._get_particle_data(halo, locs, size)
-
+ fnames = locations[halo]
+ # Everything else
+ CoM = na.array([float(line[7]), float(line[8]), float(line[9])])
+ max_dens_point = na.array([float(line[3]), float(line[4]),
+ float(line[5]), float(line[6])])
+ group_total_mass = float(line[1])
+ max_radius = float(line[13])
+ bulk_vel = na.array([float(line[10]), float(line[11]),
+ float(line[12])])
+ rms_vel = float(line[14])
+ self._groups.append(LoadedHalo(halo, size, CoM, max_dens_point,
+ group_total_mass, max_radius, bulk_vel, rms_vel, fnames))
+ halo += 1
def _collect_halo_data_locations(self):
# The halos are listed in order in the file.
@@ -2158,5 +2172,7 @@
>>> halos = LoadHaloes("HopAnalysis")
"""
self.basename = basename
-
+ LoadedHaloList.__init__(self, self.basename)
+
+
\ No newline at end of file
http://bitbucket.org/yt_analysis/yt/changeset/935b6e1eacca/
changeset: r3796:935b6e1eacca
branch: yt
user: sskory
date: 2011-03-05 19:48:17
summary: Getting non-saved particle field data not quite working yet.
affected #: 1 file (1.2 KB)
--- a/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 09:25:58 2011 -0800
+++ b/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 10:48:17 2011 -0800
@@ -788,9 +788,10 @@
return self.center_of_mass()
class LoadedHalo(Halo):
- def __init__(self, id, size=None, CoM=None,
+ def __init__(self, pf, id, size=None, CoM=None,
max_dens_point=None, group_total_mass=None, max_radius=None, bulk_vel=None,
rms_vel=None, fnames=None):
+ self.pf = pf
self.id = id
self.size = size
self.CoM = CoM
@@ -804,22 +805,42 @@
self.bin_count = None
self.overdensity = None
self.saved_fields = {}
+ self.particle_mask = None
def __getitem__(self, key):
+ # This function will try to get particle data in one of three ways,
+ # in descending preference.
+ # 1. From saved_fields, e.g. we've already got it.
+ # 2. From the halo h5 files off disk.
+ # 3. Use the unique particle indexes of the halo to select a missing
+ # field from an AMR Sphere.
try:
# We've already got it.
return self.saved_fields[key]
except KeyError:
- # Gotta go get it.
+ # Gotta go get it from the halo h5 files.
field_data = self._get_particle_data(self.id, self.fnames,
self.size, key)
+ if key == 'particle_position_x': field_data = None
if field_data is not None:
self.saved_fields[key] = field_data
return self.saved_fields[key]
else:
# Dynamically create the masking array for particles, and get
# the data using standard yt methods.
- pass
+ ds = self.pf.h.sphere(self.CoM, 1.05 * self.max_radius)
+ if self.particle_mask is None:
+ pid = self.__getitem__('particle_index')
+ sp_pid = ds['particle_index']
+ print 'testing'
+ for id in pid:
+ if id not in sp_pid: print 'trouble!', id
+ # The result of searchsorted is an array with the positions
+ # of the indexes in pid as they are in sp_pid. This is
+ # because each element of pid is in sp_pid only once.
+ self.particle_mask = na.searchsorted(sp_pid, pid)
+ print 'here'
+ return ds[key][self.particle_mask]
def _get_particle_data(self, halo, fnames, size, field):
# Given a list of file names, a halo, its size, and the desired field,
@@ -1263,7 +1284,8 @@
class LoadedHaloList(HaloList):
_name = "Loaded"
- def __init__(self, basename):
+ def __init__(self, pf, basename):
+ self.pf = pf
self._groups = []
self.basename = basename
self._retrieve_halos()
@@ -1290,8 +1312,9 @@
bulk_vel = na.array([float(line[10]), float(line[11]),
float(line[12])])
rms_vel = float(line[14])
- self._groups.append(LoadedHalo(halo, size, CoM, max_dens_point,
- group_total_mass, max_radius, bulk_vel, rms_vel, fnames))
+ self._groups.append(LoadedHalo(self.pf, halo, size, CoM,
+ max_dens_point, group_total_mass, max_radius, bulk_vel,
+ rms_vel, fnames))
halo += 1
def _collect_halo_data_locations(self):
@@ -2153,7 +2176,7 @@
HaloFinder = HOPHaloFinder
class LoadHaloes(GenericHaloFinder, LoadedHaloList):
- def __init__(self, basename):
+ def __init__(self, pf, basename):
r"""Load the full halo data into memory.
This function takes the output of `GenericHaloFinder.dump` and
@@ -2172,7 +2195,7 @@
>>> halos = LoadHaloes("HopAnalysis")
"""
self.basename = basename
- LoadedHaloList.__init__(self, self.basename)
+ LoadedHaloList.__init__(self, pf, self.basename)
\ No newline at end of file
http://bitbucket.org/yt_analysis/yt/changeset/d8d56a72cda6/
changeset: r3797:d8d56a72cda6
branch: yt
user: sskory
date: 2011-03-05 23:20:22
summary: Adding periodic support for particles to AMR Spheres using ParticleIO.
affected #: 4 files (3.6 KB)
--- a/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 10:48:17 2011 -0800
+++ b/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 14:20:22 2011 -0800
@@ -821,7 +821,7 @@
# Gotta go get it from the halo h5 files.
field_data = self._get_particle_data(self.id, self.fnames,
self.size, key)
- if key == 'particle_position_x': field_data = None
+ #if key == 'particle_position_x': field_data = None
if field_data is not None:
self.saved_fields[key] = field_data
return self.saved_fields[key]
--- a/yt/data_objects/data_containers.py Sat Mar 05 10:48:17 2011 -0800
+++ b/yt/data_objects/data_containers.py Sat Mar 05 14:20:22 2011 -0800
@@ -136,13 +136,19 @@
self.real_grid = grid
self.child_mask = 1
self.ActiveDimensions = self.data['x'].shape
+ self.DW = grid.pf.domain_right_edge - grid.pf.domain_left_edge
+
def __getitem__(self, field):
if field not in self.data.keys():
if field == "RadiusCode":
center = self.field_parameters['center']
- tr = na.sqrt( (self['x'] - center[0])**2.0 +
- (self['y'] - center[1])**2.0 +
- (self['z'] - center[2])**2.0 )
+ tempx = na.abs(self['x'] - center[0])
+ tempx = na.minimum(tempx, self.DW[0] - tempx)
+ tempy = na.abs(self['y'] - center[1])
+ tempy = na.minimum(tempy, self.DW[1] - tempy)
+ tempz = na.abs(self['z'] - center[2])
+ tempz = na.minimum(tempz, self.DW[2] - tempz)
+ tr = na.sqrt( tempx**2.0 + tempy**2.0 + tempz**2.0 )
else:
raise KeyError(field)
else: tr = self.data[field]
--- a/yt/data_objects/particle_io.py Sat Mar 05 10:48:17 2011 -0800
+++ b/yt/data_objects/particle_io.py Sat Mar 05 14:20:22 2011 -0800
@@ -129,7 +129,10 @@
ParticleIOHandler.__init__(self, pf, source)
def _get_args(self):
- return (1, (na.array(self.center, dtype='float64'), self.radius))
+ DLE = na.array(self.pf.domain_left_edge, dtype='float64')
+ DRE = na.array(self.pf.domain_right_edge, dtype='float64')
+ return (1, (na.array(self.center, dtype='float64'), self.radius,
+ 1, DLE, DRE))
class ParticleIOHandlerDisk(ParticleIOHandlerImplemented):
_source_type = "disk"
--- a/yt/utilities/hdf5_light_reader.c Sat Mar 05 10:48:17 2011 -0800
+++ b/yt/utilities/hdf5_light_reader.c Sat Mar 05 14:20:22 2011 -0800
@@ -38,6 +38,7 @@
#define npy_float128 npy_longdouble
#endif
+#define MIN(a,b) ((a) <= (b) ? (a) : (b))
static PyObject *_hdf5ReadError;
herr_t iterate_dataset(hid_t loc_id, const char *name, void *nodelist);
@@ -78,6 +79,8 @@
/* These cannot contain any pointers */
npy_float64 center[3];
npy_float64 radius;
+ npy_float64 period[3];
+ int periodic;
} sphere_validation;
typedef struct cylinder_validation_ {
@@ -992,6 +995,8 @@
/* These are borrowed references */
PyArrayObject *center = (PyArrayObject *) PyTuple_GetItem(InputData, 0);
PyObject *radius = (PyObject *) PyTuple_GetItem(InputData, 1);
+ PyObject *operiodic = PyTuple_GetItem(InputData, 2);
+ npy_float64 DW;
/* This will get freed in the finalization of particle validation */
sphere_validation *sv = (sphere_validation *)
@@ -1004,6 +1009,18 @@
sv->radius = (npy_float64) PyFloat_AsDouble(radius);
+ sv->periodic = PyInt_AsLong(operiodic);
+ if(sv->periodic == 1) {
+ PyArrayObject *domain_left_edge = (PyArrayObject *) PyTuple_GetItem(InputData, 3);
+ PyArrayObject *domain_right_edge = (PyArrayObject *) PyTuple_GetItem(InputData, 4);
+ for (i = 0; i < 3; i++){
+ DW = (*(npy_float64*) PyArray_GETPTR1(domain_right_edge, i))
+ - (*(npy_float64*) PyArray_GETPTR1(domain_left_edge, i));
+ sv->period[i] = DW;
+ //fprintf(stderr, "Setting period equal to %lf\n", sv->period[i]);
+ }
+ }
+
data->count_func = NULL;
data->count_func_float = count_particles_sphere_FLOAT;
data->count_func_double = count_particles_sphere_DOUBLE;
@@ -1523,6 +1540,7 @@
double pradius;
+ if (vdata->periodic == 0) {
for (ind = 0; ind < data->particles_to_check; ind++) {
pradius = 0.0;
tempr = (particle_position_x[ind] - vdata->center[0]); pradius += tempr*tempr;
@@ -1537,6 +1555,25 @@
data->mask[ind] = 0;
}
}
+ } else {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0;
+ tempr = fabs(particle_position_x[ind] - vdata->center[0]);
+ tempr = MIN(tempr, vdata->period[0] - tempr); pradius += tempr*tempr;
+ tempr = fabs(particle_position_y[ind] - vdata->center[1]);
+ tempr = MIN(tempr, vdata->period[1] - tempr); pradius += tempr*tempr;
+ tempr = fabs(particle_position_z[ind] - vdata->center[2]);
+ tempr = MIN(tempr, vdata->period[2] - tempr); 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;
}
@@ -1560,6 +1597,7 @@
double pradius;
+ if (vdata->periodic == 0) {
for (ind = 0; ind < data->particles_to_check; ind++) {
pradius = 0.0;
tempr = (particle_position_x[ind] - vdata->center[0]); pradius += tempr*tempr;
@@ -1574,6 +1612,25 @@
data->mask[ind] = 0;
}
}
+ } else {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0;
+ tempr = fabs(particle_position_x[ind] - vdata->center[0]);
+ tempr = MIN(tempr, vdata->period[0] - tempr); pradius += tempr*tempr;
+ tempr = fabs(particle_position_y[ind] - vdata->center[1]);
+ tempr = MIN(tempr, vdata->period[1] - tempr); pradius += tempr*tempr;
+ tempr = fabs(particle_position_z[ind] - vdata->center[2]);
+ tempr = MIN(tempr, vdata->period[2] - tempr); 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;
}
@@ -1597,6 +1654,7 @@
long double pradius;
+ if (vdata->periodic == 0) {
for (ind = 0; ind < data->particles_to_check; ind++) {
pradius = 0.0;
tempr = (particle_position_x[ind] - vdata->center[0]); pradius += tempr*tempr;
@@ -1611,6 +1669,25 @@
data->mask[ind] = 0;
}
}
+ } else {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0;
+ tempr = fabs(particle_position_x[ind] - vdata->center[0]);
+ tempr = MIN(tempr, vdata->period[0] - tempr); pradius += tempr*tempr;
+ tempr = fabs(particle_position_y[ind] - vdata->center[1]);
+ tempr = MIN(tempr, vdata->period[1] - tempr); pradius += tempr*tempr;
+ tempr = fabs(particle_position_z[ind] - vdata->center[2]);
+ tempr = MIN(tempr, vdata->period[2] - tempr); 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;
}
http://bitbucket.org/yt_analysis/yt/changeset/8ada22103b23/
changeset: r3798:8ada22103b23
branch: yt
user: sskory
date: 2011-03-05 23:51:39
summary: With this commit I believe that halo loading is working.
affected #: 2 files (857 bytes)
--- a/yt/analysis_modules/api.py Sat Mar 05 14:20:22 2011 -0800
+++ b/yt/analysis_modules/api.py Sat Mar 05 14:51:39 2011 -0800
@@ -44,7 +44,8 @@
parallelHF, \
HOPHaloFinder, \
FOFHaloFinder, \
- HaloFinder
+ HaloFinder, \
+ LoadHaloes
from .halo_mass_function.api import \
HaloMassFcn, \
--- a/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 14:20:22 2011 -0800
+++ b/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 14:51:39 2011 -0800
@@ -362,19 +362,30 @@
return None
self.bin_count = bins
# Cosmology
- h = self.data.pf.hubble_constant
- Om_matter = self.data.pf.omega_matter
- z = self.data.pf.current_redshift
+ try:
+ h = self.data.pf.hubble_constant
+ Om_matter = self.data.pf.omega_matter
+ z = self.data.pf.current_redshift
+ period = self.data.pf.domain_right_edge - \
+ self.data.pf.domain_left_edge
+ cm = self.data.pf["cm"]
+ thissize = self.indices.size
+ except AttributeError:
+ # For LoadedHaloes
+ h = self.pf.hubble_constant
+ Om_matter = self.pf.omega_matter
+ z = self.pf.current_redshift
+ period = self.pf.domain_right_edge - \
+ self.pf.domain_left_edge
+ cm = self.pf["cm"]
+ thissize = self.size
rho_crit_now = 1.8788e-29 * h**2.0 * Om_matter # g cm^-3
Msun2g = 1.989e33
rho_crit = rho_crit_now * ((1.0 + z)**3.0)
-
# Get some pertinent information about the halo.
self.mass_bins = na.zeros(self.bin_count+1, dtype='float64')
- dist = na.empty(self.indices.size, dtype='float64')
+ dist = na.empty(thissize, dtype='float64')
cen = self.center_of_mass()
- period = self.data.pf.domain_right_edge - \
- self.data.pf.domain_left_edge
mark = 0
# Find the distances to the particles. I don't like this much, but I
# can't see a way to eliminate a loop like this, either here or in
@@ -399,7 +410,7 @@
# Calculate the over densities in the bins.
self.overdensity = self.mass_bins * Msun2g / \
(4./3. * math.pi * rho_crit * \
- (self.radial_bins * self.data.pf["cm"])**3.0)
+ (self.radial_bins * cm)**3.0)
class HOPHalo(Halo):
@@ -806,6 +817,7 @@
self.overdensity = None
self.saved_fields = {}
self.particle_mask = None
+ self.ds_sort = None
def __getitem__(self, key):
# This function will try to get particle data in one of three ways,
@@ -827,20 +839,21 @@
return self.saved_fields[key]
else:
# Dynamically create the masking array for particles, and get
- # the data using standard yt methods.
+ # the data using standard yt methods. The 1.05 is there to
+ # account for possible silliness having to do with whether
+ # the maximum density or center of mass was used to calculate
+ # the maximum radius.
ds = self.pf.h.sphere(self.CoM, 1.05 * self.max_radius)
if self.particle_mask is None:
pid = self.__getitem__('particle_index')
sp_pid = ds['particle_index']
- print 'testing'
- for id in pid:
- if id not in sp_pid: print 'trouble!', id
+ self.ds_sort = sp_pid.argsort()
+ sp_pid = sp_pid[self.ds_sort]
# The result of searchsorted is an array with the positions
# of the indexes in pid as they are in sp_pid. This is
# because each element of pid is in sp_pid only once.
self.particle_mask = na.searchsorted(sp_pid, pid)
- print 'here'
- return ds[key][self.particle_mask]
+ return ds[key][self.ds_sort][self.particle_mask]
def _get_particle_data(self, halo, fnames, size, field):
# Given a list of file names, a halo, its size, and the desired field,
@@ -2181,7 +2194,10 @@
This function takes the output of `GenericHaloFinder.dump` and
re-establishes the list of halos in memory. This enables the full set
- of halo analysis features without running the halo finder again.
+ of halo analysis features without running the halo finder again. To
+ be precise, the particle data for each halo is only read in when
+ necessary, so examining a single halo will not require as much memory
+ as is required for halo finding.
Parameters
----------
@@ -2192,7 +2208,8 @@
Examples
--------
- >>> halos = LoadHaloes("HopAnalysis")
+ >>> pf = load("data0005")
+ >>> halos = LoadHaloes(pf, "HopAnalysis")
"""
self.basename = basename
LoadedHaloList.__init__(self, pf, self.basename)
http://bitbucket.org/yt_analysis/yt/changeset/b7ed27f9089a/
changeset: r3799:b7ed27f9089a
branch: yt
user: sskory
date: 2011-03-05 23:57:46
summary: Missed some api import lines.
affected #: 2 files (98 bytes)
--- a/yt/analysis_modules/api.py Sat Mar 05 14:51:39 2011 -0800
+++ b/yt/analysis_modules/api.py Sat Mar 05 14:57:46 2011 -0800
@@ -35,11 +35,13 @@
Halo, \
HOPHalo, \
parallelHOPHalo, \
+ LoadedHalo, \
FOFHalo, \
HaloList, \
HOPHaloList, \
FOFHaloList, \
parallelHOPHaloList, \
+ LoadedHaloList, \
GenericHaloFinder, \
parallelHF, \
HOPHaloFinder, \
--- a/yt/analysis_modules/halo_finding/api.py Sat Mar 05 14:51:39 2011 -0800
+++ b/yt/analysis_modules/halo_finding/api.py Sat Mar 05 14:57:46 2011 -0800
@@ -32,13 +32,16 @@
Halo, \
HOPHalo, \
parallelHOPHalo, \
+ LoadedHalo, \
FOFHalo, \
HaloList, \
HOPHaloList, \
FOFHaloList, \
parallelHOPHaloList, \
+ LoadedHaloList, \
GenericHaloFinder, \
parallelHF, \
HOPHaloFinder, \
FOFHaloFinder, \
- HaloFinder
+ HaloFinder, \
+ LoadHaloes
http://bitbucket.org/yt_analysis/yt/changeset/7d7f9ccfa9df/
changeset: r3800:7d7f9ccfa9df
branch: yt
user: sskory
date: 2011-03-06 01:36:59
summary: Adding a mylog.info to alert the user when disk I/O is happening.
affected #: 1 file (268 bytes)
--- a/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 14:57:46 2011 -0800
+++ b/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 16:36:59 2011 -0800
@@ -853,6 +853,9 @@
# of the indexes in pid as they are in sp_pid. This is
# because each element of pid is in sp_pid only once.
self.particle_mask = na.searchsorted(sp_pid, pid)
+ # We won't store this field below in saved_fields because
+ # that would mean keeping two copies of it, one in the yt
+ # machinery and one here.
return ds[key][self.ds_sort][self.particle_mask]
def _get_particle_data(self, halo, fnames, size, field):
@@ -860,6 +863,7 @@
# this returns the particle data for that halo.
# First get the list of fields from the first file. Not all fields
# are saved all the time (e.g. creation_time, particle_type).
+ mylog.info("Getting field %s from hdf5 halo particle files." % field)
f = h5py.File(fnames[0])
fields = f["Halo%08d" % halo].keys()
# If we dont have this field, we can give up right now.
http://bitbucket.org/yt_analysis/yt/changeset/b16467a8e107/
changeset: r3801:b16467a8e107
branch: yt
user: sskory
date: 2011-03-06 17:21:59
summary: Cleaning up the data/pf stuff.
affected #: 1 file (301 bytes)
--- a/yt/analysis_modules/halo_finding/halo_objects.py Sat Mar 05 16:36:59 2011 -0800
+++ b/yt/analysis_modules/halo_finding/halo_objects.py Sun Mar 06 08:21:59 2011 -0800
@@ -75,6 +75,7 @@
self._max_dens = halo_list._max_dens
self.id = id
self.data = halo_list._data_source
+ self.pf = self.data.pf
if indices is not None:
self.indices = halo_list._base_indices[indices]
else:
@@ -362,23 +363,13 @@
return None
self.bin_count = bins
# Cosmology
- try:
- h = self.data.pf.hubble_constant
- Om_matter = self.data.pf.omega_matter
- z = self.data.pf.current_redshift
- period = self.data.pf.domain_right_edge - \
- self.data.pf.domain_left_edge
- cm = self.data.pf["cm"]
- thissize = self.indices.size
- except AttributeError:
- # For LoadedHaloes
- h = self.pf.hubble_constant
- Om_matter = self.pf.omega_matter
- z = self.pf.current_redshift
- period = self.pf.domain_right_edge - \
- self.pf.domain_left_edge
- cm = self.pf["cm"]
- thissize = self.size
+ h = self.pf.hubble_constant
+ Om_matter = self.pf.omega_matter
+ z = self.pf.current_redshift
+ period = self.pf.domain_right_edge - \
+ self.pf.domain_left_edge
+ cm = self.pf["cm"]
+ thissize = max(self.size, self.indices.size)
rho_crit_now = 1.8788e-29 * h**2.0 * Om_matter # g cm^-3
Msun2g = 1.989e33
rho_crit = rho_crit_now * ((1.0 + z)**3.0)
@@ -818,6 +809,7 @@
self.saved_fields = {}
self.particle_mask = None
self.ds_sort = None
+ self.indices = na.array([]) # Never used for a LoadedHalo.
def __getitem__(self, key):
# This function will try to get particle data in one of three ways,
http://bitbucket.org/yt_analysis/yt/changeset/444ffe56a902/
changeset: r3802:444ffe56a902
branch: yt
user: sskory
date: 2011-03-06 17:24:26
summary: Merging
affected #: 0 files (0 bytes)
--- a/yt/analysis_modules/halo_profiler/multi_halo_profiler.py Sun Mar 06 08:21:59 2011 -0800
+++ b/yt/analysis_modules/halo_profiler/multi_halo_profiler.py Sun Mar 06 08:24:26 2011 -0800
@@ -34,8 +34,7 @@
load
from yt.data_objects.profiles import \
BinnedProfile1D, EmptyProfileData
-from yt.analysis_modules.halo_finding.api import \
- HaloFinder
+from yt.analysis_modules.halo_finding.api import *
from .halo_filters import \
VirialFilter
from yt.data_objects.field_info_container import \
@@ -54,15 +53,23 @@
class HaloProfiler(ParallelAnalysisInterface):
"Radial profiling, filtering, and projections for halos in cosmological simulations."
- def __init__(self, dataset, halos='multiple', halo_list_file='HopAnalysis.out', halo_list_format='yt_hop',
- halo_finder_function=HaloFinder, halo_finder_args=None, halo_finder_kwargs=None,
- use_density_center=False, density_center_exponent=1.0, use_field_max_center=None,
+ def __init__(self, dataset, output_dir=None,
+ halos='multiple', halo_list_file='HopAnalysis.out',
+ halo_list_format='yt_hop', halo_finder_function=parallelHF,
+ halo_finder_args=None,
+ halo_finder_kwargs=dict(threshold=160.0, safety=1.5,
+ dm_only=False, resize=True,
+ fancy_padding=True, rearrange=True),
+ use_density_center=False, density_center_exponent=1.0,
+ use_field_max_center=None,
halo_radius=0.1, radius_units='1', n_profile_bins=50,
profile_output_dir='radial_profiles', projection_output_dir='projections',
projection_width=8.0, projection_width_units='mpc', project_at_level='max',
velocity_center=['bulk', 'halo'], filter_quantities=['id','center']):
"""
Initialize a HaloProfiler object.
+ :param output_dir (str): if specified, all output will be put into this path instead of
+ in the dataset directories. Default: None.
:param halos (str): "multiple" for profiling more than one halo. In this mode halos are read in
from a list or identified with a halo finder. In "single" mode, the one and only halo
center is identified automatically as the location of the peak in the density field.
@@ -114,7 +121,7 @@
"""
self.dataset = dataset
-
+ self.output_dir = output_dir
self.profile_output_dir = profile_output_dir
self.projection_output_dir = projection_output_dir
self.n_profile_bins = n_profile_bins
@@ -132,6 +139,10 @@
self.filtered_halos = []
self._projection_halo_list = []
+ # Create output directory if specified
+ if self.output_dir is not None:
+ self.__check_directory(self.output_dir)
+
# Set halo finder function and parameters, if needed.
self.halo_finder_function = halo_finder_function
self.halo_finder_args = halo_finder_args
@@ -153,7 +164,8 @@
# dictionary: a dictionary containing fields and their corresponding columns.
self.halo_list_file = halo_list_file
if halo_list_format == 'yt_hop':
- self.halo_list_format = {'id':0, 'mass':1, 'center':[7, 8, 9], 'velocity':[10, 11, 12], 'r_max':13}
+ self.halo_list_format = {'id':0, 'mass':1, 'np': 2,
+ 'center':[7, 8, 9], 'velocity':[10, 11, 12], 'r_max':13}
elif halo_list_format == 'enzo_hop':
self.halo_list_format = {'id':0, 'center':[4, 5, 6]}
elif halo_list_format == 'p-groupfinder':
@@ -169,7 +181,8 @@
self.density_center_exponent = density_center_exponent
if self.use_density_center:
def _MatterDensityXTotalMass(field, data):
- return na.power((data['Matter_Density'] * data['TotalMassMsun']), self.density_center_exponent)
+ return na.power((data['Matter_Density'] * data['TotalMassMsun']),
+ self.density_center_exponent)
def _Convert_MatterDensityXTotalMass(data):
return 1
add_field("MatterDensityXTotalMass", units=r"",
@@ -288,8 +301,14 @@
# Add profile fields necessary for calculating virial quantities.
if virial_filter: self._check_for_needed_profile_fields()
- outputDir = "%s/%s" % (self.pf.fullpath, self.profile_output_dir)
- self.__check_directory(outputDir)
+ # Create output directory.
+ if self.output_dir is not None:
+ self.__check_directory("%s/%s" % (self.output_dir, self.pf.directory))
+ my_output_dir = "%s/%s/%s" % (self.output_dir, self.pf.directory,
+ self.profile_output_dir)
+ else:
+ my_output_dir = "%s/%s" % (self.pf.fullpath, self.profile_output_dir)
+ self.__check_directory(my_output_dir)
# Profile all halos.
for halo in self._get_objs('all_halos', round_robin=True):
@@ -305,7 +324,7 @@
if filter_result and len(self.profile_fields) > 0:
- profile_filename = "%s/Halo_%04d_profile.dat" % (outputDir, halo['id'])
+ profile_filename = "%s/Halo_%04d_profile.dat" % (my_output_dir, halo['id'])
profiledHalo = self._get_halo_profile(halo, profile_filename, virial_filter=virial_filter)
@@ -456,8 +475,14 @@
(self.pf.parameters['RefineBy']**proj_level)
projectionResolution = int(self.projection_width / proj_dx)
- outputDir = "%s/%s" % (self.pf.fullpath, self.projection_output_dir)
- self.__check_directory(outputDir)
+ # Create output directory.
+ if self.output_dir is not None:
+ self.__check_directory("%s/%s" % (self.output_dir, self.pf.directory))
+ my_output_dir = "%s/%s/%s" % (self.output_dir, self.pf.directory,
+ self.projection_output_dir)
+ else:
+ my_output_dir = "%s/%s" % (self.pf.fullpath, self.projection_output_dir)
+ self.__check_directory(my_output_dir)
center = [0.5 * (self.pf.parameters['DomainLeftEdge'][w] + self.pf.parameters['DomainRightEdge'][w])
for w in range(self.pf.parameters['TopGridRank'])]
@@ -519,7 +544,7 @@
if save_cube:
axis_labels = ['x', 'y', 'z']
dataFilename = "%s/Halo_%04d_%s_data.h5" % \
- (outputDir, halo['id'], axis_labels[w])
+ (my_output_dir, halo['id'], axis_labels[w])
mylog.info("Saving projection data to %s." % dataFilename)
output = h5py.File(dataFilename, "a")
@@ -535,7 +560,7 @@
output.close()
if save_images:
- pc.save("%s/Halo_%04d" % (outputDir, halo['id']), force_save=True)
+ pc.save("%s/Halo_%04d" % (my_output_dir, halo['id']), force_save=True)
del region
@@ -573,13 +598,17 @@
if filename is None:
filename = self.halo_list_file
- hopFile = "%s/%s" % (self.pf.fullpath, filename)
+ if self.output_dir is not None:
+ self.__check_directory("%s/%s" % (self.output_dir, self.pf.directory))
+ hop_file = "%s/%s/%s" % (self.output_dir, self.pf.directory, filename)
+ else:
+ hop_file = "%s/%s" % (self.pf.fullpath, filename)
- if not(os.path.exists(hopFile)):
- mylog.info("Hop file not found, running hop to get halos.")
- self._run_hop(hopFile)
+ if not(os.path.exists(hop_file)):
+ mylog.info("Halo finder file not found, running halo finder to get halos.")
+ self._run_hop(hop_file)
- self.all_halos = self._read_halo_list(hopFile)
+ self.all_halos = self._read_halo_list(hop_file)
def _read_halo_list(self, listFile):
"""
@@ -685,11 +714,11 @@
return None
@parallel_blocking_call
- def _run_hop(self, hopFile):
+ def _run_hop(self, hop_file):
"Run hop to get halos."
hop_results = self.halo_finder_function(self.pf, *self.halo_finder_args, **self.halo_finder_kwargs)
- hop_results.write_out(hopFile)
+ hop_results.write_out(hop_file)
del hop_results
self.pf.h.clear_all_data()
@@ -752,13 +781,13 @@
fid.close()
@parallel_root_only
- def __check_directory(self, outputDir):
- if (os.path.exists(outputDir)):
- if not(os.path.isdir(outputDir)):
- mylog.error("Output directory exists, but is not a directory: %s." % outputDir)
- raise IOError(outputDir)
+ def __check_directory(self, my_output_dir):
+ if (os.path.exists(my_output_dir)):
+ if not(os.path.isdir(my_output_dir)):
+ mylog.error("Output directory exists, but is not a directory: %s." % my_output_dir)
+ raise IOError(my_output_dir)
else:
- os.mkdir(outputDir)
+ os.mkdir(my_output_dir)
def shift_projections(pf, pc, oldCenter, newCenter, axis):
"""
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list