From 4d853bbe09cb3501b8bd7004dea88a712eeca5d9 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 23 Jan 2014 17:04:46 -0600 Subject: [PATCH 01/12] added density thresholding for halos and galaxies in preparation scripts --- python_tools/pipeline_source/defaults.py | 7 +- .../pipeline_source/prepareCatalogs.in.py | 237 +++++++++++++++--- 2 files changed, 214 insertions(+), 30 deletions(-) diff --git a/python_tools/pipeline_source/defaults.py b/python_tools/pipeline_source/defaults.py index 2a3de54..1b9a52d 100644 --- a/python_tools/pipeline_source/defaults.py +++ b/python_tools/pipeline_source/defaults.py @@ -111,6 +111,9 @@ haloFileDummy = 'NNNNN' # use "none" to get all halos minHaloMasses = [1.2e13] +# density threshold for halo catalogs +haloDenList = [0.0002] + # locations of data in the halo catalog haloFileMCol = 6 haloFileXCol = 0 @@ -135,7 +138,6 @@ lbox = 999.983 # Mpc/h omegaM = 0.2847979853038958 hubble = 0.6962 -#galDens = 0.000225 hodParmList = [ {'name' : "dr9mid", #BOSS: Manera et al. 2012, eq. 26 'Mmin' : 0.0, @@ -143,7 +145,8 @@ hodParmList = [ 'sigma_logM' : 0.596, 'alpha' : 1.0127, 'Mcut' : 1.19399e13, - 'galDens' : 0.0002, + 'galDens' : 0.0002, # density passed to HOD code + 'galDensFinal' : 0.0002, # subsample galaxies to reach this density }, ] diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index e5e6c87..b6d62b0 100644 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/python_tools/pipeline_source/prepareCatalogs.in.py @@ -549,7 +549,7 @@ for iSubSample in xrange(len(subSamples)): prevSubSample = thisSubSample # ----------------------------------------------------------------------------- -# now halos +# now halos - filter on mass if (args.script or args.all) and haloFileBase != "": for minHaloMass in minHaloMasses: @@ -612,8 +612,8 @@ if (args.script or args.all) and haloFileBase != "": numSubvolumes, numSlices, True, lbox, minRadius, omegaM, dataFileNameList = fileList) -if (args.halos or args.all) and haloFileBase != "": - print " Doing halos" +if (args.halos or args.all) and haloFileBase != "" and len(minHaloMasses) > 0: + print " Doing halos - mass" sys.stdout.flush() for minHaloMass in minHaloMasses: @@ -664,11 +664,20 @@ if (args.halos or args.all) and haloFileBase != "": outFile.close() if dataFormat == "sdf": + idTag = "id" + for line in open(dataFile): + if "int64_t id" in line: + idTag = "id" + break + if "int64_t ident" in line: + idTag = "ident" + break + SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" if minHaloMass == "none": - command = "%s -a 200000 %s mass ident x y z vz vy vx parent_id | awk '{if ($9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, outFileName ) + command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | awk '{if ($9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, idTag, dataFile, outFileName ) else: - command = "%s -a 200000 %s mass ident x y z vz vy vx parent_id | awk '{if ($1>%g && $9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, minHaloMass, outFileName ) + command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | awk '{if ($1>%g && $9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, idTag, dataFile, minHaloMass, outFileName ) #os.system(command) subprocess.call(command, shell=True) outFile = open(outFileName, 'a') @@ -697,6 +706,154 @@ if (args.halos or args.all) and haloFileBase != "": outFile.close() inFile.close() +# ----------------------------------------------------------------------------- +# now halos - filter on density +if (args.script or args.all) and haloFileBase != "": + + for haloDen in haloDenList: + print " Doing halo script", haloDen + sys.stdout.flush() + + # estimate number of halos to get density + if haloFileDummy == '': + dataFile = catalogDir+haloFileBase+fileNums[0] + else: + dataFile = catalogDir+haloFileBase.replace(haloFileDummy, + fileNums[0]) + + numPart = haloDen * lbox**3 + meanDen = lbox/numPart**(1./3.) + if meanDen < 1: + minRadius = 4 + else: + minRadius = int(np.ceil(meanDen)) + + setName = prefix+"halos_den"+str(haloDen) + fileList = [] + for (iRedshift, redshift) in enumerate(redshifts): + sampleName = getSampleName(setName, redshift, False) + outFileName = sampleName+".dat" + fileList.append(outFileName) + + writeScript(setName, prefix+"halos_den"+str(haloDen)+"_z", "multidark", + scriptDir, catalogDir, fileNums, + redshifts, + numSubvolumes, numSlices, False, lbox, minRadius, omegaM, + dataFileNameList = fileList) + + if doPecVel: + writeScript(setName, prefix+"halos_den"+str(haloDen)+"_z", + "multidark", + scriptDir, catalogDir, fileNums, + redshifts, + numSubvolumes, numSlices, True, lbox, minRadius, omegaM, + dataFileNameList = fileList) + +if (args.halos or args.all) and haloFileBase != "": + print " Doing halos - density" + sys.stdout.flush() + + for haloDen in haloDenList: + print " halo density = ", haloDen + sys.stdout.flush() + + for (iRedshift, redshift) in enumerate(redshifts): + print " z = ", redshift + sys.stdout.flush() + + if haloFileDummy == '': + dataFile = catalogDir+haloFileBase+fileNums[iRedshift] + else: + dataFile = catalogDir+haloFileBase.replace(haloFileDummy, + fileNums[iRedshift]) + # check the actual density of halos + inFile = open(dataFile, 'r') + numPart = 0 + if dataFormat == "sdf": + for line in inFile: + if "nhalos" in line: + numPart = int(line.split()[3].strip(';')) + break + if "npart" in line: + numPart = int(line.split()[3].strip(';')) + break + inFile.close() + else: + for (iLine, line) in enumerate(inFile): + if iLine < haloFileNumComLines: continue + numPart += 1 + inFile.close() + + numPartExpect = int(np.ceil(haloDen * lbox**3)) + if numPart < numPartExpect: + print " ERROR: not enough halos to support that density! Maximum is %g" % (1.*numPart / lbox**3) + exit(-1) + + sampleName = prefix+"halos_den"+str(haloDen)+"_z"+redshifts[iRedshift] + outFileName = catalogDir+"/"+sampleName+".dat" + outFile = open(outFileName, 'w') + outFile.write("%f\n" %(lbox)) + outFile.write("%s\n" %(omegaM)) + outFile.write("%s\n" %(hubble)) + outFile.write("%s\n" %(redshift)) + outFile.write("%d\n" %(numPartExpect)) + outFile.close() + + if dataFormat == "sdf": + idTag = "id" + for line in open(dataFile): + if "int64_t id" in line: + idTag = "id" + break + if "int64_t ident" in line: + idTag = "ident" + break + + SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" + command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | sort -rg -k 1 | awk '{if ($9==-1 && NR < %d) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, idTag, numPartExpect, outFileName ) + #os.system(command) + subprocess.call(command, shell=True) + outFile = open(outFileName, 'a') + outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n") + outFile.close() + else: + inFile = open(dataFile, 'r') + haloList = [] + for (iHalo,line) in enumerate(inFile): + if iHalo < haloFileNumComLines: continue + line = line.split(haloFileColSep) + x = float(line[haloFileXCol]) * haloFilePosRescale + y = float(line[haloFileYCol]) * haloFilePosRescale + z = float(line[haloFileZCol]) * haloFilePosRescale + vz = float(line[haloFileVZCol]) + vy = float(line[haloFileVYCol]) + vx = float(line[haloFileVXCol]) + mass = float(line[haloFileMCol]) + haloList.append((x,y,z,vz,vy,vx,mass)) + inFile.close() + + haloList = np.array(haloList) + haloMasses = haloList[:,6] + haloList = haloList[np.argsort(haloMasses)[::-1]] + + outFile = open(outFileName, 'a') + for iHalo in xrange(numPartExpect): + x = haloList[iHalo,0] + y = haloList[iHalo,1] + z = haloList[iHalo,2] + vz = haloList[iHalo,3] + vy = haloList[iHalo,4] + vx = haloList[iHalo,5] + mass = haloList[iHalo,6] + + # write to output file + outFile.write("%d %e %e %e %e %e %e %e\n" %(iHalo,x,y,z, + vz,vy,vx,mass)) + + outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n") + outFile.close() + + # ----------------------------------------------------------------------------- # now the HOD parFileText = """ @@ -754,27 +911,13 @@ if (args.script or args.all) and haloFileBase != "": print " ", thisHod['name'] # estimate number of halos to get density - if haloFileDummy == '': - dataFile = catalogDir+haloFileBase+fileNums[0] + numPart = thisHod['galDensFinal'] * lbox**3 + meanDen = lbox/numPart**(1./3.) + if meanDen < 1: + minRadius = 4 else: - dataFile = catalogDir+haloFileBase.replace(haloFileDummy, - fileNums[0]) - if dataFormat == "sdf": - numPart = 0 - SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" - command = "%s -a 200000 %s mass parent_id | awk '{if ($1>%g && $2==-1) print $1}' | wc" % (SDFcvt_PATH, dataFile, thisHod['Mcut']) - numPart = subprocess.check_output(command, shell=True) - numPart = int(numPart.split()[0]) - else: - inFile = open(dataFile, 'r') - for (iHalo,line) in enumerate(inFile): - if iHalo < haloFileNumComLines: continue - line = line.split(haloFileColSep) - if float(line[haloFileMCol]) > thisHod['Mcut']: numPart += 1 - inFile.close() - - minRadius = int(np.ceil(float(lbox)/numPart**(1./3.))) - + minRadius = int(np.ceil(meanDen)) + setName = prefix+"hod_"+thisHod['name'] writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark", scriptDir, catalogDir, fileNums, redshifts, @@ -824,6 +967,7 @@ if (args.hod or args.all) and haloFileBase != "": sigma_logM=thisHod['sigma_logM'], alpha=thisHod['alpha'], Mcut=thisHod['Mcut'], + #galden=-1, galden=thisHod['galDens'], haloFile=haloFile, haloFileFormat=dataFormat, @@ -857,9 +1001,46 @@ if (args.hod or args.all) and haloFileBase != "": print line exit(-1) - outFileName = catalogDir+"/"+sampleName+".dat" - os.system("mv %s/hod_%s.mock %s" % (catalogDir, sampleName, outFileName)) - os.system("rm %s/hod.*" % catalogDir) + # now randomly subsample the galaxies to get desired density + inFile = open(catalogDir+"/hod_"+sampleName+".mock") + numPartExpect = thisHod['galDensFinal'] * lbox**3 + numPartActual = 0 + actualDen = 0 + for (iLine,line) in enumerate(inFile): + numPartActual += 1 + inFile.close() + + if numPartActual < numPartExpect: + print " ERROR: not enough halos to support that density! Maximum is %g" % (1.*numPartActual / lbox**3) + exit(-1) + + actualDen = 1.*numPartActual / lbox**3 + keepFraction = float(thisHod['galDensFinal']) / actualDen + + inFile = open(catalogDir+"/hod_"+sampleName+".mock") + outFile = open(catalogDir+"/"+sampleName+".dat", 'w') + + outFile.write("%f\n" %(lbox)) + outFile.write("%s\n" %(omegaM)) + outFile.write("%s\n" %(hubble)) + outFile.write("%s\n" %(redshift)) + outFile.write("%d\n" %(numPartExpect)) + + numKept = 0 + for (iLine,line) in enumerate(inFile): + if np.random.uniform() > keepFraction: continue + outFile.write(line) + numKept += 1 + + inFile.close() + os.unlink(catalogDir+"/hod_"+sampleName+".mock") + + outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n") + outFile.close() + + + #os.system("mv %s/hod_%s.mock %s" % (catalogDir, sampleName, outFileName)) + #os.system("rm %s/hod.*" % catalogDir) os.system("rm ./hod.par") os.system("rm ./hod-usedvalues") From f39351bd3c157108a9d098b0f876f0e489db0b78 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 23 Jan 2014 20:11:29 -0600 Subject: [PATCH 02/12] fixed overflow with long filenames --- c_tools/hod/populate_simulation.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/c_tools/hod/populate_simulation.c b/c_tools/hod/populate_simulation.c index 20438a7..43aa51b 100644 --- a/c_tools/hod/populate_simulation.c +++ b/c_tools/hod/populate_simulation.c @@ -63,7 +63,7 @@ void populate_simulation() FILE *fp,*fpa[9],*fp2,*fpb[9],*fpc[9],*fps[9],*fpt; int i,j,k,n,imass,n1,j_start=0,i1,galcnt[1000],halocnt[1000],imag; double mass,xg[3],vg[3],nsat,nc[10],ncen,mlo,mag,err1,err2,r,fac,sigv; - char aa[1000]; + char aa[5000]; float x1,xh[3],vh[3],vgf[3]; long IDUM3 = -445; @@ -74,7 +74,7 @@ void populate_simulation() int ngal,nsati[9],ALL_FILES=0,TRACK_GALAXIES=0,WARREN_MASS_CORRECTION=0,haloid; int nadded; FILE *fpBuh; - char buhFile[1000]; + char buhFile[5000]; nadded = 0; @@ -565,7 +565,7 @@ void calc_nbody_two_halo(float **gal, int *id, int ngal) static float *rad, *xi; static int nr, flag = 0; - char fname[100]; + char fname[5000]; FILE *fp; sprintf(fname,"%s.nbody_2halo",Task.root_filename); From f91b565d6fc83fd393cf0c10e3b1ecc035cf79bd Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 23 Jan 2014 21:05:07 -0600 Subject: [PATCH 03/12] fixed some bugs in density handling; now more flexible approach to SDF files --- python_tools/pipeline_source/defaults.py | 4 +- .../pipeline_source/prepareCatalogs.in.py | 100 ++++++++++++++++-- 2 files changed, 95 insertions(+), 9 deletions(-) diff --git a/python_tools/pipeline_source/defaults.py b/python_tools/pipeline_source/defaults.py index 1b9a52d..422797c 100644 --- a/python_tools/pipeline_source/defaults.py +++ b/python_tools/pipeline_source/defaults.py @@ -109,10 +109,10 @@ haloFileDummy = 'NNNNN' # minimum halo mass cuts to apply for the halo catalog # use "none" to get all halos -minHaloMasses = [1.2e13] +minHaloMasses = [] # density threshold for halo catalogs -haloDenList = [0.0002] +haloDenList = [] # locations of data in the halo catalog haloFileMCol = 6 diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index b6d62b0..e80dfd6 100644 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/python_tools/pipeline_source/prepareCatalogs.in.py @@ -465,6 +465,7 @@ for iSubSample in xrange(len(subSamples)): rescale_position = hubble/1000./scale shift = lbox/2. rescale_velocity = 3.08567802e16/3.1558149984e16 + command = "%s -a 200000 %s x y z vz vy vx mass | awk '{print $1*%g+%g, $2*%g+%g, $3*%g+%g, $4*%g, $5*%g, $6*%g, $7}' > %s" % (SDFcvt_PATH, dataFile, rescale_position, shift, @@ -565,10 +566,31 @@ if (args.script or args.all) and haloFileBase != "": numPart = 0 if dataFormat == "sdf": SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" + massTag = "mass" + for line in open(dataFile): + if "float m200b" in line: + massTag = "m200b" + break + if "float mass" in line: + massTag = "mass" + break + + pidTag = "" + iLine = 0 + for line in open(dataFile): + iLine += 1 + if iLine > 100: break + if "parent_id" in line: + pidTag = "parent_id" + break + if "pid" in line: + pidTag = "pid" + break + if minHaloMass == "none": - command = "%s -a 200000 %s mass parent_id | awk '{if ($2==-1) print $1}' | wc" % (SDFcvt_PATH, dataFile) + command = "%s -a 200000 %s %s %s | awk '{if ($2==-1) print $1}' | wc" % (SDFcvt_PATH, dataFile, massTag, pidTag) else: - command = "%s -a 200000 %s mass parent_id | awk '{if ($1>%g && $2==-1) print $1}' | wc" % (SDFcvt_PATH, dataFile, minHaloMass) + command = "%s -a 200000 %s %s %s | awk '{if ($1>%g && $2==-1) print $1}' | wc" % (SDFcvt_PATH, dataFile, massTag, pidTag, minHaloMass) numPart = subprocess.check_output(command, shell=True) numPart = int(numPart.split()[0]) else: @@ -673,11 +695,32 @@ if (args.halos or args.all) and haloFileBase != "" and len(minHaloMasses) > 0: idTag = "ident" break + massTag = "mass" + for line in open(dataFile): + if "float m200b" in line: + massTag = "m200b" + break + if "float mass" in line: + massTag = "mass" + break + + pidTag = "" + iLine = 0 + for line in open(dataFile): + iLine += 1 + if iLine > 100: break + if "parent_id" in line: + pidTag = "parent_id" + break + if "pid" in line: + pidTag = "pid" + break + SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" if minHaloMass == "none": - command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | awk '{if ($9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, idTag, dataFile, outFileName ) + command = "%s -a 200000 %s %s %s x y z vz vy vx %s | awk '{if ($9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, massTag, idTag, pidTag, outFileName ) else: - command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | awk '{if ($1>%g && $9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, idTag, dataFile, minHaloMass, outFileName ) + command = "%s -a 200000 %s %s %s x y z vz vy vx %s | awk '{if ($1>%g && $9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, massTag, idTag, minHaloMass, pidTag, outFileName ) #os.system(command) subprocess.call(command, shell=True) outFile = open(outFileName, 'a') @@ -809,8 +852,29 @@ if (args.halos or args.all) and haloFileBase != "": idTag = "ident" break + massTag = "mass" + for line in open(dataFile): + if "float m200b" in line: + massTag = "m200b" + break + if "float mass" in line: + massTag = "mass" + break + + pidTag = "" + iLine = 0 + for line in open(dataFile): + iLine += 1 + if iLine > 100: break + if "parent_id" in line: + pidTag = "parent_id" + break + if "pid" in line: + pidTag = "pid" + break + print "TEST", massTag, pidTag SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" - command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | sort -rg -k 1 | awk '{if ($9==-1 && NR < %d) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, idTag, numPartExpect, outFileName ) + command = "%s -a 200000 %s %s %s x y z vz vy vx %s | sort -rg -k 1 | awk '{if ($9==-1 && NR < %d) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, massTag, idTag, pidTag, numPartExpect, outFileName ) #os.system(command) subprocess.call(command, shell=True) outFile = open(outFileName, 'a') @@ -945,8 +1009,30 @@ if (args.hod or args.all) and haloFileBase != "": if dataFormat == "sdf": inFile = haloFile outFile = haloFile+"_temp" + + massTag = "mass" + for line in open(inFile): + if "float m200b" in line: + massTag = "m200b" + break + if "float mass" in line: + massTag = "mass" + break + + pidTag = "" + iLine = 0 + for line in open(inFile): + iLine += 1 + if iLine > 100: break + if "parent_id" in line: + pidTag = "parent_id" + break + if "pid" in line: + pidTag = "pid" + break + SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" - command = "%s -a 200000 %s mass x y z vx vy vz parent_id | awk '{if ($8 ==-1) print $1, $2, $3, $4, $5, $6, $7}'>>%s" % (SDFcvt_PATH, inFile, outFile) + command = "%s -a 200000 %s %s x y z vx vy vz %s | awk '{if ($8 ==-1) print $1, $2, $3, $4, $5, $6, $7}'>>%s" % (SDFcvt_PATH, inFile, massTag, pidTag, outFile) #os.system(command) subprocess.call(command, shell=True) haloFile = outFile @@ -1011,7 +1097,7 @@ if (args.hod or args.all) and haloFileBase != "": inFile.close() if numPartActual < numPartExpect: - print " ERROR: not enough halos to support that density! Maximum is %g" % (1.*numPartActual / lbox**3) + print " ERROR: not enough galaxies to support that density! Maximum is %g" % (1.*numPartActual / lbox**3) exit(-1) actualDen = 1.*numPartActual / lbox**3 From d3554b07324f7f0a51f3329cdb9d0a5190904cb4 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Mon, 27 Jan 2014 12:23:59 -0600 Subject: [PATCH 04/12] storing velocity components in the netCDF file --- c_tools/mock/generateMock.cpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index 8192f42..ddf7b19 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -364,11 +364,11 @@ void createBox(SimuData *simu, vector& targets, vector& snapshot_spl // PMS FILE *fp = fopen("mask_index.txt", "w"); - fprintf(fp, "%d", boxed->NumPart); + fprintf(fp, "%ld", boxed->NumPart); fclose(fp); fp = fopen("total_particles.txt", "w"); - fprintf(fp, "%d", boxed->NumPart); + fprintf(fp, "%ld", boxed->NumPart); fclose(fp); printf("Done!\n"); // END PMS @@ -430,6 +430,9 @@ void saveBox(SimuData *&boxed, const std::string& outbox, generateMock_info& arg long *snapshot_split = boxed->as("snapshot_split"); int num_snapshots = *boxed->as("num_snapshots"); long *uniqueID = boxed->as("uniqueID"); + float *velX = boxed->Vel[0]; + float *velY = boxed->Vel[1]; + float *velZ = boxed->Vel[2]; if (!f.is_valid()) { @@ -469,6 +472,13 @@ void saveBox(SimuData *&boxed, const std::string& outbox, generateMock_info& arg v5->put(tmp_int, boxed->NumPart); delete[] tmp_int; } + + NcVar *v6 = f.add_var("vel_x", ncFloat, NumSnap_dim); + NcVar *v7 = f.add_var("vel_y", ncFloat, NumSnap_dim); + NcVar *v8 = f.add_var("vel_z", ncFloat, NumSnap_dim); + v6->put(velX, boxed->NumPart); + v7->put(velY, boxed->NumPart); + v8->put(velZ, boxed->NumPart); } void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& args_info) From b466730c426838cf440290a0a21113a242c0399b Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Mon, 27 Jan 2014 20:26:27 -0600 Subject: [PATCH 05/12] added utility to read off particle velocities --- .../void_python_tools/partUtil/partUtil.py | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/python_tools/void_python_tools/partUtil/partUtil.py b/python_tools/void_python_tools/partUtil/partUtil.py index 7444939..ba16b8a 100644 --- a/python_tools/void_python_tools/partUtil/partUtil.py +++ b/python_tools/void_python_tools/partUtil/partUtil.py @@ -105,6 +105,30 @@ def loadPart(workDir, sampleDir, sample): return partData, boxLen, volNorm, isObservationData +# ----------------------------------------------------------------------------- +def loadPartVel(workDir, sampleDir, sample): + #print " Loading particle velocities..." + sys.stdout.flush() + + infoFile = workDir+"/"+sampleDir+"/zobov_slice_"+sample.fullName+".par" + File = NetCDFFile(infoFile, 'r') + isObservation = getattr(File, 'is_observation') + + if isObservation: + print "No velocities for observations!" + return -1 + + vx = File.variables['vel_x'][0:] + vy = File.variables['vel_y'][0:] + vx = File.variables['vel_z'][0:] + + File.close() + + partVel = np.column_stack((vx,vy,vz)) + + return partVel + + # ----------------------------------------------------------------------------- def shiftPart(inPart, newCenter, periodicLine, boxLen): From 479b9b84227f3f8de9c9c7085de17354ff732bff Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Wed, 29 Jan 2014 10:08:57 -0600 Subject: [PATCH 06/12] fixed bug in velocity output --- c_tools/mock/generateMock.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index ddf7b19..5e93c3f 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -473,9 +473,9 @@ void saveBox(SimuData *&boxed, const std::string& outbox, generateMock_info& arg delete[] tmp_int; } - NcVar *v6 = f.add_var("vel_x", ncFloat, NumSnap_dim); - NcVar *v7 = f.add_var("vel_y", ncFloat, NumSnap_dim); - NcVar *v8 = f.add_var("vel_z", ncFloat, NumSnap_dim); + NcVar *v6 = f.add_var("vel_x", ncFloat, NumPart_dim); + NcVar *v7 = f.add_var("vel_y", ncFloat, NumPart_dim); + NcVar *v8 = f.add_var("vel_z", ncFloat, NumPart_dim); v6->put(velX, boxed->NumPart); v7->put(velY, boxed->NumPart); v8->put(velZ, boxed->NumPart); From c27b09b08defd4b37d4e6a1c5fa1a706dd0b6f49 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 30 Jan 2014 21:39:44 -0600 Subject: [PATCH 07/12] multidark loader now take all velocity components; velocity output bug fixed --- c_tools/mock/generateMock.cpp | 7 +++++-- c_tools/mock/loaders/multidark_loader.cpp | 7 ++++--- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index 5e93c3f..7f44b99 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -354,7 +354,8 @@ void createBox(SimuData *simu, vector& targets, vector& snapshot_spl for (int j = 0; j < 3; j++) { boxed->Pos[j] = new float[boxed->NumPart]; - boxed->Vel[j] = 0; + boxed->Vel[j] = new float[boxed->NumPart]; + //boxed->Vel[j] = 0; mul[j] = 1.0/(ranges[2*j+1] - ranges[2*j+0]); } cout << "Min range = " << ranges[0] << " " << ranges[2] << " " << ranges[4] << endl; @@ -413,6 +414,7 @@ void buildBox(SimuData *simu, long num_targets, long loaded, for (int j = 0; j < 3; j++) { boxed->Pos[j][loaded] = max(min((simu->Pos[j][pid]-ranges[j*2])*mul[j], double(1)), double(0)); + boxed->Vel[j][loaded] = simu->Vel[j][pid]; assert(boxed->Pos[j][loaded] >= 0); assert(boxed->Pos[j][loaded] <= 1); } @@ -551,7 +553,8 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a for (int j = 0; j < 3; j++) { boxed->Pos[j] = new float[boxed->NumPart]; - boxed->Vel[j] = 0; + boxed->Vel[j] = new float[boxed->NumPart]; + //boxed->Vel[j] = 0; mul[j] = 1.0/(ranges[2*j+1] - ranges[2*j+0]); } diff --git a/c_tools/mock/loaders/multidark_loader.cpp b/c_tools/mock/loaders/multidark_loader.cpp index 3d9c8a8..b9d0118 100644 --- a/c_tools/mock/loaders/multidark_loader.cpp +++ b/c_tools/mock/loaders/multidark_loader.cpp @@ -76,9 +76,10 @@ public: long estimated = (preproc == 0) ? simu->NumPart : preproc->getEstimatedPostprocessed(simu->NumPart); long allocated = estimated; - for (int k = 0; k < 3; k++) + for (int k = 0; k < 3; k++) { simu->Pos[k] = new float[allocated]; - simu->Vel[2] = new float[allocated]; + simu->Vel[k] = new float[allocated]; + } simu->Id = new long[allocated]; long *uniqueID = new long[allocated]; long *index = new long[allocated]; @@ -98,7 +99,7 @@ public: SingleParticle p; fp >> p.ID >> p.Pos[0] >> p.Pos[1] - >> p.Pos[2] >> p.Vel[2] >> tempData >> tempData >> tempData; + >> p.Pos[2] >> p.Vel[2] >> p.Vel[1] >> p.Vel[0] >> tempData; if (p.ID == -99 && p.Pos[0] == -99 && p.Pos[1] == -99 && From 303650fa4bb6f049bd99436af7aa1c270e684470 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 31 Jan 2014 10:44:42 -0600 Subject: [PATCH 08/12] updated velocity reader --- python_tools/void_python_tools/partUtil/partUtil.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python_tools/void_python_tools/partUtil/partUtil.py b/python_tools/void_python_tools/partUtil/partUtil.py index ba16b8a..aeec74c 100644 --- a/python_tools/void_python_tools/partUtil/partUtil.py +++ b/python_tools/void_python_tools/partUtil/partUtil.py @@ -120,7 +120,7 @@ def loadPartVel(workDir, sampleDir, sample): vx = File.variables['vel_x'][0:] vy = File.variables['vel_y'][0:] - vx = File.variables['vel_z'][0:] + vz = File.variables['vel_z'][0:] File.close() From eb7446e4f579b44e86a6de54c28e6547e1431633 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 31 Jan 2014 10:50:32 -0600 Subject: [PATCH 09/12] make particle loading more flexible --- .../void_python_tools/partUtil/partUtil.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/python_tools/void_python_tools/partUtil/partUtil.py b/python_tools/void_python_tools/partUtil/partUtil.py index aeec74c..46cca76 100644 --- a/python_tools/void_python_tools/partUtil/partUtil.py +++ b/python_tools/void_python_tools/partUtil/partUtil.py @@ -25,16 +25,20 @@ from netCDF4 import Dataset import sys from void_python_tools.backend import * import void_python_tools.apTools as vp +import pickle NetCDFFile = Dataset ncFloat = 'f8' # ----------------------------------------------------------------------------- -def loadPart(workDir, sampleDir, sample): +def loadPart(sampleDir): #print " Loading particle data..." sys.stdout.flush() - infoFile = workDir+"/"+sampleDir+"/zobov_slice_"+sample.fullName+".par" + with open(sampleDir+"/sample_info.dat", 'rb') as input: + sample = pickle.load(input) + + infoFile = sampleDir+"/zobov_slice_"+sample.fullName+".par" File = NetCDFFile(infoFile, 'r') ranges = np.zeros((3,2)) ranges[0][0] = getattr(File, 'range_x_min') @@ -49,7 +53,7 @@ def loadPart(workDir, sampleDir, sample): mul = np.zeros((3)) mul[:] = ranges[:,1] - ranges[:,0] - partFile = workDir+"/"+sampleDir+"/zobov_slice_"+sample.fullName + partFile = sampleDir+"/zobov_slice_"+sample.fullName iLine = 0 partData = [] part = np.zeros((3)) @@ -106,11 +110,14 @@ def loadPart(workDir, sampleDir, sample): return partData, boxLen, volNorm, isObservationData # ----------------------------------------------------------------------------- -def loadPartVel(workDir, sampleDir, sample): +def loadPartVel(sampleDir): #print " Loading particle velocities..." sys.stdout.flush() - infoFile = workDir+"/"+sampleDir+"/zobov_slice_"+sample.fullName+".par" + with open(sampleDir+"/sample_info.dat", 'rb') as input: + sample = pickle.load(input) + + infoFile = sampleDir+"/zobov_slice_"+sample.fullName+".par" File = NetCDFFile(infoFile, 'r') isObservation = getattr(File, 'is_observation') From b19fbca2ada5cf916deeccab8e5799fc1db14810 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Mon, 3 Feb 2014 04:40:50 -0600 Subject: [PATCH 10/12] improved handling of SDF files --- .../pipeline_source/prepareCatalogs.in.py | 203 +++++++----------- .../void_python_tools/backend/classes.py | 2 + .../void_python_tools/backend/launchers.py | 21 +- 3 files changed, 91 insertions(+), 135 deletions(-) diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index e80dfd6..46cd652 100644 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/python_tools/pipeline_source/prepareCatalogs.in.py @@ -69,6 +69,39 @@ if not os.access(filename, os.F_OK): parms = imp.load_source("name", filename) globals().update(vars(parms)) +#------------------------------------------------------------------------------ +def getSDFTags(dataFile): + idTag = "id" + for line in open(dataFile): + if "int64_t id" in line: + idTag = "id" + break + if "int64_t ident" in line: + idTag = "ident" + break + + massTag = "mass" + for line in open(dataFile): + if "float m200b" in line: + massTag = "m200b" + break + if "float mass" in line: + massTag = "mass" + break + + pidTag = "" + iLine = 0 + for line in open(dataFile): + iLine += 1 + if iLine > 100: break + if "parent_id" in line: + pidTag = "parent_id" + break + if "pid" in line: + pidTag = "pid" + break + + return idTag, massTag, pidTag #------------------------------------------------------------------------------ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): @@ -566,26 +599,7 @@ if (args.script or args.all) and haloFileBase != "": numPart = 0 if dataFormat == "sdf": SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" - massTag = "mass" - for line in open(dataFile): - if "float m200b" in line: - massTag = "m200b" - break - if "float mass" in line: - massTag = "mass" - break - - pidTag = "" - iLine = 0 - for line in open(dataFile): - iLine += 1 - if iLine > 100: break - if "parent_id" in line: - pidTag = "parent_id" - break - if "pid" in line: - pidTag = "pid" - break + idTag, massTag, pidTag = getSDFTags(dataFile) if minHaloMass == "none": command = "%s -a 200000 %s %s %s | awk '{if ($2==-1) print $1}' | wc" % (SDFcvt_PATH, dataFile, massTag, pidTag) @@ -686,35 +700,7 @@ if (args.halos or args.all) and haloFileBase != "" and len(minHaloMasses) > 0: outFile.close() if dataFormat == "sdf": - idTag = "id" - for line in open(dataFile): - if "int64_t id" in line: - idTag = "id" - break - if "int64_t ident" in line: - idTag = "ident" - break - - massTag = "mass" - for line in open(dataFile): - if "float m200b" in line: - massTag = "m200b" - break - if "float mass" in line: - massTag = "mass" - break - - pidTag = "" - iLine = 0 - for line in open(dataFile): - iLine += 1 - if iLine > 100: break - if "parent_id" in line: - pidTag = "parent_id" - break - if "pid" in line: - pidTag = "pid" - break + idTag, massTag, pidTag = getSDFTags(dataFile) SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" if minHaloMass == "none": @@ -832,6 +818,7 @@ if (args.halos or args.all) and haloFileBase != "": print " ERROR: not enough halos to support that density! Maximum is %g" % (1.*numPart / lbox**3) exit(-1) + sampleName = prefix+"halos_den"+str(haloDen)+"_z"+redshifts[iRedshift] outFileName = catalogDir+"/"+sampleName+".dat" outFile = open(outFileName, 'w') @@ -843,48 +830,42 @@ if (args.halos or args.all) and haloFileBase != "": outFile.close() if dataFormat == "sdf": - idTag = "id" - for line in open(dataFile): - if "int64_t id" in line: - idTag = "id" - break - if "int64_t ident" in line: - idTag = "ident" - break - - massTag = "mass" - for line in open(dataFile): - if "float m200b" in line: - massTag = "m200b" - break - if "float mass" in line: - massTag = "mass" - break - - pidTag = "" - iLine = 0 - for line in open(dataFile): - iLine += 1 - if iLine > 100: break - if "parent_id" in line: - pidTag = "parent_id" - break - if "pid" in line: - pidTag = "pid" - break - print "TEST", massTag, pidTag + idTag, massTag, pidTag = getSDFTags(dataFile) + + tempFile = catalogDir+"/temp_"+sampleName+".dat" SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" - command = "%s -a 200000 %s %s %s x y z vz vy vx %s | sort -rg -k 1 | awk '{if ($9==-1 && NR < %d) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, massTag, idTag, pidTag, numPartExpect, outFileName ) + command = "%s -a 200000 %s %s %s x y z vz vy vx %s | awk '{if ($9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, massTag, idTag, pidTag, tempFile ) #os.system(command) subprocess.call(command, shell=True) - outFile = open(outFileName, 'a') + + # do the check again since now we've filtered out subhalos + numPart = 0 + for line in enumerate(tempFile): + numPart += 1 + + actualDen = 1.*numPart / lbox**3 + keepFraction = haloDen / actualDen + if numPart < numPartExpect: + print " ERROR: not enough galaxies to support that density! Maximum is %g" % (1.*numPart / lbox**3) + exit(-1) + + numKept = 0 + for (iLine,line) in enumerate(inFile): + if np.random.uniform() > keepFraction: continue + outFile.write(line) + numKept += 1 + outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n") outFile.close() else: + actualDen = 1.*numPart / lbox**3 + keepFraction = haloDen / actualDen + inFile = open(dataFile, 'r') - haloList = [] + outFile = open(outFileName, 'a') for (iHalo,line) in enumerate(inFile): if iHalo < haloFileNumComLines: continue + if np.random.uniform() > keepFraction: continue line = line.split(haloFileColSep) x = float(line[haloFileXCol]) * haloFilePosRescale y = float(line[haloFileYCol]) * haloFilePosRescale @@ -893,31 +874,13 @@ if (args.halos or args.all) and haloFileBase != "": vy = float(line[haloFileVYCol]) vx = float(line[haloFileVXCol]) mass = float(line[haloFileMCol]) - haloList.append((x,y,z,vz,vy,vx,mass)) - inFile.close() - - haloList = np.array(haloList) - haloMasses = haloList[:,6] - haloList = haloList[np.argsort(haloMasses)[::-1]] - - outFile = open(outFileName, 'a') - for iHalo in xrange(numPartExpect): - x = haloList[iHalo,0] - y = haloList[iHalo,1] - z = haloList[iHalo,2] - vz = haloList[iHalo,3] - vy = haloList[iHalo,4] - vx = haloList[iHalo,5] - mass = haloList[iHalo,6] - - # write to output file - outFile.write("%d %e %e %e %e %e %e %e\n" %(iHalo,x,y,z, - vz,vy,vx,mass)) + outFile.write("%d %e %e %e %e %e %e %e\n" %(iHalo,x,y,z, + vz,vy,vx,mass)) outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n") + inFile.close() outFile.close() - # ----------------------------------------------------------------------------- # now the HOD parFileText = """ @@ -1010,26 +973,7 @@ if (args.hod or args.all) and haloFileBase != "": inFile = haloFile outFile = haloFile+"_temp" - massTag = "mass" - for line in open(inFile): - if "float m200b" in line: - massTag = "m200b" - break - if "float mass" in line: - massTag = "mass" - break - - pidTag = "" - iLine = 0 - for line in open(inFile): - iLine += 1 - if iLine > 100: break - if "parent_id" in line: - pidTag = "parent_id" - break - if "pid" in line: - pidTag = "pid" - break + idTag, massTag, pidTag = getSDFTags(inFile) SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" command = "%s -a 200000 %s %s x y z vx vy vz %s | awk '{if ($8 ==-1) print $1, $2, $3, $4, $5, $6, $7}'>>%s" % (SDFcvt_PATH, inFile, massTag, pidTag, outFile) @@ -1106,14 +1050,17 @@ if (args.hod or args.all) and haloFileBase != "": inFile = open(catalogDir+"/hod_"+sampleName+".mock") outFile = open(catalogDir+"/"+sampleName+".dat", 'w') - outFile.write("%f\n" %(lbox)) - outFile.write("%s\n" %(omegaM)) - outFile.write("%s\n" %(hubble)) - outFile.write("%s\n" %(redshift)) - outFile.write("%d\n" %(numPartExpect)) + #outFile.write("%f\n" %(lbox)) + #outFile.write("%s\n" %(omegaM)) + #outFile.write("%s\n" %(hubble)) + #outFile.write("%s\n" %(redshift)) + #outFile.write("%d\n" %(numPartExpect)) numKept = 0 for (iLine,line) in enumerate(inFile): + if iLine < 5: + outFile.write(line) + continue if np.random.uniform() > keepFraction: continue outFile.write(line) numKept += 1 diff --git a/python_tools/void_python_tools/backend/classes.py b/python_tools/void_python_tools/backend/classes.py index 7f9c84f..cd0af1b 100644 --- a/python_tools/void_python_tools/backend/classes.py +++ b/python_tools/void_python_tools/backend/classes.py @@ -21,6 +21,8 @@ import os from numpy import abs +import matplotlib +matplotlib.use('Agg') LIGHT_SPEED = 299792.458 diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 07836c5..530fb01 100644 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -325,6 +325,9 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, else: maskIndex = -1 maxDen = 0.2 + # TEST + #maxDen = 1.e99 + # END TEST if not (continueRun and jobSuccessful(logFile, "Done!\n")): for fileName in glob.glob(zobovDir+"/part._"+sampleName+".*"): @@ -413,6 +416,9 @@ def launchPrune(sample, binPath, mockIndex = -1 maxDen = 0.2 observationLine = "" + # TEST + #maxDen = 1.e99 + # END TEST #periodicLine = " --periodic='" #if sample.numSubvolumes == 1: periodicLine += "xy" @@ -577,7 +583,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, centralRadius = stack.rMin * 0.25 - # restrict to relavent ranges of sample + # restrict to relevant ranges of sample zMin = max(sample.zRange[0],stack.zMin) * 3000 zMax = min(sample.zRange[1],stack.zMax) * 3000 @@ -596,6 +602,9 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, else: maskIndex = 999999999 maxDen = 0.2 + # TEST + #maxDen = 1.e99 + # END TEST if INCOHERENT: nullTestFlag = "INCOHERENT" @@ -1219,8 +1228,8 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None, while badChain: ntries += 1 - #Rexpect = (stack.rMin+stack.rMax)/2 - #Rtruncate = stack.rMin*3. + 1 # TEST + Rexpect = (stack.rMin+stack.rMax)/2 + Rtruncate = stack.rMin*3. + 1 # TEST #if sample.dataType == "observation": # ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, # Niter=300000, @@ -1231,11 +1240,9 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None, # Niter=300000, # Nburn=100000, # Rextracut=Rtruncate) - #badChain = (args[0][0] > 0.5 or args[0][1] > stack.rMax or \ - # args[0][2] > stack.rMax) and \ + #badChain = (args[0][0] > 0.5 or args[0][1] > 2.*stack.rMax or \ + # args[0][2] > 2.*stack.rMax) and \ # (ntries < maxtries) - #ret,fits,args = vp.compute_radial_inertia(voidDir, stack.rMax, mode="symmetric", nBootstraps=5) - #ret,fits,args = vp.compute_inertia(voidDir, stack.rMax, nBootstraps=100, rMaxInertia=1.0) ret,fits,args = vp.compute_inertia(voidDir, stack.rMax, mode="2d", nBootstraps=500, rMaxInertia=0.7) #ret,fits,args = vp.compute_inertia(voidDir, stack.rMax, mode="symmetric", nBootstraps=500, rMaxInertia=100) badChain = False From 8dbf8f8b54713098aa383b79e6a36b14782e78c1 Mon Sep 17 00:00:00 2001 From: Paul Matthew Sutter Date: Mon, 3 Feb 2014 15:48:15 +0100 Subject: [PATCH 11/12] fixed bug when loading halos --- python_tools/pipeline_source/prepareCatalogs.in.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index 46cd652..ba7cd42 100644 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/python_tools/pipeline_source/prepareCatalogs.in.py @@ -706,7 +706,7 @@ if (args.halos or args.all) and haloFileBase != "" and len(minHaloMasses) > 0: if minHaloMass == "none": command = "%s -a 200000 %s %s %s x y z vz vy vx %s | awk '{if ($9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, massTag, idTag, pidTag, outFileName ) else: - command = "%s -a 200000 %s %s %s x y z vz vy vx %s | awk '{if ($1>%g && $9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, massTag, idTag, minHaloMass, pidTag, outFileName ) + command = "%s -a 200000 %s %s %s x y z vz vy vx %s | awk '{if ($1>%g && $9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, massTag, idTag, pidTag, minHaloMass, outFileName ) #os.system(command) subprocess.call(command, shell=True) outFile = open(outFileName, 'a') From 264d55f52b62d11216d41356980cdb16dfeea54a Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Mon, 3 Feb 2014 10:08:00 -0600 Subject: [PATCH 12/12] utility to pull void member particles --- .../void_python_tools/partUtil/partUtil.py | 143 ++++++++++++++++++ 1 file changed, 143 insertions(+) diff --git a/python_tools/void_python_tools/partUtil/partUtil.py b/python_tools/void_python_tools/partUtil/partUtil.py index 46cca76..0f8a1f7 100644 --- a/python_tools/void_python_tools/partUtil/partUtil.py +++ b/python_tools/void_python_tools/partUtil/partUtil.py @@ -156,3 +156,146 @@ def shiftPart(inPart, newCenter, periodicLine, boxLen): np.copysign(boxLen[2], part[shiftUs,2]) return part + +# ----------------------------------------------------------------------------- + +class Bunch: + def __init__(self, **kwds): + self.__dict__.update(kwds) + +class Catalog: + numVoids = 0 + numPartTot = 0 + numZonesTot = 0 + boxLen = np.zeros((3)) + part = None + zones2Parts = None + void2Zones = None + voids = None + +# ----------------------------------------------------------------------------- +def loadVoidCatalog(sampleDir): + #print " Loading particle data..." + sys.stdout.flush() + + catalog = Catalog() + with open(sampleDir+"/sample_info.dat", 'rb') as input: + sample = pickle.load(input) + + print "Loading info..." + infoFile = sampleDir+"/zobov_slice_"+sample.fullName+".par" + File = NetCDFFile(infoFile, 'r') + ranges = np.zeros((3,2)) + ranges[0][0] = getattr(File, 'range_x_min') + ranges[0][1] = getattr(File, 'range_x_max') + ranges[1][0] = getattr(File, 'range_y_min') + ranges[1][1] = getattr(File, 'range_y_max') + ranges[2][0] = getattr(File, 'range_z_min') + ranges[2][1] = getattr(File, 'range_z_max') + catalog.boxLen[0] = ranges[0][1] - ranges[0][0] + catalog.boxLen[1] = ranges[1][1] - ranges[1][0] + catalog.boxLen[2] = ranges[2][1] - ranges[2][0] + File.close() + + print "Loading all particles..." + partData, boxLen, volNorm, isObservationData = loadPart(sampleDir) + numPartTot = len(partData) + catalog.numPartTot = numPartTot + catalog.part = [] + for i in xrange(len(partData)): + catalog.part.append(Bunch(x = partData[i][0], + y = partData[i][1], + z = partData[i][2], + volume = 0, + ra = 0, + dec = 0, + redshift = 0, + uniqueID = 0)) + + + print "Loading volumes..." + volFile = sampleDir+"/vol_"+sample.fullName+".dat" + File = file(volFile) + chk = np.fromfile(File, dtype=np.int32,count=1) + vols = np.fromfile(File, dtype=np.float32,count=numPartTot) + for ivol in xrange(len(vols)): + catalog.part[ivol].volume = vols[ivol] / volNorm + + print "Loading voids..." + fileName = sampleDir+"/voidDesc_central_"+sample.fullName+".out" + catData = np.loadtxt(fileName, comments="#", skiprows=2) + catalog.voids = [] + for line in catData: + catalog.voids.append(Bunch(iVoid = int(line[0]), + voidID = int(line[1]), + coreParticle = line[2], + coreDens = line[3], + zoneVol = line[4], + zoneNumPart = line[5], + numZones = int(line[6]), + voidVol = line[7], + numPart = int(line[8]), + densCon = line[9], + voidProp = line[10], + radius = pow(line[7]/volNorm*3./4./np.pi, 1./3.), + barycenter = np.zeros((3)))) + + print "Read %d voids" % len(catalog.voids) + + print "Loading barycenters..." + iLine = 0 + for line in open(sampleDir+"/barycenters_central_"+sample.fullName+".out"): + line = line.split() + catalog.voids[iLine].barycenter[0] = float(line[1]) + catalog.voids[iLine].barycenter[1] = float(line[2]) + catalog.voids[iLine].barycenter[2] = float(line[3]) + iLine += 1 + + + print "Loading zone-void membership info..." + zoneFile = sampleDir+"/voidZone_"+sample.fullName+".dat" + catalog.void2Zones = [] + File = file(zoneFile) + numZonesTot = np.fromfile(File, dtype=np.int32,count=1) + catalog.numZonesTot = numZonesTot + for iZ in xrange(numZonesTot): + numZones = np.fromfile(File, dtype=np.int32,count=1) + catalog.void2Zones.append(Bunch(numZones = numZones, + zoneIDs = [])) + + for p in xrange(numZones): + zoneID = np.fromfile(File, dtype=np.int32,count=1) + catalog.void2Zones[iZ].zoneIDs.append(zoneID) + + + print "Loading particle-zone membership info..." + zonePartFile = sampleDir+"/voidPart_"+sample.fullName+".dat" + catalog.zones2Parts = [] + File = file(zonePartFile) + chk = np.fromfile(File, dtype=np.int32,count=1) + numZonesTot = np.fromfile(File, dtype=np.int32,count=1) + for iZ in xrange(numZonesTot): + numPart = np.fromfile(File, dtype=np.int32,count=1) + catalog.zones2Parts.append(Bunch(numPart = numPart, + partIDs = [])) + + for p in xrange(numPart): + partID = np.fromfile(File, dtype=np.int32,count=1) + catalog.zones2Parts[iZ].partIDs.append(partID) + + return catalog + + + +# ----------------------------------------------------------------------------- +def getVoidPart(catalog, voidID): + + partOut = [] + for iZ in xrange(catalog.void2Zones[voidID].numZones): + zoneID = catalog.void2Zones[voidID].zoneIDs[iZ] + for p in xrange(catalog.zones2Parts[zoneID].numPart): + partID = catalog.zones2Parts[zoneID].partIDs[p] + partOut.append(catalog.part[partID]) + + return partOut +