[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