Merge branch 'master' of ssh://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;
-        closestMatchDist = rdist;
-      }
   
-    } 
+      newMatch.matchID = iVoid2;
+      newMatch.commonVol = 0;
+      newMatch.dist = rdist; 
 
-    MATCHPROPS newMatch;
+      if (catalog1.voids[iVoid1].matches.size() < MAX_MATCHES) {
-    newMatch.matchID = closestMatchID;
-    newMatch.commonVol = 1;
-    newMatch.dist = closestMatchDist; 
         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;
+      }
+    }
   } 
-*/
 
+  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;
@@ -205,33 +214,10 @@ int main(int argc, char **argv) {
                 if (rdist <= 0.1*r1 || rdist <= 0.1*r2) match = true;
               }
        
-             
               if (match) {
-                bool alreadyMatched = false;
-                for (iMatch = 0; 
-                     iMatch < catalog1.voids[iVoid1].matches.size(); iMatch++) {
-                   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;
+                  catalog2.part[partID2].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/particleInfo.cpp

@@ -10,6 +10,7 @@ bool loadParticleInfo(ParticleInfo& info,
 		      const std::string& extra_info)
 {
       int numpart;
+      int isObservation;
 
   NcFile f_info(extra_info.c_str());
  
@@ -23,6 +24,7 @@ bool loadParticleInfo(ParticleInfo& info,
   info.ranges[2][0] = f_info.get_att("range_z_min")->as_double(0);
   info.ranges[2][1] = f_info.get_att("range_z_max")->as_double(0);
   info.mask_index = f_info.get_att("mask_index")->as_int(0); //PMS
+  isObservation = f_info.get_att("is_observation")->as_int(0); //PMS
 
   for (int i = 0; i < 3; i++)
     info.length[i] = info.ranges[i][1] - info.ranges[i][0];
@@ -59,6 +61,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 },
@@ -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);

c_tools/stacking/pruneVoids.cpp

@@ -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;
@@ -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] < 1)
+        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] < 2)
+        voids[iVoid].barycenter[2] = boxLen[2] - voids[iVoid].barycenter[2];
+    }
+
 
     // compute central density
     centralRad = voids[iVoid].radius/args_info.centralRadFrac_arg;
@@ -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");
@@ -531,6 +566,10 @@ 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]);
@@ -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,

python_tools/pipeline_source/prepareCatalogs.in.py

@@ -79,7 +79,7 @@ startCatalogStage = 1
 endCatalogStage   = 4
                
 startAPStage = 1
-endAPStage = 7
+endAPStage = 1
 
 ZOBOV_PATH = "@CMAKE_BINARY_DIR@/zobov/"
 CTOOLS_PATH = "@CMAKE_BINARY_DIR@/c_tools/"
@@ -124,7 +124,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 +135,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

@@ -314,7 +314,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 +359,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
@@ -600,13 +601,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 +905,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+\