From a45eca0b6e71347cbf2d4dfdd0c0ae74cbcdbb6e Mon Sep 17 00:00:00 2001 From: "Paul M. Sutter" Date: Thu, 6 Jun 2024 20:26:47 +0200 Subject: [PATCH] carried boundary flags all the way through to pruning; cleaned up pruning routines --- c_source/pruning/pruneVoids.cpp | 224 +++++++----------- c_source/pruning/pruneVoids.ggo | 2 + .../example_observation.py | 2 +- python_source/backend/launchers.py | 26 ++ python_source/backend/surveyTools.py | 9 +- 5 files changed, 124 insertions(+), 139 deletions(-) diff --git a/c_source/pruning/pruneVoids.cpp b/c_source/pruning/pruneVoids.cpp index 34be238..adad968 100644 --- a/c_source/pruning/pruneVoids.cpp +++ b/c_source/pruning/pruneVoids.cpp @@ -59,6 +59,7 @@ typedef struct partStruct { float x, y, z, vol; int nadj, ncnt; int *adj; + int edgeFlag; } PART; typedef struct zoneStruct { @@ -390,7 +391,7 @@ int main(int argc, char **argv) { } } - // and finally volumes + // and volumes printf(" Loading particle volumes...\n"); fp = fopen(args.partVol_arg, "r"); fread(&mask_index, 1, 4, fp); @@ -404,7 +405,16 @@ int main(int argc, char **argv) { } fclose(fp); -/* + // and finally edge flag info + printf(" Loading particle edge flags...\n"); + fp = fopen(args.partEdge_arg, "r"); + for (p = 0; p < mask_index; p++) { + fscanf(fp, "%d", &part[p].edgeFlag); + //printf("EDGE VALUE %d\n", part[p].edgeFlag); + } + fclose(fp); + +/* this was used for testing at one point // and finally finally adjacencies printf(" Loading particle adjacencies...\n"); fp = fopen(args.partAdj_arg, "r"); @@ -458,61 +468,59 @@ int main(int argc, char **argv) { interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; printf(" Read voids (%.2f sec)...\n", interval); - // load voids *again* using Guilhem's code so we can get tree + // load voids *again* using Guilhem's code so we can get tree information clock3 = clock(); - //if (!args.isObservation_flag) { - printf(" Re-loading voids and building tree..\n"); - ZobovRep zobovCat; - if (!loadZobov(args.voidDesc_arg, args.zone2Part_arg, - args.void2Zone_arg, - 0, zobovCat)) { + printf(" Re-loading voids and building tree..\n"); + ZobovRep zobovCat; + if (!loadZobov(args.voidDesc_arg, args.zone2Part_arg, args.void2Zone_arg, + 0, zobovCat)) { printf("Error loading catalog!\n"); return -1; - } - VoidTree *tree; - tree = new VoidTree(zobovCat); - zobovCat.allZones.erase(zobovCat.allZones.begin(), zobovCat.allZones.end()); + } + VoidTree *tree; + tree = new VoidTree(zobovCat); + zobovCat.allZones.erase(zobovCat.allZones.begin(), zobovCat.allZones.end()); - // copy tree information to our own data structures - for (iVoid = 0; iVoid < numVoids; iVoid++) { - voidID = voids[iVoid].voidID; - voids[iVoid].parentID = tree->getParent(voidID); - voids[iVoid].numChildren = tree->getChildren(voidID).size(); + // copy tree information to our own data structures + for (iVoid = 0; iVoid < numVoids; iVoid++) { + voidID = voids[iVoid].voidID; + voids[iVoid].parentID = tree->getParent(voidID); + voids[iVoid].numChildren = tree->getChildren(voidID).size(); - // compute level in tree - int level = 0; - int parentID = tree->getParent(voidID); - while (parentID != -1) { - level++; - parentID = tree->getParent(parentID); - } - voids[iVoid].level = level; + // compute level in tree + int level = 0; + int parentID = tree->getParent(voidID); + while (parentID != -1) { + level++; + parentID = tree->getParent(parentID); } - //} // end re-load + voids[iVoid].level = level; + } clock4 = clock(); interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; printf(" Re-read voids (%.2f sec)...\n", interval); - // check boundaries printf(" Computing void properties...\n"); + // allocate space for a particle buffer maxNumPart = 0; for (iVoid = 0; iVoid < numVoids; iVoid++) { if (voids[iVoid].numPart > maxNumPart) maxNumPart = voids[iVoid].numPart; } - voidPart = (PART *) malloc(maxNumPart * sizeof(PART)); + // main processing of each void for (iVoid = 0; iVoid < numVoids; iVoid++) { - voidID = voids[iVoid].voidID; - printf(" DOING %d (of %d) %d %d %f\n", iVoid+1, numVoids, voidID, + printf(" Working on void %d (of %d) %d %d %f\n",iVoid+1, numVoids, voidID, voids[iVoid].numPart, voids[iVoid].radius); voids[iVoid].center[0] = part[voids[iVoid].coreParticle].x; voids[iVoid].center[1] = part[voids[iVoid].coreParticle].y; voids[iVoid].center[2] = part[voids[iVoid].coreParticle].z; + + voids[iVoid].voidType = CENTRAL_VOID; // first load up particles into a buffer clock3 = clock(); @@ -523,6 +531,7 @@ int main(int argc, char **argv) { for (p = 0; p < zones2Parts[zoneID].numPart; p++) { partID = zones2Parts[zoneID].partIDs[p]; + // something went haywire if (partID > mask_index || (part[partID].vol < 1.e-27 && part[partID].vol > 0.)) { printf("BAD PART!? %d %d %e", partID, mask_index, part[partID].vol); @@ -533,6 +542,12 @@ int main(int argc, char **argv) { voidPart[i].y = part[partID].y; voidPart[i].z = part[partID].z; voidPart[i].vol = part[partID].vol; + voidPart[i].edgeFlag = part[partID].edgeFlag; + + // check for edge contamination + if (voidPart[i].edgeFlag > 0) { + voids[iVoid].voidType = EDGE_VOID; + } /* // testing for edge contamination @@ -549,7 +564,7 @@ int main(int argc, char **argv) { */ i++; } - } + } // loading particles clock4 = clock(); interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; @@ -632,75 +647,30 @@ int main(int argc, char **argv) { interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; //printf(" %.2f for central density\n", interval); - //coreParticle = voids[iVoid].coreParticle; - //voids[iVoid].rescaledCoreDens = voids[iVoid].coreDens*(pow(1.*mockIndex/numPartTot,3)); - // // compute distance from core to nearest mock - // minDist = 1.e99; - // for (p = mockIndex; p < numPartTot; p++) { - // dist[0] = part[coreParticle].x - part[p].x; - // dist[1] = part[coreParticle].y - part[p].y; - // dist[2] = part[coreParticle].z - part[p].z; - // - // dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); - // if (dist2 < minDist) minDist = dist2; - // } - // voids[iVoid].nearestMockFromCore = sqrt(minDist); - // - // // compute distance from core to nearest mock - // minDist = 1.e99; - // for (p = 0; p < mockIndex; p++) { - // dist[0] = part[coreParticle].x - part[p].x; - // dist[1] = part[coreParticle].y - part[p].y; - // dist[2] = part[coreParticle].z - part[p].z; - // - // dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); - // if (dist2 < minDist && dist2 > 1.e-10) minDist = dist2; - // } - // voids[iVoid].nearestGalFromCore = sqrt(minDist); + // compute maximum extent of void + clock3 = clock(); + maxDist = 0.; + for (p = 0; p < voids[iVoid].numPart; p++) { + dist[0] = fabs(voidPart[p].x - voids[iVoid].macrocenter[0]); + dist[1] = fabs(voidPart[p].y - voids[iVoid].macrocenter[1]); + dist[2] = fabs(voidPart[p].z - voids[iVoid].macrocenter[2]); - // compute maximum extent -/* - if (args.isObservation_flag) { - maxDist = 0.; - for (p = 0; p < voids[iVoid].numPart; p++) { - for (p2 = p; p2 < voids[iVoid].numPart; p2++) { - - dist[0] = voidPart[p].x - voidPart[p2].x; - dist[1] = voidPart[p].y - voidPart[p2].y; - dist[2] = voidPart[p].z - voidPart[p2].z; + if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]); + if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]); + if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]); - dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); - if (dist2 > maxDist) maxDist = dist2; - } - } - voids[iVoid].maxRadius = sqrt(maxDist)/2.; - } else { -*/ - - clock3 = clock(); - maxDist = 0.; - for (p = 0; p < voids[iVoid].numPart; p++) { - - dist[0] = fabs(voidPart[p].x - voids[iVoid].macrocenter[0]); - dist[1] = fabs(voidPart[p].y - voids[iVoid].macrocenter[1]); - dist[2] = fabs(voidPart[p].z - voids[iVoid].macrocenter[2]); - - if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]); - if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]); - if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]); - - dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); - if (dist2 > maxDist) maxDist = dist2; - } - voids[iVoid].maxRadius = sqrt(maxDist); - clock4 = clock(); - interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; - //printf(" %.2f for maximum extent\n", interval); -// } + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + if (dist2 > maxDist) maxDist = dist2; + } + voids[iVoid].maxRadius = sqrt(maxDist); + 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(); if (args.isObservation_flag) { - // compute distance from center to nearest mock minDist = 1.e99; for (p = mockIndex; p < numPartTot; p++) { dist[0] = voids[iVoid].macrocenter[0] - part[p].x; @@ -727,7 +697,6 @@ int main(int argc, char **argv) { if (args.useComoving_flag) { redshift = gsl_interp_eval(interp, dL, redshifts, voids[iVoid].redshiftInMpc, acc); - //printf("HELLO %e %e\n", redshift, args.zMax_arg); nearestEdge = fabs((redshift-args.zMax_arg)*LIGHT_SPEED/100.); voids[iVoid].redshift = redshift; } else { @@ -796,7 +765,7 @@ int main(int argc, char **argv) { for (int i = 0; i < 3; i++) { for (int j = i; j < 3; j++) { if (i == j) inertia[i*3+j] += dist[0]*dist[0] + dist[1]*dist[1] + - dist[2]*dist[2]; + dist[2]*dist[2]; inertia[i*3+j] -= dist[i]*dist[j]; } } @@ -828,24 +797,16 @@ int main(int argc, char **argv) { if (gsl_vector_get(voids[iVoid].eval,i) > largest) largest = gsl_vector_get(voids[iVoid].eval,i); } - // TEST voids[iVoid].ellip = 1.0 - sqrt(sqrt(fabs(smallest/largest))); - //if (a < c) ca = a/c; - //if (a >= c) ca = c/a; - //voids[iVoid].ellip = fabs(1.0 - ca); - - //if (a < c) ca = a*a/(c*c); - //if (a >= c) ca = (c*c)/(a*a); - //voids[iVoid].ellip = sqrt(fabs(1.0 - ca)); - clock4 = clock(); interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; - //printf(" %.2f for ellipticity\n", interval); } // iVoid gsl_eigen_symmv_free(eigw); - + + + // now filter and categorize the voids based on various criteria int numWrong = 0; int numHighDen = 0; int numCentral = 0; @@ -861,15 +822,6 @@ int main(int argc, char **argv) { voids[iVoid].accepted = 1; } -/* - int j = 0; - for (iVoid = 0; iVoid < voids.size(); iVoid++) { - if (voids[iVoid].densCon < 1.5) { -// voids[iVoid].accepted = -4; - } - } -*/ - // toss out voids that are obviously wrong int iGood = 0; for (iVoid = 0; iVoid < voids.size(); iVoid++) { @@ -883,6 +835,7 @@ int main(int argc, char **argv) { voids.resize(iGood); printf(" 1st filter: rejected %d obviously bad\n", numWrong); + // toss out voids that are too small iGood = 0; for (iVoid = 0; iVoid < voids.size(); iVoid++) { if (voids[iVoid].radius < args.rMin_arg) { @@ -895,9 +848,9 @@ int main(int argc, char **argv) { printf(" 2nd filter: rejected %d too small\n", numTooSmall); + // toss out voids near non-periodic box edges iGood = 0; for (iVoid = 0; iVoid < voids.size(); iVoid++) { - // *always* clean out near edges since there are no mocks there if (tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestEdge || tolerance*voids[iVoid].radius > voids[iVoid].nearestEdge) { numNearZ++; @@ -908,6 +861,7 @@ int main(int argc, char **argv) { voids.resize(iGood); printf(" 3rd filter: rejected %d too close to high redshift boundaries\n", numNearZ); + // toss out voids that are beyond redshift boundaries numNearZ = 0; iGood = 0; for (iVoid = 0; iVoid < voids.size(); iVoid++) { @@ -921,24 +875,18 @@ int main(int argc, char **argv) { } voids.resize(iGood); - //Maubert - Uncommented this part : to be sure that voids do not cross maximum redshift asked for in zrange iGood = 0; for (iVoid = 0; iVoid < voids.size(); iVoid++) { - // just in case - if (args.isObservation_flag && - voids[iVoid].redshift > args.zMax_arg) { + if (args.isObservation_flag && voids[iVoid].redshift > args.zMax_arg) { numNearZ++; } else { voids[iGood++] = voids[iVoid]; } } voids.resize(iGood); - // Maubert - End of Uncommented part - - printf(" 4th filter: rejected %d outside redshift boundaries\n", numNearZ); - // take only top-level voids + // find top-level voids numAreParents = 0; iGood = 0; for (iVoid = 0; iVoid < voids.size(); iVoid++) { @@ -950,7 +898,7 @@ int main(int argc, char **argv) { } } - + // mark high-density voids for (iVoid = 0; iVoid < voids.size(); iVoid++) { if (voids[iVoid].centralDen > args.maxCentralDen_arg) { voids[iVoid].accepted = -1; @@ -961,6 +909,17 @@ int main(int argc, char **argv) { } } + + // count voids near survey edges + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + if (voids[iVoid].voidType == CENTRAL_VOID) { + numCentral++; + } else { + 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; @@ -970,6 +929,7 @@ int main(int argc, char **argv) { numEdge++; } } +*/ printf(" Number kept: %d (out of %d)\n", (int) voids.size(), numVoids); printf(" We have %d edge voids\n", numEdge); @@ -988,7 +948,6 @@ int main(int argc, char **argv) { prefix = ""; for (int i = 0; i < 2; i++) { dataPortion = dataPortions[i]; - outputVoids(outputDir, sampleName, prefix, dataPortion, mockIndex, voids, @@ -999,7 +958,6 @@ int main(int argc, char **argv) { prefix = "untrimmed_"; for (int i = 0; i < 2; i++) { dataPortion = dataPortions[i]; - outputVoids(outputDir, sampleName, prefix, dataPortion, mockIndex, voids, @@ -1010,7 +968,6 @@ int main(int argc, char **argv) { prefix = "untrimmed_dencut_"; for (int i = 0; i < 2; i++) { dataPortion = dataPortions[i]; - outputVoids(outputDir, sampleName, prefix, dataPortion, mockIndex, voids, @@ -1020,7 +977,6 @@ int main(int argc, char **argv) { prefix = "trimmed_nodencut_"; for (int i = 0; i < 2; i++) { dataPortion = dataPortions[i]; - outputVoids(outputDir, sampleName, prefix, dataPortion, mockIndex, voids, @@ -1207,7 +1163,7 @@ void outputVoids(string outputDir, string sampleName, string prefix, } // end iVoid - closeFiles(fpZobov, fpCenters, fpBarycenter, - fpDistances, fpShapes, fpSkyPositions); + closeFiles(fpZobov, fpCenters, fpBarycenter, + fpDistances, fpShapes, fpSkyPositions); } // end outputVoids diff --git a/c_source/pruning/pruneVoids.ggo b/c_source/pruning/pruneVoids.ggo index 77f4521..1f02fc3 100644 --- a/c_source/pruning/pruneVoids.ggo +++ b/c_source/pruning/pruneVoids.ggo @@ -16,6 +16,8 @@ option "partVol" - "Particle volume file from ZOBOV" string required option "partAdj" - "Adjacency file from ZOBOV" string required +option "partEdge" - "Boundary flag file" string required + option "zone2Part" - "Particle file from ZOBOV" string required option "mockIndex" - "Beginning index of mock particles" int required diff --git a/examples/example_observation/example_observation.py b/examples/example_observation/example_observation.py index 2c7da8d..34ac582 100644 --- a/examples/example_observation/example_observation.py +++ b/examples/example_observation/example_observation.py @@ -30,7 +30,7 @@ continueRun = False # 1 : extract redshift slices from data # 2 : void extraction using zobov # 3 : removal of small voids and voids near the edge -startCatalogStage = 1 +startCatalogStage = 3 endCatalogStage = 3 basePath = os.path.dirname(os.path.abspath(__file__)) diff --git a/python_source/backend/launchers.py b/python_source/backend/launchers.py index 5465878..10f7436 100644 --- a/python_source/backend/launchers.py +++ b/python_source/backend/launchers.py @@ -446,6 +446,31 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None, else: volFileToUse = outputDir+"/vol_"+sampleName+".dat" + # re-weight the volumes of any edge galaxies to prevent watershed + # from spilling outside of survey region + if sample.dataType == "observation": + + # read in the appropriate volume file + with open(volFileToUse, mode="rb") as File: + numPartTot = np.fromfile(File, dtype=np.int32,count=1)[0] + vols = np.fromfile(File, dtype=np.float32, count=numPartTot) + + # read in the edge flag information + edgeFile = outputDir+"/galaxy_edge_flags.txt" + edgeFlags = np.loadtxt(edgeFile, dtype=np.int32) + + # set edge galaxy volumes to nearly 0 (implying very high density) + vols[ edgeFlags>0 ] = 1.e-4 + + volFile = outputDir+"/vol_weighted_"+sampleName+".dat" + with open(volFile, mode='wb') as File: + numPartTot.astype(np.int32).tofile(File) + vols.astype(np.float32).tofile(File) + + volFileToUse = outputDir+"/vol_weighted_"+sampleName+".dat" + else: + volFileToUse = outputDir+"/vol_"+sampleName+".dat" + cmd = [binPath+"/jozov2", \ outputDir+"/adj_"+sampleName+".dat", \ @@ -525,6 +550,7 @@ def launchPrune(sample, binPath, cmd += " --zone2Part=" + outputDir+"/voidPart_"+str(sampleName)+".dat" cmd += " --partVol=" + outputDir+"/vol_"+str(sampleName)+".dat" cmd += " --partAdj=" + outputDir+"/adj_"+str(sampleName)+".dat" + cmd += " --partEdge=" + outputDir+"galaxy_edge_flags.txt" cmd += " --extraInfo=" + outputDir+"/zobov_slice_"+str(sampleName)+\ ".par" cmd += " --tolerance=" + str(boundaryTolerance) diff --git a/python_source/backend/surveyTools.py b/python_source/backend/surveyTools.py index b1e32e5..798b143 100644 --- a/python_source/backend/surveyTools.py +++ b/python_source/backend/surveyTools.py @@ -161,18 +161,19 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, edgeMaskFile, phi, theta = convertAngle(RA, Dec) # check the mask edges - neighbors = healpy.get_all_neighbours(nside, theta, phi) - isOnMaskEdge = any(p == 0 for p in neighbors) + ipix = healpy.ang2pix(nside, theta, phi) + neighbors = healpy.get_all_neighbours(nside, ipix) + isOnMaskEdge = any(mask[p] == 0 for p in neighbors) # check the redshift boundaries - tol = 0.01 # TODO: made this user-adjustable + tol = 0.05 # TODO: made this user-adjustable zbuffer = (zmax-zmin)*tol isOnHighZEdge = (z >= zmax-zbuffer) isOnLowZEdge = (z <= zmin+zbuffer) - print("DOING %f %f %f %f\n" % (zbuffer, z, zmax, zmin) ) if isOnMaskEdge: edgeFile.write("1\n") + edgeMask[ipix] = 1 elif isOnHighZEdge: edgeFile.write("2\n") elif isOnLowZEdge: