Merge branch 'master' of https://bitbucket.org/cosmicvoids/void_identification

c_tools/analysis/voidOverlap.cpp

@@ -88,6 +88,8 @@ int main(int argc, char **argv) {
   int closestMatchID;
   float closestMatchDist;
   float commonVolRatio;
+  MATCHPROPS newMatch;
+  int MAX_MATCHES = 20;
 
   CATALOG catalog1, catalog2;
 
@@ -134,38 +136,44 @@ int main(int argc, char **argv) {
       printf("Will assume z-direction is periodic.\n");
     }
 
-  printf(" Determining overlap...\n");
+  // find closest voids
-
+  printf(" Finding nearest matches...\n");
-/*
-  // just use centers
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
-    closestMatchDist = 1.e99;
     for (iVoid2 = 0; iVoid2 < catalog2.numVoids; iVoid2++) {
       rdist = getDist(catalog1, catalog2, iVoid1, iVoid2, 
                       periodicX, periodicY, periodicZ);
 
-      if (rdist < closestMatchDist) {
+  
-        closestMatchID = iVoid2;
+      newMatch.matchID = iVoid2;
-        closestMatchDist = rdist;
+      newMatch.commonVol = 0;
+      newMatch.dist = rdist; 
+
+      if (catalog1.voids[iVoid1].matches.size() < MAX_MATCHES) {
+        catalog1.voids[iVoid1].matches.push_back(newMatch);
+      } else {
+        // find the farthest match
+        float farthestMatchDist = 0;
+        int   farthestMatchID = 0;
+        for (iMatch = 0; iMatch < MAX_MATCHES; iMatch++) {
+          if (catalog1.voids[iVoid1].matches[iMatch].dist > farthestMatchDist){
+            farthestMatchDist = catalog1.voids[iVoid1].matches[iMatch].dist;
+            farthestMatchID   = iMatch;
+          }           
+        } 
+        if (rdist < farthestMatchDist)
+          catalog1.voids[iVoid1].matches[farthestMatchID] = newMatch;
       }
-
     }
-    
-    MATCHPROPS newMatch;
-    newMatch.matchID = closestMatchID;
-    newMatch.commonVol = 1;
-    newMatch.dist = closestMatchDist; 
-    catalog1.voids[iVoid1].matches.push_back(newMatch);
   } 
-*/
 
+  printf(" Determining overlap...\n");
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
-    printf("  Working on void %d of %d...\n", iVoid1, catalog1.numVoids); 
+    printf("  Working on void %d of %d...\n", iVoid1+1, catalog1.numVoids); 
     voidID1 = catalog1.voids[iVoid1].voidID;
 
-    for (iVoid2 = 0; iVoid2 < catalog2.numVoids; iVoid2++) {
+    for (iMatch = 0; iMatch < catalog1.voids[iVoid1].matches.size();iMatch++) {
+      iVoid2 = catalog1.voids[iVoid1].matches[iMatch].matchID;
       voidID2 = catalog2.voids[iVoid2].voidID;
-      match = false;
 
       for (iZ1 = 0; iZ1 < catalog1.void2Zones[voidID1].numZones; iZ1++) {
         zoneID1 = catalog1.void2Zones[voidID1].zoneIDs[iZ1];
@@ -180,6 +188,7 @@ int main(int argc, char **argv) {
               partID2 = catalog2.zones2Parts[zoneID2].partIDs[p2];     
 
               match = false;
+
               if (args.useID_flag) {
                 if (catalog1.part[partID1].uniqueID == 
                     catalog2.part[partID2].uniqueID) match = true;
@@ -203,35 +212,12 @@ int main(int argc, char **argv) {
                 r2 = pow(3./4./M_PI*catalog2.part[partID2].volume / 
                          catalog2.numPartTot, 1./3.);
                 if (rdist <= 0.1*r1 || rdist <= 0.1*r2) match = true;
-             }
+              }
-       
        
               if (match) {
-                bool alreadyMatched = false;
+                catalog1.voids[iVoid1].matches[iMatch].commonVol += 
-                for (iMatch = 0; 
+                  catalog2.part[partID2].volume;
-                     iMatch < catalog1.voids[iVoid1].matches.size(); iMatch++) {
+              } // end if match
-                   if (catalog1.voids[iVoid1].matches[iMatch].matchID == 
-                       iVoid2) {
-                     alreadyMatched = true;
-                     break;
-                   }
-                 } 
-
-                 if (alreadyMatched) {
-                   //catalog1.voids[iVoid1].matches[iMatch].commonVol += 1;
-                   catalog1.voids[iVoid1].matches[iMatch].commonVol += 
-                     catalog1.part[partID1].volume;
-                 } else {
-                    MATCHPROPS newMatch;
-                    newMatch.matchID = iVoid2;
-                    //newMatch.commonVol = 1;
-                    newMatch.commonVol = catalog1.part[partID1].volume;
-    
-                    newMatch.dist = getDist(catalog1, catalog2, iVoid1, iVoid2,
-                                            periodicX, periodicY, periodicZ);
-                    catalog1.voids[iVoid1].matches.push_back(newMatch);
-                } // end if match
-              } 
 
             } // end p2
           } // end iZ2
@@ -240,11 +226,13 @@ int main(int argc, char **argv) {
     } // end iVoid2
   } // end match finding
 
+  printf(" Sorting matches...\n");
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
     sortMatches(catalog1.voids[iVoid1].matches);
   }
 
   // count up significant matches
+  printf(" Categorizing matches...\n");
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
     closestMatchDist = 0.;
     for (iMatch = 0; iMatch < catalog1.voids[iVoid1].matches.size(); iMatch++) {
@@ -257,12 +245,12 @@ int main(int argc, char **argv) {
   }
 
   // output summary
+  printf(" Output...\n");
   std::string filename;
   filename = string(args.outfile_arg);
-  filename = filename.append("_summary.out");
+  filename = filename.append("summary.out");
   fp = fopen(filename.c_str(), "w");
-  //fp = fopen(args.outfile_arg, "w");
+  fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, common volume ratio 2, relative dist, num matches, num significant matches\n");
-  fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, relative dist, num matches, num significant matches\n");
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
     int voidID = catalog1.voids[iVoid1].voidID;
     if (catalog1.voids[iVoid1].numMatches > 0) {
@@ -272,27 +260,31 @@ int main(int argc, char **argv) {
       commonVolRatio = catalog1.voids[iVoid1].matches[0].commonVol / 
                        //catalog1.voids[iVoid1].numPart;
                        catalog1.voids[iVoid1].vol;
+      float volRatio = catalog1.voids[iVoid1].matches[0].commonVol / 
+                       catalog2.voids[iVoid2].vol;
       rdist = catalog1.voids[iVoid1].matches[0].dist;
       rdist /= catalog1.voids[iVoid1].radius;
 
-      fprintf(fp, "%d %.2f %.2f %.2f %.2f %d %d\n", voidID, 
+      fprintf(fp, "%d %.2f %.2f %.2f %.2f %.2f %d %d\n", voidID, 
                                    catalog1.voids[iVoid1].radiusMpc,
                                    rRatio, 
                                    commonVolRatio,
+                                   volRatio,
                                    rdist, 
                                    catalog1.voids[iVoid1].numMatches,
                                    catalog1.voids[iVoid1].numBigMatches);
 
     } else {
-      fprintf(fp, "%d %f 0.0 0.0 0.0 0 0\n", voidID, 
+      fprintf(fp, "%d %.2f 0.0 0.0 0.0 0.0 0 0\n", voidID, 
-                                           catalog1.voids[iVoid1].radius);
+                                           catalog1.voids[iVoid1].radiusMpc);
     }
   } // end printing
   fclose(fp);
 
   // output detail
+  printf(" Output detail...\n");
   filename = string(args.outfile_arg);
-  filename = filename.append("_detail.out");
+  filename = filename.append("detail.out");
   fp = fopen(filename.c_str(), "w");
   int MAX_OUT = 10;
   fprintf(fp, "# void ID, match common vol\n");
@@ -300,14 +292,13 @@ int main(int argc, char **argv) {
     int voidID = catalog1.voids[iVoid1].voidID;
     fprintf(fp,"%d: ", voidID);
     for (iMatch = 0; iMatch < MAX_OUT; iMatch++) {
-
+      if (iMatch < catalog1.voids[iVoid1].matches.size()) {
-      if (iMatch < catalog1.voids[iVoid].matches.size()) {
         commonVolRatio = catalog1.voids[iVoid1].matches[iMatch].commonVol / 
                          //catalog1.voids[iVoid1].numPart;
                          catalog1.voids[iVoid1].vol;
 
-        fprintf(fp, "%.2f ", commonVolRatio);
+        fprintf(fp, "%.3f ", catalog1.voids[iVoid1].matches[iMatch].dist);
-
+        //fprintf(fp, "%.2f ", commonVolRatio);
       } else {
         fprintf(fp, "0.00 ");
       }
@@ -443,10 +434,6 @@ void loadCatalog(const char *partFile, const char *volFile,
     catalog.voids[i].voidProb = voidProb;
 
     catalog.voids[i].radius = pow(voidVol/catalog.numPartTot*3./4./M_PI, 1./3.);
-    //catalog.voids[i].biggestMatchID = -1;
-    //catalog.voids[i].biggestMatchVol = 0;
-    //catalog.voids[i].closestMatchID = -1;
-    //catalog.voids[i].closestMatchDist = 1.e99;
     catalog.voids[i].numMatches = 0;
     catalog.voids[i].numBigMatches = 0;
   
@@ -548,14 +535,15 @@ void sortMatches(std::vector<MATCHPROPS>& matches) {
   MATCHPROPS tempMatch;
   bool swapped;
 
-  if (matches.size() == 1) return;
+  if (matches.size() <= 1) return;
 
-  swapped = false;
+  swapped = true;
-  while (!swapped) {
+  while (swapped) {
     swapped = false;
     for (int iMatch = 0; iMatch < matches.size() - 1; iMatch++) {
-      if (matches[iMatch].commonVol > matches[iMatch+1].commonVol) {
+      if (matches[iMatch].dist > matches[iMatch+1].dist) {
-        tempMatch = matches[iMatch];
+      //if (matches[iMatch].commonVol < matches[iMatch+1].commonVol) {
+        tempMatch = matches[iMatch+1];
         matches[iMatch+1] = matches[iMatch];
         matches[iMatch] = tempMatch;
         swapped = true;

c_tools/libzobov/loadZobov.cpp

@@ -1,3 +1,4 @@
+#include <cassert>
 #include <string>
 #include <fstream>
 #include <iostream>
@@ -152,6 +153,7 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName,
       for (int j = 0; j < z.allVoids[volId].zId.size(); j++)
         {
           int zzid = z.allVoids[volId].zId[j];
+          assert(zzid < z.allZones.size());
 	  actualNumber += z.allZones[zzid].pId.size();
         }
       

c_tools/libzobov/particleInfo.cpp

@@ -1,3 +1,4 @@
+#include <cstdlib>
 #include <netcdfcpp.h>
 #include <CosmoTool/fortran.hpp>
 #include "particleInfo.hpp"
@@ -5,24 +6,55 @@
 using namespace std;
 using namespace CosmoTool;
 
+template<bool failure>
+double retrieve_attr_safe_double(NcFile& f, const char *name, double defval)
+{
+  NcAtt *a = f.get_att(name);
+  if (a == 0)
+    {
+      if (failure)
+        abort();
+      return defval;
+    }
+  return a->as_double(0);
+}
+
+template<bool failure>
+int retrieve_attr_safe_int(NcFile& f, const char *name, int defval)
+{
+  NcAtt *a = f.get_att(name);
+  if (a == 0)
+    {
+      if (failure)
+        abort();
+      return defval;
+    }
+  return a->as_int(0);
+}
+
+  
+
 bool loadParticleInfo(ParticleInfo& info,
 		      const std::string& particles, 
 		      const std::string& extra_info)
 {
       int numpart;
+      int isObservation;
 
   NcFile f_info(extra_info.c_str());
+  NcError nerr(NcError::verbose_nonfatal);
  
   if (!f_info.is_valid())
     return false;
 
-  info.ranges[0][0] = f_info.get_att("range_x_min")->as_double(0);
+  info.ranges[0][0] = retrieve_attr_safe_double<true>(f_info, "range_x_min", 0);
-  info.ranges[0][1] = f_info.get_att("range_x_max")->as_double(0);
+  info.ranges[0][1] = retrieve_attr_safe_double<true>(f_info, "range_x_max", 0);
-  info.ranges[1][0] = f_info.get_att("range_y_min")->as_double(0);
+  info.ranges[1][0] = retrieve_attr_safe_double<true>(f_info, "range_y_min", 0);
-  info.ranges[1][1] = f_info.get_att("range_y_max")->as_double(0);
+  info.ranges[1][1] = retrieve_attr_safe_double<true>(f_info, "range_y_max", 0);
-  info.ranges[2][0] = f_info.get_att("range_z_min")->as_double(0);
+  info.ranges[2][0] = retrieve_attr_safe_double<true>(f_info, "range_z_min", 0);
-  info.ranges[2][1] = f_info.get_att("range_z_max")->as_double(0);
+  info.ranges[2][1] = retrieve_attr_safe_double<true>(f_info, "range_z_max", 0);
-  info.mask_index = f_info.get_att("mask_index")->as_int(0); //PMS
+  info.mask_index = retrieve_attr_safe_int<true>(f_info, "mask_index", 0);
+  isObservation = retrieve_attr_safe_int<false>(f_info, "is_observation", 0);
 
   for (int i = 0; i < 3; i++)
     info.length[i] = info.ranges[i][1] - info.ranges[i][0];
@@ -59,6 +91,14 @@ bool loadParticleInfo(ParticleInfo& info,
       for (int i = 0; i < numpart; i++)
 	info.particles[i].z = mul*f.readReal32();
       f.endCheckpoint();
+
+      if (!isObservation) {
+        for (int i = 0; i < numpart; i++) {
+	        info.particles[i].x += info.ranges[0][0];
+	        info.particles[i].y += info.ranges[1][0];
+	        info.particles[i].z += info.ranges[2][0];
+        }
+      }
     }
   catch (const NoSuchFileException& e)
     {

c_tools/mock/generateFromCatalog.cpp

@@ -466,6 +466,7 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn
   fp.add_att("range_z_min", -Lmax/100.);
   fp.add_att("range_z_max", Lmax/100.);
   fp.add_att("mask_index", pdata.mask_index); // PMS
+  fp.add_att("is_observation", 1); // PMS
 
   int nOutputPart = pdata.mask_index;
   //int nOutputPart = pdata.pos.size();

c_tools/mock/generateMock.cpp

@@ -221,7 +221,12 @@ void generateOutput(SimuData *data, int axis,
 // cleared. That way new particles can be appended if this is a multi-file snapshot.
 void selectBox(SimuData *simu, std::vector<long>& targets, generateMock_info& args_info)
 {
-  float subsample = args_info.subsample_given ? args_info.subsample_arg : 1.0;
+  float subsample;
+  if (args_info.subsample_given) {
+    subsample = args_info.subsample_arg;
+  } else {
+    subsample = 1.0;
+  }
   double ranges[3][2] = {
     { args_info.rangeX_min_arg, args_info.rangeX_max_arg },
     { args_info.rangeY_min_arg, args_info.rangeY_max_arg },
@@ -312,9 +317,9 @@ void buildBox(SimuData *simu, long num_targets, long loaded,
       
       for (int j = 0; j < 3; j++)
 	{
-	  boxed->Pos[j][loaded] = (simu->Pos[j][pid]-ranges[j*2])*mul[j];
+	  boxed->Pos[j][loaded] = max(min((simu->Pos[j][pid]-ranges[j*2])*mul[j], double(1)), double(0));
-	  assert(boxed->Pos[j][loaded] > 0);
+	  assert(boxed->Pos[j][loaded] >= 0);
-	  assert(boxed->Pos[j][loaded] < 1);
+	  assert(boxed->Pos[j][loaded] <= 1);
 	}
       uniqueID[loaded] = (simu_uniqueID != 0) ? simu_uniqueID[pid] : 0;
       expansion_fac[loaded] = efac[pid];
@@ -338,6 +343,7 @@ void saveBox(SimuData *&boxed, const std::string& outbox)
   f.add_att("range_z_min", ranges[4]);
   f.add_att("range_z_max", ranges[5]);
   f.add_att("mask_index", -1);
+  f.add_att("is_observation", 0);
 
   NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
   NcDim *NumSnap_dim = f.add_dim("numsnap_dim", num_snapshots);
@@ -394,6 +400,7 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a
   mul = new double[3];
   ranges = new double[6];
   snapshot_split = new long[*num_snapshots];
+  expansion_fac = new double[boxed->NumPart];
 
 
   boxed->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
@@ -402,6 +409,7 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a
   boxed->new_attribute("particle_id", particle_id, delete_adaptor<long>);
   boxed->new_attribute("num_snapshots", num_snapshots, delete_adaptor<int>);
   boxed->new_attribute("snapshot_split", snapshot_split, delete_adaptor<long>);
+  boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor<double>);
 
   v_id->get(particle_id, boxed->NumPart);
   v_snap->get(snapshot_split, *num_snapshots);

c_tools/mock/loaders/gadget_loader.cpp

@@ -77,7 +77,7 @@ public:
 
 SimulationLoader *gadgetLoader(const std::string& snapshot, double Mpc_unitLength, int flags)
 {
-  bool singleFile;
+  bool singleFile = false;
   int num_files;
   SimuData *d;
 

c_tools/stacking/pruneVoids.cpp

@@ -87,7 +87,7 @@ int main(int argc, char **argv) {
   double nearestEdge, redshift;
   char line[500], junkStr[10];
   int mask_index;
-  double ranges[2][3], boxLen[3], mul; 
+  double ranges[3][2], boxLen[3], mul; 
   double volNorm, radius;
   int clock1, clock2;
   int periodicX=0, periodicY=0, periodicZ=0;
@@ -103,6 +103,9 @@ int main(int argc, char **argv) {
                                              args_info.periodic_arg);
 
   // check for periodic box
+  periodicX = 0;
+  periodicY = 0;
+  periodicZ = 0;
   if (!args_info.isObservation_flag) {
     if ( strchr(args_info.periodic_arg, 'x') != NULL) {
       periodicX = 1;
@@ -143,9 +146,9 @@ int main(int argc, char **argv) {
   temp = (float *) malloc(numPartTot * sizeof(float));
  
   volNorm = numPartTot/(boxLen[0]*boxLen[1]*boxLen[2]);
-  printf("VOL NORM = %f\n", volNorm);
+  printf(" VOL NORM = %f\n", volNorm);
 
-  printf("CENTRAL DEN = %f\n", args_info.maxCentralDen_arg);
+  printf(" CENTRAL DEN = %f\n", args_info.maxCentralDen_arg);
 
   fread(&dummy, 1, 4, fp);
   fread(temp, numPartTot, 4, fp);
@@ -311,13 +314,16 @@ int main(int argc, char **argv) {
     voids[iVoid].barycenter[2] = 0.;
 
     for (p = 0; p < voids[iVoid].numPart; p++) {
-      dist[0] = fabs(voidPart[p].x - voids[iVoid].center[0]);
+      dist[0] = voidPart[p].x - voids[iVoid].center[0];
-      dist[1] = fabs(voidPart[p].y - voids[iVoid].center[1]);
+      dist[1] = voidPart[p].y - voids[iVoid].center[1];
-      dist[2] = fabs(voidPart[p].z - voids[iVoid].center[2]);
+      dist[2] = voidPart[p].z - voids[iVoid].center[2];
 
-      if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]);
+      if (periodicX && fabs(dist[0]) > boxLen[0]/2.) 
-      if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]);
+          dist[0] = dist[0] - copysign(boxLen[0], dist[0]);
-      if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]);
+      if (periodicY && fabs(dist[1]) > boxLen[1]/2.) 
+          dist[1] = dist[1] - copysign(boxLen[1], dist[1]);
+      if (periodicZ && fabs(dist[2]) > boxLen[2]/2.) 
+          dist[2] = dist[2] - copysign(boxLen[2], dist[2]);
 
       voids[iVoid].barycenter[0] += voidPart[p].vol*(dist[0]);
       voids[iVoid].barycenter[1] += voidPart[p].vol*(dist[1]);
@@ -331,12 +337,25 @@ int main(int argc, char **argv) {
     voids[iVoid].barycenter[1] += voids[iVoid].center[1];
     voids[iVoid].barycenter[2] += voids[iVoid].center[2];
 
-    if (periodicX) 
+    if (periodicX) {
-      voids[iVoid].barycenter[0] = fmod(voids[iVoid].barycenter[0], boxLen[0]);
+      if (voids[iVoid].barycenter[0] > boxLen[0])
-    if (periodicY) 
+        voids[iVoid].barycenter[0] = voids[iVoid].barycenter[0] - boxLen[0];
-      voids[iVoid].barycenter[1] = fmod(voids[iVoid].barycenter[1], boxLen[1]);
+      if (voids[iVoid].barycenter[0] < 0)
-    if (periodicZ) 
+        voids[iVoid].barycenter[0] = boxLen[0] + voids[iVoid].barycenter[0];
-      voids[iVoid].barycenter[2] = fmod(voids[iVoid].barycenter[2], boxLen[2]);
+    }
+    if (periodicY) {
+      if (voids[iVoid].barycenter[1] > boxLen[1])
+        voids[iVoid].barycenter[1] = voids[iVoid].barycenter[1] - boxLen[1];
+      if (voids[iVoid].barycenter[1] < 0)
+        voids[iVoid].barycenter[1] = boxLen[1] + voids[iVoid].barycenter[1];
+    }
+    if (periodicZ) {
+      if (voids[iVoid].barycenter[2] > boxLen[2])
+        voids[iVoid].barycenter[2] = voids[iVoid].barycenter[2] - boxLen[2];
+      if (voids[iVoid].barycenter[2] < 0)
+        voids[iVoid].barycenter[2] = boxLen[2] + voids[iVoid].barycenter[2];
+    }
+
 
     // compute central density
     centralRad = voids[iVoid].radius/args_info.centralRadFrac_arg;
@@ -378,8 +397,8 @@ int main(int argc, char **argv) {
       for (p = 0; p < voids[iVoid].numPart; p++) {
   
         dist[0] = fabs(voidPart[p].x - voids[iVoid].barycenter[0]);
-        dist[0] = fabs(voidPart[p].y - voids[iVoid].barycenter[1]);
+        dist[1] = fabs(voidPart[p].y - voids[iVoid].barycenter[1]);
-        dist[0] = fabs(voidPart[p].z - voids[iVoid].barycenter[2]);
+        dist[2] = fabs(voidPart[p].z - voids[iVoid].barycenter[2]);
 
         if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]);
         if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]);
@@ -450,6 +469,11 @@ int main(int argc, char **argv) {
     voids[iVoid].nearestEdge = nearestEdge;
   } // iVoid
 
+  int numWrong = 0;
+  int numHighDen = 0;
+  int numEdge = 0;
+  int numTooSmall = 0;
+
   printf(" Picking winners and losers...\n");
   for (iVoid = 0; iVoid < numVoids; iVoid++) {
     voids[iVoid].accepted = 1;
@@ -462,36 +486,43 @@ int main(int argc, char **argv) {
  
     if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) {
       voids[iVoid].accepted = -1;
+      numHighDen++;
     }
 
     // toss out voids that are obviously wrong
-    if (voids[iVoid].densCon > 1.e3) {
+    if (voids[iVoid].densCon > 1.e4) {
       voids[iVoid].accepted = -4;
+      numWrong++;
     }
 
     if (strcmp(args_info.dataPortion_arg, "edge")  == 0 &&
         tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) {
       voids[iVoid].accepted = -3;
+      numEdge++;
     }
 
     if (strcmp(args_info.dataPortion_arg, "central") == 0 && 
         tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestMock) {
       voids[iVoid].accepted = -3;
+      numEdge++;
     }
 
     if (voids[iVoid].radius < args_info.rMin_arg) {
       voids[iVoid].accepted = -2;
+      numTooSmall++;
     }
 
     // *always* clean out near edges since there are no mocks there
     if (tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestEdge) {
       voids[iVoid].accepted = -3;
+      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++;
     }
   }
 
@@ -501,6 +532,10 @@ int main(int argc, char **argv) {
   }
 
   printf(" Number kept: %d (out of %d)\n", numKept, numVoids);
+  printf("   Rejected %d near the edge\n", numEdge); 
+  printf("   Rejected %d too small\n", numTooSmall); 
+  printf("   Rejected %d obviously bad\n", numWrong); 
+  printf("   Rejected %d too high central density\n", numHighDen); 
   
   printf(" Output...\n");
   fp = fopen(args_info.output_arg, "w");
@@ -529,8 +564,12 @@ int main(int argc, char **argv) {
              voids[iVoid].densCon, 
              voids[iVoid].voidProb);
 
-     fprintf(fpBarycenter, "%d %e %e %e\n", 
+     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]);
@@ -540,11 +579,19 @@ 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.;

crossCompare/analysis/datasetsToAnalyze.py

@@ -2,14 +2,14 @@
 
 
 workDir = "/home/psutter2/workspace/Voids/"
-figDir  = "./figs"
+dataDir = "/home/psutter2/workspace/Voids/multidark/crossCompare/mergerTree/"
-dataDir = "./data"
  
 
 CTOOLS_PATH = "../../c_tools/"
 
+baseSampleDir = "multidark/md_ss0.000175/sample_md_ss0.000175_z0.56_d00/"
+
 sampleDirList = [ 
-                  "multidark/md_ss0.000175/sample_md_ss0.000175_z0.56_d00/",
                   "multidark/md_ss0.0004/sample_md_ss0.0004_z0.56_d00/",
                 ]
 

crossCompare/analysis/mergerTree.py

@@ -30,30 +30,30 @@ if not os.access(filename, os.F_OK):
 parms = imp.load_source("name", filename)
 globals().update(vars(parms))
 
-if not os.access(figDir, os.F_OK):
+if not os.access(dataDir, os.F_OK):
-  os.makedirs(figDir)
+  os.makedirs(dataDir)
 
 outFileName = dataDir + "/" + dataNameBase #+ ".dat"
 
+with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
+  baseSample = pickle.load(input)
+
 for (iSample, sampleDir) in enumerate(sampleDirList):
-  if iSample == 0: continue
 
   with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
     sample = pickle.load(input)
 
-  with open(workDir+sampleDirList[0]+"/sample_info.dat", 'rb') as input:
+  print " Working with", sample.fullName, "...",
-    baseSample = pickle.load(input)
-
-  print " Working with", sample.fullName, "..."
 
   sampleName = sample.fullName
 
   binPath = CTOOLS_PATH+"/analysis/voidOverlap"
   logFile = os.getcwd()+"/mergerTree.log"
-  stepOutputFileName = outFileName + "_" + sampleName + "_"
+  stepOutputFileName = outFileName + "_" + baseSample.fullName + "_" + \
+                       sampleName + "_"
   #stepOutputFileName = os.getcwd()+"/data/overlap_"
 
-  launchVoidOverlap(baseSample, sample, workDir+sampleDirList[0], 
+  launchVoidOverlap(baseSample, sample, workDir+baseSampleDir, 
                     workDir+sampleDir, binPath, 
                     thisDataPortion="central", logFile=logFile,
                     continueRun=False, workDir=workDir,

pipeline/generateCatalog.py

@@ -20,6 +20,7 @@ if (len(sys.argv) > 1):
     print "  Cannot find parameter file %s!" % filename
     exit(-1)
   parms = imp.load_source("name", filename)
+  regenerateFlag = False
   globals().update(vars(parms))
 else:
   print " Using default parameters"
@@ -90,7 +91,7 @@ for sample in dataSampleList:
     launchGenerate(sample, GENERATE_PATH, workDir=workDir, 
                    inputDataDir=inputDataDir, zobovDir=zobovDir,
                    figDir=figDir, logFile=logFile, useLCDM=useLCDM,
-                   continueRun=continueRun)
+                   continueRun=continueRun, regenerate=regenerateFlag)
 
   # --------------------------------------------------------------------------
   if (startCatalogStage <= 2) and (endCatalogStage >= 2) and not sample.isCombo:

python_tools/pipeline_source/prepareCatalogs.in.py

@@ -79,8 +79,9 @@ startCatalogStage = 1
 endCatalogStage   = 4
                
 startAPStage = 1
-endAPStage = 7
+endAPStage = 1
 
+regenerateFlag = False
 ZOBOV_PATH = "@CMAKE_BINARY_DIR@/zobov/"
 CTOOLS_PATH = "@CMAKE_BINARY_DIR@/c_tools/"
 freshStack = True
@@ -124,7 +125,7 @@ newSample = Sample(dataFile = "{dataFile}",
                    zBoundaryMpc = ({zMinMpc}, {zMaxMpc}),
                    omegaM    = {omegaM},
                    minVoidRadius = {minRadius},
-                   profileBinSize = 1.0,
+                   profileBinSize = "auto",
                    includeInHubble = True,
                    partOfCombo = False,
                    isCombo = False,
@@ -135,12 +136,20 @@ newSample = Sample(dataFile = "{dataFile}",
                    useLightCone = {useLightCone},
                    subsample = {subsample})
 dataSampleList.append(newSample)
-newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+2, True, False)
+newSample.addStack(0.0, 5.0, 20, 25, False, False)
-newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+4, True, False)
+newSample.addStack(0.0, 5.0, 30, 35, False, False)
-newSample.addStack({zMin}, {zMax}, 2*{minRadius}+2, 2*{minRadius}+6, True, False)
+newSample.addStack(0.0, 5.0, 40, 45, False, False)
-newSample.addStack({zMin}, {zMax}, 2*{minRadius}+6, 2*{minRadius}+10, True, False)
+newSample.addStack(0.0, 5.0, 50, 55, False, False)
-newSample.addStack({zMin}, {zMax}, 2*{minRadius}+10, 2*{minRadius}+18, True, False)
+newSample.addStack(0.0, 5.0, 60, 65, False, False)
-newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, False)
+newSample.addStack(0.0, 5.0, 70, 75, False, False)
+newSample.addStack(0.0, 5.0, 80, 85, False, False)
+newSample.addStack(0.0, 5.0, 90, 95, False, False)
+#newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+2, True, False)
+#newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+4, True, False)
+#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+2, 2*{minRadius}+6, True, False)
+#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+6, 2*{minRadius}+10, True, False)
+#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+10, 2*{minRadius}+18, True, False)
+#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, False)
                """
   for (iFile, redshift) in enumerate(redshifts):
     fileNum = fileNums[iFile]

python_tools/void_python_tools/apTools/profiles/__init__.py

@@ -1,6 +1,7 @@
 from build import *
 from draw import *
 from fit import *
+from inertia import *
 from mcmc import *
 from generateExpFigure import *
 from getSurveyProps import *

python_tools/void_python_tools/backend/launchers.py

@@ -24,11 +24,18 @@ ncFloat = 'f8' # Double precision
 # -----------------------------------------------------------------------------
 def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, 
                    zobovDir=None, figDir=None, logFile=None, useLCDM=False,
-                   continueRun=None):
+                   continueRun=None,regenerate=False):
 
   if sample.dataType == "observation":
     sampleName = sample.fullName
 
+    if regenerate:
+      inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par"
+      outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName
+    else:
+      inputParameterFlag = ""
+      outputFile = zobovDir+"/zobov_slice_" + sampleName
+
     if sample.dataFile == "":
       datafile = inputDataDir+"/"+sampleName
     else:
@@ -50,16 +57,17 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
       zMax %g
       density_fake %g
       %s
-      """ % (datafile, maskFile, zobovDir+"/zobov_slice_"+sampleName,
+      %s
+      """ % (datafile, maskFile, outputFile,
              zobovDir+"/zobov_slice_"+sampleName+".par",
              sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity,
-             useLCDMFlag)
+             useLCDMFlag, inputParameterFlag)
 
     parmFile = os.getcwd()+"/generate_"+sample.fullName+".par"
 
     file(parmFile, mode="w").write(conf)
 
-    if not (continueRun and jobSuccessful(logFile, "Done!\n")):
+    if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")):
       cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile)
       os.system(cmd)
       if jobSuccessful(logFile, "Done!\n"):
@@ -97,6 +105,13 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
 
     datafile = inputDataDir+"/"+sample.dataFile
 
+    if regenerate:
+      inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par"
+      outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName 
+    else:
+      inputParameterFlag = ""
+      outputFile = zobovDir+"/zobov_slice_" + sampleName
+
     if sample.usePecVel:
       includePecVelString = "peculiarVelocities"
     else:
@@ -134,20 +149,21 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
       rangeZ_min %g
       rangeZ_max %g
       subsample %g
-      """ % (dataFileLine, zobovDir+"/zobov_slice_"+sampleName,
+      %s
+      """ % (dataFileLine, outputFile,
              zobovDir+"/zobov_slice_"+sampleName+".par",
              includePecVelString,
              useLightConeString,
              sample.dataUnit,
              xMin, xMax, yMin, yMax,
              sample.zBoundaryMpc[0], sample.zBoundaryMpc[1],
-             sample.subsample)
+             sample.subsample,inputParameterFlag)
 
     parmFile = os.getcwd()+"/generate_"+sample.fullName+".par"
 
     file(parmFile, mode="w").write(conf)
 
-    if not (continueRun and jobSuccessful(logFile, "Done!\n")):
+    if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")):
       cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile)
       os.system(cmd)
       if jobSuccessful(logFile, "Done!\n"):
@@ -258,6 +274,7 @@ def launchPrune(sample, binPath, thisDataPortion=None,
        sample.boxLen <= 1.e-1: 
       periodicLine += "z"
     periodicLine += "' "
+    periodicLine = ""
 
   if not (continueRun and jobSuccessful(logFile, "NetCDF: Not a valid ID\n")):
     cmd = binPath
@@ -297,6 +314,9 @@ def launchPrune(sample, binPath, thisDataPortion=None,
                           str(thisDataPortion)+"_"+\
                           str(sampleName)+".out"
     cmd += " &> " + logFile
+    f=file("run_prune.sh",mode="w")
+    f.write(cmd)
+    f.close()
     os.system(cmd)
 
     if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or \
@@ -314,7 +334,8 @@ def launchPrune(sample, binPath, thisDataPortion=None,
 def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir, 
                       binPath, thisDataPortion=None, 
                       logFile=None, workDir=None,
-                      continueRun=None, outputFile=None):
+                      continueRun=None, outputFile=None, 
+                      matchMethod=None):
 
   sampleName1 = sample1.fullName
   sampleName2 = sample2.fullName
@@ -358,7 +379,7 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir,
     cmd += " --zonePartFile2=" + sample2Dir+"/voidPart_" + \
            str(sampleName2)+".dat"
 
-    cmd += " --useID"
+    if matchMethod == "useID": cmd += " --useID"
     cmd += periodicLine
     cmd += " --outfile=" + outputFile
     cmd += " &> " + logFile
@@ -379,14 +400,14 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
                 voidDir=None, freshStack=True, runSuffix=None,
                 zobovDir=None,
                 INCOHERENT=False, ranSeed=None, summaryFile=None, 
-                continueRun=None, dataType=None):
+                continueRun=None, dataType=None, prefixRun=""):
 
   sampleName = sample.fullName
 
   runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin,
                              stack.rMax, thisDataPortion)
 
-  logFile = logDir+"/stack_"+sampleName+"_"+runSuffix+".out"
+  logFile = logDir+("/%sstack_"%prefixRun)+sampleName+"_"+runSuffix+".out"
  
   treeFile = voidDir+"/tree.data"
 
@@ -413,7 +434,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
     maxDen = 0.2*float(maskIndex)/float(totalPart)
   else:
     maskIndex = 999999999
-    maxDen = 0.2
+    maxDen = -0.8
 
   if INCOHERENT:
     nullTestFlag = "INCOHERENT"
@@ -454,7 +475,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
    zobovDir+"/vol_"+sampleName+".dat",
    stack.rMin,
    stack.rMax,
-   zobovDir+"/zobov_slice_"+sampleName,
+   zobovDir+("/%szobov_slice_"%prefixRun)+sampleName,
    zobovDir+"/zobov_slice_"+sampleName+".par",
    maxDen,
    centralRadius,
@@ -470,7 +491,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
    zobovDir+"/boundaryDistances_"+thisDataPortion+"_"+sampleName+".out",
    rescaleFlag)
 
-  parmFile = os.getcwd()+"/stack_"+sample.fullName+".par"
+  parmFile = os.getcwd()+("/%sstack_"%prefixRun)+sample.fullName+".par"
 
   fp = file(parmFile, mode="w")
   fp.write(conf)
@@ -600,13 +621,18 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
 
     os.system("mv %s %s" % ("tree.data", treeFile))
     os.system("mv %s %s" % ("void_indexes.txt", voidDir+"/"))
-    os.system("mv %s %s" % ("posx.nc", voidDir+"/posx.nc"))
+    os.system("mv %s %s" % ("posx.nc", voidDir+"/"))
-    os.system("mv %s %s" % ("posy.nc", voidDir+"/posy.nc"))
+    os.system("mv %s %s" % ("posy.nc", voidDir+"/"))
-    os.system("mv %s %s" % ("posz.nc", voidDir+"/posz.nc"))
+    os.system("mv %s %s" % ("posz.nc", voidDir+"/"))
-    os.system("mv %s %s" % ("redshifts.nc", voidDir+"/redshifts.nc"))
+    os.system("mv %s %s" % ("z_void_indexes.txt", voidDir+"/"))
+    os.system("mv %s %s" % ("z_posx.nc", voidDir+"/"))
+    os.system("mv %s %s" % ("z_posy.nc", voidDir+"/"))
+    os.system("mv %s %s" % ("z_posz.nc", voidDir+"/"))
+    os.system("mv %s %s" % ("redshifts.nc", voidDir+"/"))
     os.system("mv %s %s" % ("indexes.nc", voidDir+"/"))
     os.system("mv %s %s" % ("kdtree_stackvoids.dat", voidDir+"/"))
     os.system("mv %s %s" % ("centers.txt", voidDir+"/"))
+    os.system("mv %s %s" % ("z_centers.txt", voidDir+"/"))
     os.system("mv %s %s" % ("sky_positions.txt", voidDir+"/"))
     os.system("mv %s %s" % ("check.txt", voidDir+"/"))
     os.system("mv %s %s" % ("tracer.txt", voidDir+"/"))
@@ -899,22 +925,24 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
     while badChain:
       Rexpect = (stack.rMin+stack.rMax)/2
       Rtruncate = stack.rMin*3. + 1 # TEST
-      if sample.dataType == "observation":
+      #if sample.dataType == "observation":
-        ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
+      #  ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
-                                      Niter=300000,
+      #                                Niter=300000,
-                                      Nburn=100000,
+      #                                Nburn=100000,
-                                      Rextracut=Rtruncate)
+      #                                Rextracut=Rtruncate)
-      else:
+      #else:
-        ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
+      #  ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
-                                      Niter=300000,
+      #                                Niter=300000,
-                                      Nburn=100000,
+      #                                Nburn=100000,
-                                      Rextracut=Rtruncate)
+      #                                Rextracut=Rtruncate)
-      badChain = (args[0][0] > 0.5 or args[0][1] > stack.rMax or \
+      #badChain = (args[0][0] > 0.5 or args[0][1] > stack.rMax or \
-                  args[0][2] > stack.rMax) and \
+      #            args[0][2] > stack.rMax) and \
-                 (ntries < maxtries)
+      #           (ntries < maxtries)
+      ret,fits,args = vp.compute_inertia(voidDir, stack.rMax)
+      badChain = False
       ntries += 1
 
-    np.save(voidDir+"/chain.npy", ret)
+    #np.save(voidDir+"/chain.npy", ret)
     np.savetxt(voidDir+"/fits.out", fits)
 
     plotTitle = "Sample: "+sample.nickName+\

zobov/voztie.c

@@ -8,6 +8,7 @@ int main(int argc, char *argv[]) {
   FILE *part, *adj, *vol;
   char partfile[200], *suffix, adjfile[200], volfile[200], *outDir;
   float *vols, volstemp;
+  int *cnt_adj;
   
   PARTADJ *adjs;
 
@@ -78,13 +79,17 @@ int main(int argc, char *argv[]) {
   if (mockIndex == -1) mockIndex = np;
   // END PMS
 
+  cnt_adj = (int *)malloc(np*sizeof(int));
+  if (cnt_adj == NULL)
+    printf("Could not allocate cnt_adj.\n");
+
   adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ));
   if (adjs == NULL) printf("Couldn't allocate adjs.\n");
   vols = (float *)malloc(np*sizeof(float));
   if (vols == NULL) printf("Couldn't allocate vols.\n");
   orig = (int *)malloc(nvpmax*sizeof(int));
   if (orig == NULL) printf("Couldn't allocate orig.\n");
-  if ((vols == NULL) || (orig == NULL) || (adjs == NULL)) {
+  if ((cnt_adj==NULL) || (vols == NULL) || (orig == NULL) || (adjs == NULL)) {
     printf("Not enough memory to allocate. Exiting.\n");
     exit(0);
   }
@@ -238,11 +243,24 @@ printf("\n");
 // END OMS
   /* Adjacencies: first the numbers of adjacencies, 
      and the number we're actually going to write per particle */
+
+  // Recount the number of adjacencies after merge
+  for(i=0;i<mockIndex;i++)
+    cnt_adj[i] = 0;
+  for(i=0;i<mockIndex;i++)
+    {
+      for (j=0;j<adjs[i].nadj;j++)
+        {
+          cnt_adj[adjs[i].adj[j]]++;
+          cnt_adj[i]++;
+        }
+    }
+
 // PMS
   for (i=0;i<mockIndex;i++)
   //for (i=0;i<np;i++)
 // END PMS
-    fwrite(&adjs[i].nadj,1,sizeof(int),adj);
+    fwrite(&cnt_adj[i],1,sizeof(int),adj);
     
   /* Now the lists of adjacencies (without double counting) */
   // PMS
@@ -253,9 +271,11 @@ printf("\n");
       nout = 0;
       for (j=0;j<adjs[i].nadj; j++) if (adjs[i].adj[j] > i) nout++;
       fwrite(&nout,1,sizeof(int),adj);      
-      for (j=0;j<adjs[i].nadj; j++) 
+      for (j=0;j<adjs[i].nadj; j++) {
-	if (adjs[i].adj[j] > i) 
+        int id = adjs[i].adj[j];
+	if (id > i) 
 	  fwrite(&(adjs[i].adj[j]),1,sizeof(int),adj);
+      }
     }
 
   fclose(adj);
@@ -275,5 +295,10 @@ printf("\n");
 
   fclose(vol);
 
+  free(vols);
+  free(cnt_adj);
+  free(adjs);
+  free(orig);
+
   return(0);
 }