fixed and calrified zobov vs survey normalization schemes

c_source/pruning/pruneVoids.cpp

@@ -271,6 +271,9 @@ int main(int argc, char **argv) {
   part = (PART *) malloc(numPartTot * sizeof(PART));
   temp = (float *) malloc(numPartTot * sizeof(float));
 
+  // this is the volume normalization understood by zobov,
+  // which is all particles (including mock boundary galaxies used in old
+  // versions of VIDE) divided by the total volume of the bounding box
   volNorm = numPartTot/(boxLen[0]*boxLen[1]*boxLen[2]);
   printf(" VOL NORM = %f\n", volNorm);
 

python_source/backend/surveyTools.py

@@ -236,7 +236,7 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
 
   # paint galaxy flags onto healpix map for diagnostics
   # TODO - drop this when done testing
-  log.write("  Saving diagnostic map to file...\n")
+  log.write("  Saving diagnostic maps to file...\n")
   flagMap = np.zeros(len(contourMap))
   justEdgeRA = RA[flagList == 1]
   justEdgeDec = Dec[flagList == 1]
@@ -244,13 +244,28 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
   phi, theta = convertAngle(justEdgeRA, justEdgeDec)
 
   ipix = healpy.ang2pix(nside, theta, phi)
-  np.put(flagMap, ipix, 1)
+  for i in ipix:
+    flagMap[i] += 1
+  #np.put(flagMap, ipix, 1)
 
   healpy.write_map(outputDir+"/flagged_galaxies.fits", flagMap, 
                    overwrite=True,
                    dtype=np.dtype('float64'))
 
 
+  allGalMap = np.zeros(len(contourMap))
+  phi, theta = convertAngle(RA, Dec)
+  
+  ipix = healpy.ang2pix(nside, theta, phi)
+  for i in ipix:
+    allGalMap[i] += 1
+  #np.put(allGalMap, ipix, 1)
+
+  healpy.write_map(outputDir+"/all_galaxies.fits", allGalMap, 
+                   overwrite=True,
+                   dtype=np.dtype('float64'))
+
+
 #  # output galaxy edge flags
 #  edgeFile = open(edgeGalFile, "w")
 #

python_source/voidUtil/catalogUtil.py

@@ -49,6 +49,12 @@ def loadPart(sampleDir):
     ranges[2][0] = getattr(File, 'range_z_min')
     ranges[2][1] = getattr(File, 'range_z_max')
     isObservation = getattr(File, 'is_observation')
+
+    # old verison of VIDe includes the boundary tracers in the file
+    nGal = 0
+    if hasattr(File, 'mask_index'):
+      nGal = getattr(File, 'mask_index')
+
     File.close()
     mul = np.zeros((3))
     mul[:] = ranges[:,1] - ranges[:,0]
@@ -99,31 +105,22 @@ def loadPart(sampleDir):
       uniqueID = np.fromfile(File, dtype=np.int64,count=Np)
       chk = np.fromfile(File, dtype=np.int32,count=1)
 
+    # if it's an old catalog, trim it
+    if nGal > 0:
+      x = x[0:nGal]
+      y = y[0:nGal]
+      z = z[0:nGal]
+      RA = RA[0:nGal]
+      Dec = Dec[0:nGal]
+      redshift = redshift[0:nGal]
+      uniqueID = uniqueID[0:nGal]
+
     partData = np.column_stack((x,y,z))
     extraData = np.column_stack((RA,Dec,redshift,uniqueID))
 
     boxLen = mul
-
-    # TODO - should use built-in volNorm function instead
-
-    #if isObservation == 1:
-    #  # look for the mask file
-    #  if os.access(sample.maskFile, os.F_OK):
-    #    maskFile = sample.maskFile
-    #  else:
-    #    maskFile = sampleDir+"/"+os.path.basename(sample.maskFile)
-    #    print "Using maskfile found in:", maskFile
-    #  props = getSurveyProps(maskFile, sample.zBoundary[0],
-    #                            sample.zBoundary[1],
-    #                            sample.zBoundary[0],
-    #                            sample.zBoundary[1], "all",
-    #                            selectionFuncFile=sample.selFunFile,
-    #                            useComoving=sample.useComoving)
-    #  boxVol = props[0]
-    #  volNorm = maskIndex/boxVol
-    #else:
     boxVol = np.prod(boxLen)
-    volNorm = len(x)/boxVol
+    volNorm = Np/boxVol # this is the zobov normalization
 
     isObservationData = isObservation == 1
 
@@ -148,19 +145,25 @@ def getVolNorm(sampleDir):
     mul = np.zeros((3))
     mul[:] = ranges[:,1] - ranges[:,0]
 
-    partFile = sampleDir+"/zobov_slice_"+sample.fullName
-    with open(partFile, mode="rb") as File:
-      chk = np.fromfile(File, dtype=np.int32,count=1)
-      Np = np.fromfile(File, dtype=np.int32,count=1)[0]
+    # getting the total number of tracers is an adventure depending
+    # on which verison of VIDE we're pulling up
+    maskIndexFile = sampleDir+"/mask_index.txt"
+    totalPartFile = sampleDir+"/total_particles.txt"
+
+    if os.path.exists(maskIndexFile):
+      nTot = float(open(totalPartFile, "r").read())
+      nGal = float(open(maskIndexFile, "r").read())
+    else:
+      nTot = float(open(totalPartFile, "r").read())
+      nGal = nTot
 
     boxLen = mul
 
     boxVol = np.prod(boxLen)
-    volNorm = Np/boxVol
-    volNormObs = volNorm
+    volNormZobov = nTot/boxVol
+    volNormObs = volNormZobov
 
     if isObservation == 1:
-      # look for the mask file
       if os.access(sample.maskFile, os.F_OK):
         maskFile = sample.maskFile
       else:
@@ -174,10 +177,10 @@ def getVolNorm(sampleDir):
                                 sample.omegaM,
                                 selectionFuncFile=sample.selFunFile,
                                 useComoving=sample.useComoving)
-      boxVol = props[0]
-      volNormObs = Np/boxVol
+      surveyVol = props[0]
+      volNormObs = nGal/surveyVol
 
-    return volNorm, volNormObs
+    return volNormZobov, volNormObs
 
 
 # -----------------------------------------------------------------------------
@@ -278,7 +281,7 @@ class Catalog:
   numVoids = 0
   numPartTot = 0
   numZonesTot = 0
-  volNorm = 0 # normalization used by zobov across entire volumne
+  volNormZobov = 0 # normalization used by zobov across entire volumne
   volNormObs = 0 # normalization based on average density of survey volume
   boxLen = np.zeros((3))
   ranges = np.zeros((3,2))
@@ -324,8 +327,8 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True,
   catalog.ranges = ranges
   File.close()
 
-  volNorm, volNormObs = getVolNorm(sampleDir)
-  catalog.volNorm = volNorm
+  volNormZobov, volNormObs = getVolNorm(sampleDir)
+  catalog.volNormZobov = volNormZobov
   catalog.volNormObs = volNormObs
 
   if untrimmed: