start of shape measurement in void outputs; some improvements and bug fixes to AP analysis

c_tools/stacking/pruneVoids.cpp

@@ -10,6 +10,7 @@
 //   for producing a "clean" void catalog for other purposes.
 
 #include "gsl/gsl_math.h"
+#include "gsl/gsl_eigen.h"
 #include "string.h"
 #include "ctype.h"
 #include "stdlib.h"
@@ -43,6 +44,8 @@ typedef struct voidStruct {
   float nearestEdge;
   float center[3], barycenter[3];
   int accepted;
+  gsl_vector *eval;
+  gsl_matrix *evec;
 } VOID;
 
 int main(int argc, char **argv) {
@@ -75,6 +78,7 @@ int main(int argc, char **argv) {
   int numVoids, mockIndex, numKept;
   double tolerance;
   FILE *fp, *fpBarycenter, *fpDistances, *fpSkyPositions, *fpInfo;
+  FILE *fpShapes;
   PART *part, *voidPart;
   ZONE2PART *zones2Parts;
   VOID2ZONE *void2Zones;
@@ -92,6 +96,8 @@ int main(int argc, char **argv) {
   int clock1, clock2;
   int periodicX=0, periodicY=0, periodicZ=0;
 
+  gsl_eigen_symmv_workspace *eigw = gsl_eigen_symmv_alloc(3);
+
   numVoids = args_info.numVoids_arg;
   mockIndex = args_info.mockIndex_arg;
   tolerance = args_info.tolerance_arg;
@@ -209,6 +215,9 @@ int main(int argc, char **argv) {
 
     voids[i-1].radius = pow(voidVol/volNorm*3./4./M_PI, 1./3.);
     voids[i-1].accepted = 1;
+  
+    voids[i-1].eval = gsl_vector_alloc(3);
+    voids[i-1].evec = gsl_matrix_alloc(3,3);
   }
   fclose(fp);
 
@@ -356,11 +365,10 @@ int main(int argc, char **argv) {
         voids[iVoid].barycenter[2] = boxLen[2] - voids[iVoid].barycenter[2];
     }
 
-
     // compute central density
     centralRad = voids[iVoid].radius/args_info.centralRadFrac_arg;
-    centralRad *= centralRad;
     centralDen = 0.;
+    int numCentral = 0;
     for (p = 0; p < voids[iVoid].numPart; p++) {
       dist[0] = fabs(voidPart[p].x - voids[iVoid].barycenter[0]);
       dist[1] = fabs(voidPart[p].y - voids[iVoid].barycenter[1]);
@@ -371,9 +379,10 @@ int main(int argc, char **argv) {
       if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]);
 
       dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2);
-      if (dist2 < centralRad) centralDen += 1;
+      if (sqrt(dist2) < centralRad) numCentral += 1;
     }
-    voids[iVoid].centralDen = centralDen / (4./3. * M_PI * pow(centralRad, 3./2.));
+    voids[iVoid].centralDen = numCentral / (volNorm*4./3. * M_PI * 
+                                            pow(centralRad, 3.));
 
     // compute maximum extent
 /*
@@ -467,8 +476,32 @@ int main(int argc, char **argv) {
     }
 
     voids[iVoid].nearestEdge = nearestEdge;
+
+    // compute eigenvalues and vectors for orientation and shape
+    double inertia[9];
+    for (int p = 0; p < voids[iVoid].numPart; p++) {
+      dist[0] = voidPart[p].x - voids[iVoid].barycenter[0];
+      dist[1] = voidPart[p].y - voids[iVoid].barycenter[1];
+      dist[2] = voidPart[p].z - voids[iVoid].barycenter[2];
+
+      for (int i = 0; i < 3; i++) {
+      for (int j = i; j < 3; j++) {
+        if (i == j) inertia[i*3+j] += dist[0]*dist[0] + dist[1]*dist[1] + 
+                                    dist[2]*dist[2];
+        inertia[i*3+j] -= dist[i]*dist[j];      
+      }
+      }
+    }
+    inertia[1*3+0] = inertia[0*3+1];
+    inertia[2*3+0] = inertia[0*3+2];
+    inertia[2*3+1] = inertia[1*3+2];
+
+    gsl_matrix_view m = gsl_matrix_view_array(inertia, 3, 3);
+    gsl_eigen_symmv(&m.matrix, voids[iVoid].eval, voids[iVoid].evec, eigw);
   } // iVoid
 
+  gsl_eigen_symmv_free(eigw);
+    
   int numWrong = 0;
   int numHighDen = 0;
   int numEdge = 0;
@@ -515,14 +548,14 @@ int main(int argc, char **argv) {
     // *always* clean out near edges since there are no mocks there
     if (tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestEdge) {
       voids[iVoid].accepted = -3;
-      numEdge++;
+      if (voids[iVoid].accepted == 1) numEdge++;
     }
 
     // assume the lower z-boundary is "soft" in observations
     if (args_info.isObservation_flag && 
         voids[iVoid].redshift < args_info.zMin_arg) {
       voids[iVoid].accepted = -3;
-      numEdge++;
+      if (voids[iVoid].accepted == 1) numEdge++;
     }
   }
 
@@ -543,10 +576,12 @@ int main(int argc, char **argv) {
   fpInfo = fopen(args_info.outInfo_arg, "w");
   fpDistances = fopen(args_info.outDistances_arg, "w");
   fpSkyPositions = fopen(args_info.outSkyPositions_arg, "w");
+  fpShapes = fopen(args_info.outShapes_arg, "w");
   fprintf(fp, "%d particles, %d voids.\n", mockIndex, numKept);
   fprintf(fp, "see column in master void file\n");
   fprintf(fpInfo, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast\n");
   fprintf(fpSkyPositions, "# RA, dec, redshift, radius (Mpc/h), void ID\n");
+  fprintf(fpShapes, "# void ID, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n");
   for (iVoid = 0; iVoid < numVoids; iVoid++) {
 
     if (voids[iVoid].accepted != 1) continue;
@@ -566,10 +601,6 @@ int main(int argc, char **argv) {
 
      fprintf(fpBarycenter, "%d %e %e %e\n", 
              voids[iVoid].voidID,
-// TEST
-             //voids[iVoid].center[0],
-             //voids[iVoid].center[1],
-             //voids[iVoid].center[2]);
              voids[iVoid].barycenter[0],
              voids[iVoid].barycenter[1],
              voids[iVoid].barycenter[2]);
@@ -579,19 +610,11 @@ int main(int argc, char **argv) {
              voids[iVoid].nearestMock);
 
      double outCenter[3];
-// TEST
-     //outCenter[0] = voids[iVoid].center[0];
-     //outCenter[1] = voids[iVoid].center[1];
-     //outCenter[2] = voids[iVoid].center[2];
      outCenter[0] = voids[iVoid].barycenter[0];
      outCenter[1] = voids[iVoid].barycenter[1];
      outCenter[2] = voids[iVoid].barycenter[2];
 
      if (args_info.isObservation_flag) {
-// TEST
-       //outCenter[0] = (voids[iVoid].center[0]-boxLen[0]/2.)*100.;
-       //outCenter[1] = (voids[iVoid].center[1]-boxLen[1]/2.)*100.;
-       //outCenter[2] = (voids[iVoid].center[2]-boxLen[2]/2.)*100.;
        outCenter[0] = (voids[iVoid].barycenter[0]-boxLen[0]/2.)*100.;
        outCenter[1] = (voids[iVoid].barycenter[1]-boxLen[1]/2.)*100.;
        outCenter[2] = (voids[iVoid].barycenter[2]-boxLen[2]/2.)*100.;
@@ -609,12 +632,29 @@ int main(int argc, char **argv) {
              voids[iVoid].densCon);
 
      fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n",
-             atan((voids[iVoid].barycenter[1]-boxLen[1]/2.)/(voids[iVoid].barycenter[0]-boxLen[0]/2.)) * 180/M_PI + 180,  
+             atan((voids[iVoid].barycenter[1]-boxLen[1]/2.) / 
-             asin((voids[iVoid].barycenter[2]-boxLen[2]/2.)/voids[iVoid].redshiftInMpc) * 180/M_PI,  
+                  (voids[iVoid].barycenter[0]-boxLen[0]/2.)) * 180/M_PI + 180,  
+             asin((voids[iVoid].barycenter[2]-boxLen[2]/2.) / 
+                   voids[iVoid].redshiftInMpc) * 180/M_PI,  
              voids[iVoid].redshift,
              voids[iVoid].radius,
              voids[iVoid].voidID);
 
+     fprintf(fpShapes, "%d %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n", 
+             voids[iVoid].voidID,
+             gsl_vector_get(voids[iVoid].eval, 0),
+             gsl_vector_get(voids[iVoid].eval, 1),
+             gsl_vector_get(voids[iVoid].eval, 2),
+             gsl_matrix_get(voids[iVoid].evec, 0 ,0),
+             gsl_matrix_get(voids[iVoid].evec, 0 ,1),
+             gsl_matrix_get(voids[iVoid].evec, 0 ,2),
+             gsl_matrix_get(voids[iVoid].evec, 1 ,0),
+             gsl_matrix_get(voids[iVoid].evec, 1 ,1),
+             gsl_matrix_get(voids[iVoid].evec, 1 ,2),
+             gsl_matrix_get(voids[iVoid].evec, 2 ,0),
+             gsl_matrix_get(voids[iVoid].evec, 2 ,1),
+             gsl_matrix_get(voids[iVoid].evec, 2 ,2)
+            );
   }
   fclose(fp);
   fclose(fpInfo);

python_tools/void_python_tools/backend/launchers.py

@@ -293,6 +293,9 @@ def launchPrune(sample, binPath, thisDataPortion=None,
     cmd += " --outSkyPositions=" + zobovDir+"/sky_positions_"+\
                           str(thisDataPortion)+"_"+\
                           str(sampleName)+".out"
+    cmd += " --outShapes=" + zobovDir+"/shapes_"+\
+                          str(thisDataPortion)+"_"+\
+                          str(sampleName)+".out"
     cmd += " --outDistances=" + zobovDir+"/boundaryDistances_"+\
                           str(thisDataPortion)+"_"+\
                           str(sampleName)+".out"
@@ -872,6 +875,14 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
       print "no voids here; skipping!"
       return
 
+    numVoids = int(open(voidDir+"/num_voids.txt", "r").readline())
+    if numVoids < 15:
+      print "not enough voids to fit; skipping!"
+      fp = open(voidDir+"/NOFIT", "w")
+      fp.write("not enough voids: %d \n" % numVoids)
+      fp.close()
+      return
+
     if stack.zMin < sample.zRange[0] or stack.zMax > sample.zRange[1]:
       print "outside sample redshift range; skipping!"
       return
@@ -918,7 +929,7 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
       #badChain = (args[0][0] > 0.5 or args[0][1] > stack.rMax or \
       #            args[0][2] > stack.rMax) and \
       #           (ntries < maxtries)
-      ret,fits,args = vp.compute_inertia(voidDir, stack.rMax)
+      ret,fits,args = vp.compute_inertia(voidDir, stack.rMax, mode="symmetric")
       badChain = False
       ntries += 1
 
@@ -983,7 +994,7 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
 def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
                  INCOHERENT=None, workDir=None, figDir=None, errorBars=None, 
                  ZOBOV_PATH=None, continueRun=None, voidDir=None, 
-                 doPlot = True):
+                 doPlot = True, setName=None):
 
   for thisDataPortion in dataPortions:
     print "    For data portion", thisDataPortion
@@ -1225,7 +1236,8 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
           plotTitle = ''
         else:
           #plotTitle = "all samples, "+thisDataPortion+" voids"
-          plotTitle = ''
+          plotTitle = setName
+          #plotTitle = ''
         vp.do_all_obs(zbase, allExpList, aveDistList,
                       rlist, plotTitle=plotTitle, sampleNames=shortSampleNames,
                       plotAve=True, mulfac = 1.0, biasLine = 1.16,