guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-05 07:41:11 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

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

Browse source
This commit is contained in:
Guilhem Lavaux 2012-11-06 15:42:18 -05:00
commit a59b9aae84
11 changed files with 223 additions and 1526 deletions
Download patch file Download diff file Expand all files Collapse all files

18
GetQhull.cmake
View file

@ -1,18 +0,0 @@
SET(QHULL_BASE_PATH CACHE PATH "Qhull base path")
find_path(QHULL_INCLUDE_PATH qhull_a.h HINTS ${QHULL_BASE_PATH}/src/libqhull)
find_path(QHULL_CPP_INCLUDE_PATH Qhull.h HINTS ${QHULL_BASE_PATH}/src/libqhullcpp)
find_library(QHULL_LIBRARY qhull_p HINTS ${QHULL_BASE_PATH}/lib)
find_library(QHULL_CPP_LIBRARY qhullcpp HINTS ${QHULL_BASE_PATH}/lib)
find_library(QHULL_P_LIBRARY qhullstatic_p HINTS ${QHULL_BASE_PATH}/lib)
if ((NOT QHULL_INCLUDE_PATH) OR (NOT QHULL_CPP_LIBRARY))
message(SEND_ERROR "Qhull library not found")
endif((NOT QHULL_INCLUDE_PATH) OR (NOT QHULL_CPP_LIBRARY))
set(QHULL_INCLUDES ${QHULL_INCLUDE_PATH} ${QHULL_INCLUDE_PATH}/.. ${QHULL_CPP_INCLUDE_PATH} ${QHULL_BASE_PATH}/src)
set(QHULL_LIBRARIES ${QHULL_CPP_LIBRARY} ${QHULL_P_LIBRARY})
add_definitions(-Dqh_QHpointer)
mark_as_advanced(QHULL_INCLUDE_PATH QHULL_CPP_INCLUDE_PATH QHULL_LIBRARY QHULL_CPP_LIBRARY QHULL_P_LIBRARY)

4
c_tools/libzobov/voidTree.hpp
View file

@ -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;
}

1
c_tools/mock/generateFromCatalog.cpp
View file

@ -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();
}

1
c_tools/mock/generateMock.cpp
View file

@ -391,6 +391,7 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info&
f.add_att("range_y_max", ranges[1][1]);
f.add_att("range_z_min", ranges[2][0]);
f.add_att("range_z_max", ranges[2][1]);
f.add_att("mask_index", -1);
NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim);

29
c_tools/stacking/pruneVoids.cpp
View file

@ -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,

8
pipeline/generateCatalog.py
View file

@ -23,9 +23,8 @@ if (len(sys.argv) > 1):
filename = sys.argv[1]
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file!"
print " Cannot find parameter file %s!" % filename
exit(-1)
#parms = __import__(filename[:-3], globals(), locals(), ['*'])
parms = imp.load_source("name", filename)
globals().update(vars(parms))
else:
@ -108,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!"

13
pipeline/prepareCatalogs.py
View file

@ -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/"
@ -196,8 +195,10 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False)
mySubvolume = "%d%d" % (iX, iY)
sampleName = getSampleName(prefix, base, redshift, useVel,
iSlice=iSlice, iVol=mySubvolume)
sampleName = getSampleName(prefix, base, sliceMin, useVel,
iSlice=-1, iVol=mySubvolume)
#sampleName = getSampleName(prefix, base, redshift, useVel,
# iSlice=iSlice, iVol=mySubvolume)
scriptFile.write(sampleInfo.format(dataFile=dataFileName,
dataFormat=dataFormat,
@ -295,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,

5
pipeline/prepareGadgetCatalog.py
View file

@ -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/"

1453
python_tools/void_python_tools/backend/backend.py
View file

File diff suppressed because it is too large Load diff

63
python_tools/void_python_tools/backend/launchers.py
View file

@ -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))
@ -150,6 +149,9 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
else:
print "already done!"
if os.access("comoving_distance.txt", os.F_OK):
os.system("mv %s %s" % ("comoving_distance.txt", zobovDir))
if os.access(parmFile, os.F_OK):
os.unlink(parmFile)
@ -460,34 +462,37 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
return
# figure out box volume and average density
maskFile = sample.maskFile
sulFunFile = sample.selFunFile
if sample.dataType == "observation":
maskFile = sample.maskFile
sulFunFile = sample.selFunFile
if not os.access(sample.selFunFile, os.F_OK) and not volumeLimited:
print " Cannot find", selFunFile, "!"
exit(-1)
if not os.access(sample.selFunFile, os.F_OK) and not sample.volumeLimited:
print " Cannot find", selFunFile, "!"
exit(-1)
sys.stdout = open(logFile, 'a')
sys.stderr = open(logFile, 'a')
zMin = sample.zRange[0]
zMax = sample.zRange[1]
if not sample.volumeLimited:
props = vp.getSurveyProps(maskFile, stack.zMin,
stack.zMax, zMin, zMax, "all",
selectionFuncFile=sample.selFunFile)
sys.stdout = open(logFile, 'a')
sys.stderr = open(logFile, 'a')
zMin = sample.zRange[0]
zMax = sample.zRange[1]
if not sample.volumeLimited:
props = vp.getSurveyProps(maskFile, stack.zMin,
stack.zMax, zMin, zMax, "all",
selectionFuncFile=sample.selFunFile)
else:
zMinForVol = sample.zBoundary[0]
zMaxForVol = sample.zBoundary[1]
props = vp.getSurveyProps(maskFile, zMinForVol,
zMaxForVol, zMin, zMax, "all")
sys.stdout = sys.__stdout__
sys.stderr = sys.__stderr__
boxVol = props[0]
nbar = props[1]
if sample.volumeLimited:
nbar = 1.0
else:
zMinForVol = sample.zBoundary[0]
zMaxForVol = sample.zBoundary[1]
props = vp.getSurveyProps(maskFile, zMinForVol,
zMaxForVol, zMin, zMax, "all")
sys.stdout = sys.__stdout__
sys.stderr = sys.__stderr__
boxVol = props[0]
nbar = props[1]
if sample.volumeLimited:
nbar = 1.0
boxVol = sample.boxLen**3
summaryLine = runSuffix + " " + \
thisDataPortion + " " + \
@ -809,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,
@ -1173,7 +1178,11 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
voidDir = sample.zobovDir+"/stacks_" + runSuffix
centersFile = voidDir+"/centers.txt"
if os.access(centersFile, os.F_OK):
voidRedshifts = np.loadtxt(centersFile)[:,5]
voidRedshifts = np.loadtxt(centersFile)
if voidRedshifts.ndim > 1:
voidRedshifts = voidRedshifts[:,5]
else:
voidRedshifts = voidRedshifts[5]
#fp.write(str(len(voidRedshifts))+" ")
np.savetxt(fp, voidRedshifts[None])
else:

154
python_tools/void_python_tools/plotting/plotTools.py
View file

@ -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")