extra info in void catalog output (void type, max radius, nearest flagged galaxy), added filtering routines on these properties for post-analysis

This commit is contained in:
Paul M. Sutter 2025-04-21 15:02:05 -04:00
parent f9fc9e8990
commit dadfd8d40a
2 changed files with 70 additions and 45 deletions

View file

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

View file

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