setting the stage for on-the-fly mask and edge calculation

c_source/pruning/pruneVoids.cpp

@@ -506,7 +506,7 @@ int main(int argc, char **argv) {
   for (iVoid = 0; iVoid < numVoids; iVoid++) {
 
     voidID = voids[iVoid].voidID;
-    printf("  DOING %d (of %d) %d %d %f\n", iVoid, numVoids, voidID, 
+    printf("  DOING %d (of %d) %d %d %f\n", iVoid+1, numVoids, voidID, 
                                            voids[iVoid].numPart, 
                                            voids[iVoid].radius);
 

examples/example_observation/example_observation.py

@@ -82,8 +82,13 @@ newSample = Sample(
                    # assume sample is volume-limited?
                    volumeLimited = True,
 
-                   # HEALpix mask file
-                   maskFile = inputDataDir+"/example_observation_mask.fits",
+                   # HEALpix mask file - set to None to auto-compute
+                   maskFile = None,
+                   #maskFile = inputDataDir+"/example_observation_mask.fits",
+
+                   # if maskFile blank, desired resolution for HEALpix
+                   # mask mapping, otherwise pulled from maskFile
+                   nsideForMask = 128,
 
                    # radial selection function (if not volume limited)
                    selFunFile = None,
@@ -100,12 +105,12 @@ newSample = Sample(
 
                    # density of mock particles in cubic Mpc/h
                    # (make this as high as you can afford)
-                   fakeDensity = 0.0000001,
+                   fakeDensity = 0.05,
 
                    # if true, convert to comoving space using LCDM cosmology
                    useComoving = True,
 
-                   # cosmology
+                   # cosmology assuming flat universe
                    omegaM = 0.3,
 
                    )

python_source/backend/__init__.py

@@ -20,3 +20,4 @@
 from .classes import *
 from .launchers import *
 from .cosmologyTools import *
+from .surveyTools import *

python_source/backend/cosmologyTools.py

@@ -22,11 +22,10 @@
 
 import numpy as np
 import scipy.integrate as integrate
-import healpy as healpy
 import os
 from backend import *
 
-__all__=['expansion', 'angularDiameter', 'aveExpansion', 'getSurveyProps']
+__all__=['expansion', 'angularDiameter', 'aveExpansion']
 
 # returns 1/E(z) for the given cosmology
 def expansion(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
@@ -51,62 +50,3 @@ def aveExpansion(zStart, zEnd, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
   ave = integrate.quad(expansion, zStart, zEnd, args=(Om, Ot, w0, wa))[0]
   ave = (zEnd-zStart)/ave
   return ave
-
-# -----------------------------------------------------------------------------
-# returns the volume and galaxy density for a given redshit slice
-def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectionFuncFile=None, useComoving=False):
-
-  LIGHT_SPEED = 299792.458
-
-  mask = healpy.read_map(maskFile)
-  area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi
-
-  if useComoving:
-    zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=0.27)
-    zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=0.27)
-    selFunMin = LIGHT_SPEED/100.*angularDiameter(selFunMin, Om=0.27)
-    selFunMax = LIGHT_SPEED/100.*angularDiameter(selFunMax, Om=0.27)
-  else:
-    zmin = zmin * 3000
-    zmax = zmax * 3000
-    selFunMin *= 3000
-    selFunMax *= 3000
-  volume = area * (zmax**3  - zmin**3) / 3
-
-  if selectionFuncFile == None:
-    nbar = 1.0
-
-  elif not os.access(selectionFuncFile, os.F_OK):
-    print("   Warning, selection function file %s not found, using default of uniform selection." % selectionFuncFile)
-    nbar = 1.0
-
-  else:
-    selfunc = np.genfromtxt(selectionFuncFile)
-    selfunc = np.array(selfunc)
-    selfunc[:,0] = selfunc[:,0]/100.
-    selfuncUnity = selfunc
-    selfuncUnity[:,1] = 1.0
-    selfuncMin = selfunc[0,0]
-    selfuncMax = selfunc[-1,0]
-    selfuncDx = selfunc[1,0] - selfunc[0,0]
-    selfuncN = np.size(selfunc[:,0])
-  
-    selFunMin = max(selFunMin, selfuncMin)
-    selFunMax = min(selFunMax, selfuncMax)
-  
-  
-    def f(z): return selfunc[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 
-    def fTotal(z): return selfuncUnity[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 
-  
-    zrange = np.linspace(selFunMin, selFunMax)
-  
-    nbar = scipy.integrate.quad(f, selFunMin, selFunMax)
-    nbar = nbar[0]
-  
-    ntotal = scipy.integrate.quad(fTotal, 0.0, max(selfuncUnity[:,0]))
-    ntotal = ntotal[0]
-  
-    nbar = ntotal / area / nbar
-
-
-  return (volume, nbar)

python_source/backend/launchers.py

@@ -36,6 +36,7 @@ from   pylab import figure
 from   netCDF4 import Dataset
 from backend.classes import *
 from backend.cosmologyTools import *
+from backend.surveyTools import *
 import pickle
 import scipy.interpolate as interpolate
 
@@ -531,7 +532,7 @@ def launchPrune(sample, binPath,
     cmd += " --outputDir=" + zobovDir
     cmd += " --sampleName=" + str(sampleName)
     log = open(logFile, 'w')
-    log.write(f"Command is {cmd}\n")
+    #log.write(f"Command is {cmd}\n")
     subprocess.call(cmd, stdout=log, stderr=log, shell=True)
     log.close()
 

python_source/backend/surveyTools.py

@@ -0,0 +1,129 @@
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/apTools/chi2/cosmologyTools.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+# 
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+# a suite of functions to compute expansion rates, angular diameter 
+# distances, and expected void stretching
+
+import numpy as np
+import healpy as healpy
+import os
+from backend import *
+
+__all__=['getSurveyProps', 'getNside', 'figureOutMask', 'findEdgeGalaxies']
+
+# -----------------------------------------------------------------------------
+# returns the volume and galaxy density for a given redshit slice
+def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectionFuncFile=None, useComoving=False):
+
+  LIGHT_SPEED = 299792.458
+
+  mask = healpy.read_map(maskFile)
+  area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi
+
+  if useComoving:
+    zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=0.27)
+    zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=0.27)
+    selFunMin = LIGHT_SPEED/100.*angularDiameter(selFunMin, Om=0.27)
+    selFunMax = LIGHT_SPEED/100.*angularDiameter(selFunMax, Om=0.27)
+  else:
+    zmin = zmin * 3000
+    zmax = zmax * 3000
+    selFunMin *= 3000
+    selFunMax *= 3000
+  volume = area * (zmax**3  - zmin**3) / 3
+
+  if selectionFuncFile == None:
+    nbar = 1.0
+
+  elif not os.access(selectionFuncFile, os.F_OK):
+    print("   Warning, selection function file %s not found, using default of uniform selection." % selectionFuncFile)
+    nbar = 1.0
+
+  else:
+    selfunc = np.genfromtxt(selectionFuncFile)
+    selfunc = np.array(selfunc)
+    selfunc[:,0] = selfunc[:,0]/100.
+    selfuncUnity = selfunc
+    selfuncUnity[:,1] = 1.0
+    selfuncMin = selfunc[0,0]
+    selfuncMax = selfunc[-1,0]
+    selfuncDx = selfunc[1,0] - selfunc[0,0]
+    selfuncN = np.size(selfunc[:,0])
+  
+    selFunMin = max(selFunMin, selfuncMin)
+    selFunMax = min(selFunMax, selfuncMax)
+  
+  
+    def f(z): return selfunc[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 
+    def fTotal(z): return selfuncUnity[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 
+  
+    zrange = np.linspace(selFunMin, selFunMax)
+  
+    nbar = scipy.integrate.quad(f, selFunMin, selFunMax)
+    nbar = nbar[0]
+  
+    ntotal = scipy.integrate.quad(fTotal, 0.0, max(selfuncUnity[:,0]))
+    ntotal = ntotal[0]
+  
+    nbar = ntotal / area / nbar
+
+
+  return (volume, nbar)
+
+
+# -----------------------------------------------------------------------------
+# returns the nside resolution from the given maskfile
+def getNside(maskFile):
+  nside = 1.0
+
+  return nside
+
+
+# -----------------------------------------------------------------------------
+# computes the mask from a given datafile and writes it to a file
+def figureOutMask(galFile, nside, outMaskFile):
+
+  npix = healpy.nside2npix(nside)
+  mask = np.zeros((npix))
+
+for line in open(galFile):
+  line = line.split()
+  RA  = np.float(line[3])
+  Dec = np.float(line[4])
+  z   = np.float(line[5])
+
+  phi   = np.pi/180.*RA
+  theta = Dec*np.pi/180.
+  theta = np.pi/2. - Dec*np.pi/180.
+  pos = np.zeros((3))
+
+  pix = healpy.ang2pix(nside, theta, phi)
+  mask[pix] = 1.
+ 
+  healpy.write_map(outMaskFile, mask)
+
+  return mask
+
+# -----------------------------------------------------------------------------
+# figures out which galaxies live on a mask edge, and also writes the edge
+#   map to an auxillary file
+def findEdgeGalaxies(galFile, maskFile):
+
+
+  return