[Yt-svn] yt-commit r769 - trunk/yt/lagos/hop

[mturk at wrangler.dreamhost.com]

Author: mturk
Date: Thu Sep 11 07:00:21 2008
New Revision: 769
URL: http://yt.spacepope.org/changeset/769

Log:
By making a couple more numpy views, as well as reducing the number of calls to
expensive re-orderings of the particles, I have sped up the parts of HopOutput that
deal with particle data -- this will likely be visible as a global speedup in
dealing with Halo collections.

Also, I fixed a couple typos in the C code, added a debugger for free-ing (that
is currently commented out) and added a tiny bit more output information.



Modified:
   trunk/yt/lagos/hop/EnzoHop.c
   trunk/yt/lagos/hop/HopOutput.py
   trunk/yt/lagos/hop/hop.h

Modified: trunk/yt/lagos/hop/EnzoHop.c
==============================================================================
--- trunk/yt/lagos/hop/EnzoHop.c	(original)
+++ trunk/yt/lagos/hop/EnzoHop.c	Thu Sep 11 07:00:21 2008
@@ -90,7 +90,7 @@
     mass    = (PyArrayObject *) PyArray_FromAny(omass,
                     PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
                     NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
-    if((!zpos)||(PyArray_SIZE(mass) != num_particles)) {
+    if((!mass)||(PyArray_SIZE(mass) != num_particles)) {
     PyErr_Format(_HOPerror,
              "EnzoHop: xpos and mass must be the same length.");
     goto _fail;
@@ -113,19 +113,24 @@
   
  	/* Copy positions into kd structure. */
 
+    fprintf(stdout, "Filling in %d particles\n", num_particles);
 	for (i = 0; i < num_particles; i++) {
-	  kd->p[i].r[0] = (float)*(npy_float64*) PyArray_GETPTR1(xpos, i);
-	  kd->p[i].r[1] = (float)*(npy_float64*) PyArray_GETPTR1(ypos, i);
-	  kd->p[i].r[2] = (float)*(npy_float64*) PyArray_GETPTR1(zpos, i);
+	  kd->p[i].r[0] = (float)(*(npy_float64*) PyArray_GETPTR1(xpos, i));
+	  kd->p[i].r[1] = (float)(*(npy_float64*) PyArray_GETPTR1(ypos, i));
+	  kd->p[i].r[2] = (float)(*(npy_float64*) PyArray_GETPTR1(zpos, i));
 	  kd->p[i].fMass = (float)(*(npy_float64*) PyArray_GETPTR1(mass, i)/totalmass);
 	}
 
     HC my_comm;
     my_comm.s = newslice();
     my_comm.gl = (Grouplist*)malloc(sizeof(Grouplist));
+    if(my_comm.gl == NULL) {
+        fprintf(stderr, "failed allocating Grouplist\n");
+        goto _fail;
+    }
     initgrouplist(my_comm.gl);
 
-    fprintf(stderr, "Calling hop... %d\n",num_particles);
+    fprintf(stderr, "Calling hop... %d %0.3e\n",num_particles,thresh);
     hop_main(kd, &my_comm, thresh);
 
     fprintf(stderr, "Calling regroup...\n");

Modified: trunk/yt/lagos/hop/HopOutput.py
==============================================================================
--- trunk/yt/lagos/hop/HopOutput.py	(original)
+++ trunk/yt/lagos/hop/HopOutput.py	Thu Sep 11 07:00:21 2008
@@ -43,16 +43,17 @@
         self.__run_hop()
         mylog.info("Parsing outputs")
         self.__parse_output()
-        mylog.debug("Finished.")
+        mylog.debug("Finished. (%s)", len(self))
 
     def __obtain_particles(self):
         if self.dm_only: ii = self.__get_dm_indices()
         else: ii = slice(None)
-        self._base_indices = ii
         self.particle_fields = {}
         for field in ["particle_position_%s" % ax for ax in 'xyz'] + \
                      ["ParticleMassMsun"]:
+            tot_part = self.data_source[field].size
             self.particle_fields[field] = self.data_source[field][ii]
+        self._base_indices = na.arange(tot_part)[ii]
 
     def __run_hop(self):
         self.densities, self.tags = \
@@ -80,6 +81,7 @@
         counts = na.bincount(self.tags+1)
         sort_indices = na.argsort(self.tags)
         grab_indices = na.indices(self.tags.shape).ravel()[sort_indices]
+        dens = self.densities[sort_indices]
         cp = 0
         for i in unique_ids:
             cp_c = cp + counts[i+1]
@@ -88,10 +90,10 @@
                 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])
+            md_i = na.argmax(dens[cp:cp_c])
             px, py, pz = [self.particle_fields['particle_position_%s'%ax][group_indices]
                                             for ax in 'xyz']
-            self._max_dens[i] = (self.densities[sort_indices][cp:cp_c][md_i],
+            self._max_dens[i] = (dens[cp:cp_c][md_i],
                                  px[md_i], py[md_i], pz[md_i])
             cp += counts[i+1]
 
@@ -147,8 +149,7 @@
         self.hop_output = hop_output
         self.id = id
         self.data = hop_output.data_source
-        self.indices = indices
-        self._base_indices = hop_output._base_indices
+        self.indices = hop_output._base_indices[indices]
         
     def center_of_mass(self):
         """
@@ -210,7 +211,7 @@
         return r.max()
 
     def __getitem__(self, key):
-        return self.data[key][self._base_indices][self.indices]
+        return self.data[key][self.indices]
 
     def get_sphere(self, center_of_mass=True):
         """

Modified: trunk/yt/lagos/hop/hop.h
==============================================================================
--- trunk/yt/lagos/hop/hop.h	(original)
+++ trunk/yt/lagos/hop/hop.h	Thu Sep 11 07:00:21 2008
@@ -1,4 +1,5 @@
 #include "slice.h"
+//#define free(A) if(A==NULL)fprintf(stderr,"FREEING DOUBLE\n");fprintf(stderr,"Freeing "#A" ("__FILE__":%d)\n",__LINE__);free(A);
 
 /* ----------------------------------------------------------------------- */
 /* The following structures track all the information about the groups */