From 326756b2bc2c2d542180cabbb8bd6da17c5bd47c Mon Sep 17 00:00:00 2001 From: "Paul M. Sutter" Date: Tue, 4 Mar 2025 14:07:36 -0500 Subject: [PATCH] + radial profiles now full actual density from survey volume, not zobov normalization + getVolNorm provides both zobov normalization and average density from survey volume for observations + significant update and cleanup to plotting routines --- c_source/pruning/pruneVoids.cpp | 3 -- python_source/backend/surveyTools.py | 10 +++-- python_source/voidUtil/catalogUtil.py | 45 +++++++++++++---------- python_source/voidUtil/periodic_kdtree.py | 8 ++-- python_source/voidUtil/plotUtil.py | 30 +++++++-------- python_source/voidUtil/profileUtil.py | 2 +- python_source/voidUtil/xcorUtil.py | 34 +++++++++-------- python_source/voidUtil/xcorlib.py | 2 +- 8 files changed, 73 insertions(+), 61 deletions(-) diff --git a/c_source/pruning/pruneVoids.cpp b/c_source/pruning/pruneVoids.cpp index 82c6273..32ec124 100644 --- a/c_source/pruning/pruneVoids.cpp +++ b/c_source/pruning/pruneVoids.cpp @@ -667,7 +667,6 @@ int main(int argc, char **argv) { interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; //printf(" %.2f for maximum extent\n", interval); -/* // compute distance from center to nearest mock boundary particle // (with new boundary handling this will go away) clock3 = clock(); @@ -686,7 +685,6 @@ int main(int argc, char **argv) { } else { voids[iVoid].nearestMock = 1.e99; } - if (args.isObservation_flag) { voids[iVoid].redshiftInMpc = sqrt(pow(voids[iVoid].macrocenter[0] - boxLen[0]/2.,2) + @@ -745,7 +743,6 @@ int main(int argc, char **argv) { clock4 = clock(); interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; //printf(" %.2f for nearest edge\n", interval); -*/ // compute eigenvalues and vectors for orientation and shape clock3 = clock(); diff --git a/python_source/backend/surveyTools.py b/python_source/backend/surveyTools.py index e331105..44ac8f0 100644 --- a/python_source/backend/surveyTools.py +++ b/python_source/backend/surveyTools.py @@ -194,11 +194,9 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile, rayDist = meanPartSep raySteps = np.arange(zmin, zmax, rayDist) log.write(" Will take %d steps per ray with %.2e distance between steps\n" % (len(raySteps), rayDist)) - #print(meanPartSep, len(raySteps)) contourPixels = np.nonzero(contourMap)[0] log.write(" We have %d rays to work with\n" % (len(contourPixels))) - #print(contourPixels) for pixel in contourPixels: #print("Working with pixel %d" % pixel) vec = healpy.pix2vec(nside,pixel) @@ -227,6 +225,11 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile, dist, nearest = galTree.query(minEdge) flagList[nearest] = 3 + log.write(" Flagging statistics:\n") + log.write(" %d original galaxies\n" % len(flagList)) + log.write(" %d near edges\n" % len(flagList[flagList==1])) + log.write(" %d near redshift bounds\n" % len(flagList[flagList==2])) + log.write(" %d remaining\n" % len(flagList[flagList==0])) # output flag information log.write(" Saving galaxy flags to file...\n") np.savetxt(edgeGalFile, flagList, fmt="%d") @@ -243,7 +246,8 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile, ipix = healpy.ang2pix(nside, theta, phi) np.put(flagMap, ipix, 1) - healpy.write_map(outputDir+"/flagged_galaxies.fits", flagMap, overwrite=True, + healpy.write_map(outputDir+"/flagged_galaxies.fits", flagMap, + overwrite=True, dtype=np.dtype('float64')) diff --git a/python_source/voidUtil/catalogUtil.py b/python_source/voidUtil/catalogUtil.py index 0aca19a..ead60a2 100644 --- a/python_source/voidUtil/catalogUtil.py +++ b/python_source/voidUtil/catalogUtil.py @@ -115,6 +115,8 @@ def loadPart(sampleDir): boxLen = mul + # TODO - should use built-in volNorm function instead + #if isObservation == 1: # # look for the mask file # if os.access(sample.maskFile, os.F_OK): @@ -165,26 +167,29 @@ def getVolNorm(sampleDir): boxLen = mul - #if isObservation == 1: - # # look for the mask file - # if os.access(sample.maskFile, os.F_OK): - # maskFile = sample.maskFile - # else: - # maskFile = sampleDir+"/"+os.path.basename(sample.maskFile) - # print "Using maskfile found in:", maskFile - # props = getSurveyProps(maskFile, sample.zBoundary[0], - # sample.zBoundary[1], - # sample.zBoundary[0], - # sample.zBoundary[1], "all", - # selectionFuncFile=sample.selFunFile, - # useComoving=sample.useComoving) - # boxVol = props[0] - # volNorm = maskIndex/boxVol - #else: boxVol = np.prod(boxLen) volNorm = Np/boxVol + volNormObs = volNorm - return volNorm + if isObservation == 1: + # look for the mask file + if os.access(sample.maskFile, os.F_OK): + maskFile = sample.maskFile + else: + maskFile = sampleDir+"/"+os.path.basename(sample.maskFile) + print("Using maskfile found in: " + maskFile) + props = getSurveyProps(maskFile, sample.zBoundary[0], + sample.zBoundary[1], + sample.zBoundary[0], + sample.zBoundary[1], + "all", + sample.omegaM, + selectionFuncFile=sample.selFunFile, + useComoving=sample.useComoving) + boxVol = props[0] + volNormObs = maskIndex/boxVol + + return volNorm, volNormObs # ----------------------------------------------------------------------------- @@ -285,7 +290,8 @@ class Catalog: numVoids = 0 numPartTot = 0 numZonesTot = 0 - volNorm = 0 + volNorm = 0 # normalization used by zobov across entire volumne + volNormObs = 0 # normalization based on average density of survey volume boxLen = np.zeros((3)) ranges = np.zeros((3,2)) part = None @@ -328,8 +334,9 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True, catalog.ranges = ranges File.close() - volNorm = getVolNorm(sampleDir) + volNorm, volNormObs = getVolNorm(sampleDir) catalog.volNorm = volNorm + catalog.volNormObs = volNormObs if untrimmed: prefix = "untrimmed_" diff --git a/python_source/voidUtil/periodic_kdtree.py b/python_source/voidUtil/periodic_kdtree.py index 3a6aba6..be7d7f2 100644 --- a/python_source/voidUtil/periodic_kdtree.py +++ b/python_source/voidUtil/periodic_kdtree.py @@ -279,12 +279,12 @@ class PeriodicCKDTree(cKDTree): retshape = np.shape(x)[:-1] if retshape!=(): if k>1: - dd = np.empty(retshape+(k,),dtype=np.float) + dd = np.empty(retshape+(k,),dtype=float) dd.fill(np.inf) ii = np.empty(retshape+(k,),dtype=np.int) ii.fill(self.n) elif k==1: - dd = np.empty(retshape,dtype=np.float) + dd = np.empty(retshape,dtype=float) dd.fill(np.inf) ii = np.empty(retshape,dtype=np.int) ii.fill(self.n) @@ -310,7 +310,7 @@ class PeriodicCKDTree(cKDTree): else: return np.inf, self.n elif k>1: - dd = np.empty(k,dtype=np.float) + dd = np.empty(k,dtype=float) dd.fill(np.inf) ii = np.empty(k,dtype=np.int) ii.fill(self.n) @@ -368,7 +368,7 @@ class PeriodicCKDTree(cKDTree): save substantial amounts of time by putting them in a PeriodicCKDTree and using query_ball_tree. """ - x = np.asarray(x).astype(np.float) + x = np.asarray(x).astype(float) if x.shape[-1] != self.m: raise ValueError("Searching for a %d-dimensional point in a " \ "%d-dimensional KDTree" % (x.shape[-1], self.m)) diff --git a/python_source/voidUtil/plotUtil.py b/python_source/voidUtil/plotUtil.py index 18674a1..f01a6c9 100644 --- a/python_source/voidUtil/plotUtil.py +++ b/python_source/voidUtil/plotUtil.py @@ -57,8 +57,8 @@ def plotNumberFunction(catalogList, # cumulative: if True, plots cumulative number function # binWidth: width of histogram bins in Mpc/h # returns: -# ellipDistList: array of len(catalogList), -# each element has array of size bins, number, +/- 1 sigma +# distList: array of len(catalogList), +# each element has array of size bins, number, +/- 1 sigma print("Plotting number function") @@ -72,7 +72,7 @@ def plotNumberFunction(catalogList, else: plt.ylabel(r"log ($dn/dR$ [$h^3$ Gpc$^{-3}$])", fontsize=14) - ellipDistList = [] + distList = [] for (iSample,catalog) in enumerate(catalogList): sample = catalog.sampleInfo @@ -126,24 +126,25 @@ def plotNumberFunction(catalogList, alpha = 0.55 fill_between(binCentersToUse, lower, upper, - label=lineTitle, color=lineColor, + label='', color=lineColor, alpha=alpha, ) lineStyle = '-' plt.plot(binCentersToUse, mean, lineStyle, color=lineColor, + label=lineTitle, linewidth=3) - ellipDistList.append((binCentersToUse, mean, lower, upper)) + distList.append((binCentersToUse, mean, lower, upper)) plt.legend(loc = "upper right", fancybox=True, prop={'size':14}) plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight") - plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") + #plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight") - return ellipDistList + return distList # ----------------------------------------------------------------------------- @@ -186,10 +187,10 @@ def plotEllipDist(catalogList, plt.legend(loc = "upper right", fancybox=True, prop={'size':14}) plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight") - plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") + #plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight") - return + return ellipDistList # ----------------------------------------------------------------------------- @@ -297,7 +298,7 @@ def plotVoidCells(catalog, plotName="cells"+str(voidID) plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight") - plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") + #plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight") return @@ -305,7 +306,7 @@ def plotVoidCells(catalog, # ----------------------------------------------------------------------------- def plotVoronoi(catalog, - voidIndex=-1, + voidID=-1, figDir="./", plotName="voronoi", plotCenter = (0,0,0), @@ -314,7 +315,7 @@ def plotVoronoi(catalog, # plots the voronoi cells belonging to a void (or a slice if no void given) # catalog: void catalog -# voidIndec: index of void to plot (-1 to use all particles) +# voidID: ID of void to plot (-1 to use all particles) # figDir: output directory for figures # plotName: name to prefix to all outputs @@ -327,7 +328,7 @@ def plotVoronoi(catalog, print(" Selecting particles...") - if (voidIndex == -1): + if (voidID == -1): xmin = plotCenter[0]-plotWidth[0]/2. xmax = plotCenter[0]+plotWidth[0]/2. ymin = plotCenter[1]-plotWidth[1]/2. @@ -348,7 +349,6 @@ def plotVoronoi(catalog, #ax.set_ylim(ymin,ymax) #ax_set.zlim(zmin,zmax) else: - voidID = catalog.voids[voidIndex].voidID partList = getVoidPart(catalog, voidID) plotCenter = catalog.voids[voidIndex].macrocenter plotRadius = catalog.voids[voidIndex].radius @@ -381,7 +381,7 @@ def plotVoronoi(catalog, print(" Saving...") plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight") - plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") + #plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight") print(" Done!") diff --git a/python_source/voidUtil/profileUtil.py b/python_source/voidUtil/profileUtil.py index b5ba21d..8802ae3 100644 --- a/python_source/voidUtil/profileUtil.py +++ b/python_source/voidUtil/profileUtil.py @@ -75,7 +75,7 @@ def buildProfile(catalog, rMin, rMax, nBins=10): deltaV = 4*np.pi/3*(radii[1:]**3-radii[0:(radii.size-1)]**3) thisProfile = np.float32(thisProfile) thisProfile /= deltaV - thisProfile /= catalog.volNorm + thisProfile /= catalog.volNormObs allProfiles.append(thisProfile) binCenters = 0.5*(radii[1:] + radii[:-1]) diff --git a/python_source/voidUtil/xcorUtil.py b/python_source/voidUtil/xcorUtil.py index 6b92850..4f178c3 100644 --- a/python_source/voidUtil/xcorUtil.py +++ b/python_source/voidUtil/xcorUtil.py @@ -27,13 +27,15 @@ from voidUtil import getArray def computeXcor(catalog, figDir="./", + namePrefix="", Nmesh = 256, - Nbin = 100 + Nbin = 100, ): # Computes and plots void-void and void-matter(galaxy) correlations # catalog: catalog to analyze # figDir: where to place plots +# namePrefix: prefix to add to file names # Nmesh: number of grid cells in cic mesh-interpolation # Nbin: number of bins in final plots @@ -84,18 +86,18 @@ def computeXcor(catalog, plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$') plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$') plt.title(r'Dark matter') - plt.savefig(figDir+'/dm.eps', bbox_inches="tight") - plt.savefig(figDir+'/dm.pdf', bbox_inches="tight") - plt.savefig(figDir+'/dm.png', bbox_inches="tight") + #plt.savefig(figDir+'/dm.eps', bbox_inches="tight") + plt.savefig(figDir+'/'+namePrefix+'dm.pdf', bbox_inches="tight") + plt.savefig(figDir+'/'+namePrefix+'dm.png', bbox_inches="tight") plt.clf() plt.imshow(np.sum(dv[:,:,:]+1,2)/Nmesh,extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian') plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$') plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$') plt.title(r'Voids') - plt.savefig(figDir+'/dv.eps', bbox_inches="tight") #, dpi=300 - plt.savefig(figDir+'/dv.pdf', bbox_inches="tight") #, dpi=300 - plt.savefig(figDir+'/dv.png', bbox_inches="tight") #, dpi=300 + #plt.savefig(figDir+'/dv.eps', bbox_inches="tight") #, dpi=300 + plt.savefig(figDir+'/'+namePrefix+'dv.pdf', bbox_inches="tight") #, dpi=300 + plt.savefig(figDir+'/'+namePrefix+'dv.png', bbox_inches="tight") #, dpi=300 plt.clf() @@ -116,10 +118,11 @@ def computeXcor(catalog, plt.yscale('log') plt.xlim(kmin,kmax) plt.ylim(10**np.floor(np.log10(abs(Pvm[1:]).min()))/margin, max(10**np.ceil(np.log10(Pmm.max())),10**np.ceil(np.log10(Pvv.max())))*margin) - plt.legend([pa, pb, pc],['tt', 'vt', 'vv'],'best',prop={'size':12}) - plt.savefig(figDir+'/power.eps', bbox_inches="tight") - plt.savefig(figDir+'/power.pdf', bbox_inches="tight") - plt.savefig(figDir+'/power.png', bbox_inches="tight") + plt.legend([pa, pb, pc],['tt', 'vt', 'vv']) + #plt.legend([pa, pb, pc],['tt', 'vt', 'vv'],'best',prop={'size':12}) + #plt.savefig(figDir+'/power.eps', bbox_inches="tight") + plt.savefig(figDir+'/'+namePrefix+'power.pdf', bbox_inches="tight") + plt.savefig(figDir+'/'+namePrefix+'power.png', bbox_inches="tight") plt.clf() pa ,= plt.plot(rm, Xmm, 'k-o', ms=0.8*ms, mew=mew) @@ -137,10 +140,11 @@ def computeXcor(catalog, plt.yscale('log') plt.xlim(rmin,rmax) plt.ylim(min(10**np.floor(np.log10(abs(Xvm).min())),10**np.floor(np.log10(abs(Xmm).min())))/margin, max(10**np.ceil(np.log10(Xmm.max())),10**np.ceil(np.log10(Xvv.max())))*margin) - plt.legend([pa, pb, pc],['tt', 'vt', 'vv'],'best',prop={'size':12}) - plt.savefig(figDir+'/correlation.eps', bbox_inches="tight") - plt.savefig(figDir+'/correlation.pdf', bbox_inches="tight") - plt.savefig(figDir+'/correlation.png', bbox_inches="tight") + plt.legend([pa, pb, pc],['tt', 'vt', 'vv']) + #plt.legend([pa, pb, pc],['tt', 'vt', 'vv'],'best',prop={'size':12}) + #plt.savefig(figDir+'/correlation.eps', bbox_inches="tight") + plt.savefig(figDir+'/'+namePrefix+'correlation.pdf', bbox_inches="tight") + plt.savefig(figDir+'/'+namePrefix+'correlation.png', bbox_inches="tight") plt.clf() diff --git a/python_source/voidUtil/xcorlib.py b/python_source/voidUtil/xcorlib.py index 2ef6ff0..f300fad 100644 --- a/python_source/voidUtil/xcorlib.py +++ b/python_source/voidUtil/xcorlib.py @@ -71,7 +71,7 @@ def powcor(d1, d2, Lbox, Nbin = 10, scale = 'lin', cor = False, cic = True, dim dx = Lbox/Nmesh Xr = np.fft.irfftn(Pk)*(Nmesh/Lbox)**3 - (Nmx, rm, Xrm, SXrm) = shellavg(np.real(Xr), dx, Nmesh, Nbin = Nbin/2, xmin = dx, xmax = 140., scale = scale, dim = dim) + (Nmx, rm, Xrm, SXrm) = shellavg(np.real(Xr), dx, Nmesh, Nbin = Nbin//2, xmin = dx, xmax = 140., scale = scale, dim = dim) return ((Nm, km, Pkm, SPkm),(Nmx, rm, Xrm, SXrm))