more work on merger tree analysis code

c_tools/analysis/voidOverlap.cpp

@@ -38,15 +38,22 @@ typedef struct voidZoneStruct {
   int *zoneIDs;
 } VOID2ZONE;
 
+typedef struct matchProps {
+  int matchID;
+  float commonVol;
+  float dist;
+} MATCHPROPS;
+
 typedef struct voidStruct {
   float vol, coreDens, zoneVol, densCon, voidProb, radius;
   int voidID, numPart, numZones, coreParticle, zoneNumPart;
   float maxRadius, nearestMock, centralDen, redshift, redshiftInMpc;
   float nearestEdge;
   float barycenter[3];
-  std::vector<int> matches;
+  std::vector<MATCHPROPS> matches;
-  int biggestMatchID, numMatches;
+  int numMatches;
-  float biggestMatchVol;
+  int numBigMatches;
+  float radiusMpc;
 } VOID;
 
 typedef struct catalog {
@@ -63,16 +70,24 @@ void loadCatalog(const char *partFile, const char *volFile,
                  const char *infoFile, const char *barycenterFile,
                  const char *zonePartFile, CATALOG& catalog);
 
-float getDist(CATALOG& catalog1, CATALOG& catalog2, int& iVoid1, int& iVoid2);
+float getDist(CATALOG& catalog1, CATALOG& catalog2, int& iVoid1, int& iVoid2,
+               bool periodicX, bool periodicY, bool operiodicZ);
+
+void sortMatches(std::vector<MATCHPROPS>& matches);
+
+// ----------------------------------------------------------------------------
 
 int main(int argc, char **argv) {
 
-  int p1, p2, iZ1, iZ2, iVoid1, iVoid2, iVoid, zoneID1, zoneID2;
+  int p1, p2, iZ1, iZ2, iVoid1, iVoid2, iVoid, zoneID1, zoneID2, iMatch;
   int partID1, partID2;
   int voidID1, voidID2;
   bool periodicX=false, periodicY=false, periodicZ=false, match;
   float dist[3], rdist, r1, r2;
   FILE *fp;
+  int closestMatchID;
+  float closestMatchDist;
+  float commonVolRatio;
 
   CATALOG catalog1, catalog2;
 
@@ -120,6 +135,30 @@ int main(int argc, char **argv) {
     }
 
   printf(" Determining overlap...\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;
+      }
+
+    } 
+    
+    MATCHPROPS newMatch;
+    newMatch.matchID = closestMatchID;
+    newMatch.commonVol = 1;
+    newMatch.dist = closestMatchDist; 
+    catalog1.voids[iVoid1].matches.push_back(newMatch);
+  } 
+*/
+
   //for (iVoid1 = 0; iVoid1 < 1; iVoid1++) {
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
     printf("  Working on void %d of %d...\n", iVoid1, catalog1.numVoids); 
@@ -130,20 +169,21 @@ int main(int argc, char **argv) {
       match = false;
 
       for (iZ1 = 0; iZ1 < catalog1.void2Zones[voidID1].numZones; iZ1++) {
-        if (match) break;
+  //      if (match) break;
         zoneID1 = catalog1.void2Zones[voidID1].zoneIDs[iZ1];
 
         for (p1 = 0; p1 < catalog1.zones2Parts[zoneID1].numPart; p1++) {
-          if (match) break;
+ //         if (match) break;
           partID1 = catalog1.zones2Parts[zoneID1].partIDs[p1];
 
           for (iZ2 = 0; iZ2 < catalog2.void2Zones[voidID2].numZones; iZ2++) {
-            if (match) break;
+ //           if (match) break;
             zoneID2 = catalog2.void2Zones[voidID2].zoneIDs[iZ2];
 
             for (p2 = 0; p2 < catalog2.zones2Parts[zoneID2].numPart; p2++) {
               partID2 = catalog2.zones2Parts[zoneID2].partIDs[p2];     
 
+              match = false;
               if (args.useID_flag) {
                 if (catalog1.part[partID1].uniqueID == 
                     catalog2.part[partID2].uniqueID) match = true;
@@ -169,11 +209,36 @@ int main(int argc, char **argv) {
                 if (rdist <= 0.1*r1 || rdist <= 0.1*r2) match = true;
              }
        
+             
               if (match) {
-                //printf("MATCHED TO %d\n", iVoid2);
+                bool alreadyMatched = false;
-                catalog1.voids[iVoid1].matches.push_back(iVoid2);
+                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;
+                 } 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);
+                    //printf("MATCHED TO %d\n", iVoid2);
+                    //catalog1.voids[iVoid1].matches.push_back(iVoid2);
+                } // end if match
+//                break;
+              } 
 
             } // end p2
           } // end iZ2
@@ -182,53 +247,79 @@ int main(int argc, char **argv) {
     } // end iVoid2
   } // end iVoid1
 
+  for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
+    sortMatches(catalog1.voids[iVoid1].matches);
+  }
 
   // output properties of matches
   fp = fopen(args.outfile_arg, "w");
+  fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, relative dist, num matches, num significant matches\n");
   //for (iVoid1 = 0; iVoid1 < 1; iVoid1++) {
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
-    for (iVoid = 0; iVoid < catalog1.voids[iVoid1].matches.size(); iVoid++) {
-      iVoid2 = catalog1.voids[iVoid1].matches[iVoid];
-       rdist = getDist(catalog1, catalog2, iVoid1, iVoid2);
 
-      if (rdist/catalog1.voids[iVoid1].radius > 1 || 
+    // find closest match
-          rdist/catalog2.voids[iVoid2].radius > 1) {
+    closestMatchDist = 0.;
-        catalog1.voids[iVoid].matches[iVoid] = -1;
+    for (iMatch = 0; iMatch < catalog1.voids[iVoid1].matches.size(); iMatch++) {
-        continue;
+      rdist = catalog1.voids[iVoid1].matches[iMatch].dist;
-      }
+      //iVoid2 = catalog1.voids[iVoid1].matches[iVoid];
+       //rdist = getDist(catalog1, catalog2, iVoid1, iVoid2,
+       //                periodicX, periodicY, periodicZ);
 
-      if (catalog2.voids[iVoid2].vol > 
+      commonVolRatio = catalog1.voids[iVoid1].matches[iMatch].commonVol / 
-          catalog1.voids[iVoid1].biggestMatchVol) {
+                     //  catalog1.voids[iVoid1].numPart;
-        catalog1.voids[iVoid1].biggestMatchID = iVoid2;
+                     catalog1.voids[iVoid1].vol;
-        catalog1.voids[iVoid1].biggestMatchVol = catalog2.voids[iVoid2].vol;
+ 
+      if (commonVolRatio > 0.1) catalog1.voids[iVoid1].numBigMatches++;
+
+      if (commonVolRatio > closestMatchDist) {
+        closestMatchID = iMatch;
+        closestMatchDist = commonVolRatio;
       }
-      catalog1.voids[iVoid1].numMatches++;
+      //if (rdist < closestMatchDist) {
+      //  closestMatchID = iMatch;
+      //  closestMatchDist = rdist;
+      //}
+
       //printf("CENTER %d %d %e %e %e\n", iVoid1, iVoid2, rdist, rdist/catalog1.voids[iVoid1].radius, rdist/catalog2.voids[iVoid2].radius);
-    } 
       
+    } 
+    catalog1.voids[iVoid1].numMatches = catalog1.voids[iVoid1].matches.size();
+
+    // actually print out
     int voidID = catalog1.voids[iVoid1].voidID;
     if (catalog1.voids[iVoid1].numMatches > 0) {
-      iVoid2 = catalog1.voids[iVoid1].biggestMatchID;
+      iVoid2 = catalog1.voids[iVoid1].matches[closestMatchID].matchID;
-      float volRatio = catalog2.voids[iVoid2].vol / 
+      //iVoid2 = catalog1.voids[iVoid1].biggestMatchID;
+      float rRatio = catalog2.voids[iVoid2].radius / 
+                     catalog1.voids[iVoid1].radius;
+      float commonVolRatio = 
+                     catalog1.voids[iVoid1].matches[closestMatchID].commonVol / 
+                     //catalog1.voids[iVoid1].numPart;
                      catalog1.voids[iVoid1].vol;
-      rdist = getDist(catalog1, catalog2, iVoid1, iVoid2);
+      rdist = getDist(catalog1, catalog2, iVoid1, iVoid2,
+                      periodicX, periodicY, periodicZ);
       rdist /= catalog1.voids[iVoid1].radius;
 
-      fprintf(fp, "%d %e %e %d\n", voidID, 
+      fprintf(fp, "%d %.2f %.2f %.2f %.2f %d %d\n", voidID, 
-                                   volRatio, 
+                                   catalog1.voids[iVoid1].radiusMpc,
+                                   rRatio, 
+                                   commonVolRatio,
                                    rdist, 
-                                   catalog1.voids[iVoid1].numMatches);
+                                   catalog1.voids[iVoid1].numMatches,
+                                   catalog1.voids[iVoid1].numBigMatches);
 
     } else {
-      fprintf(fp, "%d 0.0 0.0 0\n", voidID);
+      fprintf(fp, "%d %f 0.0 0.0 0.0 0 0\n", voidID, 
+                                           catalog1.voids[iVoid1].radius);
     }
   }
   fclose(fp);
 
-  printf("\n Well, bye.\n");
+  printf("\nDone!\n");
   return 0;
 } // end main
 
+// ----------------------------------------------------------------------------
 // ----------------------------------------------------------------------------
 void loadCatalog(const char *partFile, const char *volFile, 
                  const char *voidFile, const char *zoneFile,
@@ -245,7 +336,6 @@ void loadCatalog(const char *partFile, const char *volFile,
   float ranges[3][2];
 
   printf("Loading info...\n");
-printf("HELLO %s\n", infoFile);
   NcFile f_info(infoFile);
   ranges[0][0] = f_info.get_att("range_x_min")->as_double(0);
   ranges[0][1] = f_info.get_att("range_x_max")->as_double(0);
@@ -269,6 +359,8 @@ printf("HELLO %s\n", infoFile);
   catalog.part.resize(numPartTot);
   catalog.numPartTot = numPartTot;
 
+  volNorm = numPartTot/(catalog.boxLen[0]*catalog.boxLen[1]*catalog.boxLen[2]);
+
   temp = (float *) malloc(numPartTot * sizeof(float));
   temp2 = (long *) malloc(numPartTot * sizeof(long));
 
@@ -348,9 +440,14 @@ printf("HELLO %s\n", infoFile);
     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].biggestMatchID = -1;
-    catalog.voids[i].biggestMatchVol = 0;
+    //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;
+  
+    catalog.voids[i].radiusMpc = pow(voidVol/volNorm*3./4./M_PI, 1./3.);
 
     i++;
   }
@@ -418,7 +515,8 @@ printf("HELLO %s\n", infoFile);
 } // end loadCatalog
 
 // ----------------------------------------------------------------------------
-float getDist(CATALOG& catalog1, CATALOG& catalog2, int& iVoid1, int& iVoid2) {
+float getDist(CATALOG& catalog1, CATALOG& catalog2, int& iVoid1, int& iVoid2,
+              bool periodicX, bool periodicY, bool periodicZ) {
 
   float rdist, dist[3];
 
@@ -429,8 +527,38 @@ float getDist(CATALOG& catalog1, CATALOG& catalog2, int& iVoid1, int& iVoid2) {
   dist[2] = fabs(catalog1.voids[iVoid1].barycenter[2] - 
                  catalog2.voids[iVoid2].barycenter[2]);
 
+  if (periodicX) dist[0] = fmin(dist[0], 1.0 - dist[0]);
+  if (periodicY) dist[1] = fmin(dist[1], 1.0 - dist[1]);
+  if (periodicZ) dist[2] = fmin(dist[2], 1.0 - dist[2]);
+
   rdist = sqrt(dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2]);
 
   return rdist;
 
 } // end getDist
+
+
+// ----------------------------------------------------------------------------
+void sortMatches(std::vector<MATCHPROPS>& matches) {
+
+  //printf("SORTING %d\n", matches.size());
+  MATCHPROPS tempMatch;
+  bool swapped;
+
+  if (matches.size() == 1) return;
+
+  swapped = false;
+  while (!swapped) {
+    swapped = false;
+    for (int iMatch = 0; iMatch < matches.size() - 1; iMatch++) {
+      if (matches[iMatch].commonVol > matches[iMatch+1].commonVol) {
+        tempMatch = matches[iMatch];
+        matches[iMatch+1] = matches[iMatch];
+        matches[iMatch] = tempMatch;
+        swapped = true;
+      }
+    }
+  }
+
+  return;
+} // end sortMatches

crossCompare/analysis/mergerTree.py

@@ -62,16 +62,20 @@ for (iSample, sampleDir) in enumerate(sampleDirList):
   if iSample == 1:
     os.system("cp %s %s" % (stepOutputFileName, outFileName))
   else:
-    outFile = fp.open("temp.out", "w")
+    outFile = open("temp.out", "w")
-    stepOutFile1 = fp.open(outFileName, "r")
+    stepOutFile1 = open(outFileName, "r")
-    stepOutFile2 = fp.open(stepOutputFileName, "r")
+    stepOutFile2 = open(stepOutputFileName, "r")
 
     for line1 in stepOutFile1:
+      line1 = line1.rstrip()
       line2 = stepOutFile2.readline()
       outFile.write(line1+" "+line2)
 
+    os.system("cp %s %s" % ("temp.out", outFileName))
     outFile.close()
     stepOutFile1.close()
     stepOutFile2.close()
 
-if os.access("mergerTree.log", os.F_OK): os.unlink("mergerTree.log")
+#if os.access("mergerTree.log", os.F_OK): os.unlink("mergerTree.log")
+if os.access("temp.out", os.F_OK): os.unlink("temp.out")
+#if os.access("thisStep.out", os.F_OK): os.unlink("thisStep.out")

crossCompare/plotting/datasetsToPlot.py

@@ -5,14 +5,14 @@ workDir = "/home/psutter2/workspace/Voids/"
 figDir  = "./figs"
 
 sampleDirList = [ 
-                  "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/",
+#                  "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/",
-                  "multidark/md_halos_min1.393e12_pv/sample_md_halos_min1.393e12_pv_z0.56_d00/",
+#                  "multidark/md_halos_min1.393e12_pv/sample_md_halos_min1.393e12_pv_z0.56_d00/",
-                  "random/ran_ss0.000175/sample_ran_ss0.000175_z0.56_d00/",
+#                  "random/ran_ss0.000175/sample_ran_ss0.000175_z0.56_d00/",
-                  "random/ran_ss0.1/sample_ran_ss0.1_z0.56_d00/",
+#                  "random/ran_ss0.1/sample_ran_ss0.1_z0.56_d00/",
-                  "multidark/md_hod_dr9mid_pv/sample_md_hod_dr9mid_pv_z0.56_d00/",
+#                  "multidark/md_hod_dr9mid_pv/sample_md_hod_dr9mid_pv_z0.56_d00/",
-                  "multidark/md_ss0.000175_pv/sample_md_ss0.000175_pv_z0.56_d00/",
+#                  "multidark/md_ss0.000175_pv/sample_md_ss0.000175_pv_z0.56_d00/",
                   "sdss_dr9/sample_lss.dr9cmassmid.dat/",
-                  "lanl/masked/masked_lanl_hod_dr9mid_pv/sample_masked_lanl_hod_dr9mid_pv_z0.0/" ]
+                  "lanl/masked/masked_lanl_hod_dr9mid_pv/sample_masked_lanl_hod_dr9mid_pv_z0.5/" ]
 
 dataPortion = "central"