[Yt-svn] yt-commit r555 - in trunk/yt: lagos/hop raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri Jun 13 16:07:19 PDT 2008
Author: mturk
Date: Fri Jun 13 16:07:17 2008
New Revision: 555
URL: http://yt.spacepope.org/changeset/555
Log:
Couple fixes to the HOP stuff -- first off, previously we were not correctly
correcting for the possible presence of stars. Now the particles you are
returned are definitely *just* the ones that HOP was run on.
I removed the superfluous .next(self) method in the hop list.
Calculation of the radius now accounts for periodicity.
The center of mass is NOT accounted for with periodicity.
I fixed some SHOW STOPPING BUGS in the center_of_mass calculation and the
bulk_velocity calculation. Oops.
I've also added a HaloMassesPositionPlot object. (This should probably get a
new name.) It accepts a hop_results list and plots the halo positions and
radii, true-to-life. I'll add labels another day.
I've also restricted tags to a 32 bit integer on the backend, which fixed some
issues with bincount. (If you have more than INT_MAX halos, yowza!)
Modified:
trunk/yt/lagos/hop/EnzoHop.c
trunk/yt/lagos/hop/HopOutput.py
trunk/yt/raven/Plot3DInterface.py
Modified: trunk/yt/lagos/hop/EnzoHop.c
==============================================================================
--- trunk/yt/lagos/hop/EnzoHop.c (original)
+++ trunk/yt/lagos/hop/EnzoHop.c Fri Jun 13 16:07:17 2008
@@ -153,15 +153,15 @@
PyArray_DescrFromType(NPY_FLOAT64));
PyArrayObject *particle_group_id = (PyArrayObject *)
PyArray_SimpleNewFromDescr(1, PyArray_DIMS(xpos),
- PyArray_DescrFromType(NPY_INT64));
+ PyArray_DescrFromType(NPY_INT32));
for (i = 0; i < num_particles; i++) {
// Density is in kd->p[i].fDensity
*(npy_float64*)(PyArray_GETPTR1(particle_density, i)) =
(npy_float64) kd->p[i].fDensity;
// tag is in gl.s->ntag[i+1]
- *(npy_int64*)(PyArray_GETPTR1(particle_group_id, i)) =
- (npy_int64) my_comm.s->ntag[i+1];
+ *(npy_int32*)(PyArray_GETPTR1(particle_group_id, i)) =
+ (npy_int32) my_comm.s->ntag[i+1];
}
kdFinish(kd);
Modified: trunk/yt/lagos/hop/HopOutput.py
==============================================================================
--- trunk/yt/lagos/hop/HopOutput.py (original)
+++ trunk/yt/lagos/hop/HopOutput.py Fri Jun 13 16:07:17 2008
@@ -42,6 +42,7 @@
def __run_hop(self):
if self.dm_only: ii = self.__get_dm_indices()
else: ii = slice(None)
+ self._base_indices = ii
self.densities, self.tags = \
RunHOP(self.data_source["particle_position_x"][ii],
self.data_source["particle_position_y"][ii],
@@ -62,17 +63,19 @@
def __parse_output(self):
unique_ids = na.unique(self.tags)
- bins = na.digitize(self.tags, unique_ids)
- counts = na.bincount(self.tags+1)
+ counts = na.bincount((self.tags+1)
sort_indices = na.argsort(self.tags)
grab_indices = na.indices(self.tags.shape).ravel()[sort_indices]
cp = 0
for i in unique_ids:
cp_c = cp + counts[i+1]
+ if i == -1:
+ cp += counts[i+1]
+ continue
group_indices = grab_indices[cp:cp_c]
self._groups.append(HopGroup(self, i, group_indices))
md_i = na.argmax(self.densities[sort_indices][cp:cp_c])
- px, py, pz = [self.data_source['particle_position_%s'%ax][group_indices]
+ px, py, pz = [self.data_source['particle_position_%s'%ax][self._base_indices][group_indices]
for ax in 'xyz']
self._max_dens[i] = (self.densities[sort_indices][cp:cp_c][md_i],
px[md_i], py[md_i], pz[md_i])
@@ -84,17 +87,12 @@
def __iter__(self):
return HopIterator(self)
- def next(self):
- if self.__index == len(self): raise StopIteration
- tr = self[self.__index]
- self.__index += 1
- return tr
-
def __getitem__(self, key):
return self._groups[key]
def write_out(self, filename="HopAnalysis.out"):
f = open(filename,"w")
+ f.write("# Center of mass does NOT account for periodicity!\n")
f.write("\t".join(["# Group","Mass","# part","max dens"
"x","y","z", "center-of-mass",
"x","y","z",
@@ -115,7 +113,7 @@
class HopIterator(object):
def __init__(self, hop):
self.hop = hop
- self.index = 0
+ self.index = -1
def next(self):
self.index += 1
@@ -129,12 +127,14 @@
self.id = id
self.data = hop_output.data_source
self.indices = indices
+ self._base_indices = hop_output._base_indices
def center_of_mass(self):
+ # Center of mass does not account for periodicity and is likely wrong!
pm = self["ParticleMassMsun"]
cx = (self["particle_position_x"] * pm).sum()
- cy = (self["particle_position_x"] * pm).sum()
- cz = (self["particle_position_x"] * pm).sum()
+ cy = (self["particle_position_y"] * pm).sum()
+ cz = (self["particle_position_z"] * pm).sum()
return na.array([cx,cy,cz])/pm.sum()
def maximum_density(self):
@@ -151,19 +151,23 @@
def bulk_velocity(self):
pm = self["ParticleMassMsun"]
vx = (self["particle_velocity_x"] * pm).sum()
- vy = (self["particle_velocity_x"] * pm).sum()
- vz = (self["particle_velocity_x"] * pm).sum()
+ vy = (self["particle_velocity_y"] * pm).sum()
+ vz = (self["particle_velocity_z"] * pm).sum()
return na.array([vx,vy,vz])/pm.sum()
def maximum_radius(self):
- center = self.center_of_mass()
- r = na.sqrt((self["particle_position_x"]-center[0])**2.0
- + (self["particle_position_y"]-center[1])**2.0
- + (self["particle_position_z"]-center[2])**2.0)
+ center = self.maximum_density_location()
+ #center = self.center_of_mass()
+ rx = na.abs(self["particle_position_x"]-center[0])
+ ry = na.abs(self["particle_position_y"]-center[1])
+ rz = na.abs(self["particle_position_z"]-center[2])
+ r = na.sqrt(na.minimum(rx, 1.0-rx)**2.0
+ + na.minimum(ry, 1.0-ry)**2.0
+ + na.minimum(rz, 1.0-rz)**2.0)
return r.max()
def __getitem__(self, key):
- return self.data[key][self.indices]
+ return self.data[key][self._base_indices][self.indices]
def get_sphere(self, center_of_mass=True):
if center_of_mass: center = self.center_of_mass()
Modified: trunk/yt/raven/Plot3DInterface.py
==============================================================================
--- trunk/yt/raven/Plot3DInterface.py (original)
+++ trunk/yt/raven/Plot3DInterface.py Fri Jun 13 16:07:17 2008
@@ -141,7 +141,7 @@
def __register_callbacks(self):
# More should go here for functional changes to the object
- s2plot.cs2scb(self.__my_callback) # Install a dynamic callback
+ s2plot.cs2scb(self.__my_callback) # Install a dynamic callback
def __my_callback(self, t, kc):
s2plot.ds2dvr(self.vrid, 0)
@@ -203,3 +203,60 @@
def setup_plot_points(self):
xyz = [f[0](self.profile._data_source[f[1]]) for f in self.bf]
self.xyz = [xyz[0].size] + xyz + [1]
+
+class HaloMassesPositionPlot(object):
+ def __init__(self, hop_results, window_opts="/S2MONO",
+ cmap="jet"):
+ self.hop_results = hop_results
+ self.window_opts = window_opts
+ self.cmap = cmap
+ self.started = False
+ self.__setup_data()
+
+ def __setup_data(self):
+ c, r, m = [], [], []
+ for g in self.hop_results:
+ c.append(g.center_of_mass())
+ r.append(g.maximum_radius())
+ m.append(g.total_mass())
+ self.centers = na.array(c)
+ self.radii = na.array(r)
+ self.masses = na.log10(na.array(m))
+ self._dmax = self.masses.max()
+ self._dmin = self.masses.min()
+
+ def __setup_spheres(self):
+ cm = matplotlib.cm.get_cmap(self.cmap)
+ scaled_masses = ((self.masses-self._dmin)/(self._dmax-self._dmin))
+ for i in xrange(self.radii.size):
+ r,g,b,a = cm(scaled_masses[i])
+ s2plot.ns2sphere(self.centers[i,0],
+ self.centers[i,1],
+ self.centers[i,2],
+ self.radii[i], r,g,b)
+
+ def run(self, pre_call=None):
+ self.__setup_s2plot()
+ self.__setup_spheres()
+ self.__setup_labels()
+ self.__register_callbacks()
+ if pre_call is not None: pre_call(self)
+ s2plot.s2disp(-1, 1)
+
+ def __setup_labels(self):
+ pass
+
+ def __setup_s2plot(self):
+ s2plot.s2opendo(self.window_opts)
+ s2plot.s2swin(0.,1., 0.,1., 0.,1.)
+ opts = "BCDE"
+ s2plot.s2box(opts,0,0, opts,0,0, opts,0,0)
+
+ amb = {'r':0.8, 'g':0.8, 'b':0.8} # ambient light
+ s2plot.ss2srm(s2plot.SHADE_FLAT); # Set shading type to FLAT
+ s2plot.ss2sl(amb, 0, None, None, 0) # Ambient lighting only
+
+ self.started = True
+
+ def __register_callbacks(self):
+ pass
More information about the yt-svn
mailing list