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

- fixed bug where omegaM was always 0.3 regardless of user input

Browse source 
- wrote out flagging of boundary and edge galaxies
- turned off placing of boundary particles
This commit is contained in:
Paul M. Sutter 2024-06-05 17:28:16 +02:00
parent aded7a7c2c
commit f59fee9bf8
4 changed files with 164 additions and 102 deletions
Download patch file Download diff file Expand all files Collapse all files

16
c_source/prep/prepObservation.cpp
View file

@ -113,7 +113,7 @@ void loadData(const string& fname, NYU_VData & data)
// Maubert - avoid double counting of last particle in array if EOF is after a blank line // Maubert - avoid double counting of last particle in array if EOF is after a blank line
if (!f) if (!f)
{ {
continue; continue;
} }
// End Maubert // End Maubert
d.uniqueID = d.index; d.uniqueID = d.index;
@ -273,7 +273,7 @@ void generateSurfaceMask(prepObservation_info& args ,
output_data.Lmax = Rmax; output_data.Lmax = Rmax;
// PMS - write a small text file with galaxy position // PMS - write a small text file with galaxy position (for diagnostic purposes)
FILE *fp; FILE *fp;
fp = fopen("galaxies.txt", "w"); fp = fopen("galaxies.txt", "w");
for (int i = 0; i < data.size(); i++) { for (int i = 0; i < data.size(); i++) {
@ -289,6 +289,9 @@ void generateSurfaceMask(prepObservation_info& args ,
cout << format("Rmax is %g, surface volume is %g") % (Rmax/100) % (volume/(4*M_PI)) << endl; cout << format("Rmax is %g, surface volume is %g") % (Rmax/100) % (volume/(4*M_PI)) << endl;
volume *= Rmax*Rmax*Rmax/3/1e6; volume *= Rmax*Rmax*Rmax/3/1e6;
numToInsert = (int)floor(volume*args.density_fake_arg); numToInsert = (int)floor(volume*args.density_fake_arg);
// TEST NOT USING MOCK PARTICLES
numToIntsert = 0
// END TEST
cout << format("3d volume to fill: %g (Mpc/h)^3") % volume << endl; cout << format("3d volume to fill: %g (Mpc/h)^3") % volume << endl;
cout << format("Will insert %d particles") % numToInsert << endl; cout << format("Will insert %d particles") % numToInsert << endl;
@ -325,7 +328,7 @@ void generateSurfaceMask(prepObservation_info& args ,
} }
while (!stop_here); while (!stop_here);
// PMS : write mock galaxies to a small file for plotting // PMS : write mock galaxies to a small file for diagnostic purposes
fprintf(fp, "%e %e %e\n", fprintf(fp, "%e %e %e\n",
(p.xyz[0]), (p.xyz[0]),
(p.xyz[1]), (p.xyz[1]),
@ -344,7 +347,7 @@ void generateSurfaceMask(prepObservation_info& args ,
fclose(fp); fclose(fp);
// PMS // PMS
// TEST - insert mock galaxies along box edge // TEST - insert mock galaxies along box edge - this is for tesselation safety
fp = fopen("mock_boundary.txt", "w"); fp = fopen("mock_boundary.txt", "w");
double dx[3]; double dx[3];
dx[0] = output_data.box[0][1] - output_data.box[0][0]; dx[0] = output_data.box[0][1] - output_data.box[0][0];
@ -410,8 +413,9 @@ void generateSurfaceMask(prepObservation_info& args ,
rmax = args.zMax_arg * LIGHT_SPEED; rmax = args.zMax_arg * LIGHT_SPEED;
} }
//for (int q = 0; q < 0; q++) { // TEST NOT USING BOUNDARY PARTICLES
for (int q = 0; q < full_mask_list.size(); q++) { for (int q = 0; q < 0; q++) {
//for (int q = 0; q < full_mask_list.size(); q++) {
vec3 v = mask.pix2vec(full_mask_list[q]); vec3 v = mask.pix2vec(full_mask_list[q]);
Position p; Position p;

2
examples/example_observation/example_observation.py
View file

@ -83,7 +83,7 @@ newSample = Sample(
volumeLimited = True, volumeLimited = True,
# HEALpix mask file - set to None to auto-compute # HEALpix mask file - set to None to auto-compute
maskFile = None, maskFile = "",
#maskFile = inputDataDir+"/example_observation_mask.fits", #maskFile = inputDataDir+"/example_observation_mask.fits",
# if maskFile blank, desired resolution for HEALpix # if maskFile blank, desired resolution for HEALpix

147
python_source/backend/launchers.py
View file

@ -47,25 +47,36 @@ LIGHT_SPEED = 299792.458
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
def launchPrep(sample, binPath, workDir=None, inputDataDir=None, def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
zobovDir=None, figDir=None, logFile=None, useComoving=False, outputDir=None, figDir=None, logFile=None, useComoving=False,
continueRun=None, regenerate=False): continueRun=None, regenerate=False):
if sample.dataType == "observation": if sample.dataType == "observation":
sampleName = sample.fullName sampleName = sample.fullName
if regenerate: if regenerate:
inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par" inputParameterFlag = "inputParameter " + outputDir +
outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName "/zobov_slice_"+sampleName+".par"
outputFile = outputDir + "/regenerated_zobov_slice_" + sampleName
else: else:
inputParameterFlag = "" inputParameterFlag = ""
outputFile = zobovDir+"/zobov_slice_" + sampleName outputFile = outputDir + "/zobov_slice_" + sampleName
if sample.dataFile == "": if sample.dataFile == "":
datafile = inputDataDir+"/"+sampleName datafile = inputDataDir+"/"+sampleName
else: else:
datafile = inputDataDir+"/"+sample.dataFile datafile = inputDataDir+"/"+sample.dataFile
maskFile = sample.maskFile if sample.maskFile == "":
maskFile = outputDir + "/constructed_mask.fits"
figureOutMask(dataFile, sample.nsideForMask, maskFile)
else:
maskFile = sample.maskFile
edgeGalFile = outputDir + "/galaxy_edge_flags.txt"
edgeMaskFile = outputDir + "/mask_edge_map.fits"
findEdgeGalaxies(dataFile, maskFile, edgeGalFile, edgeMaskFile,
sample.zBoundary[0], sample.zBoundary[1], sample.omegaM,
useComoving)
if useComoving: if useComoving:
useComovingFlag = "useComoving" useComovingFlag = "useComoving"
@ -84,7 +95,7 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
%s %s
omegaM %g omegaM %g
""" % (datafile, maskFile, outputFile, """ % (datafile, maskFile, outputFile,
zobovDir+"/zobov_slice_"+sampleName+".par", outputDir+"/zobov_slice_"+sampleName+".par",
sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity, sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity,
useComovingFlag, inputParameterFlag, sample.omegaM) useComovingFlag, inputParameterFlag, sample.omegaM)
@ -108,21 +119,21 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
if os.access(parmFile, os.F_OK): os.unlink(parmFile) if os.access(parmFile, os.F_OK): os.unlink(parmFile)
if os.access("contour_map.fits", os.F_OK): if os.access("contour_map.fits", os.F_OK):
os.system("mv %s %s" % ("contour_map.fits", zobovDir)) os.system("mv %s %s" % ("contour_map.fits", outputDir))
os.system("mv %s %s" % ("mask_map.fits", zobovDir)) os.system("mv %s %s" % ("mask_map.fits", outputDir))
if os.access("comoving_distance.txt", os.F_OK): if os.access("comoving_distance.txt", os.F_OK):
os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) os.system("mv %s %s" % ("comoving_distance.txt", outputDir))
if os.access("mask_index.txt", os.F_OK): if os.access("mask_index.txt", os.F_OK):
os.system("mv %s %s" % ("mask_index.txt", zobovDir)) os.system("mv %s %s" % ("mask_index.txt", outputDir))
os.system("mv %s %s" % ("total_particles.txt", zobovDir)) os.system("mv %s %s" % ("total_particles.txt", outputDir))
if os.access("galaxies.txt", os.F_OK): if os.access("galaxies.txt", os.F_OK):
os.system("mv %s %s" % ("galaxies.txt", zobovDir)) os.system("mv %s %s" % ("galaxies.txt", outputDir))
os.system("mv %s %s" % ("mock_galaxies.txt", zobovDir)) os.system("mv %s %s" % ("mock_galaxies.txt", outputDir))
os.system("mv %s %s" % ("mock_boundary.txt", zobovDir)) os.system("mv %s %s" % ("mock_boundary.txt", outputDir))
os.system("mv %s %s" % ("mock_sphere.txt", zobovDir)) os.system("mv %s %s" % ("mock_sphere.txt", outputDir))
else: # simulation else: # simulation
sampleName = sample.fullName sampleName = sample.fullName
@ -142,16 +153,16 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
if prevSubSample == -1: if prevSubSample == -1:
inputParameterFlag = "" inputParameterFlag = ""
outputFile = zobovDir+"/zobov_slice_" + sampleName + "_ss" + \ outputFile = outputDir+"/zobov_slice_" + sampleName + "_ss" + \
thisSubSample thisSubSample
keepFraction = float(thisSubSample) keepFraction = float(thisSubSample)
subSampleLine = "subsample %g" % keepFraction subSampleLine = "subsample %g" % keepFraction
resubSampleLine = "" resubSampleLine = ""
firstSubSample = keepFraction firstSubSample = keepFraction
else: else:
inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+\ inputParameterFlag = "inputParameter " + outputDir+"/zobov_slice_"+\
sampleName+"_ss"+prevSubSample+".par" sampleName+"_ss"+prevSubSample+".par"
outputFile = zobovDir+"/zobov_slice_" + sampleName + "_ss" + \ outputFile = outputDir+"/zobov_slice_" + sampleName + "_ss" + \
thisSubSample thisSubSample
keepFraction = float(thisSubSample)/float(prevSubSample) keepFraction = float(thisSubSample)/float(prevSubSample)
subSampleLine = "subsample %s" % firstSubSample subSampleLine = "subsample %s" % firstSubSample
@ -243,8 +254,8 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
# remove intermediate files # remove intermediate files
if (prevSubSample != -1): if (prevSubSample != -1):
os.unlink(zobovDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample+".par") os.unlink(outputDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample+".par")
os.unlink(zobovDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample) os.unlink(outputDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample)
doneLine = "Done! %5.2e\n" % keepFraction doneLine = "Done! %5.2e\n" % keepFraction
if not jobSuccessful(logFile, doneLine): if not jobSuccessful(logFile, doneLine):
@ -256,27 +267,27 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
if jobSuccessful(logFile, doneLine): print("done") if jobSuccessful(logFile, doneLine): print("done")
# place the final subsample # place the final subsample
os.system("mv %s %s" % (zobovDir+"/zobov_slice_"+sampleName+"_ss"+\ os.system("mv %s %s" % (outputDir+"/zobov_slice_"+sampleName+"_ss"+\
prevSubSample, zobovDir+"/zobov_slice_"+sampleName)) prevSubSample, outputDir+"/zobov_slice_"+sampleName))
os.system("mv %s %s" % (zobovDir+"/zobov_slice_"+sampleName+"_ss"+\ os.system("mv %s %s" % (outputDir+"/zobov_slice_"+sampleName+"_ss"+\
prevSubSample+".par", zobovDir+"/zobov_slice_"+sampleName+".par")) prevSubSample+".par", outputDir+"/zobov_slice_"+sampleName+".par"))
if os.access("comoving_distance.txt", os.F_OK): if os.access("comoving_distance.txt", os.F_OK):
os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) os.system("mv %s %s" % ("comoving_distance.txt", outputDir))
if os.access(parmFile, os.F_OK): if os.access(parmFile, os.F_OK):
os.unlink(parmFile) os.unlink(parmFile)
if os.access("mask_index.txt", os.F_OK): if os.access("mask_index.txt", os.F_OK):
os.system("mv %s %s" % ("mask_index.txt", zobovDir)) os.system("mv %s %s" % ("mask_index.txt", outputDir))
os.system("mv %s %s" % ("total_particles.txt", zobovDir)) os.system("mv %s %s" % ("total_particles.txt", outputDir))
#os.system("mv %s %s" % ("sample_info.txt", zobovDir)) #os.system("mv %s %s" % ("sample_info.txt", outputDir))
# add to sample info file # add to sample info file
if sample.dataType == "observation": if sample.dataType == "observation":
(boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0], (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0],
sample.zRange[1], sample.zRange[0], sample.zRange[1], "all", sample.zRange[1], sample.zRange[0], sample.zRange[1], "all",
useComoving=useComoving) sample.omegaM, useComoving=useComoving)
else: else:
iX = float(sample.mySubvolume[0]) iX = float(sample.mySubvolume[0])
iY = float(sample.mySubvolume[1]) iY = float(sample.mySubvolume[1])
@ -290,16 +301,16 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
boxVol = (xMax-xMin)*(yMax-yMin)*(zMax-zMin) boxVol = (xMax-xMin)*(yMax-yMin)*(zMax-zMin)
nbar = 1.0 nbar = 1.0
numTracers = int(open(zobovDir+"/mask_index.txt", "r").read()) numTracers = int(open(outputDir+"/mask_index.txt", "r").read())
numTotal = int(open(zobovDir+"/total_particles.txt", "r").read()) numTotal = int(open(outputDir+"/total_particles.txt", "r").read())
meanSep = (1.*numTracers/boxVol/nbar)**(-1/3.) meanSep = (1.*numTracers/boxVol/nbar)**(-1/3.)
# save this sample's information # save this sample's information
with open(zobovDir+"/sample_info.dat", mode='wb') as output: with open(outputDir+"/sample_info.dat", mode='wb') as output:
pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL) pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL)
fp = open(zobovDir+"/sample_info.txt", 'w') fp = open(outputDir+"/sample_info.txt", 'w')
fp.write("Sample name: %s\n" % sample.fullName) fp.write("Sample name: %s\n" % sample.fullName)
fp.write("Sample nickname: %s\n" % sample.nickName) fp.write("Sample nickname: %s\n" % sample.nickName)
fp.write("Data type: %s\n" % sample.dataType) fp.write("Data type: %s\n" % sample.dataType)
@ -323,13 +334,13 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
fp.close() fp.close()
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
numZobovDivisions=None, numZobovThreads=None, numZobovDivisions=None, numZobovThreads=None,
mergingThreshold=0.2): mergingThreshold=0.2):
sampleName = sample.fullName sampleName = sample.fullName
datafile = zobovDir+"zobov_slice_"+sampleName datafile = outputDir+"zobov_slice_"+sampleName
logFile = logDir+"/zobov_"+sampleName+".out" logFile = logDir+"/zobov_"+sampleName+".out"
@ -339,8 +350,8 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
os.unlink(vozScript) os.unlink(vozScript)
if sample.dataType == "observation": if sample.dataType == "observation":
maskIndex = open(zobovDir+"/mask_index.txt", "r").read() maskIndex = open(outputDir+"/mask_index.txt", "r").read()
totalPart = open(zobovDir+"/total_particles.txt", "r").read() totalPart = open(outputDir+"/total_particles.txt", "r").read()
maxDen = mergingThreshold*float(maskIndex)/float(totalPart) maxDen = mergingThreshold*float(maskIndex)/float(totalPart)
else: else:
maskIndex = -1 maskIndex = -1
@ -349,18 +360,18 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
print(" WARNING! You are using a single ZOBOV division with a simulation. Periodic boundaries will not be respected!") print(" WARNING! You are using a single ZOBOV division with a simulation. Periodic boundaries will not be respected!")
if not (continueRun and jobSuccessful(logFile, "Done!\n")): if not (continueRun and jobSuccessful(logFile, "Done!\n")):
for fileName in glob.glob(zobovDir+"/part._"+sampleName+".*"): for fileName in glob.glob(outputDir+"/part._"+sampleName+".*"):
os.unlink(fileName) os.unlink(fileName)
if os.access(zobovDir+"/adj_"+sampleName+".dat", os.F_OK): if os.access(outputDir+"/adj_"+sampleName+".dat", os.F_OK):
os.unlink(zobovDir+"adj_"+sampleName+".dat") os.unlink(outputDir+"adj_"+sampleName+".dat")
if os.access(zobovDir+"/voidDesc_"+sampleName+".out", os.F_OK): if os.access(outputDir+"/voidDesc_"+sampleName+".out", os.F_OK):
os.unlink(zobovDir+"/voidDesc_"+sampleName+".out") os.unlink(outputDir+"/voidDesc_"+sampleName+".out")
cmd = [binPath+"/vozinit", datafile, "0.1", "1.0", str(numZobovDivisions), \ cmd = [binPath+"/vozinit", datafile, "0.1", "1.0", str(numZobovDivisions), \
"_"+sampleName, str(numZobovThreads), \ "_"+sampleName, str(numZobovThreads), \
binPath, zobovDir, str(maskIndex)] binPath, outputDir, str(maskIndex)]
log = open(logFile, 'w') log = open(logFile, 'w')
subprocess.call(cmd, stdout=log, stderr=log) subprocess.call(cmd, stdout=log, stderr=log)
log.close() log.close()
@ -375,13 +386,13 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
sample.selFunFile != None: sample.selFunFile != None:
# load volumes # load volumes
volFile = zobovDir+"/vol_"+sampleName+".dat" volFile = outputDir+"/vol_"+sampleName+".dat"
with open(volFile, mode="rb") as File: with open(volFile, mode="rb") as File:
numPartTot = np.fromfile(File, dtype=np.int32,count=1) numPartTot = np.fromfile(File, dtype=np.int32,count=1)
vols = np.fromfile(File, dtype=np.float32,count=numPartTot) vols = np.fromfile(File, dtype=np.float32,count=numPartTot)
# load redshifts # load redshifts
partFile = zobovDir+"/zobov_slice_"+sample.fullName partFile = outputDir+"/zobov_slice_"+sample.fullName
with open(partFile, mode="rb") as File: with open(partFile, mode="rb") as File:
chk = np.fromfile(File, dtype=np.int32,count=1) chk = np.fromfile(File, dtype=np.int32,count=1)
Np = np.fromfile(File, dtype=np.int32,count=1) Np = np.fromfile(File, dtype=np.int32,count=1)
@ -412,7 +423,7 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
redshifts = np.fromfile(File, dtype=np.float32,count=Np) redshifts = np.fromfile(File, dtype=np.float32,count=Np)
chk = np.fromfile(File, dtype=np.int32,count=1) chk = np.fromfile(File, dtype=np.int32,count=1)
# z # redshift
chk = np.fromfile(File, dtype=np.int32,count=1) chk = np.fromfile(File, dtype=np.int32,count=1)
redshifts = np.fromfile(File, dtype=np.float32,count=Np) redshifts = np.fromfile(File, dtype=np.float32,count=Np)
chk = np.fromfile(File, dtype=np.int32,count=1) chk = np.fromfile(File, dtype=np.int32,count=1)
@ -428,29 +439,29 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
for i in range(len(vols)): for i in range(len(vols)):
vols[i] *= selfunc(redshifts[i]) vols[i] *= selfunc(redshifts[i])
volFile = zobovDir+"/vol_weighted_"+sampleName+".dat" volFile = outputDir+"/vol_weighted_"+sampleName+".dat"
with open(volFile, mode='wb') as File: with open(volFile, mode='wb') as File:
numPartTot.astype(np.int32).tofile(File) numPartTot.astype(np.int32).tofile(File)
vols.astype(np.float32).tofile(File) vols.astype(np.float32).tofile(File)
volFileToUse = zobovDir+"/vol_weighted_"+sampleName+".dat" volFileToUse = outputDir+"/vol_weighted_"+sampleName+".dat"
else: else:
volFileToUse = zobovDir+"/vol_"+sampleName+".dat" volFileToUse = outputDir+"/vol_"+sampleName+".dat"
cmd = [binPath+"/jozov2", \ cmd = [binPath+"/jozov2", \
zobovDir+"/adj_"+sampleName+".dat", \ outputDir+"/adj_"+sampleName+".dat", \
volFileToUse, \ volFileToUse, \
zobovDir+"/voidPart_"+sampleName+".dat", \ outputDir+"/voidPart_"+sampleName+".dat", \
zobovDir+"/voidZone_"+sampleName+".dat", \ outputDir+"/voidZone_"+sampleName+".dat", \
zobovDir+"/voidDesc_"+sampleName+".out", \ outputDir+"/voidDesc_"+sampleName+".out", \
str(maxDen), str(maskIndex)] str(maxDen), str(maskIndex)]
log = open(logFile, 'a') log = open(logFile, 'a')
subprocess.call(cmd, stdout=log, stderr=log) subprocess.call(cmd, stdout=log, stderr=log)
log.close() log.close()
# don't need the subbox files # don't need the subbox files
for fileName in glob.glob(zobovDir+"/part._"+sampleName+".*"): for fileName in glob.glob(outputDir+"/part._"+sampleName+".*"):
os.unlink(fileName) os.unlink(fileName)
if jobSuccessful(logFile, "Done!\n"): if jobSuccessful(logFile, "Done!\n"):
@ -467,19 +478,19 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
def launchPrune(sample, binPath, def launchPrune(sample, binPath,
summaryFile=None, logFile=None, zobovDir=None, summaryFile=None, logFile=None, outputDir=None,
continueRun=None, useComoving=False, mergingThreshold=0.2, continueRun=None, useComoving=False, mergingThreshold=0.2,
boundaryTolerance=1.0): boundaryTolerance=1.0):
sampleName = sample.fullName sampleName = sample.fullName
numVoids = sum(1 for line in \ numVoids = sum(1 for line in \
open(zobovDir+"/voidDesc_"+sampleName+".out")) open(outputDir+"/voidDesc_"+sampleName+".out"))
numVoids -= 2 numVoids -= 2
if sample.dataType == "observation": if sample.dataType == "observation":
mockIndex = open(zobovDir+"/mask_index.txt", "r").read() mockIndex = open(outputDir+"/mask_index.txt", "r").read()
totalPart = open(zobovDir+"/total_particles.txt", "r").read() totalPart = open(outputDir+"/total_particles.txt", "r").read()
maxDen = mergingThreshold*float(mockIndex)/float(totalPart) maxDen = mergingThreshold*float(mockIndex)/float(totalPart)
observationLine = " --isObservation" observationLine = " --isObservation"
#periodicLine = " --periodic=''" #periodicLine = " --periodic=''"
@ -497,7 +508,7 @@ def launchPrune(sample, binPath,
if sample.minVoidRadius == -1: if sample.minVoidRadius == -1:
minRadius = -1 minRadius = -1
for line in open(zobovDir+"/sample_info.txt"): for line in open(outputDir+"/sample_info.txt"):
if "Estimated mean tracer separation" in line: if "Estimated mean tracer separation" in line:
minRadius = float(line.split()[5]) minRadius = float(line.split()[5])
break break
@ -510,13 +521,13 @@ def launchPrune(sample, binPath,
if not (continueRun and (jobSuccessful(logFile, "NetCDF: Not a valid ID\n") \ if not (continueRun and (jobSuccessful(logFile, "NetCDF: Not a valid ID\n") \
or jobSuccessful(logFile, "Done!\n"))): or jobSuccessful(logFile, "Done!\n"))):
cmd = binPath cmd = binPath
cmd += " --partFile=" + zobovDir+"/zobov_slice_"+str(sampleName) cmd += " --partFile=" + outputDir+"/zobov_slice_"+str(sampleName)
cmd += " --voidDesc=" + zobovDir+"/voidDesc_"+str(sampleName)+".out" cmd += " --voidDesc=" + outputDir+"/voidDesc_"+str(sampleName)+".out"
cmd += " --void2Zone="+zobovDir+"/voidZone_"+str(sampleName)+".dat" cmd += " --void2Zone="+outputDir+"/voidZone_"+str(sampleName)+".dat"
cmd += " --zone2Part=" + zobovDir+"/voidPart_"+str(sampleName)+".dat" cmd += " --zone2Part=" + outputDir+"/voidPart_"+str(sampleName)+".dat"
cmd += " --partVol=" + zobovDir+"/vol_"+str(sampleName)+".dat" cmd += " --partVol=" + outputDir+"/vol_"+str(sampleName)+".dat"
cmd += " --partAdj=" + zobovDir+"/adj_"+str(sampleName)+".dat" cmd += " --partAdj=" + outputDir+"/adj_"+str(sampleName)+".dat"
cmd += " --extraInfo=" + zobovDir+"/zobov_slice_"+str(sampleName)+\ cmd += " --extraInfo=" + outputDir+"/zobov_slice_"+str(sampleName)+\
".par" ".par"
cmd += " --tolerance=" + str(boundaryTolerance) cmd += " --tolerance=" + str(boundaryTolerance)
cmd += " --mockIndex=" + str(mockIndex) cmd += " --mockIndex=" + str(mockIndex)
@ -529,7 +540,7 @@ def launchPrune(sample, binPath,
cmd += periodicLine cmd += periodicLine
cmd += useComovingFlag cmd += useComovingFlag
cmd += " --omegaM=" + str(sample.omegaM) cmd += " --omegaM=" + str(sample.omegaM)
cmd += " --outputDir=" + zobovDir cmd += " --outputDir=" + outputDir
cmd += " --sampleName=" + str(sampleName) cmd += " --sampleName=" + str(sampleName)
log = open(logFile, 'w') log = open(logFile, 'w')
#log.write(f"Command is {cmd}\n") #log.write(f"Command is {cmd}\n")
@ -635,7 +646,7 @@ def launchVelocityStack(sample, stack, binPath,
velField_file, velField_file,
thisDataPortion=None, logDir=None, thisDataPortion=None, logDir=None,
voidDir=None, runSuffix=None, voidDir=None, runSuffix=None,
zobovDir=None, outputDir=None,
summaryFile=None, summaryFile=None,
continueRun=None, dataType=None, prefixRun=""): continueRun=None, dataType=None, prefixRun=""):

101
python_source/backend/surveyTools.py
View file

@ -29,24 +29,20 @@ __all__=['getSurveyProps', 'getNside', 'figureOutMask', 'findEdgeGalaxies']
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
# returns the volume and galaxy density for a given redshit slice # returns the volume and galaxy density for a given redshit slice
def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectionFuncFile=None, useComoving=False): def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion,
omegaM, selectionFuncFile=None, useComoving=False):
LIGHT_SPEED = 299792.458 #LIGHT_SPEED = 299792.458
mask = healpy.read_map(maskFile) mask = healpy.read_map(maskFile)
area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi
if useComoving: if useComoving:
zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=0.27) zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=omegaM)
zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=0.27) zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=omegaM)
selFunMin = LIGHT_SPEED/100.*angularDiameter(selFunMin, Om=0.27)
selFunMax = LIGHT_SPEED/100.*angularDiameter(selFunMax, Om=0.27)
else: else:
zmin = zmin * 3000 zmin = zmin * LIGHT_SPEED/100.
zmax = zmax * 3000 zmax = zmax * LIGHT_SPEED/100.
selFunMin *= 3000
selFunMax *= 3000
volume = area * (zmax**3 - zmin**3) / 3
if selectionFuncFile == None: if selectionFuncFile == None:
nbar = 1.0 nbar = 1.0
@ -90,40 +86,91 @@ def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectio
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
# returns the nside resolution from the given maskfile # returns the nside resolution from the given maskfile
def getNside(maskFile): def getNside(maskFile):
nside = 1.0
return nside mask = healpy.read_map(maskFile)
return healpy.get_nside(mask)
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
# computes the mask from a given datafile and writes it to a file # helper function to convert RA,dec to phi,theta
def convertAngle(RA, Dec):
phi = np.pi/180.*RA
theta = Dec*np.pi/180.
theta = np.pi/2. - Dec*np.pi/180.
return phi, theta
# -----------------------------------------------------------------------------
# computes the mask from a given galaxy datafile and writes it to a file
def figureOutMask(galFile, nside, outMaskFile): def figureOutMask(galFile, nside, outMaskFile):
npix = healpy.nside2npix(nside) npix = healpy.nside2npix(nside)
mask = np.zeros((npix)) mask = np.zeros((npix))
for line in open(galFile): for line in open(galFile):
line = line.split() line = line.split()
RA = np.float(line[3]) RA = np.float(line[3])
Dec = np.float(line[4]) Dec = np.float(line[4])
z = np.float(line[5]) z = np.float(line[5])
phi = np.pi/180.*RA phi, theta = convertAngle(RA, Dec)
theta = Dec*np.pi/180.
theta = np.pi/2. - Dec*np.pi/180.
pos = np.zeros((3))
pix = healpy.ang2pix(nside, theta, phi) pix = healpy.ang2pix(nside, theta, phi)
mask[pix] = 1. mask[pix] = 1.
healpy.write_map(outMaskFile, mask) healpy.write_map(outMaskFile, mask)
return mask return mask
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
# figures out which galaxies live on a mask edge, and also writes the edge # figures out which galaxies live on a mask edge, and also writes the edge
# map to an auxillary file # map to an auxillary file
def findEdgeGalaxies(galFile, maskFile): def findEdgeGalaxies(galFile, maskFile, edgeGalFile, edgeMaskFile,
zmin, zmax, omegaM, useComoving=False):
if useComoving:
zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=omegaM)
zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=omegaM)
else:
zmin *= LIGHT_SPEED/100.
zmax *= LIGHT_SPEED/100.
mask = healpy.read_map(maskFile)
nside = healpy.get_nside(mask)
npix = healpy.nside2npix(nside)
edgeMask = np.zeros((npix))
edgeFile = open(edgeGalFile, "w")
for line in open(galFile):
line = line.split()
RA = np.float(line[3])
Dec = np.float(line[4])
z = np.float(line[5])
if useComoving:
z = LIGHT_SPEED/100.*angularDiameter(z, Om=omegaM)
else:
z *= LIGHT_SPEED/100.
phi, theta = convertAngle(RA, Dec)
isOnEdge = False
# check the mask edges
neighbors = healpy.get_all_neighbors(nside, theta, phi)
isOnEdge = any(p == 0 for p in neighbors):
# check the redshift boundaries
if isOnEdge:
edgeFile.write("1")
else:
edgeFile.write("0")
edgeFile.close()
healpy.write_map(edgeMaskFile, edgeMask)
return return