futher fixes for handling sub-volumes in analysis

c_tools/libzobov/voidTree.hpp

@@ -126,8 +126,8 @@ public:
         // Go bigger, though I would say we should not to.
       }
     while (iter_candidate != candidateList.begin()) ;
-    if (vid_good_candidate < 0)
-      std::cout << "Failure to lookup parent (2)" << std::endl;
+    //if (vid_good_candidate < 0)
+    //  std::cout << "Failure to lookup parent (2)" << std::endl;
     return vid_good_candidate;
   }
 

c_tools/mock/generateFromCatalog.cpp

@@ -414,6 +414,7 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn
       for (uint32_t i = 0; i < pdata.pos.size(); i++)
 	{
 	  f.writeReal32((pdata.pos[i].xyz[j]+Lmax)/(2*Lmax));
+if (i < 10) printf("TEST WRITE %d %e\n", (pdata.pos[i].xyz[j]+Lmax)/(2*Lmax));
 	}
       f.endCheckpoint();
     }

c_tools/stacking/pruneVoids.cpp

@@ -254,8 +254,9 @@ int main(int argc, char **argv) {
   for (iVoid = 0; iVoid < numVoids; iVoid++) {
 
     voidID = voids[iVoid].voidID;
-    //printf("  DOING %d (of %d) %d %d\n", iVoid, numVoids, voidID, 
-     //                                    voids[iVoid].numPart);
+      printf("  DOING %d (of %d) %d %d %f\n", iVoid, numVoids, voidID, 
+                                           voids[iVoid].numPart, 
+                                           voids[iVoid].radius);
 
     voids[iVoid].center[0] = part[voids[iVoid].coreParticle].x;
     voids[iVoid].center[1] = part[voids[iVoid].coreParticle].y;
@@ -411,30 +412,36 @@ int main(int argc, char **argv) {
   } // iVoid
 
   printf(" Picking winners and losers...\n");
-  numKept = numVoids;
+  for (iVoid = 0; iVoid < numVoids; iVoid++) {
+    voids[iVoid].accepted = 1;
+  }
+
   for (iVoid = 0; iVoid < numVoids; iVoid++) {
     if (strcmp(args_info.dataPortion_arg, "edge")  == 0 &&
         tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) {
-      numKept--;
       voids[iVoid].accepted = 0;
     }
 
     if (strcmp(args_info.dataPortion_arg, "central") == 0 && 
         tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestMock) {
-      numKept--;
+      voids[iVoid].accepted = 0;
+    }
+
+    if (voids[iVoid].radius < args_info.rMin_arg) {
       voids[iVoid].accepted = 0;
     }
 
     if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) {
-      numKept--;
       voids[iVoid].accepted = -1;
     }
 
-    if (voids[iVoid].radius < args_info.rMin_arg) {
-      numKept--;
-      voids[iVoid].accepted = 0;
-    }
   }
+
+  numKept = 0;
+  for (iVoid = 0; iVoid < numVoids; iVoid++) {
+    if (voids[iVoid].accepted == 1) numKept++;
+  }
+
   printf(" Number kept: %d (out of %d)\n", numKept, numVoids);
   
   printf(" Output...\n");
@@ -452,7 +459,7 @@ int main(int argc, char **argv) {
     if (voids[iVoid].accepted != 1) continue;
 
      fprintf(fp, "%d %d %d %f %f %d %d %f %d %f %f\n", 
-             i, 
+             iVoid, 
              voids[iVoid].voidID, 
              voids[iVoid].coreParticle, 
              voids[iVoid].coreDens, 

pipeline/generateCatalog.py

@@ -107,6 +107,9 @@ if (startCatalogStage <= 4) and (endCatalogStage >= 4):
   sys.stdout.flush()
 
   for thisDataPortion in dataPortions:
-    plotNumberCounts(workDir, dataSampleList, figDir, showPlot=True, dataPortion=thisDataPortion, setName=setName)
+    plotRedshiftDistribution(workDir, dataSampleList, figDir, showPlot=False, 
+                             dataPortion=thisDataPortion, setName=catalogName)
+    plotSizeDistribution(workDir, dataSampleList, figDir, showPlot=False, 
+                             dataPortion=thisDataPortion, setName=catalogName)
 
 print "\n Done!"

pipeline/prepareCatalogs.py

@@ -7,7 +7,6 @@
 import numpy as np
 import os
 import sys
-from Scientific.IO.NetCDF import NetCDFFile
 import void_python_tools as vp
 import argparse
 
@@ -95,10 +94,10 @@ from void_python_tools.backend.classes import *
 
 continueRun = False
 startCatalogStage = 1
-endCatalogStage   = 3
+endCatalogStage   = 4
                
 startAPStage = 1
-endAPStage = 6
+endAPStage = 7
 
 ZOBOV_PATH = os.getenv("PWD")+"/../zobov/"
 CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/"
@@ -297,7 +296,7 @@ if args.script or args.all:
   for line in inFile: numPart += 1
   inFile.close()
 
-  minRadius = int(np.ceil(lbox/numPart**(1./3.)))
+  minRadius = 2*int(np.ceil(lbox/numPart**(1./3.)))
 
   if dataFormat == "multidark":
     writeScript(prefix, "halos", scriptDir, catalogDir, redshifts, 

pipeline/prepareGadgetCatalog.py

@@ -5,7 +5,6 @@
 import numpy as np
 import os
 import sys
-from Scientific.IO.NetCDF import NetCDFFile
 import void_python_tools as vp
 import argparse
 
@@ -96,10 +95,10 @@ from void_python_tools.backend.classes import *
 
 continueRun = False
 startCatalogStage = 1
-endCatalogStage   = 3
+endCatalogStage   = 4
                
 startAPStage = 1
-endAPStage = 6
+endAPStage = 7
 
 ZOBOV_PATH = os.getenv("PWD")+"/../zobov/"
 CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/"

python_tools/void_python_tools/backend/launchers.py

@@ -76,7 +76,6 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
 
     if os.access("contour_map.fits", os.F_OK):
       os.system("mv %s %s" % ("contour_map.fits", zobovDir))
-      os.system("mv %s %s" % ("mask_map.fits", zobovDir))
 
     if os.access("comoving_distance.txt", os.F_OK):
       os.system("mv %s %s" % ("comoving_distance.txt", zobovDir))
@@ -815,7 +814,7 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
     maxtries = 5
     while badChain:
       Rexpect = (stack.rMin+stack.rMax)/2
-      Rtruncate = stack.rMax*3. + 1 # TEST
+      Rtruncate = stack.rMin*3. + 1 # TEST
       ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
                                     Niter=300000,
                                     Nburn=100000,

python_tools/void_python_tools/plotting/plotTools.py

@@ -1,4 +1,5 @@
-__all__=['plotNumberCounts']
+__all__=['plotRedshiftDistribution', 'plotSizeDistribution', 'plot1dProfiles',
+         'plotMarg1d']
 
 from void_python_tools.backend.classes import *
 from plotDefs import *
@@ -7,11 +8,12 @@ import os
 import pylab as plt
 
 # -----------------------------------------------------------------------------
-def plotNumberCounts(workDir=None, sampleList=None, figDir=None, plotNameBase="numbercount", 
+def plotRedshiftDistribution(workDir=None, sampleList=None, figDir=None, 
+                     plotNameBase="zdist", 
                      showPlot=False, dataPortion=None, setName=None):
 
   plt.clf()
-  plt.xlabel("Comoving Distance (Mpc/h)")
+  plt.xlabel("Redshift")
   plt.ylabel("Number of Voids")
 
   plotTitle = setName
@@ -67,3 +69,149 @@ def plotNumberCounts(workDir=None, sampleList=None, figDir=None, plotNameBase="n
     os.system("display %s" % figDir+"/fig_"+plotName+".png")
 
 
+# -----------------------------------------------------------------------------
+def plotSizeDistribution(workDir=None, sampleList=None, figDir=None, 
+                     plotNameBase="sizedist", 
+                     showPlot=False, dataPortion=None, setName=None):
+
+  plt.clf()
+  plt.xlabel("Void Radius (Mpc/h)")
+  plt.ylabel("Number of Voids")
+
+  plotTitle = setName
+
+  plotName = plotNameBase
+
+  xMin = 1.e00
+  xMax = 0
+
+  for (iSample,sample) in enumerate(sampleList):
+
+    sampleName = sample.fullName
+    lineTitle = sampleName
+
+    filename = workDir+"/sample_"+sampleName+"/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
+
+    xMin = 5
+    xMax = 140
+
+    range = (xMin, xMax)
+    nbins = np.ceil((xMax-xMin)/5)
+
+    #thisMax = np.max(data[:,5])
+    #thisMin = np.min(data[:,5])
+    #if thisMax > xMax: xMax = thisMax
+    #if thisMin < xMin: xMin = thisMin
+
+    plt.hist(data[:,4], bins=nbins,
+             label=lineTitle, color=colorList[iSample],
+             histtype = "step", range=range,
+             linewidth=linewidth)
+
+  plt.legend(title = "Samples", loc = "upper right")
+  plt.title(plotTitle)
+
+  plt.xlim(xMin, xMax)
+  #plt.xlim(xMin, xMax*1.4) # make room for legend
+
+  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 showPlot:
+    os.system("display %s" % figDir+"/fig_"+plotName+".png")
+
+
+# -----------------------------------------------------------------------------
+def plot1dProfiles(workDir=None, sampleList=None, figDir=None, 
+                   plotNameBase="1dprofile", 
+                   showPlot=False, dataPortion=None, setName=None):
+
+  plt.clf()
+  plt.xlabel(r"$R/R_{v,\mathrm{max}}$")
+  plt.ylabel(r"$n / \bar n$")
+
+  for (iSample,sample) in enumerate(sampleList):
+
+    sampleName = sample.fullName
+
+    for (iStack,stack) in enumerate(sample.stacks):
+
+      plotTitle = setName
+      plotName = plotNameBase
+
+      runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, 
+                                 stack.rMax, dataPortion)
+
+      plotTitle = sampleName + ", z = "+str(stack.zMin)+"-"+str(stack.zMax)+", R = "+str(stack.rMin)+"-"+str(stack.rMax)+ r" $h^{-1}$ Mpc"
+
+      filename = workDir+"/sample_"+sampleName+"/stacks_"+runSuffix+\
+                 "profile_1d.txt"
+
+      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
+
+      data[:,1] /= stack.rMax
+      plt.ylim(ymin=0.0, ymax=np.amax(data[:,2])+0.1)
+      plt.xlim(xmin=0.0, xmax=2.1)
+      plt.plot(data[:,1], data[:,2], label=lineTitle, color=colorList[0],
+               linewidth=linewidth)
+
+      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 showPlot:
+        os.system("display %s" % figDir+"/fig_"+plotName+".png")
+
+
+# -----------------------------------------------------------------------------
+def plotMarg1d(workDir=None, sampleList=None, figDir=None, 
+               plotNameBase="marg1d", 
+               showPlot=False, dataPortion=None, setName=None):
+
+    plotNames = ("Om", "w0", "wa")
+    plotTitles = ("$\Omega_M$", "$w_0$", "$w_a$")
+    files = ("Om", "w0", "wa")
+
+    for iPlot in range(len(plotNames)):
+
+      plt.clf()
+
+      plotName = plotNameBase+"_"+plotNames[iPlot]+"_"+dataPortion
+      plotTitle = plotTitles[iPlot]
+      dataFile = workDir + "/likelihoods_"+dataPortion+"_"+files[iPlot]+".dat"
+
+      plt.xlabel(plotTitle, fontsize="20")
+      plt.ylabel("Likelihood", fontsize="20")
+      plt.ylim(0.0, 1.0)
+
+      data = np.loadtxt(dataFile, comments="#")
+      plt.plot(data[:,0], data[:,1], color='k', linewidth=linewidth)
+
+      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 showPlot:
+        os.system("display %s" % figDir+"/fig_"+plotName+".png")
+
+
+