From dadfd8d40ae34226b45477cf0e0a822e5e3eff20 Mon Sep 17 00:00:00 2001 From: "Paul M. Sutter" Date: Mon, 21 Apr 2025 15:02:05 -0400 Subject: [PATCH] extra info in void catalog output (void type, max radius, nearest flagged galaxy), added filtering routines on these properties for post-analysis --- c_source/pruning/pruneVoids.cpp | 78 +++++++++++++-------------- python_source/voidUtil/catalogUtil.py | 37 ++++++++++--- 2 files changed, 70 insertions(+), 45 deletions(-) diff --git a/c_source/pruning/pruneVoids.cpp b/c_source/pruning/pruneVoids.cpp index d84cef6..a3e4be2 100644 --- a/c_source/pruning/pruneVoids.cpp +++ b/c_source/pruning/pruneVoids.cpp @@ -76,8 +76,7 @@ typedef struct voidStruct { float vol, coreDens, zoneVol, densCon, voidProb, radius; float rescaledCoreDens; int voidID, numPart, numZones, coreParticle, zoneNumPart; - float maxRadius, nearestMock, centralDen, redshift, redshiftInMpc; - float nearestMockFromCore, nearestGalFromCore; + float maxRadius, nearestFlag, centralDen, redshift, redshiftInMpc; float nearestEdge; float center[3], macrocenter[3]; int accepted; @@ -191,7 +190,7 @@ int main(int argc, char **argv) { int numVoids, mockIndex; double tolerance; FILE *fp; - PART *part, *voidPart; + PART *part, *voidPart, *flaggedPart; ZONE2PART *zones2Parts; VOID2ZONE *void2Zones; vector voids; @@ -407,12 +406,24 @@ int main(int argc, char **argv) { // and finally edge flag info printf(" Loading particle edge flags...\n"); + int numFlagged = 0; fp = fopen(args.partEdge_arg, "r"); for (p = 0; p < mask_index; p++) { fscanf(fp, "%d", &part[p].edgeFlag); + if (part[p].edgeFlag > 0) numFlagged++; } fclose(fp); + // copy flagged galaxies to separate array for faster processing later + flaggedPart = (PART *) malloc(numFlagged * sizeof(PART)); + int iFlag = 0; + for (p = 0; p < mask_index; p++) { + if (part[p].edgeFlag > 0) { + flaggedPart[iFlag] = part[p]; + iFlag++; + } + } + /* this was used for testing at one point // and finally finally adjacencies printf(" Loading particle adjacencies...\n"); @@ -521,7 +532,7 @@ int main(int argc, char **argv) { voids[iVoid].voidType = CENTRAL_VOID; - // first load up particles into a buffer + // first load up this void's particles into a buffer clock3 = clock(); i = 0; for (iZ = 0; iZ < void2Zones[voidID].numZones; iZ++) { @@ -665,24 +676,24 @@ int main(int argc, char **argv) { clock4 = clock(); 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(); + + // compute distance from macrocenter to nearest flagged particle + clock3 = clock(); if (args.isObservation_flag) { minDist = 1.e99; - for (p = mockIndex; p < numPartTot; p++) { - dist[0] = voids[iVoid].macrocenter[0] - part[p].x; - dist[1] = voids[iVoid].macrocenter[1] - part[p].y; - dist[2] = voids[iVoid].macrocenter[2] - part[p].z; + + for (p = 0; p < numFlagged; p++) { + dist[0] = voids[iVoid].macrocenter[0] - flaggedPart[p].x; + dist[1] = voids[iVoid].macrocenter[1] - flaggedPart[p].y; + dist[2] = voids[iVoid].macrocenter[2] - flaggedPart[p].z; dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); if (dist2 < minDist) minDist = dist2; } - voids[iVoid].nearestMock = sqrt(minDist); + voids[iVoid].nearestFlag = sqrt(minDist); } else { - voids[iVoid].nearestMock = 1.e99; + voids[iVoid].nearestFlag = 1.e99; } if (args.isObservation_flag) { voids[iVoid].redshiftInMpc = @@ -919,18 +930,6 @@ int main(int argc, char **argv) { numEdge++; } } -/* - // mark voids near survey edges - for (iVoid = 0; iVoid < voids.size(); iVoid++) { - if (tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) { - voids[iVoid].voidType = CENTRAL_VOID; - numCentral++; - } else { - voids[iVoid].voidType = EDGE_VOID; - numEdge++; - } - } -*/ printf(" Number kept: %d (out of %d)\n", (int) voids.size(), numVoids); printf(" We have %d edge voids\n", numEdge); @@ -981,6 +980,7 @@ void openFiles(string outputDir, string sampleName, int mockIndex, int numKept, FILE** fpZobov, FILE** fpCenters, FILE** fpBarycenter, FILE** fpShapes, + FILE** fpExtra, FILE** fpSkyPositions) { *fpZobov = fopen((outputDir+"/"+prefix+"voidDesc_"+dataPortion+"_"+sampleName).c_str(), "w"); @@ -997,6 +997,9 @@ void openFiles(string outputDir, string sampleName, *fpShapes = fopen((outputDir+"/"+prefix+"shapes_"+dataPortion+"_"+sampleName).c_str(), "w"); fprintf(*fpShapes, "# void ID, ellip, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n"); + + *fpExtra = fopen((outputDir+"/"+prefix+"extraInfo_"+dataPortion+"_"+sampleName).c_str(), "w"); + fprintf(*fpExtra, "# void type, max radius, nearest edge\n"); } // end openFiles @@ -1004,6 +1007,7 @@ void openFiles(string outputDir, string sampleName, // ---------------------------------------------------------------------------- void closeFiles(FILE* fpZobov, FILE* fpCenters, FILE* fpBarycenter, FILE* fpShapes, + FILE* fpExtra, FILE* fpSkyPositions) { fclose(fpZobov); @@ -1011,6 +1015,7 @@ void closeFiles(FILE* fpZobov, FILE* fpCenters, fclose(fpBarycenter); //fclose(fpDistances); fclose(fpShapes); + fclose(fpExtra); fclose(fpSkyPositions); } // end closeFile @@ -1025,13 +1030,13 @@ void outputVoids(string outputDir, string sampleName, string prefix, int iVoid; VOID outVoid; FILE *fp, *fpZobov, *fpCenters, *fpCentersNoCut, *fpBarycenter, - *fpShapes, *fpSkyPositions; + *fpShapes, *fpExtra, *fpSkyPositions; openFiles(outputDir, sampleName, prefix, dataPortion, mockIndex, voids.size(), &fpZobov, &fpCenters, &fpBarycenter, - &fpShapes, &fpSkyPositions); + &fpShapes, &fpExtra, &fpSkyPositions); for (iVoid = 0; iVoid < voids.size(); iVoid++) { @@ -1080,16 +1085,6 @@ void outputVoids(string outputDir, string sampleName, string prefix, outVoid.macrocenter[1], outVoid.macrocenter[2]); - /* - fprintf(fpDistances, "%d %e %e %e %e %e\n", - outVoid.voidID, - outVoid.nearestMock, - outVoid.radius, - outVoid.rescaledCoreDens, - outVoid.nearestMockFromCore, - outVoid.nearestGalFromCore); - */ - fprintf(fpCenters, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d %d %d %d %.2f\n", outCenter[0], outCenter[1], @@ -1139,9 +1134,14 @@ void outputVoids(string outputDir, string sampleName, string prefix, gsl_matrix_get(outVoid.evec, 2 ,2) ); + fprintf(fpExtra, "%d %.5f %.5f\n", + outVoid.voidType, + outVoid.maxRadius, + outVoid.nearestFlag + ); } // end iVoid closeFiles(fpZobov, fpCenters, fpBarycenter, - fpShapes, fpSkyPositions); + fpShapes, fpExtra, fpSkyPositions); } // end outputVoids diff --git a/python_source/voidUtil/catalogUtil.py b/python_source/voidUtil/catalogUtil.py index c189bb3..1fa3d74 100644 --- a/python_source/voidUtil/catalogUtil.py +++ b/python_source/voidUtil/catalogUtil.py @@ -367,7 +367,8 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True, numPart = int(line[8]), densCon = line[9], voidProb = line[10], - radius = 0., # this is read in later + # below values to be read in or computed later + radius = 0., macrocenter = np.zeros((3)), redshift = 0, RA = 0, @@ -379,6 +380,9 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True, ellipticity = 0., eigenVals = np.zeros((3)), eigenVecs = np.zeros((3,3)), + voidType = CENTRAL_VOID, + maxRadius = 0., + nearestEdge = 0. )) catalog.numVoids = len(catalog.voids) @@ -438,6 +442,19 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True, iLine += 1 + fileName = sampleDir+"/"+prefix+"extraInfo_"+dataPortion+"_"+sample.fullName+".out" + if os.path.exists(fileName): + catData = np.loadtxt(filename, comments="#") + for (iLine,line) in enumerate(cataData): + catalog.voids[iLine].voidType = int(line[0]) + catalog.voids[iLine].maxRadius = float(line[1]) + catalog.voids[iLine].nearestEdge = float(line[2]) + + iLine += 1 + + else: + print(" Extra info file not found!") + if loadParticles: print("Loading all particles...") partData, boxLen, volNorm, isObservationData, ranges, extraData, mockPart = loadPart(sampleDir) @@ -473,7 +490,7 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True, #catalog.part[:].edgeFlags = 0 print(" Loading volumes...") - if sample.hasWeightedVolumes: + if hasattr(sample, "hasWeightedVolumes") and sample.hasWeightedVolumes: volFile = sampleDir+"/vol_weighted_"+sample.fullName+".dat" else: volFile = sampleDir+"/vol_"+sample.fullName+".dat" @@ -589,12 +606,13 @@ def getVoidPart(catalog, voidID): return partOut # ----------------------------------------------------------------------------- +# various handy catalog filtering routines + def filterVoidsOnSize(catalog, rMin): catalog.voids = [v for v in catalog.voids if v.radius >= rMin] catalog.numVoids = len(catalog.voids) return catalog -# ----------------------------------------------------------------------------- def filterVoidsOnTreeLevel(catalog, level): catalog.voids = [v for v in catalog.voids if v.treeLevel == level] if level == -1: @@ -602,24 +620,31 @@ def filterVoidsOnTreeLevel(catalog, level): catalog.numVoids = len(catalog.voids) return catalog -# ----------------------------------------------------------------------------- def filterVoidsOnCentralDen(catalog, maxCentralDen): catalog.voids = [v for v in catalog.voids if v.centralDen <= maxCentralDen] catalog.numVoids = len(catalog.voids) return catalog -# ----------------------------------------------------------------------------- def filterVoidsOnCoreDen(catalog, maxCoreDen): catalog.voids = [v for v in catalog.voids if v.coreDens <= maxCoreDen] catalog.numVoids = len(catalog.voids) return catalog -# ----------------------------------------------------------------------------- def filterVoidsOnDenCon(catalog, minDenCon): catalog.voids = [v for v in catalog.voids if v.densCon >= minDenCon] catalog.numVoids = len(catalog.voids) return catalog +def filterVoidsOnType(catalog, voidType): + catalog.voids = [v for v in catalog.voids if v.voidType == voidType] + catalog.numVoids = len(catalog.voids) + return catalog + +def filterVoidsOnNearestEdge(catalog, factor=1.0): + catalog.voids = [v for v in catalog.voids if \ + factor*v.maxRadius <= v.nearestEdge] + catalog.numVoids = len(catalog.voids) + return catalog # -----------------------------------------------------------------------------