added script to compare stacks and radial profiles of matched voids

c_tools/analysis/voidOverlap.cpp

@@ -271,7 +271,7 @@ int main(int argc, char **argv) {
   filename = string(args.outfile_arg);
   filename = filename.append("summary.out");
   fp = fopen(filename.c_str(), "w");
-  fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, common volume ratio 2, relative dist, num matches, num significant matches\n");
+  fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, common volume ratio 2, relative dist, num matches, num significant matches, match ID\n");
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
     int voidID = catalog1.voids[iVoid1].voidID;
     if (catalog1.voids[iVoid1].numMatches > 0) {
@@ -286,14 +286,15 @@ int main(int argc, char **argv) {
       rdist = catalog1.voids[iVoid1].matches[0].dist;
       rdist /= catalog1.voids[iVoid1].radius;
 
-      fprintf(fp, "%d %.4f %.4f %.4f %.4f %.4f %d %d\n", voidID, 
+      fprintf(fp, "%d %.4f %.4f %.4f %.4f %.4f %d %d %d\n", voidID, 
                                    catalog1.voids[iVoid1].radiusMpc,
                                    rRatio, 
                                    commonVolRatio,
                                    volRatio,
                                    rdist, 
                                    catalog1.voids[iVoid1].numMatches,
-                                   catalog1.voids[iVoid1].numBigMatches);
+                                   catalog1.voids[iVoid1].numBigMatches,
+                                   catalog2.voids[iVoid2].voidID);
 
     } else {
       fprintf(fp, "%d %.2f 0.0 0.0 0.0 0.0 0 0\n", voidID, 

crossCompare/analysis/compareProfiles.py

@@ -0,0 +1,121 @@
+#!/usr/bin/env python
+#+
+#   VIDE -- Void IDEntification pipeline -- ./pipeline/apAnalysis.py
+#   Copyright (C) 2010-2013 Guilhem Lavaux
+#   Copyright (C) 2011-2013 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.
+#+
+
+# takes a set of voids with given radii, and stacks and plots matches in 
+# other catalogs
+
+from void_python_tools.backend import *
+import imp
+import pickle
+import os
+import numpy as np
+import argparse
+
+# ------------------------------------------------------------------------------
+
+mergerNameBase = "mergerTree"
+
+parser = argparse.ArgumentParser(description='Analyze.')
+parser.add_argument('--parm', dest='parm', default='datasetsToAnalyze.py', help='path to parameter file')
+args = parser.parse_args()
+
+# ------------------------------------------------------------------------------
+
+filename = args.parm
+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(dataDir, os.F_OK):
+  os.makedirs(dataDir)
+
+mergerFileBase = dataDir + "/" + mergerNameBase
+
+baseIDList = []
+
+# get list of base voids
+with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
+  baseSample = pickle.load(input)
+  baseSampleName = baseSample.fullName
+  baseVoidList = np.loadtxt(workDir+baseSampleDir+"/centers_nocut_central_"+\
+                    baseSampleName+".out")
+
+for stack in baseSample.stacks:
+  print " Stack:", stack.rMin    
+
+  accepted = (baseVoidList[:,4] > stack.rMin) & (baseVoidList[:,4] < stack.rMax)
+  baseIDList = baseVoidList[accepted][:,7]
+
+  for (iSample, sampleDir) in enumerate(sampleDirList):
+    with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
+      sample = pickle.load(input)
+
+    print "  Working with", sample.fullName, "...",
+    sys.stdout.flush()
+    sampleName = sample.fullName
+
+    # get list of appropriate voids
+    if sample == baseSample:
+      idList = baseIDList
+    else:
+      idList = 2
+      matchList = np.loadtxt(mergerFileBase+"_"+baseSampleName+"_"+sampleName+\
+                            "_summary.out")
+      accepted = (matchList[:,0] == baseIDList).any()
+      idList = matchList[accepted][:,8]
+
+    voidBaseDir = workDir+"/"+sampleDir+"stacks"
+        
+    runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, 
+                               stack.rMax, thisDataPortion, 
+                               customLine="selected")
+    stack.rMin = 0.
+    stack.rMax = 1000.
+
+    voidDir = voidBaseDir+"_"+runSuffix
+
+    if not os.access(voidDir,os.F_OK): os.makedirs(voidDir)
+
+    print "    Stacking voids...",
+    sys.stdout.flush()
+        
+    STACK_PATH = CTOOLS_PATH+"/stacking/stackVoidsZero"
+
+    launchStack(sample, stack, STACK_PATH, 
+                   thisDataPortion=thisDataPortion,
+                   logDir=logDir, voidDir=voidDir, 
+                   zobovDir=workDir+"/"+sampleDir,
+                   freshStack=freshStack, INCOHERENT=False,
+                   ranSeed=101010, summaryFile=None, 
+                   continueRun=continueRun,
+                   dataType=sample.dataType,
+                   idList=idList)
+ 
+    print "   Profiling stacks...",
+    sys.stdout.flush()
+    logFile = logDir+"/profile_"+sampleName+"_"+runSuffix+".out"
+    launchProfile(sample, stack, voidDir=voidDir, 
+                  logFile=logFile, continueRun=False)
+
+print " Done!"

crossCompare/analysis/mergerTree.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python
 #+
 #   VIDE -- Void IDEntification pipeline -- ./crossCompare/analysis/mergerTree.py
 #   Copyright (C) 2010-2013 Guilhem Lavaux
@@ -17,9 +18,6 @@
 #   with this program; if not, write to the Free Software Foundation, Inc.,
 #   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 #+
-#!/usr/bin/env python
-
-# plots cumulative distributions of number counts
 
 from void_python_tools.backend import *
 from void_python_tools.plotting import *
@@ -35,13 +33,13 @@ import argparse
 dataNameBase = "mergerTree"
 
 parser = argparse.ArgumentParser(description='Analyze.')
-parser.add_argument('--parmFile', dest='parmFile', default='datasetsToAnalyze.py',
+parser.add_argument('--parm', dest='parm', default='datasetsToAnalyze.py',
                     help='path to parameter file')
 args = parser.parse_args()
 
 # ------------------------------------------------------------------------------
 
-filename = args.parmFile
+filename = args.parm
 print " Loading parameters from", filename
 if not os.access(filename, os.F_OK):
   print "  Cannot find parameter file %s!" % filename
@@ -78,7 +76,7 @@ for (iSample, sampleDir) in enumerate(sampleDirList):
                     thisDataPortion="central", logFile=logFile,
                     continueRun=False, workDir=workDir,
                     outputFile=stepOutputFileName,
-                    matchMethod="proximity")
+                    matchMethod="useID")
 
   # attach columns to summary file
   #if iSample == 1:
@@ -101,3 +99,5 @@ for (iSample, sampleDir) in enumerate(sampleDirList):
 #if os.access("mergerTree.log", os.F_OK): os.unlink("mergerTree.log")
 #if os.access("temp.out", os.F_OK): os.unlink("temp.out")
 #if os.access("thisStep.out", os.F_OK): os.unlink("thisStep.out")
+
+print " Done!"

crossCompare/plotting/plotNumberFunc.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python
 #+
 #   VIDE -- Void IDEntification pipeline -- ./crossCompare/plotting/plotNumberFunc.py
 #   Copyright (C) 2010-2013 Guilhem Lavaux
@@ -17,7 +18,6 @@
 #   with this program; if not, write to the Free Software Foundation, Inc.,
 #   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 #+
-#!/usr/bin/env python
 
 # plots cumulative distributions of number counts
 
@@ -36,8 +36,8 @@ from scipy.stats import ks_2samp
 
 plotNameBase = "compdist"
 
-obsFudgeFactor = 1.0 # what fraction of the volume are we *reall* capturing?
-#obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing?
+#obsFudgeFactor = 1.0 # what fraction of the volume are we *reall* capturing?
+obsFudgeFactor = .15 # what fraction of the volume are we *reall* capturing?
 
 linewidth = 1
 
@@ -45,7 +45,7 @@ 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',
+parser.add_argument('--parm', dest='parm', default='datasetsToPlot.py',
                     help='path to parameter file')
 args = parser.parse_args()
 
@@ -55,7 +55,7 @@ errorBarsX = np.linspace(0, plotMax, num=nErrorBars)
 
 # ------------------------------------------------------------------------------
 
-filename = args.parmFile
+filename = args.parm
 print " Loading parameters from", filename
 if not os.access(filename, os.F_OK):
   print "  Cannot find parameter file %s!" % filename
@@ -102,7 +102,7 @@ for dataPortion in dataPortions:
 
     boxVol *= 1.e-9 # Mpc->Gpc
 
-    filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\
+    filename = workDir+"/"+sampleDirList[iSample]+"/centers_nocut_"+dataPortion+"_"+\
                sampleName+".out"
     if not os.access(filename, os.F_OK):
       print "File not found: ", filename
@@ -137,8 +137,13 @@ for dataPortion in dataPortions:
              ecolor=colorList[iSample],
              fmt=None, label='_nolegend_', capsize=0)
 
-    hist, bin_edges = np.histogram(data, bins=100, range=(0,100)) 
-    allData.append(hist)
+    hist, bin_edges = np.histogram(data, bins=40, range=(0,100))
+    #allData.append(hist)
+
+    binCenters = 0.5*(bin_edges[1:] + bin_edges[:-1])
+    #plt.plot(binCenters, hist, '-',
+    #         label=lineTitle, color=colorList[iSample],
+    #         linewidth=linewidth)
 
   plt.legend(title = "Samples", loc = "upper right", prop={'size':8})
   #plt.title(plotTitle)
@@ -155,19 +160,19 @@ 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")
 
-dataFile = figDir+"/data_"+plotName+".dat"
-fp = open(dataFile, 'w')
-fp.write("# R [Mpc/h], N [h^3 Gpc^-3]\n")
-fp.write("# ")
-for sample in dataSampleList:
-  fp.write(sample.fullName+" ")
-fp.write("\n")
-for i in xrange(100):
-  fp.write(str(bin_edges[i]) + " ")
-  for iSample in xrange(len(dataSampleList)):
-    fp.write(str(allData[iSample][i])+" ")
-  fp.write("\n")
-fp.close()
+#dataFile = figDir+"/data_"+plotName+".dat"
+#fp = open(dataFile, 'w')
+#fp.write("# R [Mpc/h], N [h^3 Gpc^-3]\n")
+#fp.write("# ")
+#for sample in dataSampleList:
+#  fp.write(sample.fullName+" ")
+#fp.write("\n")
+#for i in xrange(100):
+#  fp.write(str(bin_edges[i]) + " ")
+#  for iSample in xrange(len(dataSampleList)):
+#    fp.write(str(allData[iSample][i])+" ")
+#  fp.write("\n")
+#fp.close()
 
 if args.showPlot:
   os.system("display %s" % figDir+"/fig_"+plotName+".png")

python_tools/void_python_tools/backend/classes.py

@@ -184,9 +184,9 @@ def jobSuccessful(logFile, doneString):
       if doneString in line: jobDone = True
   return jobDone
 
-def getStackSuffix(zMin, zMax, rMin, rMax, dataPortion):
+def getStackSuffix(zMin, zMax, rMin, rMax, dataPortion, customLine=""):
   return "z"+str(zMin)+"-"+str(zMax)+"_"+str(rMin)+"-"+str(rMax)+\
-         "Mpc"+"_"+dataPortion
+           "Mpc"+customLine+"_"+dataPortion
 
 def my_import(name):
     mod = __import__(name)

python_tools/void_python_tools/backend/launchers.py

@@ -437,12 +437,19 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
                 voidDir=None, freshStack=True, runSuffix=None,
                 zobovDir=None,
                 INCOHERENT=False, ranSeed=None, summaryFile=None, 
-                continueRun=None, dataType=None, prefixRun=""):
+                continueRun=None, dataType=None, prefixRun="",
+                idList=None):
 
   sampleName = sample.fullName
 
+  if idList != None:
+    customLine = "selected"
+  else:
+    customLine = ""
+
   runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin,
-                             stack.rMax, thisDataPortion)
+                             stack.rMax, thisDataPortion, 
+                             customLine=customLine)
 
   logFile = logDir+("/%sstack_"%prefixRun)+sampleName+"_"+runSuffix+".out"
  
@@ -483,6 +490,16 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
   else:
     rescaleFlag = ""
 
+  if idList != None:
+    idListFlag = "idList '" + str(idList).strip('[]') + "'"
+    rMinToUse = 0.
+    rMaxToUse = 100000.
+    centralRadius = rMaxToUse
+  else:
+    idListFlag = ""
+    rMinToUse = stack.rMin
+    rMaxToUse = stack.rMax
+
   conf="""
   desc %s
   partzone %s
@@ -505,13 +522,14 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
   barycenters %s
   boundaryDistances %s
   %s
+  %s
   """ % \
   (zobovDir+"/voidDesc_"+thisDataPortion+"_"+sampleName+".out",
    zobovDir+"/voidPart_"+sampleName+".dat",
    zobovDir+"/voidZone_"+sampleName+".dat",
    zobovDir+"/vol_"+sampleName+".dat",
-   stack.rMin,
-   stack.rMax,
+   rMinToUse,
+   rMaxToUse,
    zobovDir+("/%szobov_slice_"%prefixRun)+sampleName,
    zobovDir+"/zobov_slice_"+sampleName+".par",
    maxDen,
@@ -526,7 +544,9 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
    thisDataPortion,
    zobovDir+"/barycenters_"+thisDataPortion+"_"+sampleName+".out",
    zobovDir+"/boundaryDistances_"+thisDataPortion+"_"+sampleName+".out",
-   rescaleFlag)
+   rescaleFlag,
+   idListFlag 
+   )
 
   parmFile = os.getcwd()+("/%sstack_"%prefixRun)+sample.fullName+".par"