Continuing development. Boundary handling seems to be stable and working. Now testing void finding.

c_source/prep/prepObservation.cpp

@@ -229,6 +229,7 @@ void flagEdgeGalaxies(prepObservation_info& args ,
   output_data.mask_index = output_data.id_gal.size();
 
   // write a small text file with galaxy position (for diagnostic purposes)
+  // TODO - remove this
   FILE *fp;
   fp = fopen("galaxies.txt", "w");
   for (int i = 0; i < data.size(); i++) {  
@@ -240,7 +241,7 @@ void flagEdgeGalaxies(prepObservation_info& args ,
   }
   fclose(fp);
 
-/* NOTE: temporarily moved to python for quick debugging. Will move back to
+/* NOTE: temporarily moved to python for quick prototyping. Will move back to
          here once it's all sorted
 
   // convert redshift boundaries to covmoving if necessary
@@ -463,8 +464,8 @@ int main(int argc, char **argv)
   fp = fopen("total_particles.txt", "w");
   fprintf(fp, "%d", output_data.pos.size());
   fclose(fp);
-  printf("Done!\n");
   // END PMS - TODO REMOVE THIS
+  printf("Done!\n");
 
   return 0;
 }

c_source/prep/prepObservation.ggo

@@ -19,5 +19,3 @@ option "useComoving" - "Convert to real space using LCDM cosmology" flag off
 option "omegaM" - "Omega Matter for fiducial cosmology" double optional default="0.27"
 
 option "nsideForContour" - "HEALPix NSIDE resolution for figuring out mask contours" int optional default="-1"
-
-option "meanPartSep" - "Estimated mean tracer seperation in h^3 / Mpc^3" double optional default="1"

examples/example_observation/example_observation.py

@@ -30,7 +30,7 @@ continueRun = False
 #   1 : extract redshift slices from data
 #   2 : void extraction using zobov
 #   3 : removal of small voids and voids near the edge
-startCatalogStage = 2
+startCatalogStage = 1
 endCatalogStage   = 3
 
 basePath = os.path.dirname(os.path.abspath(__file__))
@@ -51,7 +51,7 @@ figDir = os.path.join(workDir,"figs","example_observation")
 numZobovThreads = 2
 
 # optimization: number of subdivisions of the volume
-numZobovDivisions = 2
+numZobovDivisions = 1
 
 # Maximum density for merging voids
 #   0 (equivalent to infinitely large value) -> Merge everything (no threshold)

python_source/backend/launchers.py

@@ -76,24 +76,6 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
       sample.maskFile = outputDir + "/constructed_mask.fits"
       figureOutMask(datafile, sample.nsideForContour, sample.maskFile)
 
-    # 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)
-  
-    numTracers = int(open(outputDir+"/mask_index.txt", "r").read())
-    sample.meanPartSep = (1.*numTracers/boxVol/nbar)**(-1/3.)
-
-
-    # flag edge galaxies
-    galFile = outputDir + "galaxies.txt"
-    edgeGalFile = outputDir + "/galaxy_edge_flags.txt"
-    #edgeMaskFile = outputDir + "/mask_edge_map.fits"
-    contourFile = outputDir + "/contour_map.fits"
-    findEdgeGalaxies(galFile, sample.maskFile, edgeGalFile, contourFile,
-                     sample.zBoundary[0], sample.zBoundary[1], sample.omegaM,
-                     useComoving, sample.boundaryWidth, sample.meanPartSep)
-
     if useComoving:
       useComovingFlag = "useComoving"
     else:
@@ -111,12 +93,11 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
       %s
       omegaM %g
       nsideForContour %g
-      meanPartSep %g
       """ % (datafile, sample.maskFile, outputFile,
              outputDir+"/zobov_slice_"+sampleName+".par",
              sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity,
              useComovingFlag, inputParameterFlag, sample.omegaM,
-             sample.nsideForContour, sample.meanPartSep)
+             sample.nsideForContour)
 
     parmFile = os.getcwd()+"/prep_"+sample.fullName+".par"
 
@@ -126,6 +107,24 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
       arg1 = "--configFile=%s" % parmFile
       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)
+     
+        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
+        #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)
+
       if jobSuccessful(logFile, "Done!\n"):
         endTime = time.time()
         walltime = endTime - startTime
@@ -137,6 +136,16 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
     else:
       print("already done!")
 
+
+    # compute mean particle separation
+    # TODO - clean up this duplicate to cover all continueRun options
+    (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.)
+
     if os.access(parmFile, os.F_OK): os.unlink(parmFile)
 
     if os.access("contour_map.fits", os.F_OK):

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):
+                     meanPartSep, outputDir):
 
   if useComoving:
     zmin = comovingDistance(zmin, Om=omegaM)*LIGHT_SPEED
@@ -145,15 +145,50 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
   npix = healpy.nside2npix(nside)
 
   # load in galaxies
-  galPos = np.genfromtxt(galFile)
+  #galPos = np.genfromtxt(galFile)
+  galPos = np.genfromtxt(outputDir+"/galaxies.txt")
+  with open(galFile, 'rb') as File:
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+   Np = np.fromfile(File, dtype=np.int32,count=1)[0]
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+   x = np.fromfile(File, dtype=np.float32,count=Np)
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+   y = np.fromfile(File, dtype=np.float32,count=Np)
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+   z = np.fromfile(File, dtype=np.float32,count=Np)
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+   RA = np.fromfile(File, dtype=np.float32,count=Np)
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+   Dec = np.fromfile(File, dtype=np.float32,count=Np)
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+   redshift = np.fromfile(File, dtype=np.float32,count=Np)
+   chk = np.fromfile(File, dtype=np.int32,count=1)
+
+  #print(galPos.shape)
+  #galPos = np.column_stack((x,y,z))
+  #print(galPos.shape)
+
   flagList = np.zeros(len(galPos[:,0]))
   galTree = scipy.spatial.cKDTree(galPos)
 
   # flag galaxies near mask edges
   # 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)
+  #print(meanPartSep, len(raySteps))
  
   contourPixels = np.nonzero(contourMap)[0]
   #print(contourPixels)
@@ -164,34 +199,42 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
     y = raySteps * vec[1]
     z = raySteps * vec[2]
     ray = np.array((x,y,z)).T
-    #print(ray)
     
     dist, nearest = galTree.query(ray)
     flagList[nearest] = 1
-    #print(nearest)
 
   # flag galaxies near redsfhit boundaries
-  # TODO - save time by only covering portion of sphere with data
+  # TODO - save time by only covering portion of sphere that has data
-  ds = np.sqrt(healpy.nside2pixarea(nside)) / 1000.
-  phi = np.arange(0, 2*np.pi, ds*2)
-  theta = np.arange(0, np.pi, ds)
-  vec = healpy.ang2vec(theta, phi)
 
+  sphereIndices = np.arange(len(contourMap))
+  vec = healpy.pix2vec(nside, sphereIndices)
+  vec = np.asarray(vec)
   maxEdge = zmax * vec
+  maxEdge = maxEdge.T
   dist, nearest = galTree.query(maxEdge)
-  #print(nearest)
-  #print(galPos[nearest])
   flagList[nearest] = 2
 
   minEdge = zmin * vec
+  minEdge = minEdge.T
   dist, nearest = galTree.query(minEdge)
-  #print(nearest)
-  #print(galPos[nearest])
   flagList[nearest] = 3
 
   # output flag information
   np.savetxt(edgeGalFile, flagList, fmt="%d")
 
+  # paint galaxy flags onto healpix map for diagnostics
+  # TODO - drop this when done testing
+  flagMap = np.zeros(len(contourMap))
+  justEdgeRA = RA[flagList == 1]
+  justEdgeDec = Dec[flagList == 1]
+
+  phi, theta = convertAngle(justEdgeRA, justEdgeDec)
+
+  ipix = healpy.ang2pix(nside, theta, phi)
+  np.put(flagMap, ipix, 1)
+
+  healpy.write_map("./flagged_galaxies.fits", flagMap, overwrite=True,
+                   dtype=np.dtype('float64'))
 
 
 #  # output galaxy edge flags

python_source/vide_pipeline/__main__.py

@@ -112,7 +112,7 @@ for sample in dataSampleList:
   # -------------------------------------------------------------------------
   if (startCatalogStage <= 3) and (endCatalogStage >= 3) and not sample.isCombo:
 
-    print("  Pruning void catalogs", "...", end='')
+    print("  Pruning void catalogs...", end='')
     sys.stdout.flush()
         
     logFile = logDir+"/pruneVoids_"+sampleName+".out"