some bug fixes to the a-p analysis scripts; fixed major bug in prunevoids; updated catalog release scripts

c_tools/mock/generateMock.cpp

@@ -188,12 +188,13 @@ SimuData *doLoadMultidark(const char *multidarkname)
                                 &outd->Pos[0][i], &outd->Pos[1][i], 
                                 &outd->Pos[2][i], &outd->Vel[2][i]);
 
-    outd->uniqueID[i] = 1.0;
+    outd->uniqueID[i] = 1.0 * outd->Id[i];
-    //outd->uniqueID[i] = 1.0 * outd->Id[i];
 
+    if (i < 10) printf("TEST %d %d\n", i, outd->Id[i]);
     if (outd->Id[i] == -99 && 
         outd->Pos[0][i] == -99 && outd->Pos[1][i] == -99 && 
         outd->Pos[2][i] == -99 && outd->Vel[2][i] == -99) {
+   printf("FOUND END\n");
       break;
     } else {
       actualNumPart++;
@@ -396,8 +397,6 @@ void generateOutput(SimuData *data, int axis,
   f.beginCheckpoint();
   for (uint32_t i = 0; i < data->NumPart; i++)
     {
-      //printf("HELLO %d %d\n", i, data->Id[i]);
-      //f.writeReal32(data->Id[i]);
       f.writeReal32(data->uniqueID[i]);
     }
   f.endCheckpoint();

c_tools/stacking/pruneVoids.cpp

@@ -441,9 +441,6 @@ int main(int argc, char **argv) {
     }
 
     voids[iVoid].nearestEdge = nearestEdge;
-    //if (nearestEdge < voids[iVoid].nearestMock) {
-    //  voids[iVoid].nearestMock = nearestEdge;
-    //}
   } // iVoid
 
   printf(" Picking winners and losers...\n");
@@ -453,7 +450,7 @@ int main(int argc, char **argv) {
 
   for (iVoid = 0; iVoid < numVoids; iVoid++) {
     if (voids[iVoid].densCon < 1.5) {
-      //voids[iVoid].accepted = 0;
+//      voids[iVoid].accepted = -4;
     }
  
     if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) {

pipeline/applyMaskToMock.py

@@ -1,8 +1,6 @@
 #!/usr/bin/env python
 
-# prepares input catalogs based on multidark simulations
+# applies a mask to a given dataset
-#   (borrows heavily from generateMock, but doesn't hold much in memory)
-# also creates necessary analyzeVoids input files
 
 import numpy as np
 import os
@@ -32,6 +30,18 @@ parser.add_argument('--scripts', dest='script', action='store_const',
 parser.add_argument('--parmFile', dest='parmFile',
                    default="",
                    help='path to parameter file')
+parser.add_argument('--subsamples', dest='subsample', action='store_const',
+                   const=True, default=False,
+                   help='write subsamples')
+parser.add_argument('--halos', dest='halos', action='store_const',
+                   const=True, default=False,
+                   help='write halos')
+parser.add_argument('--hod', dest='hod', action='store_const',
+                   const=True, default=False,
+                   help='write hod')
+parser.add_argument('--all', dest='all', action='store_const',
+                   const=True, default=False,
+                   help='write everything')
 args = parser.parse_args()
 
 
@@ -57,9 +67,9 @@ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1):
 
 #------------------------------------------------------------------------------
 # for given dataset parameters, outputs a script for use with analyzeVoids
-def writeScript(setName, dataFileNameBase,
+def writeObservationScript(setName, dataFileNameBase, maskFileName,
-                scriptDir, catalogDir, fileNums, redshifts, numSubvolumes,
+                scriptDir, catalogDir, fileNums, redshifts, 
-                numSlices, useVel, lbox, minRadius, omegaM, subsample=1.0, 
+                useVel, minRadius, omegaM,
                 suffix=".dat"):
 
 
@@ -88,7 +98,7 @@ ranSeed = 101010
 useLCDM = False
 bias = 1.16
 
-dataPortions = ["central"]
+dataPortions = ["central", "all"]
 dataSampleList = []
 """)
 
@@ -112,25 +122,17 @@ numZobovThreads = {numZobovThreads}
 
   sampleInfo = """
 newSample = Sample(dataFile = "{dataFile}",
-                   dataFormat = "{dataFormat}",
-                   dataUnit = {dataUnit},
                    fullName = "{sampleName}",
                    nickName = "{sampleName}",
-                   dataType = "simulation",
+                   dataType = "observation",
+                   maskFile = "{maskFileName}",
                    zBoundary = ({zMin}, {zMax}),
                    zRange    = ({zMin}, {zMax}),
-                   zBoundaryMpc = ({zMinMpc}, {zMaxMpc}),
-                   omegaM    = {omegaM},
                    minVoidRadius = {minRadius},
                    includeInHubble = True,
                    partOfCombo = False,
                    isCombo = False,
-                   boxLen = {boxLen},
+                   usePecVel = {usePecVel})
-                   usePecVel = {usePecVel},
-                   numSubvolumes = {numSubvolumes},
-                   mySubvolume = "{mySubvolume}",
-                   useLightCone = {useLightCone},
-                   subsample = {subsample})
 dataSampleList.append(newSample)
 newSample.addStack({zMin}, {zMax}, {minRadius}  , {minRadius}+2, True, False)
 newSample.addStack({zMin}, {zMax}, {minRadius}  , {minRadius}+4, True, False)
@@ -142,180 +144,249 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False)
   for (iFile, redshift) in enumerate(redshifts):
     fileNum = fileNums[iFile]
     zBox = float(redshift)
-    Om = float(omegaM)
-    zBoxMpc = LIGHT_SPEED/100.*vp.angularDiameter(zBox, Om=Om)
-    boxMaxMpc = zBoxMpc + lbox
    
-    # converter from redshift to comoving distance   
+    zMin = zBox
-    zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000)
+    zMax = zBox + 0.1 
-    zVsDX = np.zeros(len(zVsDY))
+    zMin, zMax = getMaskedZRange(lbox, zBox, zMin, zMax, omegaM)
-    for i in xrange(len(zVsDY)):
+    dataFileName = dataFileNameBase + fileNum + suffix
-      zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
 
-    if useLightCone:
+    sampleName = getSampleName(setName, redshift, useVel,
-      boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \
+                                iSlice=-1, iVol=-1)
-                  LIGHT_SPEED*lbox, zVsDX, zVsDY)-zBox
-      dzSafe = 0.03
-    else:
-      boxWidthZ = lbox*100./LIGHT_SPEED
-      dzSafe = 0.0
 
-    for iSlice in xrange(numSlices):
+    scriptFile.write(sampleInfo.format(dataFile=dataFileName,
-      sliceMin = zBox + dzSafe + iSlice*(boxWidthZ-dzSafe)/numSlices
-      sliceMax = zBox + dzSafe + (iSlice+1)*(boxWidthZ-dzSafe)/numSlices
-
-      sliceMinMpc = sliceMin*LIGHT_SPEED/100.
-      sliceMaxMpc = sliceMax*LIGHT_SPEED/100.
-
-      sliceMin = "%0.2f" % sliceMin
-      sliceMax = "%0.2f" % sliceMax
-      sliceMinMpc = "%0.1f" % sliceMinMpc
-      sliceMaxMpc = "%0.1f" % sliceMaxMpc
-
-      dataFileName = dataFileNameBase + fileNum + suffix
-
-      for iX in xrange(numSubvolumes):
-        for iY in xrange(numSubvolumes):
-
-          mySubvolume = "%d%d" % (iX, iY)
-
-          sampleName = getSampleName(setName, sliceMin, useVel,
-                                     iSlice=iSlice, iVol=mySubvolume)
-
-          scriptFile.write(sampleInfo.format(dataFile=dataFileName,
                                          dataFormat=dataFormat,
-                                         dataUnit=dataUnit,
                                          sampleName=sampleName,
-                                         zMin=sliceMin,
+                                         maskFileName=maskFileName,
-                                         zMax=sliceMax,
+                                         zMin=zMin,
-                                         zMinMpc=sliceMinMpc,
+                                         zMax=zMax,
-                                         zMaxMpc=sliceMaxMpc,
-                                         omegaM=Om,
-                                         boxLen=lbox,
                                          usePecVel=useVel,
-                                         minRadius=minRadius,
+                                         minRadius=minRadius))
-                                         numSubvolumes=numSubvolumes,
-                                         mySubvolume=mySubvolume,
-                                         useLightCone=useLightCone,
-                                         subsample=subsample))
 
   scriptFile.close()
   return
 
 
 #------------------------------------------------------------------------------
-#------------------------------------------------------------------------------
-if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir)
-if not os.access(catalogDir, os.F_OK): os.mkdir(catalogDir)
-
-# -----------------------------------------------------------------------------
-# now the SDSS HOD
-parFileText = """
-% cosmology
-OMEGA_M {omegaM}
-HUBBLE  {hubble}
-OMEGA_B 0.0469
-SIGMA_8 0.82
-SPECTRAL_INDX 0.95
-ITRANS 5
-REDSHIFT {redshift}
-
-% halo definition
-%DELTA_HALO 200
-DELTA_HALO 740.74 % 200/Om_m
-M_max 1.00E+16
-
-% fit function types
-pdfs 11
-pdfc 2
-EXCLUSION 4
-
-% hod parameters 
-M_min {Mmin}
-GALAXY_DENSITY 0.0111134 % computed automatically if M_min set, use for sanity
-M1 {M1}   
-sigma_logM {sigma_logM}
-alpha {alpha}
-M_cut {Mcut}
-
-% simulation info
-real_space_xi 1
-HOD 1
-populate_sim 1
-HaloFile {haloFile}
-RESOLUTION {numPartPerSide}
-BOX_SIZE   {boxSize}
-
-% output
-root_filename hod
-               """
-
-print " Doing DR9 HOD"
-
-    # these parameters come from Manera et al 2012, eq. 26
-parFileName = "./hod.par"
-parFile = open(parFileName, 'w')
-haloFile = catalogDir+haloFileBase+fileNums[iRedshift]
-parFile.write(parFileText.format(omegaM=omegaM,
-                                     hubble=hubble,
-                                     redshift=redshift,
-                                     Mmin=1.23e13,
-                                     M1=1.e14,
-                                     sigma_logM=0.596,
-                                     alpha=1.0127,
-                                     Mcut=1.19399e13,
-                                     haloFile=haloFile,
-                                     numPartPerSide=numPart**(1/3.),
-                                     boxSize=lbox))
-parFile.close()
-
-os.system(hodPath+" "+parFileName+">& /dev/null")
-
 # now place these particles on a lightcone, restrict redshift range, apply mask
-mask = hp.read_map(maskFile)
+def applyMask(inFileName, outFileName, maskFileName, lbox, zBox, zMin, zMax, omegaM):
-nside = hp.get_nside(mask)
+  mask = hp.read_map(maskFileName)
+  nside = hp.get_nside(mask)
 
-inFile  = open('hod.mock', 'r') 
+  inFile  = open(inFileName, 'r') 
-outFile = open(catalogDir+"/mock.out"))
+  outFile = open(outFileName, 'w')
 
-zBox = float(redshiftRange[0])
+  Om = float(omegaM)
-Om = float(omegaM)
 
-# converter from redshift to comoving distance   
+  # converter from redshift to comoving distance   
-zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000)
+  zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000)
-zVsDX = np.zeros(len(zVsDY))
+  zVsDX = np.zeros(len(zVsDY))
-for i in xrange(len(zVsDY)):
+  for i in xrange(len(zVsDY)):
-  zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
+    zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
-
-for line in inFile:
-  line = line.split(',')
-  x  = float(line[0]) - lbox/2.
-  y  = float(line[1]) - lbox/2.
-  z  = float(line[2]) - lbox/2.
-  vz = float(line[5])
 
   zBoxInMpc = vp.angularDiameter(zBox, Om=Om)
 
-  redshift = np.sqrt(x*x + y*y + z*z) 
+  for (iLine,line) in enumerate(inFile):
-  redshift = np.interp(zBoxInMpc+100./LIGHT_SPEED*redshift, zVsDX, zVsDY)
+    if iLine < 5:
+      print >> outFile, line.rstrip()
+      continue
 
-  if redshift < redshiftRange[0] or redshift > redshiftRange[1]: continue
+    line = line.split(' ')
+    uniqueID = int(line[0])
+    x  = float(line[1]) - lbox/2.
+    y  = float(line[2]) - lbox/2.
+    z  = float(line[3]) - lbox/2.
+    vz = float(line[4])
 
-  RA = np.atan((y-lbox/2.)/(x-lbox/2.)) * 100/np.pi + 180.
+    # TODO if usePecvel
-  Dec = np.asin((z-lboc/2.)/(redshift*LIGHT_SPEED/100.)) * 180/np.pi
 
-  phi = np.pi/180. * RA
+    redshift = np.sqrt(x*x + y*y + z*z) 
-  theta = np.pi/2. - Dec*np.pi/180.
+    redshift = np.interp(zBoxInMpc+100./LIGHT_SPEED*redshift, zVsDX, zVsDY)
-  pos = np.zeros((3))
 
-  pix = hp.ang2pix(nside, theta, phi)
+    if redshift < zMin or redshift > zMax: continue
-  if mask[pix] <= 0: continue   
 
-  print >> outFile, RA, Dec, redshift*LIGHT_SPEED, 0., x, y, z
+    vec = np.array((x,y,z))
+    theta, phi = hp.vec2ang(vec)
+    theta = theta[0]
+    phi = phi[0]
+    RA = phi*180./np.pi
+    Dec = 180./np.pi*(np.pi/2.-theta)
 
-inFile.close()
+    pix = hp.vec2pix(nside, x, y, z)
-outFile.close()
+    if mask[pix] <= 0.2: continue   
   
-os.system("rm ./hod.*")
+    print >> outFile, RA, Dec, redshift*LIGHT_SPEED, uniqueID, x, y, z
+ 
+  inFile.close()
+  outFile.close()
+
+#------------------------------------------------------------------------------
+def getMaskedZRange(lbox, zBox, zMin, zMax, omegaM):
+
+  if zMin < zBox: zMin = zBox 
+  Om = float(omegaM)
+
+  # converter from redshift to comoving distance   
+  zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000)
+  zVsDX = np.zeros(len(zVsDY))
+  for i in xrange(len(zVsDY)):
+    zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
+
+  zBoxInMpc = vp.angularDiameter(zBox, Om=Om)
+
+  boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \
+                  LIGHT_SPEED*lbox/2., zVsDX, zVsDY)-zBox
+
+  print "RANGE", zBox, zBox+boxWidthZ
+  if zMax > zBox+boxWidthZ: zMax = zBox+boxWidthZ
+   
+  return zMin, zMax
+
+#------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
+
+inCatalogDir = catalogDir
+
+prefix = "masked_" + prefix
+voidOutputDir += "masked/"
+figDir += "masked/"
+logDir += "masked/"
+scriptDir += "masked/"
+catalogDir += "masked/"
+dataType = "observation"
+
+if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir)
+if not os.access(catalogDir, os.F_OK): os.mkdir(catalogDir)
+
+#------------------------------------------------------------------------------
+# first the directly downsampled runs
+# Note: ss0.002   ~ SDSS DR7 dim2
+#       ss0.0004 ~ SDSS DR9 mid 
+baseResolution = float(numPart)/lbox/lbox/lbox # particles/Mpc^3
+for thisSubSample in subSamples:
+
+  keepFraction = float(thisSubSample) / baseResolution
+  maxKeep = keepFraction * numPart
+  minRadius = int(np.ceil(lbox/maxKeep**(1./3)))
+
+  if args.script or args.all:
+    print " Doing subsample", thisSubSample, " scripts"
+    setName = prefix+"ss"+str(thisSubSample)
+    writeObservationScript(setName, "masked_md.ss"+str(thisSubSample)+"_z", 
+                  maskFileName, 
+                  scriptDir, catalogDir, fileNums, 
+                  redshifts, False, minRadius, omegaM)
+    writeObservationScript(setName, "masked_md.ss"+str(thisSubSample)+"_z", 
+                  maskFileName, 
+                  scriptDir, catalogDir, fileNums, 
+                  redshifts, True, minRadius, omegaM)
+  
+  if args.subsample or args.all:
+    print " Doing subsample", thisSubSample
+
+    for (iRedshift, redshift) in enumerate(redshifts):
+      print "   redshift", redshift
+  
+      zMin = float(redshift)
+      zMax = zMin + 0.1 
+      zMin, zMax = getMaskedZRange(lbox, float(redshift), zMin, zMax, omegaM)
+
+      if dataFormat == "multidark" or dataFormat == "random":
+        inFileName = inCatalogDir+"/md.ss"+str(thisSubSample)+"_z"+redshift+".dat"
+        outFileName = catalogDir+"/masked_md.ss"+str(thisSubSample)+"_z"+redshift+".dat"
+        
+        applyMask(inFileName, outFileName, maskFileName, lbox, float(redshift), 
+                  zMin, zMax, omegaM)
 
 
+# -----------------------------------------------------------------------------
+# now halos
+if (args.script or args.all) and dataFormat == "multidark":
+  print " Doing halo scripts"
+
+  for minHaloMass in minHaloMasses:
+  
+    if dataFormat == "multidark":
+      setName = prefix+"halos_min"+str(minHaloMass)
+      writeScript(setName, "md.halos_min"+str(minHaloMass)+"_z", 
+                scriptDir, catalogDir, fileNums, 
+                redshifts, 
+                numSubvolumes, numSlices, False, lbox, minRadius, omegaM)
+      writeScript(setName, "md.halos_min"+str(minHaloMass)+"_z", 
+                scriptDir, catalogDir, fileNums, 
+                redshifts, 
+                numSubvolumes, numSlices, True, lbox, minRadius, omegaM)
+
+if args.halos or args.all: 
+  print " Doing halos"
+
+  for minHaloMass in minHaloMasses:
+    print "  min halo mass = ", minHaloMass
+
+    for (iRedshift, redshift) in enumerate(redshifts):
+      print "   z = ", redshift
+
+      dataFile = catalogDir+haloFileBase+fileNums[iRedshift]
+      inFile = open(dataFile, 'r')
+      numPart = 0
+      for line in inFile: 
+        line = line.split(',')
+        if minHaloMass == "none" or float(line[6]) > minHaloMass:
+          numPart += 1
+      inFile.close()
+
+      sampleName = "md.halos_min"+str(minHaloMass)+"_z"+redshift
+      outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
+
+# -----------------------------------------------------------------------------
+# now the SDSS HOD
+
+if (args.script or args.all) and dataFormat == "multidark":
+  print " Doing DR7 HOD scripts"
+  if dataFormat == "multidark":
+    setName = prefix+"hod_dr72dim2"
+    writeScript(setName, "md.hod_dr72dim2_z",
+              scriptDir, catalogDir, fileNums, redshifts, 
+              numSubvolumes, numSlices, False, lbox, 5, omegaM)
+    writeScript(setName, "md.hod_dr72dim2_z",
+              scriptDir, catalogDir, fileNums, redshifts, 
+              numSubvolumes, numSlices, True, lbox, 5, omegaM)
+
+if args.hod or args.all:
+  print " Doing DR7 HOD"
+  for (iRedshift, redshift) in enumerate(redshifts):
+    print "  z = ", redshift
+
+    sampleName = getSampleName("md.hod_dr72dim2", redshift, False)
+    outFileName = catalogDir+"/"+sampleName+".dat"
+
+# -----------------------------------------------------------------------------
+# now the BOSS HOD
+if (args.script or args.all) and dataFormat == "multidark":
+  print " Doing DR9 HOD scripts"
+  if dataFormat == "multidark":
+    setName = prefix+"hod_dr9mid"
+    writeScript(setName, "md.hod_dr9mid_z",
+              scriptDir, catalogDir, fileNums, redshifts, 
+               numSubvolumes, numSlices, False, lbox, 15, omegaM)
+    writeScript(setName, "md.hod_dr9mid_z",
+              scriptDir, catalogDir, fileNums, redshifts, 
+               numSubvolumes, numSlices, True, lbox, 15, omegaM)
+
+if args.hod or args.all:
+  print " Doing DR9 HOD"
+  for (iRedshift, redshift) in enumerate(redshifts):
+    print "  z = ", redshift
+
+    outFileName = catalogDir+"/"+sampleName+".dat"
+
+    if dataFormat == "multidark" or dataFormat == "random":
+      sampleName = getSampleName("md.hod_dr9mid", redshift, False)
+      inFileName = inCatalogDir+"/"+sampleName+".dat"
+      outFileName = catalogDir+"/masked_"+sampleName+".dat"
+      
+      zMin = float(redshift)
+      zMax = zMin + 0.1 
+      zMin, zMax = getMaskedZRange(lbox, float(redshift), zMin, zMax, omegaM)
+        
+      applyMask(inFileName, outFileName, maskFileName, lbox, float(redshift), 
+                zMin, zMax, omegaM)

pipeline/datasets/multidark.py

@@ -49,7 +49,7 @@ numSubvolumes = 1
 prefix = "md_"
 
 # list of desired subsamples - these are in unts of h Mpc^-3!
-#subSamples = [ 1.0 ]
+#subSamples = [0.0004]
 subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002))
 
 # common filename of halo files
@@ -62,7 +62,7 @@ minHaloMasses = (("none", 2e12, 1.23e13))
 # adjust these two parameters given the memory contraints on your system:
 #   numZobovDivisions: how many sub-volumes per dimension will zobov process
 #   numZobovThreads: how many sub-volumes to process at once?   
-numZobovDivisions = 4
+numZobovDivisions = 2
 numZobovThreads = 2
 
 # simulation information
@@ -71,6 +71,9 @@ lbox = 1000 # Mpc/h
 omegaM = 0.27
 hubble = 0.70
 
+# the mask file which is used by applyMaskToMock
+maskFileName = os.getenv("HOME")+"/workspace/Voids/catalogs/boss/boss_mask_final.fits"
+
 # END CONFIGURATION
 # -----------------------------------------------------------------------------
 # -----------------------------------------------------------------------------

pipeline/generateCatalog.py

@@ -9,13 +9,6 @@ import pickle
 
 # ------------------------------------------------------------------------------
 
-def my_import(name):
-    mod = __import__(name)
-    components = name.split('.')
-    for comp in components[1:]:
-        mod = getattr(mod, comp)
-    return mod
-
 if (len(sys.argv) == 0):
   print "Usage: ./anylyzeVoids.py parameter_file.py"
   exit(-1)

pipeline/prepareCatalogs.py

@@ -11,15 +11,6 @@ import void_python_tools as vp
 import argparse
 import imp
 
-# ------------------------------------------------------------------------------
-
-def my_import(name):
-    mod = __import__(name)
-    components = name.split('.')
-    for comp in components[1:]:
-        mod = getattr(mod, comp)
-    return mod
-
 # -----------------------------------------------------------------------------
 
 LIGHT_SPEED = 299792.458
@@ -83,7 +74,7 @@ def writeScript(setName, dataFileNameBase,
 import os
 from void_python_tools.backend.classes import *
 
-continueRun = True # set to True to enable restarting aborted jobs
+continueRun = False # set to True to enable restarting aborted jobs
 startCatalogStage = 1
 endCatalogStage   = 4
                
@@ -133,6 +124,7 @@ newSample = Sample(dataFile = "{dataFile}",
                    zBoundaryMpc = ({zMinMpc}, {zMaxMpc}),
                    omegaM    = {omegaM},
                    minVoidRadius = {minRadius},
+                   profileBinSize = 1.0,
                    includeInHubble = True,
                    partOfCombo = False,
                    isCombo = False,
@@ -143,12 +135,12 @@ newSample = Sample(dataFile = "{dataFile}",
                    useLightCone = {useLightCone},
                    subsample = {subsample})
 dataSampleList.append(newSample)
-newSample.addStack({zMin}, {zMax}, {minRadius}  , {minRadius}+2, True, False)
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+2, True, False)
-newSample.addStack({zMin}, {zMax}, {minRadius}  , {minRadius}+4, True, False)
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+4, True, False)
-newSample.addStack({zMin}, {zMax}, {minRadius}+2, {minRadius}+6, True, False)
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}+2, 2*{minRadius}+6, True, False)
-newSample.addStack({zMin}, {zMax}, {minRadius}+6, {minRadius}+10, True, False)
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}+6, 2*{minRadius}+10, True, False)
-newSample.addStack({zMin}, {zMax}, {minRadius}+10, {minRadius}+18, True, False)
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}+10, 2*{minRadius}+18, True, False)
-newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False)
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, False)
                """
   for (iFile, redshift) in enumerate(redshifts):
     fileNum = fileNums[iFile]

plotting/datasetsToPlot.py

@@ -4,9 +4,9 @@
 workDir = "/home/psutter2/workspace/Voids/"
 figDir  = "./figs"
 
-sampleDirList = [ "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/",
+sampleDirList = [ 
-                  #"multidark/md_ss01.0_pv/sample_md_ss1.0_pv_z0.56_d00/",
+                  "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/",
-                  "multidark/md_halos_min1.23e13_pv/sample_md_halos_min1.23e13_pv_z0.56_d00/",
+                  "multidark/md_halos_min1.23e+13_pv/sample_md_halos_min1.23e+13_pv_z0.56_d00/",
                   "random/ran_ss0.0004/sample_ran_ss0.0004_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/",

plotting/plotCompareDensCon.py

@@ -14,8 +14,6 @@ import argparse
 
 # ------------------------------------------------------------------------------
 
-from datasetsToPlot import *
-
 plotNameBase = "compdenscon"
 
 obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing?
@@ -24,10 +22,20 @@ parser = argparse.ArgumentParser(description='Plot.')
 parser.add_argument('--show', dest='showPlot', action='store_const',
                    const=True, default=False,
                    help='display the plot (default: just write eps)')
+parser.add_argument('--parmFile', dest='parmFile', default='datasetsToPlot.py',
+                    help='path to parameter file')
 args = parser.parse_args()
 
 # ------------------------------------------------------------------------------
 
+filename = args.parmFile
+print " Loading parameters from", filename
+if not os.access(filename, os.F_OK):
+  print "  Cannot find parameter file %s!" % filename
+  exit(-1)
+parms = imp.load_source("name", filename)
+globals().update(vars(parms))
+
 if not os.access(figDir, os.F_OK):
   os.makedirs(figDir)
 
@@ -38,10 +46,10 @@ for sampleDir in sampleDirList:
     dataSampleList.append(pickle.load(input))
 
 plt.clf()
-plt.xlabel("Void Radius (Mpc/h)")
+plt.xlabel("Void Density Contrast")
 plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]")
 plt.yscale('log')
-plt.xlim(xmax=80.)
+plt.xlim(xmax=5.)
 
 plotName = plotNameBase
 
@@ -62,6 +70,7 @@ for (iSample,sample) in enumerate(dataSampleList):
 
   boxVol *= 1.e-9 # Mpc->Gpc
 
+  print " Loading", sampleName
   filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\
              sampleName+".out"
   if not os.access(filename, os.F_OK):

plotting/plotCompareDist.py

@@ -14,8 +14,6 @@ import argparse
 
 # ------------------------------------------------------------------------------
 
-from datasetsToPlot import *
-
 plotNameBase = "compdist"
 
 obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing?
@@ -24,10 +22,20 @@ parser = argparse.ArgumentParser(description='Plot.')
 parser.add_argument('--show', dest='showPlot', action='store_const',
                    const=True, default=False,
                    help='display the plot (default: just write eps)')
+parser.add_argument('--parmFile', dest='parmFile', default='datasetsToPlot.py',
+                    help='path to parameter file')
 args = parser.parse_args()
 
 # ------------------------------------------------------------------------------
 
+filename = args.parmFile
+print " Loading parameters from", filename
+if not os.access(filename, os.F_OK):
+  print "  Cannot find parameter file %s!" % filename
+  exit(-1)
+parms = imp.load_source("name", filename)
+globals().update(vars(parms))
+
 if not os.access(figDir, os.F_OK):
   os.makedirs(figDir)
 
@@ -38,10 +46,10 @@ for sampleDir in sampleDirList:
     dataSampleList.append(pickle.load(input))
 
 plt.clf()
-plt.xlabel("Void Radius (Mpc/h)")
+plt.xlabel("Void Radius [Mpc/h]")
 plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]")
 plt.yscale('log')
-plt.xlim(xmax=80.)
+plt.xlim(xmax=120.)
 
 plotName = plotNameBase
 

plotting/plotDensConVsR.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python
+
+# plots cumulative distributions of number counts versus density contrast
+
+from void_python_tools.backend import *
+from void_python_tools.plotting import *
+import void_python_tools.apTools as vp
+import imp
+import pickle
+import os
+import matplotlib.pyplot as plt
+import numpy as np
+import argparse
+
+# ------------------------------------------------------------------------------
+
+plotNameBase = "densconvsr"
+
+obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing?
+
+parser = argparse.ArgumentParser(description='Plot.')
+parser.add_argument('--show', dest='showPlot', action='store_const',
+                   const=True, default=False,
+                   help='display the plot (default: just write eps)')
+parser.add_argument('--parmFile', dest='parmFile', default='datasetsToPlot.py',
+                    help='path to parameter file')
+args = parser.parse_args()
+
+# ------------------------------------------------------------------------------
+
+filename = args.parmFile
+print " Loading parameters from", filename
+if not os.access(filename, os.F_OK):
+  print "  Cannot find parameter file %s!" % filename
+  exit(-1)
+parms = imp.load_source("name", filename)
+globals().update(vars(parms))
+
+if not os.access(figDir, os.F_OK):
+  os.makedirs(figDir)
+
+dataSampleList = []
+
+for sampleDir in sampleDirList:
+  with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
+    dataSampleList.append(pickle.load(input))
+
+plt.clf()
+plt.ylabel("Void Density Contrast")
+plt.xlabel("Void Radius [Mpc/h]")
+#plt.xlim(xmax=5.)
+plt.yscale('log')
+
+plotName = plotNameBase
+
+for (iSample,sample) in enumerate(dataSampleList):
+
+  sampleName = sample.fullName
+  lineTitle = sampleName
+
+  if sample.dataType == "observation":
+    boxVol = vp.getSurveyProps(sample.maskFile, 
+                               sample.zBoundary[0], sample.zBoundary[1], 
+                               sample.zRange[0], sample.zRange[1], "all",
+                               selectionFuncFile=sample.selFunFile)[0]
+    boxVol *= obsFudgeFactor
+  else:
+    boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] - 
+                                          sample.zBoundaryMpc[0])
+
+  boxVol *= 1.e-9 # Mpc->Gpc
+
+  print " Loading", sampleName
+  filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\
+             sampleName+".out"
+  if not os.access(filename, os.F_OK):
+    print "File not found: ", filename
+    continue
+
+  data = np.loadtxt(filename, comments="#")
+  if data.ndim == 1:
+    print " Too few!"
+    continue
+
+  plt.plot(data[:,4], data[:,8], '-',
+           label=lineTitle, color=colorList[iSample],
+           linewidth=linewidth)
+
+plt.legend(title = "Samples", loc = "upper right", prop={'size':8})
+#plt.title(plotTitle)
+
+plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
+plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
+plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
+
+if args.showPlot:
+  os.system("display %s" % figDir+"/fig_"+plotName+".png")
+
+

python_tools/void_python_tools/backend/classes.py

@@ -167,3 +167,10 @@ def jobSuccessful(logFile, doneString):
 def getStackSuffix(zMin, zMax, rMin, rMax, dataPortion):
   return "z"+str(zMin)+"-"+str(zMax)+"_"+str(rMin)+"-"+str(rMax)+\
          "Mpc"+"_"+dataPortion
+
+def my_import(name):
+    mod = __import__(name)
+    components = name.split('.')
+    for comp in components[1:]:
+        mod = getattr(mod, comp)
+    return mod

python_tools/void_python_tools/backend/launchers.py

@@ -764,12 +764,10 @@ def launchProfile(sample, stack, voidDir=None, logFile=None, continueRun=None):
       fp.close()
       return
 
-    # TEST
     if sample.profileBinSize == "auto":
-      density = 0.5 * 50 / Rcircular
+      density = 0.5 * 50 / Rcircular / 2
     else:
-      #density = 0.5 * 50 / 60
+      density = sample.profileBinSize
-      density = 0.5 * 50 / Rcircular / sample.profileBinSize
 
     sys.stdout = open(logFile, 'w')
     sys.stderr = open(logFile, 'a')
@@ -836,10 +834,16 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
     while badChain:
       Rexpect = (stack.rMin+stack.rMax)/2
       Rtruncate = stack.rMin*3. + 1 # TEST
-      ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
+      if sample.dataType == "observation":
-                                    Niter=300000,
+        ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
-                                    Nburn=100000,
+                                      Niter=300000,
-                                    Rextracut=Rtruncate)
+                                      Nburn=100000,
+                                      Rextracut=Rtruncate)
+      else:
+        ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
+                                      Niter=300000,
+                                      Nburn=100000,
+                                      Rextracut=Rtruncate)
       badChain = (args[0][0] > 0.5 or args[0][1] > stack.rMax or \
                   args[0][2] > stack.rMax) and \
                  (ntries < maxtries)
@@ -862,8 +866,8 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
       rescaleFactor = 1.0
 
     # TEST - plotting raw galaxy counts
-    vp.draw_fit(*args, delta_tick=50, vmax=1.0, delta_v=0.01,
+    #vp.draw_fit(*args, delta_tick=50, vmax=1.0, delta_v=0.01,
-    #vp.draw_fit(*args, delta_tick=0.2, vmax=1.0, delta_v=0.01,
+    vp.draw_fit(*args, delta_tick=0.2, vmax=1.0, delta_v=0.01,
                 nocontour=True, model_max=1.0, delta_model=0.1,
                 plotTitle=plotTitle, show_ratio=showRatio,
                 showContours=showContours,
@@ -965,8 +969,12 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
       #   (not not fully accurate) plots
       for (iZBin, stack) in enumerate(sample.getUniqueZBins()):
 
-        aveDist = vp.aveStretchCone(stack.zMin, stack.zMax, 
+        if sample.dataType == "observation":
-                                    skyFrac = sample.skyFraction)
+          aveDist = vp.aveStretchCone(stack.zMin, stack.zMax, 
+                                      skyFrac = sample.skyFraction)
+        else:
+          aveDist = vp.aveStretch(stack.zMin, stack.zMax)
+
         aveDistList[iZBin, 0] = stack.zMin
         aveDistList[iZBin, 1] = aveDist
         aveDistList[iZBin, 2] = 0.00125
@@ -1010,16 +1018,24 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
           #                            skyFrac = sample.skyFraction, 
           #                            voidRedshifts=voidRedshifts)
 
-          aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+          if sample.dataType == "observation":
-                                      skyFrac = sample.skyFraction)
+            aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+                                        skyFrac = sample.skyFraction)
+          else:
+            aveDist = vp.aveStretch(zBin.zMin, zBin.zMax)
+
           expList[0, iR, iZBin, 2] = aveDist
 
           for (iZCheck,zBinCheck) in enumerate(allZBins):
             for (iRCheck,rBinCheck) in enumerate(allRBins):
               if zBin.zMin == zBinCheck.zMin and zBin.zMax == zBinCheck.zMax:
                 if rBin.rMin == rBinCheck.rMin and rBin.rMax == rBinCheck.rMax:
-                  aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+                  if sample.dataType == "observation":
-                                              skyFrac = sample.skyFraction)
+                    aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+                                                skyFrac = sample.skyFraction)
+                  else:
+                    aveDist = vp.aveStretch(zBin.zMin, zBin.zMax)
+
                   allExpList[iSample, iRCheck, iZCheck, 0] = exp
                   allExpList[iSample, iRCheck, iZCheck, 1] = expError
 
@@ -1112,8 +1128,12 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
         if zBin.zMaxPlot > plotZmax: plotZmax = zBin.zMaxPlot
         if zBin.zMinPlot < plotZmin: plotZmin = zBin.zMinPlot
 
-        aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+        if sample.dataType == "observation":
-                                    skyFrac = sample.skyFraction)
+          aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+                                      skyFrac = sample.skyFraction)
+        else:
+          aveDist = vp.aveStretch(zBin.zMin, zBin.zMax)
+
         aveDistList[iZ, 0] = zBin.zMin
         aveDistList[iZ, 1] = aveDist
         aveDistList[iZ, 2] = 0.00125

python_tools/void_python_tools/plotting/plotDefs.py

@@ -1,15 +1,15 @@
 LIGHT_SPEED = 299792.458 
 
-colorList = ['r', 'b', 'g', 'y', 'c', 'm', 'y', 
+colorList = ['r', 'b', 'g', 'y', 'c', 'm', 
              'brown', 'grey',
-                          'darkred', 'orange', 'pink', 'darkblue',
+             'darkred', 'orange', 'pink', 'darkblue',
-                          'lightblue', 'chocolate',
+             'lightblue', 'chocolate',
-                          'indigo', 'lightseagreen', 'maroon', 'olive',
+             'indigo', 'lightseagreen', 'maroon', 'olive',
-                          'royalblue', 'palevioletred', 'seagreen', 'tomato',
+             'royalblue', 'palevioletred', 'seagreen', 'tomato',
-                          'aquamarine', 'darkslateblue',
+             'aquamarine', 'darkslateblue',
-                          'khaki', 'lawngreen', 'mediumorchid',
+             'khaki', 'lawngreen', 'mediumorchid',
-                           'orangered', 'thistle'
+             'orangered', 'thistle'
-                          'yellowgreen']
+             'yellowgreen']
 
 linewidth = 4
 fontsize = 12