Merge branch 'master' of bitbucket.org:cosmicvoids/void_identification

c_tools/stacking/pruneVoids.cpp

@@ -548,7 +548,7 @@ int main(int argc, char **argv) {
     }
   }
   voids.resize(iGood);
-  printf("  1st filter: reiGoodected %d obviously bad\n", numWrong); 
+  printf("  1st filter: rejected %d obviously bad\n", numWrong); 
 
   iGood = 0;
   for (iVoid = 0; iVoid < voids.size(); iVoid++) {
@@ -559,7 +559,7 @@ int main(int argc, char **argv) {
     }
   }
   voids.resize(iGood);
-  printf("  2nd filter: reiGoodected %d too small\n", numTooSmall); 
+  printf("  2nd filter: rejected %d too small\n", numTooSmall); 
 
 
   iGood = 0;
@@ -572,7 +572,7 @@ int main(int argc, char **argv) {
     }
   }
   voids.resize(iGood);
-  printf("  3rd filter: reiGoodected %d too close to high redshift boundaries\n", numNearZ); 
+  printf("  3rd filter: rejected %d too close to high redshift boundaries\n", numNearZ); 
 
   numNearZ = 0;
   iGood = 0;
@@ -586,7 +586,7 @@ int main(int argc, char **argv) {
     }
   }
   voids.resize(iGood);
-  printf("  4th filter: reiGoodected %d too close to low redshift boundaries\n", numNearZ); 
+  printf("  4th filter: rejected %d too close to low redshift boundaries\n", numNearZ); 
 
   for (iVoid = 0; iVoid < voids.size(); iVoid++) {
     if (voids[iVoid].centralDen > args.maxCentralDen_arg) {

python_tools/pipeline_source/prepareCatalogs.in.py

@@ -62,7 +62,7 @@ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1):
 def writeScript(setName, dataFileNameBase, dataFormat,
                 scriptDir, catalogDir, fileNums, redshifts, numSubvolumes,
                 numSlices, useVel, lbox, minRadius, omegaM, subsample=1.0, 
-                suffix=".dat"):
+                suffix=".dat", dataFileNameList=None):
 
 
   if useVel: setName += "_pv"
@@ -134,7 +134,7 @@ newSample = Sample(dataFile = "{dataFile}",
                    numSubvolumes = {numSubvolumes},
                    mySubvolume = "{mySubvolume}",
                    useLightCone = {useLightCone},
-                   subsample = {subsample})
+                   subsample = "{subsample}")
 dataSampleList.append(newSample)
 newSample.addStack(0.0, 5.0, 20, 25, False, False)
 newSample.addStack(0.0, 5.0, 30, 35, False, False)
@@ -165,33 +165,32 @@ newSample.addStack(0.0, 5.0, 90, 95, False, False)
     for i in xrange(len(zVsDY)):
       zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
 
-    if useLightCone:
+    boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \
-      boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \
+                LIGHT_SPEED*lbox, zVsDX, zVsDY)-zBox
-                  LIGHT_SPEED*lbox, zVsDX, zVsDY)-zBox
-      dzSafe = 0.03
-    else:
-      boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \
-                  LIGHT_SPEED*lbox, zVsDX, zVsDY)-zBox
-      #boxWidthZ = lbox*100./LIGHT_SPEED
-      dzSafe = 0.0
 
     for iSlice in xrange(numSlices):
-      sliceMin = zBox + dzSafe + iSlice*(boxWidthZ-2.*dzSafe)/numSlices
-      sliceMax = zBox + dzSafe + (iSlice+1)*(boxWidthZ-2.*dzSafe)/numSlices
 
       if useLightCone:
+        dzSafe = 0.03
+        sliceMin = zBox + dzSafe + iSlice*(boxWidthZ-2.*dzSafe)/numSlices
+        sliceMax = zBox + dzSafe + (iSlice+1)*(boxWidthZ-2.*dzSafe)/numSlices
         sliceMinMpc = sliceMin*LIGHT_SPEED/100.
         sliceMaxMpc = sliceMax*LIGHT_SPEED/100.
       else:
-        sliceMinMpc = LIGHT_SPEED/100.*vp.angularDiameter(sliceMin, Om=Om)
+        sliceMinMpc = zBoxMpc + iSlice*lbox/numSlices
-        sliceMaxMpc = LIGHT_SPEED/100.*vp.angularDiameter(sliceMax, Om=Om)
+        sliceMaxMpc = zBoxMpc + (iSlice+1)*lbox/numSlices
+        sliceMin = np.interp(sliceMinMpc*100./LIGHT_SPEED, zVsDX, zVsDY)
+        sliceMax = np.interp(sliceMaxMpc*100./LIGHT_SPEED, zVsDX, zVsDY)
 
       sliceMin = "%0.2f" % sliceMin
       sliceMax = "%0.2f" % sliceMax
       sliceMinMpc = "%0.1f" % sliceMinMpc
       sliceMaxMpc = "%0.1f" % sliceMaxMpc
 
-      dataFileName = dataFileNameBase + fileNum + suffix
+      if (dataFileNameList != None):
+        dataFileName = dataFileNameList[iFile]
+      else:
+        dataFileName = dataFileNameBase + fileNum + suffix
 
       for iX in xrange(numSubvolumes):
         for iY in xrange(numSubvolumes):
@@ -216,7 +215,7 @@ newSample.addStack(0.0, 5.0, 90, 95, False, False)
                                          numSubvolumes=numSubvolumes,
                                          mySubvolume=mySubvolume,
                                          useLightCone=useLightCone,
-                                         subsample=subsample))
+                                         subsample=str(subsample).strip('[]')))
 
   scriptFile.close()
   return
@@ -233,43 +232,59 @@ if not os.access(catalogDir, os.F_OK): os.mkdir(catalogDir)
 #       ss0.000175 ~ SDSS DR9 mid 
 baseResolution = float(numPart)/lbox/lbox/lbox # particles/Mpc^3
 prevSubSample = -1
-for thisSubSample in sorted(subSamples, reverse=True):
+subSamples = sorted(subSamples, reverse=True)
+for iSubSample in xrange(len(subSamples)):
 
-  keepFraction = float(thisSubSample) / baseResolution
+  subSampleList = subSamples[0:iSubSample+1]
-  maxKeep = keepFraction * numPart
+
+  keepFractionList = []
+  for subSample in subSampleList:
+    keepFractionList.append(float(subSample) / baseResolution)
+  thisSubSample = subSamples[iSubSample]
+  maxKeep = keepFractionList[-1] * numPart
   minRadius = int(np.ceil(lbox/maxKeep**(1./3)))
 
+  partFileList = []
+  for (iRedshift, redshift) in enumerate(redshifts):
+    if particleFileDummy == '':
+      partFileList.append(particleFileBase+fileNums[iRedshift])
+    else:
+      partFileList.append(particleFileBase.replace(particleFileDummy, 
+                                                   fileNums[iRedshift]))
+
   if args.script or args.all:
     print " Doing subsample", thisSubSample, "scripts"
     sys.stdout.flush()
     setName = prefix+"ss"+str(thisSubSample)
-    if dataFormat == "multidark":
+
-      subSampleToUse = 1.0
+    if dataFormat == "random":
-      fileToUse = prefix+"ss"+str(thisSubSample)+"_z"
-      suffix = ".dat"
-    elif dataFormat == "gadget":
-      subSampleToUse = keepFraction
-      #subSampleToUse = thisSubSample
-      fileToUse = particleFileBase
-      suffix = ""
-    elif dataFormat == "lanl":
-      subSampleToUse = keepFraction
-      #subSampleToUse = thisSubSample
-      fileToUse = particleFileBase
-      suffix = ""
-    elif dataFormat == "random":
       subSampleToUse = 1.0
       fileToUse = "ran.ss"+str(thisSubSample)+"_z"
       suffix = ".dat"
 
-    writeScript(setName, fileToUse, dataFormat, 
+      writeScript(setName, fileToUse, dataFormat, 
                   scriptDir, catalogDir, fileNums, redshifts, 
                   numSubvolumes, numSlices, False, lbox, minRadius, omegaM, 
                   subsample=subSampleToUse, suffix=suffix)
-    writeScript(setName, fileToUse, dataFormat, 
+      writeScript(setName, fileToUse, dataFormat, 
                   scriptDir, catalogDir, fileNums, redshifts, 
                   numSubvolumes, numSlices, True, lbox, minRadius, omegaM, 
                   subsample=subSampleToUse, suffix=suffix)
+    else:
+      subSampleToUse = keepFractionList
+      suffix = ""
+      fileToUse = partFileList[0]
+      writeScript(setName, fileToUse, dataFormat, 
+                  scriptDir, catalogDir, fileNums, redshifts, 
+                  numSubvolumes, numSlices, False, lbox, minRadius, omegaM, 
+                  subsample=subSampleToUse, suffix=suffix,
+                  dataFileNameList=partFileList)
+      writeScript(setName, fileToUse, dataFormat, 
+                  scriptDir, catalogDir, fileNums, redshifts, 
+                  numSubvolumes, numSlices, True, lbox, minRadius, omegaM, 
+                  subsample=subSampleToUse, suffix=suffix,
+                  dataFileNameList=partFileList)
+       
 
   if args.subsample or args.all:
     print " Doing subsample", thisSubSample
@@ -411,31 +426,37 @@ if (args.halos or args.all) and haloFileBase != "":
 
       sampleName = prefix+"halos_min"+str(minHaloMass)+"_z"+fileNums[iRedshift]
       outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
-
       outFile.write("%f\n" %(lbox))
       outFile.write("%s\n" %(omegaM))
       outFile.write("%s\n" %(hubble))
       outFile.write("%s\n" %(redshift))
       outFile.write("%d\n" %(numPart))
 
-      inFile = open(dataFile, 'r')
+      if dataFormat == "sdf":
-      for (iHalo,line) in enumerate(inFile):
+        # TODO process halo file with SDFcvt
-        if iHalo < haloFileNumComLines: continue
+        outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
-        line = line.split(haloFileColSep)
+        outFile.write("-99 -99 -99 -99 -99 -99 -99\n")
-        if minHaloMass == "none" or float(line[haloFileMCol]) > minHaloMass:
+        outFile.close()
-          x  = float(line[haloFileXCol])
+      else:
-          y  = float(line[haloFileYCol])
+        outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
-          z  = float(line[haloFileZCol])
+        inFile = open(dataFile, 'r')
-          vz = float(line[haloFileVZCol])
+        for (iHalo,line) in enumerate(inFile):
-          vy = float(line[haloFileVYCol])
+          if iHalo < haloFileNumComLines: continue
-          vx = float(line[haloFileVXCol])
+          line = line.split(haloFileColSep)
+          if minHaloMass == "none" or float(line[haloFileMCol]) > minHaloMass:
+            x  = float(line[haloFileXCol])
+            y  = float(line[haloFileYCol])
+            z  = float(line[haloFileZCol])
+            vz = float(line[haloFileVZCol])
+            vy = float(line[haloFileVYCol])
+            vx = float(line[haloFileVXCol])
 
-          # write to output file
+            # write to output file
-          outFile.write("%d %e %e %e %e %e %e\n" %(iHalo,x,y,z,vz,vy,vx))
+            outFile.write("%d %e %e %e %e %e %e\n" %(iHalo,x,y,z,vz,vy,vx))
 
-      outFile.write("-99 -99 -99 -99 -99 -99 -99\n")
+        outFile.write("-99 -99 -99 -99 -99 -99 -99\n")
-      inFile.close()
+        outFile.close()
-      outFile.close()
+        inFile.close()
 
 # -----------------------------------------------------------------------------
 # now the HOD
@@ -510,6 +531,14 @@ if (args.hod or args.all) and haloFileBase != "":
       else:
         haloFile = catalogDir+haloFileBase.replace(haloFileDummy,
                                                  fileNums[iRedshift])
+
+      if dataFormat == "sdf":
+        # TODO process halo file with SDFcvt
+        inFile = haloFile
+        outFile = haloFile+"_temp"
+
+        haloFile = outFile
+
       parFile.write(parFileText.format(omegaM=omegaM,
                                      hubble=hubble,
                                      redshift=redshift,
@@ -535,4 +564,6 @@ if (args.hod or args.all) and haloFileBase != "":
       os.system("mv %s/hod.mock %s" % (catalogDir, outFileName))
       os.system("rm %s/hod.*" % catalogDir)
 
+      os.system("rm %s" % haloFile)
+
 print " Done!"

python_tools/void_python_tools/backend/launchers.py

@@ -105,37 +105,43 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
 
     datafile = inputDataDir+"/"+sample.dataFile
 
-    if regenerate:
+    prevSubSample = -1
-      inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par"
+    for thisSubSample in sample.subsample.split(', '):
-      outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName 
+
-    else:
+      if prevSubSample == -1:
-      inputParameterFlag = ""
+        inputParameterFlag = ""
-      outputFile = zobovDir+"/zobov_slice_" + sampleName
+        outputFile = zobovDir+"/zobov_slice_" + sampleName + "_ss" + thisSubSample
+        keepFraction = float(thisSubSample)
+        subSampleLine = "subsample %g" % keepFraction
+      else:
+        inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+\
+                             sampleName+"_ss"+prevSubSample+".par"
+        outputFile = zobovDir+"/_zobov_slice_" + sampleName + "_ss" + thisSubSample
+        keepFraction = float(thisSubSample)/float(prevSubSample)
+        subSampleLine = "resubsample %g" % keepFraction
 
-    if sample.usePecVel:
-      includePecVelString = "peculiarVelocities"
-    else:
       includePecVelString = ""
+      if sample.usePecVel: includePecVelString = "peculiarVelocities"
 
-    if sample.useLightCone:
-      useLightConeString = "cosmo"
-    else:
       useLightConeString = ""
+      if sample.useLightCone: useLightConeString = "cosmo"
 
-    if sample.dataFormat == "multidark" or sample.dataFormat == "random":
+      if sample.dataFormat == "multidark" or sample.dataFormat == "random":
-      dataFileLine = "multidark " + datafile
+        dataFileLine = "multidark " + datafile
-    elif sample.dataFormat == "gadget":
+      elif sample.dataFormat == "gadget":
-      dataFileLine = "gadget " + datafile
+        dataFileLine = "gadget " + datafile
+      elif sample.dataFormat == "sdf":
+        dataFileLine = "sdf " + datafile
     
-    iX = float(sample.mySubvolume[0])
+      iX = float(sample.mySubvolume[0])
-    iY = float(sample.mySubvolume[1])
+      iY = float(sample.mySubvolume[1])
 
-    xMin = iX/sample.numSubvolumes * sample.boxLen
+      xMin = iX/sample.numSubvolumes * sample.boxLen
-    yMin = iY/sample.numSubvolumes * sample.boxLen
+      yMin = iY/sample.numSubvolumes * sample.boxLen
-    xMax = (iX+1)/sample.numSubvolumes * sample.boxLen
+      xMax = (iX+1)/sample.numSubvolumes * sample.boxLen
-    yMax = (iY+1)/sample.numSubvolumes * sample.boxLen
+      yMax = (iY+1)/sample.numSubvolumes * sample.boxLen
 
-    conf="""
+      conf="""
       %s
       output %s
       outputParameter %s
@@ -148,36 +154,52 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
       rangeY_max %g
       rangeZ_min %g
       rangeZ_max %g
-      subsample %g
+      %s
       %s
       """ % (dataFileLine, outputFile,
-             zobovDir+"/zobov_slice_"+sampleName+".par",
+             outputFile+".par",
              includePecVelString,
              useLightConeString,
              sample.dataUnit,
              xMin, xMax, yMin, yMax,
              sample.zBoundaryMpc[0], sample.zBoundaryMpc[1],
-             sample.subsample,inputParameterFlag)
+             subSampleLine,inputParameterFlag)
 
-    parmFile = os.getcwd()+"/generate_"+sample.fullName+".par"
+      parmFile = os.getcwd()+"/generate_"+sample.fullName+".par"
 
-    file(parmFile, mode="w").write(conf)
+      file(parmFile, mode="w").write(conf)
 
-    if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")):
+      doneLine = "Done! %5.2e\n" % keepFraction
-      cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile)
+      print "TEST", doneLine
-      os.system(cmd)
+      if not (continueRun and jobSuccessful(logFile, doneLine)):
-      if jobSuccessful(logFile, "Done!\n"):
+        if (prevSubSample == -1):
-        print "done"
+          cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile)
-      else:
+        else:
-        print "FAILED!"
+          cmd = "%s --configFile=%s &>> %s" % (binPath,parmFile,logFile)
-        exit(-1)
+        os.system(cmd)
       
-    else:
+        if not jobSuccessful(logFile, doneLine):
-      print "already done!"
+          print "FAILED!"
+          exit(-1)
+
+      # remove intermediate files
+      if (prevSubSample != -1):
+       os.unlink(zobovDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample+".par")
+       os.unlink(zobovDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample)
+
+      prevSubSample = thisSubSample
+
+    if (continueRun and jobSuccessful(logFile, doneLine)): print "already done!"
+    if jobSuccessful(logFile, doneLine): print "done"
+ 
+    # place the final subsample     
+    os.system("mv %s %s" % (zobovDir+"/zobov_slice_"+sampleName+"_ss"+\
+              prevSubSample, zobovDir+"/zobov_slice_"+sampleName))
+    os.system("mv %s %s" % (zobovDir+"/zobov_slice_"+sampleName+"_ss"+\
+              prevSubSample+".par", zobovDir+"/zobov_slice_"+sampleName+".par"))
 
     if os.access("comoving_distance.txt", os.F_OK):
       os.system("mv %s %s" % ("comoving_distance.txt", zobovDir))
-      #os.system("mv %s %s" % ("sample_info.txt", zobovDir))
 
     if os.access(parmFile, os.F_OK):
       os.unlink(parmFile)