mirror of
https://bitbucket.org/cosmicvoids/vide_public.git
synced 2025-07-04 23:31:12 +00:00
carried boundary flags all the way through to pruning; cleaned up pruning routines
This commit is contained in:
Paul M. Sutter
parent
dd181da42a
commit
a45eca0b6e
5 changed files with 124 additions and 139 deletions
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue