From 61acd5938330ff3f9ec6e6114f0b00022427b5ff Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Tue, 6 Nov 2012 14:41:09 -0600 Subject: [PATCH] futher fixes for handling sub-volumes in analysis --- c_tools/libzobov/voidTree.hpp | 4 +- c_tools/mock/generateFromCatalog.cpp | 1 + c_tools/stacking/pruneVoids.cpp | 29 ++-- pipeline/generateCatalog.py | 5 +- pipeline/prepareCatalogs.py | 7 +- pipeline/prepareGadgetCatalog.py | 5 +- .../void_python_tools/backend/launchers.py | 3 +- .../void_python_tools/plotting/plotTools.py | 154 +++++++++++++++++- 8 files changed, 182 insertions(+), 26 deletions(-) diff --git a/c_tools/libzobov/voidTree.hpp b/c_tools/libzobov/voidTree.hpp index a316027..bb825dc 100644 --- a/c_tools/libzobov/voidTree.hpp +++ b/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; } diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index f7efcc7..febc6d3 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/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(); } diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 4580e60..d6d5890 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/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, diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 025bbc2..7f95f66 100755 --- a/pipeline/generateCatalog.py +++ b/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!" diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index f45e31b..af9d997 100755 --- a/pipeline/prepareCatalogs.py +++ b/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, diff --git a/pipeline/prepareGadgetCatalog.py b/pipeline/prepareGadgetCatalog.py index 2fc83f3..f0548fa 100755 --- a/pipeline/prepareGadgetCatalog.py +++ b/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/" diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 92536cd..e75343e 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/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, diff --git a/python_tools/void_python_tools/plotting/plotTools.py b/python_tools/void_python_tools/plotting/plotTools.py index 4fad79a..a7c5058 100644 --- a/python_tools/void_python_tools/plotting/plotTools.py +++ b/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") + + +