From ad2d3722cc06b73eef848fe8d3a450bf59379cd5 Mon Sep 17 00:00:00 2001
From: "Paul M. Sutter" <pmsutter@gmail.com>
Date: Fri, 17 Jan 2025 21:53:25 +0800
Subject: [PATCH] Added double-checking of flagged galaxies to ensure survey
 boundaries are respected. Python components of launchers now add their status
 to log files. Rearranged some warnings for broken tessellations due to guard
 point encounters, since it's not always a bad thing.

---
 c_source/jozov2/jozov2.cpp            |  14 +++-
 c_source/zobov/voz1b1.c               |   2 +
 python_source/backend/launchers.py    | 107 ++++++++++++++++++++------
 python_source/backend/surveyTools.py  |  17 +++-
 python_source/voidUtil/catalogUtil.py |  18 +++--
 5 files changed, 122 insertions(+), 36 deletions(-)

diff --git a/c_source/jozov2/jozov2.cpp b/c_source/jozov2/jozov2.cpp
index ea14344..6e9fb78 100644
--- a/c_source/jozov2/jozov2.cpp
+++ b/c_source/jozov2/jozov2.cpp
@@ -86,13 +86,19 @@ int main(int argc,char **argv)
     }
 
   /* Check that we got all the pairs */
+  bool printedError = false;
   for (int i = 0; i < np; i++)
     {
       if (p[i].ncnt != p[i].nadj) {
-	cout
-	  << format("We didn't get all of %d's adj's; %d != %d.") 
-	  % i % p[i].ncnt % p[i].nadj
-	  << endl;
+
+      if (!printedError) {
+        cout << "WARNING! Not all particle adjancies are adding up. This is likely because some particles connected to guard particles. Ensure that your tesselation is valid." << endl;
+        printedError = true;
+      }
+	//cout
+	//  << format("We didn't get all of %d's adj's; %d != %d.") 
+	//  % i % p[i].ncnt % p[i].nadj
+	//  << endl;
 	p[i].nadj = p[i].ncnt;
       }
     }
diff --git a/c_source/zobov/voz1b1.c b/c_source/zobov/voz1b1.c
index 3daee57..23a78ab 100644
--- a/c_source/zobov/voz1b1.c
+++ b/c_source/zobov/voz1b1.c
@@ -310,12 +310,14 @@ int main(int argc, char *argv[]) {
   }
 
 // PMS - reset number of adjancies to not include links to border guards
+/*
   for (i=0; i<nvp; i++) {
     int nActual = 0;
     for (j = 0; j < adjs[i].nadj; j++)
       if (adjs[i].adj[j] < nvp) nActual++;
     adjs[i].nadj = nActual;
   }
+*/
 
   /* Get the adjacencies back to their original values */
 
diff --git a/python_source/backend/launchers.py b/python_source/backend/launchers.py
index dbdf53a..26ccdad 100644
--- a/python_source/backend/launchers.py
+++ b/python_source/backend/launchers.py
@@ -44,6 +44,10 @@ import time
 NetCDFFile = Dataset
 ncFloat = 'f8' # Double precision
 
+class Bunch:
+  def __init__(self, **kwds):
+    self.__dict__.update(kwds)
+
 #LIGHT_SPEED = 299792.458
 
 # -----------------------------------------------------------------------------
@@ -104,27 +108,37 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
     if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")):
       with open(parmFile, mode="wt") as f:
         f.write(conf)
+
+      log = open(logFile, 'w')
       arg1 = "--configFile=%s" % parmFile
-      with open(logFile, 'wt') as log:
-        subprocess.call([binPath, arg1], stdout=log, stderr=log)
+      #with open(logFile, 'wt') as log:
+      subprocess.call([binPath, arg1], stdout=log, stderr=log)
 
-        # compute mean particle separation
-        (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0],
-              sample.zRange[1], sample.zRange[0], sample.zRange[1], "all",
-              sample.omegaM, useComoving=useComoving)
+      # compute mean particle separation
+      log.write("\nFlagging edge galaxies.\n")
+      (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0],
+            sample.zRange[1], sample.zRange[0], sample.zRange[1], "all",
+            sample.omegaM, useComoving=useComoving)
      
-        numTracers = int(open(outputDir+"/mask_index.txt", "r").read())
-        sample.meanPartSep = (1.*numTracers/boxVol/nbar)**(-1/3.)
+      numTracers = int(open(outputDir+"/mask_index.txt", "r").read())
+      sample.meanPartSep = (1.*numTracers/boxVol/nbar)**(-1/3.)
 
-        # flag edge galaxies with python routine for now
-        galFile = outputDir+"/zobov_slice_"+sampleName
+      # flag edge galaxies with python routine for now
+      galFile = outputDir+"/zobov_slice_"+sampleName
+      if os.access("contour_map.fits", os.F_OK):
+        os.system("mv %s %s" % ("contour_map.fits", outputDir))
         os.system("mv %s %s" % ("galaxies.txt", outputDir))
-        #galFile = outputDir + "galaxies.txt"
-        edgeGalFile = outputDir + "/galaxy_edge_flags.txt"
-        contourFile = outputDir + "/contour_map.fits"
-        findEdgeGalaxies(galFile, sample.maskFile, edgeGalFile, contourFile,
-                        sample.zBoundary[0], sample.zBoundary[1], sample.omegaM,
-                        useComoving, sample.boundaryWidth, sample.meanPartSep, outputDir)
+      #galFile = outputDir + "galaxies.txt"
+      edgeGalFile = outputDir + "/galaxy_edge_flags.txt"
+      contourFile = outputDir + "/contour_map.fits"
+      findEdgeGalaxies(galFile, sample.maskFile, edgeGalFile, contourFile,
+                       sample.zBoundary[0], sample.zBoundary[1], sample.omegaM,
+                       useComoving, sample.boundaryWidth, sample.meanPartSep, 
+                       outputDir, log)
+
+      log.write("Done!\n")
+      log.close()
+
 
       if jobSuccessful(logFile, "Done!\n"):
         endTime = time.time()
@@ -402,7 +416,7 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
 
     # re-weight the volumes based on selection function
     if sample.dataType == "observation" and \
-       sample.selFunFile != None:
+      sample.selFunFile != None:
 
       # load volumes
       volFile = outputDir+"/vol_"+sampleName+".dat"
@@ -417,6 +431,7 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
         Np = np.fromfile(File, dtype=np.int32,count=1)
         chk = np.fromfile(File, dtype=np.int32,count=1)
 
+        # skipping through the file to get to the correct location
         # x
         chk = np.fromfile(File, dtype=np.int32,count=1)
         redshifts = np.fromfile(File, dtype=np.float32,count=Np)
@@ -468,29 +483,75 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
     else:
       volFileToUse = outputDir+"/vol_"+sampleName+".dat"
 
-    # re-weight the volumes of any edge galaxies to prevent watershed
+    # double-check to make sure galaxies were flagged correctly.
+    #   if any non-edge galaxies have bad tesselations, flag them
+    # also, re-weight the volumes of any edge galaxies to prevent watershed
     #   from spilling outside of survey region
     if sample.dataType == "observation":
 
-      # read in the appropriate volume file
-      with open(volFileToUse, mode="rb") as File:
-        numPartTot = np.fromfile(File, dtype=np.int32,count=1)[0]
-        vols = np.fromfile(File, dtype=np.float32, count=numPartTot)
-      
+      log = open(logFile, 'a')
+      log.write("\nDouble-checking all edge flags.\n\n")
+
       # read in the edge flag information
       edgeFile = outputDir+"/galaxy_edge_flags.txt"
       edgeFlags = np.loadtxt(edgeFile, dtype=np.int32)
 
+      # load in and count adjacenies
+      log.write("Loading adjancies...\n")
+      part = []
+      adjFile = outputDir+"adj_"+sampleName+".dat"
+      with open(adjFile, mode="rb") as File:
+        numPart = np.fromfile(File, dtype=np.int32,count=1)[0]
+
+        # i'm tired and can't think of a better way to do this right now
+        for p in range(numPart):
+          part.append(Bunch(nadjs = 0, adjs = []))
+
+        nadjPerPart = np.fromfile(File, dtype=np.int32,count=numPart)
+        for p in range(numPart):
+          nin = np.fromfile(File, dtype=np.int32, count=1)[0]
+          for n in range(nin):
+            pAdj = np.fromfile(File, dtype=np.int32, count=1)[0]
+            if (p < pAdj):
+              part[p].adjs.append(pAdj)
+              part[pAdj].adjs.append(p)
+ 
+      # if the number of expected adjancies != actual adjacenies, this galaxy
+      #   connected to a guard point on the boundary. This is fine for edge
+      #   galaxies, but occasionally interior galaxies slip through the
+      #   boundary handling, so we'll mark those as an edge
+      for p in range(numPart):
+        nHave = len(part[p].adjs)
+        nExpected = nadjPerPart[p]
+        #print("Working on %d: %d vs %d" % (p, nHave, nExpected))
+        if nHave != nExpected and edgeFlags[p] == 0:
+          log.write("Found an unflagged galaxy: %d. Correcting.\n" % (p))
+          edgeFlags[p] = 4
+
+      # re-write the new edge flag information
+      np.savetxt(edgeFile, edgeFlags, fmt="%d")
+      log.write("\nDone double-checking. New flags saved.\n")
+
+      # read in the appropriate volume file
+      log.write("Re-weighting edge galaxy volumes...\n")
+      with open(volFileToUse, mode="rb") as File:
+        numPartTot = np.fromfile(File, dtype=np.int32,count=1)[0]
+        vols = np.fromfile(File, dtype=np.float32, count=numPartTot)
+      
       # set edge galaxy volumes to nearly 0 (implying very high density)
       vols[ edgeFlags>0 ] = 1.e-10
 
       volFile = outputDir+"/vol_weighted_"+sampleName+".dat"
+      log.write("Saving new volume file to %s.\n" % (volFile))
       with open(volFile, mode='wb') as File:
         numPartTot.astype(np.int32).tofile(File)
         vols.astype(np.float32).tofile(File)
 
       sample.hasWeightedVolumes = True
       volFileToUse = outputDir+"/vol_weighted_"+sampleName+".dat"
+      log.write("Done re-weighting.\n\n")
+
+      log.close()
     else:
       volFileToUse = outputDir+"/vol_"+sampleName+".dat"
 
diff --git a/python_source/backend/surveyTools.py b/python_source/backend/surveyTools.py
index 0fd1d8b..e331105 100644
--- a/python_source/backend/surveyTools.py
+++ b/python_source/backend/surveyTools.py
@@ -131,7 +131,7 @@ def figureOutMask(galFile, nside, outMaskFile):
 # figures out which galaxies live on a mask or redshift edge
 def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
                      zmin, zmax, omegaM, useComoving, boundaryWidth,
-                     meanPartSep, outputDir):
+                     meanPartSep, outputDir, log):
 
   if useComoving:
     zmin = comovingDistance(zmin, Om=omegaM)*LIGHT_SPEED
@@ -140,12 +140,13 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
     zmin *= LIGHT_SPEED
     zmax *= LIGHT_SPEED
 
+  log.write("  Reading contour map...\n")
   contourMap = healpy.read_map(contourFile)
   nside = healpy.get_nside(contourMap)
   npix = healpy.nside2npix(nside)
 
   # load in galaxies
-  #galPos = np.genfromtxt(galFile)
+  log.write("  Loading galaxies...\n")
   # TODO - WHY IS THIS FASTER THAN np.column_stack???
   galPos = np.genfromtxt(outputDir+"/galaxies.txt")
   with open(galFile, 'rb') as File:
@@ -181,6 +182,7 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
   #galPos = np.column_stack((x,y,z))
   #print(galPos.shape)
 
+  log.write("  Building k-d tree...\n")
   flagList = np.zeros(len(galPos[:,0]))
   galTree = scipy.spatial.cKDTree(galPos)
 
@@ -188,10 +190,14 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
   # using the "ray marching" algorithm: follow rays along lines of sight
   #   of all mask edges, flagging nearest neighbor galaxies as we go
  
-  raySteps = np.arange(zmin, zmax, meanPartSep)
+  log.write("  Begin ray-marching...\n")
+  rayDist = meanPartSep
+  raySteps = np.arange(zmin, zmax, rayDist)
+  log.write("  Will take %d steps per ray with %.2e distance between steps\n" % (len(raySteps), rayDist))
   #print(meanPartSep, len(raySteps))
  
   contourPixels = np.nonzero(contourMap)[0]
+  log.write("  We have %d rays to work with\n" % (len(contourPixels)))
   #print(contourPixels)
   for pixel in contourPixels:
     #print("Working with pixel %d" % pixel)
@@ -207,6 +213,7 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
   # flag galaxies near redsfhit boundaries
   # TODO - save time by only covering portion of sphere that has data
 
+  log.write("  Flagging galaxies near redshift bounds...\n")
   sphereIndices = np.arange(len(contourMap))
   vec = healpy.pix2vec(nside, sphereIndices)
   vec = np.asarray(vec)
@@ -221,10 +228,12 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
   flagList[nearest] = 3
 
   # output flag information
+  log.write("  Saving galaxy flags to file...\n")
   np.savetxt(edgeGalFile, flagList, fmt="%d")
 
   # paint galaxy flags onto healpix map for diagnostics
   # TODO - drop this when done testing
+  log.write("  Saving diagnostic map to file...\n")
   flagMap = np.zeros(len(contourMap))
   justEdgeRA = RA[flagList == 1]
   justEdgeDec = Dec[flagList == 1]
@@ -234,7 +243,7 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
   ipix = healpy.ang2pix(nside, theta, phi)
   np.put(flagMap, ipix, 1)
 
-  healpy.write_map("./flagged_galaxies.fits", flagMap, overwrite=True,
+  healpy.write_map(outputDir+"/flagged_galaxies.fits", flagMap, overwrite=True,
                    dtype=np.dtype('float64'))
 
 
diff --git a/python_source/voidUtil/catalogUtil.py b/python_source/voidUtil/catalogUtil.py
index 528b55c..ef14dfa 100644
--- a/python_source/voidUtil/catalogUtil.py
+++ b/python_source/voidUtil/catalogUtil.py
@@ -448,11 +448,12 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True,
       print(" Loading edge flags...")
       edgeFlagFile = sampleDir+"/galaxy_edge_flags.txt"
       if os.path.isfile(edgeFlagFile):
-        edgeFlags = np.genfromtxt(edgeFlagFile, dtype=np.int32, skip_header=1)
+        edgeFlags = np.loadtxt(edgeFlagFile, dtype=np.int32)
         for iEdge in range(len(edgeFlags)):
           catalog.part[iEdge].edgeFlag = edgeFlags[iEdge]
       else:
         print("   Edge file not found!")
+        catalog.part[:].edgeFlag = 0
                     
     print(" Loading volumes...")
     if sample.hasWeightedVolumes:
@@ -476,22 +477,29 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True,
       # but the file only stores one half of each pair, so we need to match
       for p in range(numPart):
         nin = np.fromfile(File, dtype=np.int32, count=1)[0]
-        
         for n in range(nin):
           pAdj = np.fromfile(File, dtype=np.int32, count=1)[0]
-
           if (p < pAdj):
             catalog.part[p].adjs.append(pAdj)
             catalog.part[pAdj].adjs.append(p)
 
 
     print(" Sanity checking adjacenies...")
+    nBad = 0
+    ## TEST - writing bad galaxies to file
+    #badGalFile = open("badgal.txt", "w")
     for p in range(numPart):
       catalog.part[p].nadjs = len(catalog.part[p].adjs)
       nHave = len(catalog.part[p].adjs)
       nExpected = nadjPerPart[p]
-      if (nHave != nExpected):
-        print("   Error for particle %d: Have %d adj, expected %d" % (p, nHave, nExpected)) 
+      if (nHave != nExpected): print("   Error for particle %d: Have %d adj, expected %d (flag: %d)" % (p, nHave, nExpected, catalog.part[p].edgeFlag)) 
+        #if (catalog.part[p].edgeFlag == 0):
+        #  badGalFile.write("%e %e %e\n" % (catalog.part[p].x, catalog.part[p].y, catalog.part[p].z))
+        #  nBad += 1
+   
+    #print("NUMBER OF BAD GALAXIES: %d" % nBad)
+    #badGalFile.close()
+
 
     print(" Loading zone-void membership info...")
     zoneFile = sampleDir+"/voidZone_"+sample.fullName+".dat"