guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 15:21:11 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions
vide_public/c_source/pruning/pruneVoids.cpp
Paul M. Sutter 83d0d66d90 fixed accidental deletion of redshift calculation
2025-04-22 23:00:52 -04:00

1133 lines
37 KiB
C++
Raw Blame History

/*+
VIDE -- Void IDentification and Examination -- ./c_tools/stacking/pruneVoids.cpp
Copyright (C) 2010-2014 Guilhem Lavaux
Copyright (C) 2011-2014 P. M. Sutter
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+*/
// Reads in the void catalog and removes any void that could potentially
// be affected by a mock particle. It does this by computing the longest
// particle distance within each void and comparing it to the distance
// of the nearest mock particle. If the void could potentially by rotated
// to include this particle, we throw out the void.
// This is placed here instead of using the edgeAvoidance option in
// stackVoidsZero so that we can optionally filter the entire
// catalog at once before the stacking phase. This is useful
// for producing a "clean" void catalog for other purposes.
#include "gsl/gsl_math.h"
#include "gsl/gsl_eigen.h"
#include "string.h"
#include "ctype.h"
#include "stdlib.h"
#include <math.h>
#include <stdio.h>
#include <netcdf>
#include "pruneVoids_conf.h"
#include <vector>
#include "assert.h"
#include "voidTree.hpp"
#include "loadZobov.hpp"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#define LIGHT_SPEED 299792.458
#define MPC2Z 100./LIGHT_SPEED
#define Z2MPC LIGHT_SPEED/100.
#define CENTRAL_VOID 1
#define EDGE_VOID 2
using namespace std;
typedef struct partStruct {
float x, y, z, vol;
int nadj, ncnt;
int *adj;
int edgeFlag;
} PART;
typedef struct zoneStruct {
int numPart;
int *partIDs;
} ZONE2PART;
typedef struct voidZoneStruct {
int numZones;
int *zoneIDs;
} VOID2ZONE;
typedef struct voidStruct {
float vol, coreDens, zoneVol, densCon, voidProb, radius;
float rescaledCoreDens;
int voidID, numPart, numZones, coreParticle, zoneNumPart;
float maxRadius, nearestFlag, centralDen, redshift, redshiftInMpc;
float nearestEdge;
float center[3], macrocenter[3];
int accepted;
int voidType;
int parentIndx, numChildren, level;
bool isLeaf, hasHighCentralDen;
gsl_vector *eval;
gsl_matrix *evec;
float ellip;
} VOID;
// this defines the expansion function that we will integrate
// Laveaux & Wandelt (2012) Eq. 24
struct my_expan_params { double Om; double w0; double wa; };
double expanFun (double z, void * p) {
struct my_expan_params * params = (struct my_expan_params *)p;
double Om = (params->Om);
double w0 = (params->w0);
double wa = (params->wa);
//const double h0 = 1.0;
const double h0 = 0.71;
double ez;
double wz = w0 + wa*z/(1+z);
ez = Om*pow(1+z,3) + (1.-Om);
//ez = Om*pow(1+z,3) + pow(h0,2) * (1.-Om)*pow(1+z,3+3*wz);
ez = sqrt(ez);
//ez = sqrt(ez)/h0;
ez = 1./ez;
return ez;
}
void openFiles(string outputDir, string sampleName,
string prefix, string dataPortion,
int numPartTot, int numKept,
FILE** fpZobov, FILE** fpCenters,
FILE** fpCentersNoCut,
FILE** fpBarycenter, FILE** fpShapes,
FILE** fpSkyPositions);
void closeFiles(FILE* fpZobov, FILE* fpCenters,
FILE* fpCentersNoCut,
FILE* fpBarycenter, FILE* fpShapes,
FILE* fpSkyPositions);
void outputVoids(string outputDir, string sampleName, string prefix,
string dataPortion, int numPartTot,
vector<VOID> voids,
bool isObservation, double *boxLen,
bool doTrim, bool doCentralDenCut);
int main(int argc, char **argv) {
// initialize arguments
pruneVoids_info args;
pruneVoids_conf_params args_params;
pruneVoids_conf_init(&args);
pruneVoids_conf_params_init(&args_params);
args_params.check_required = 0;
if (pruneVoids_conf_ext (argc, argv, &args, &args_params))
return 1;
if (!args.configFile_given) {
if (pruneVoids_conf_required (&args,
PRUNEVOIDS_CONF_PACKAGE))
return 1;
} else {
args_params.check_required = 1;
args_params.initialize = 0;
if (pruneVoids_conf_config_file (args.configFile_arg,
&args,
&args_params))
return 1;
}
// initialize cosmology integrator and interpolator
gsl_function expanF;
expanF.function = &expanFun;
struct my_expan_params expanParams;
expanParams.Om = args.omegaM_arg;
expanParams.w0 = -1.0;
expanParams.wa = 0.0;
expanF.params = &expanParams;
double result, error;
size_t nEval;
int iZ, numZ = 8000;
double maxZ = 10.0, z, *dL, *redshifts;
dL = (double *) malloc(numZ * sizeof(double));
redshifts = (double *) malloc(numZ * sizeof(double));
for (iZ = 0; iZ < numZ; iZ++) {
z = iZ * maxZ/numZ;
gsl_integration_qng(&expanF, 0.0, z, 1.e-6, 1.e-6, &result, &error, &nEval);
dL[iZ] = result*LIGHT_SPEED/100.;
//printf("HERE %e %e\n", z, dL[iZ]);
redshifts[iZ] = z;
}
gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, numZ);
gsl_interp_init(interp, dL, redshifts, numZ);
gsl_interp_accel *acc = gsl_interp_accel_alloc();
int i, p, p2, numPartTot, numZonesTot, dummy, iVoid;
int numVoids;
FILE *fp;
PART *part, *voidPart, *flaggedPart;
ZONE2PART *zones2Parts;
VOID2ZONE *void2Zones;
vector<VOID> voids;
float *temp, junk, voidVol;
int junkInt, voidID, numPart, numZones, zoneID, partID, maxNumPart;
int coreParticle, zoneNumPart;
float coreDens, zoneVol, densCon, voidProb, dist[3], dist2, minDist, maxDist;
float centralRad, centralDen;
double nearestEdge, redshift;
char line[500], junkStr[10];
string outputDir, sampleName, dataPortion, prefix;
double ranges[3][2], boxLen[3], mul;
double volNorm, radius;
int clock1, clock2, clock3, clock4;
double interval;
int periodicX=0, periodicY=0, periodicZ=0;
string dataPortions[2];
gsl_eigen_symmv_workspace *eigw = gsl_eigen_symmv_alloc(3);
numVoids = args.numVoids_arg;
clock1 = clock();
// check for periodic box
periodicX = 0;
periodicY = 0;
periodicZ = 0;
if (!args.isObservation_flag) {
if ( strchr(args.periodic_arg, 'x') != NULL) {
periodicX = 1;
printf("Will assume x-direction is periodic.\n");
}
if ( strchr(args.periodic_arg, 'y') != NULL) {
periodicY = 1;
printf("Will assume y-direction is periodic.\n");
}
if ( strchr(args.periodic_arg, 'z') != NULL) {
periodicZ = 1;
printf("Will assume z-direction is periodic.\n");
}
}
// load box size
printf("\n Getting info...\n");
netCDF::NcFile f_info(args.extraInfo_arg, netCDF::NcFile::read);
f_info.getAtt("range_x_min").getValues(&ranges[0][0]);
f_info.getAtt("range_x_max").getValues(&ranges[0][1]);
f_info.getAtt("range_y_min").getValues(&ranges[1][0]);
f_info.getAtt("range_y_max").getValues(&ranges[1][1]);
f_info.getAtt("range_z_min").getValues(&ranges[2][0]);
f_info.getAtt("range_z_max").getValues(&ranges[2][1]);
printf(" Range xmin %e\n", ranges[0][0]);
printf(" Range xmax %e\n", ranges[0][1]);
printf(" Range ymin %e\n", ranges[1][0]);
printf(" Range ymax %e\n", ranges[1][1]);
printf(" Range zmin %e\n", ranges[2][0]);
printf(" Range zmax %e\n", ranges[2][1]);
boxLen[0] = ranges[0][1] - ranges[0][0];
boxLen[1] = ranges[1][1] - ranges[1][0];
boxLen[2] = ranges[2][1] - ranges[2][0];
// read in all particle positions
clock3 = clock();
printf("\n Loading particles...\n");
fp = fopen(args.partFile_arg, "r");
fread(&dummy, 1, 4, fp);
fread(&numPartTot, 1, 4, fp);
fread(&dummy, 1, 4, fp);
part = (PART *) malloc(numPartTot * sizeof(PART));
temp = (float *) malloc(numPartTot * sizeof(float));
// this is the volume normalization understood by zobov,
// which is all particles (including mock boundary galaxies used in old
// versions of VIDE) divided by the total volume of the bounding box
volNorm = numPartTot/(boxLen[0]*boxLen[1]*boxLen[2]);
fread(&dummy, 1, 4, fp);
fread(temp, numPartTot, 4, fp);
mul = ranges[0][1] - ranges[0][0];
for (p = 0; p < numPartTot; p++)
part[p].x = mul*temp[p];
fread(&dummy, 1, 4, fp);
fread(&dummy, 1, 4, fp);
fread(temp, numPartTot, 4, fp);
mul = ranges[1][1] - ranges[1][0];
for (p = 0; p < numPartTot; p++)
part[p].y = mul*temp[p];
fread(&dummy, 1, 4, fp);
fread(&dummy, 1, 4, fp);
fread(temp, numPartTot, 4, fp);
mul = ranges[2][1] - ranges[2][0];
for (p = 0; p < numPartTot; p++)
part[p].z = mul*temp[p];
if (!args.isObservation_flag) {
for (p = 0; p < numPartTot; p++) {
part[p].x += ranges[0][0];
part[p].y += ranges[1][0];
part[p].z += ranges[2][0];
}
}
fclose(fp);
clock4 = clock();
interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
printf(" Read %d particles (%.2f sec)...\n", numPartTot, interval);
// read in desired voids
clock3 = clock();
printf(" Loading voids...\n");
fp = fopen(args.voidDesc_arg ,"r");
fgets(line, sizeof(line), fp);
sscanf(line, "%d %s %d %s", &junkInt, junkStr, &junkInt, junkStr);
fgets(line, sizeof(line), fp);
voids.resize(numVoids);
i = 0;
while (fgets(line, sizeof(line), fp) != NULL) {
sscanf(line, "%d %d %d %f %f %d %d %f %d %f %f\n", &iVoid, &voidID,
&coreParticle, &coreDens, &zoneVol, &zoneNumPart, &numZones,
&voidVol, &numPart, &densCon, &voidProb);
i++;
voids[i-1].coreParticle = coreParticle;
voids[i-1].zoneNumPart = zoneNumPart;
voids[i-1].coreDens = coreDens;
voids[i-1].zoneVol = zoneVol;
voids[i-1].voidID = voidID;
voids[i-1].vol = voidVol;
voids[i-1].numPart = numPart;
voids[i-1].numZones = numZones;
voids[i-1].densCon = densCon;
voids[i-1].voidProb = voidProb;
voids[i-1].radius = pow(voidVol/volNorm*3./4./M_PI, 1./3.);
voids[i-1].accepted = 1;
voids[i-1].isLeaf = true;
voids[i-1].hasHighCentralDen = false;
voids[i-1].parentIndx = -1;
voids[i-1].numChildren = 0;
voids[i-1].level = 0;
voids[i-1].eval = gsl_vector_alloc(3);
voids[i-1].evec = gsl_matrix_alloc(3,3);
voids[i-1].ellip = 0;
}
fclose(fp);
// load up the zone membership for each void
printf(" Loading void-zone membership info...\n");
fp = fopen(args.void2Zone_arg, "r");
fread(&numZonesTot, 1, 4, fp);
void2Zones = (VOID2ZONE *) malloc(numZonesTot * sizeof(VOID2ZONE));
for (iZ = 0; iZ < numZonesTot; iZ++) {
fread(&numZones, 1, 4, fp);
void2Zones[iZ].numZones = numZones;
void2Zones[iZ].zoneIDs = (int *) malloc(numZones * sizeof(int));
for (p = 0; p < numZones; p++) {
fread(&void2Zones[iZ].zoneIDs[p], 1, 4, fp);
}
}
fclose(fp);
// now the particles-zone
printf(" Loading zone-particle membership info...\n");
fp = fopen(args.zone2Part_arg, "r");
fread(&dummy, 1, 4, fp);
fread(&numZonesTot, 1, 4, fp);
zones2Parts = (ZONE2PART *) malloc(numZonesTot * sizeof(ZONE2PART));
for (iZ = 0; iZ < numZonesTot; iZ++) {
fread(&numPart, 1, 4, fp);
zones2Parts[iZ].numPart = numPart;
zones2Parts[iZ].partIDs = (int *) malloc(numPart * sizeof(int));
for (p = 0; p < numPart; p++) {
fread(&zones2Parts[iZ].partIDs[p], 1, 4, fp);
}
}
// and volumes
printf(" Loading particle volumes...\n");
fp = fopen(args.partVol_arg, "r");
fread(&numPartTot, 1, 4, fp);
for (p = 0; p < numPartTot; p++) {
fread(&temp[0], 1, 4, fp);
part[p].vol = temp[0];
}
fclose(fp);
// and finally edge flag info
printf(" Loading particle edge flags...\n");
int numFlagged = 0;
fp = fopen(args.partEdge_arg, "r");
for (p = 0; p < numPartTot; 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 < numPartTot; 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");
fp = fopen(args.partAdj_arg, "r");
fread(&numPartTot, 1, 4, fp);
if (numPartTot != numPartTot) {
printf("NON-MATCHING MOCK INDICES!? %d %d\n", numPartTot, numPartTot);
exit(-1);
}
int tempInt;
for (p = 0; p < numPartTot; p++) {
fread(&tempInt, 1, 4, fp);
part[p].nadj = tempInt;
part[p].ncnt = 0;
if (part[p].nadj > 0)
part[p].adj = (int *) malloc(part[p].nadj * sizeof(int));
}
for (p = 0; p < numPartTot; p++) {
fread(&tempInt, 1, 4, fp);
int nin = tempInt;
if (nin > 0) {
for (int nAdj = 0; nAdj < nin; nAdj++) {
int tempAdj;
fread(&tempAdj, 1, 4, fp);
// this bit has been readjusted just in case we are
// accidentally still linking to mock particles
//if (tempAdj < numPartTot) {
assert(p < tempAdj);
//if (part[p].ncnt == part[p].nadj) {
// printf("OVERFLOW %d\n", p);
//} else if (part[tempAdj].ncnt == part[tempAdj].nadj) {
// printf("OVERFLOW %d\n", tempAdj);
//} else {
part[p].adj[part[p].ncnt] = tempAdj;
part[p].ncnt++;
if (tempAdj < numPartTot) {
part[tempAdj].adj[part[tempAdj].ncnt] = p;
part[tempAdj].ncnt++;
}
//}
//}
}
//printf("ADJ %d %d %d %d %d\n", p, nin, part[p].nadj, nAdj, tempInt);
}
}
fclose(fp);
*/
clock4 = clock();
interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
printf(" Read voids (%.2f sec)...\n", interval);
clock3 = clock();
// TODO - TEST NEW TREE BUILDING TECHNIQUE AND REMOVE THIS
/*
// load voids *again* using Guilhem's code so we can get tree information
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());
*/
// if all of a void's zones also belong to another void,
// it is a child of that void
for (iVoid = 0; iVoid < numVoids; iVoid++) {
voidID = voids[iVoid].voidID;
int numMyZones = voids[iVoid].numZones;
for (int iCheck = 0; iCheck < numVoids; iCheck++) {
int numCheckZones = voids[iCheck].numZones;
if (numMyZones > numCheckZones) continue;
if (iVoid == iCheck) continue;
int checkID = voids[iCheck].voidID;
bool allZonesMatch = true;
for (iZ = 0; iZ < void2Zones[voidID].numZones; iZ++) {
int myZone = void2Zones[voidID].zoneIDs[iZ];
bool foundMatch = false;
for (int jZ = 0; jZ < void2Zones[checkID].numZones; jZ++) {
int checkZone = void2Zones[checkID].zoneIDs[jZ];
foundMatch = myZone == checkZone;
}
if (not foundMatch) {
allZonesMatch = false;
break;
}
}
if (allZonesMatch) {
voids[iVoid].parentIndx = iCheck;
voids[iCheck].numChildren++;
}
}
} // end building tree
// compute level in tree
for (iVoid = 0; iVoid < numVoids; iVoid++) {
int level = 0;
int parentIndx = voids[iVoid].parentIndx;
while (parentIndx != -1) {
level++;
parentIndx = voids[parentIndx].parentIndx;
}
voids[iVoid].level = level;
}
// TODO - TEST NEW TREE BUILDING TECHNIQUE AND REMOVE THIS
/*
// 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;
}
*/
clock4 = clock();
interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
printf(" Building void tree (%.2f sec)...\n", interval);
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;
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 this void's particles into a buffer
clock3 = clock();
i = 0;
for (iZ = 0; iZ < void2Zones[voidID].numZones; iZ++) {
zoneID = void2Zones[voidID].zoneIDs[iZ];
for (p = 0; p < zones2Parts[zoneID].numPart; p++) {
partID = zones2Parts[zoneID].partIDs[p];
// something went haywire
if (part[partID].vol < 1.e-27 && part[partID].vol > 0.) {
printf("BAD PART!? %d %d %e", partID, numPartTot, part[partID].vol);
exit(-1);
}
voidPart[i].x = part[partID].x;
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
if (part[partID].vol < 1.e-27) {
printf("CONTAMINATED!! %d %d\n", iVoid, partID);
} else {
//printf("NORMAL!! %d %d %e\n", iVoid, partID, part[partID].vol);
}
for (int iAdj = 0; iAdj < part[partID].ncnt; iAdj++) {
if (part[partID].adj[iAdj] > numPartTot) {
printf("CONTAMINATED!! %d %d %d\n", iVoid, partID, iAdj);
}
}
*/
i++;
}
} // loading particles
clock4 = clock();
interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
//printf(" %.2f for buffer\n", interval);
// compute macrocenters
clock3 = clock();
double weight = 0.;
voids[iVoid].macrocenter[0] = 0.;
voids[iVoid].macrocenter[1] = 0.;
voids[iVoid].macrocenter[2] = 0.;
for (p = 0; p < voids[iVoid].numPart; p++) {
dist[0] = voidPart[p].x - voids[iVoid].center[0];
dist[1] = voidPart[p].y - voids[iVoid].center[1];
dist[2] = voidPart[p].z - voids[iVoid].center[2];
if (periodicX && fabs(dist[0]) > boxLen[0]/2.)
dist[0] = dist[0] - copysign(boxLen[0], dist[0]);
if (periodicY && fabs(dist[1]) > boxLen[1]/2.)
dist[1] = dist[1] - copysign(boxLen[1], dist[1]);
if (periodicZ && fabs(dist[2]) > boxLen[2]/2.)
dist[2] = dist[2] - copysign(boxLen[2], dist[2]);
voids[iVoid].macrocenter[0] += voidPart[p].vol*(dist[0]);
voids[iVoid].macrocenter[1] += voidPart[p].vol*(dist[1]);
voids[iVoid].macrocenter[2] += voidPart[p].vol*(dist[2]);
weight += voidPart[p].vol;
}
voids[iVoid].macrocenter[0] /= weight;
voids[iVoid].macrocenter[1] /= weight;
voids[iVoid].macrocenter[2] /= weight;
voids[iVoid].macrocenter[0] += voids[iVoid].center[0];
voids[iVoid].macrocenter[1] += voids[iVoid].center[1];
voids[iVoid].macrocenter[2] += voids[iVoid].center[2];
if (periodicX) {
if (voids[iVoid].macrocenter[0] > ranges[0][1])
voids[iVoid].macrocenter[0] = voids[iVoid].macrocenter[0] - boxLen[0];
if (voids[iVoid].macrocenter[0] < ranges[0][0])
voids[iVoid].macrocenter[0] = boxLen[0] + voids[iVoid].macrocenter[0];
}
if (periodicY) {
if (voids[iVoid].macrocenter[1] > ranges[1][1])
voids[iVoid].macrocenter[1] = voids[iVoid].macrocenter[1] - boxLen[1];
if (voids[iVoid].macrocenter[1] < ranges[1][0])
voids[iVoid].macrocenter[1] = boxLen[1] + voids[iVoid].macrocenter[1];
}
if (periodicZ) {
if (voids[iVoid].macrocenter[2] > ranges[2][1])
voids[iVoid].macrocenter[2] = voids[iVoid].macrocenter[2] - boxLen[2];
if (voids[iVoid].macrocenter[2] < ranges[2][0])
voids[iVoid].macrocenter[2] = boxLen[2] + voids[iVoid].macrocenter[2];
}
clock4 = clock();
interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
//printf(" %.2f for macrocenter\n", interval);
// compute central density
clock3 = clock();
centralRad = voids[iVoid].radius/args.centralRadFrac_arg;
centralDen = 0.;
int numCentral = 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 (sqrt(dist2) < centralRad) numCentral += 1;
}
voids[iVoid].centralDen = numCentral / (volNorm*4./3. * M_PI *
pow(centralRad, 3.));
clock4 = clock();
interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
//printf(" %.2f for central density\n", interval);
// 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]);
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);
// compute distance from macrocenter to nearest flagged particle
clock3 = clock();
if (args.isObservation_flag) {
minDist = 1.e99;
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].nearestFlag = sqrt(minDist);
} else {
voids[iVoid].nearestFlag = 1.e99;
}
if (args.isObservation_flag) {
voids[iVoid].redshiftInMpc =
sqrt(pow(voids[iVoid].macrocenter[0] - boxLen[0]/2.,2) +
pow(voids[iVoid].macrocenter[1] - boxLen[1]/2.,2) +
pow(voids[iVoid].macrocenter[2] - boxLen[2]/2.,2));
if (args.useComoving_flag) {
redshift = gsl_interp_eval(interp, dL, redshifts,
voids[iVoid].redshiftInMpc, acc);
voids[iVoid].redshift = redshift;
} else {
redshift = voids[iVoid].redshiftInMpc;
voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.;
}
} else {
voids[iVoid].redshiftInMpc = voids[iVoid].macrocenter[2];
if (args.useComoving_flag) {
voids[iVoid].redshift = gsl_interp_eval(interp, dL, redshifts,
voids[iVoid].redshiftInMpc, acc);
} else {
voids[iVoid].redshift = voids[iVoid].macrocenter[2]/LIGHT_SPEED*100.;
}
nearestEdge = 1.e99;
if (!periodicX) {
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].macrocenter[0] - ranges[0][0]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].macrocenter[0] - ranges[0][1]));
}
if (!periodicY) {
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].macrocenter[1] - ranges[1][0]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].macrocenter[1] - ranges[1][1]));
}
if (!periodicZ) {
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].macrocenter[2] - ranges[2][0]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].macrocenter[2] - ranges[2][1]));
}
}
voids[iVoid].nearestEdge = nearestEdge;
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();
double inertia[9];
for (int i = 0; i < 9; i++) inertia[i] = 0.;
for (int p = 0; p < voids[iVoid].numPart; p++) {
dist[0] = voidPart[p].x - voids[iVoid].macrocenter[0];
dist[1] = voidPart[p].y - voids[iVoid].macrocenter[1];
dist[2] = voidPart[p].z - voids[iVoid].macrocenter[2];
if (periodicX && fabs(dist[0]) > boxLen[0]/2.)
dist[0] = dist[0] - copysign(boxLen[0], dist[0]);
if (periodicY && fabs(dist[1]) > boxLen[1]/2.)
dist[1] = dist[1] - copysign(boxLen[1], dist[1]);
if (periodicZ && fabs(dist[2]) > boxLen[2]/2.)
dist[2] = dist[2] - copysign(boxLen[2], dist[2]);
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];
inertia[i*3+j] -= dist[i]*dist[j];
}
}
}
inertia[1*3+0] = inertia[0*3+1];
inertia[2*3+0] = inertia[0*3+2];
inertia[2*3+1] = inertia[1*3+2];
gsl_matrix_view m = gsl_matrix_view_array(inertia, 3, 3);
gsl_eigen_symmv(&m.matrix, voids[iVoid].eval, voids[iVoid].evec, eigw);
float a = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,1) +
gsl_vector_get(voids[iVoid].eval,2) -
gsl_vector_get(voids[iVoid].eval,0)));
float b = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,2) +
gsl_vector_get(voids[iVoid].eval,0) -
gsl_vector_get(voids[iVoid].eval,1)));
float c = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,0) +
gsl_vector_get(voids[iVoid].eval,1) -
gsl_vector_get(voids[iVoid].eval,2)));
float ca;
float cb = c/b;
float smallest = 1.e99;
float largest = 0.;
for (int i = 0; i < 3; i++) {
if (gsl_vector_get(voids[iVoid].eval,i) < smallest)
smallest = gsl_vector_get(voids[iVoid].eval,i);
if (gsl_vector_get(voids[iVoid].eval,i) > largest)
largest = gsl_vector_get(voids[iVoid].eval,i);
}
voids[iVoid].ellip = 1.0 - sqrt(sqrt(fabs(smallest/largest)));
clock4 = clock();
interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
} // 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;
int numEdge = 0;
int numNearZ = 0;
int numAreParents = 0;
int numTooSmall = 0;
printf(" Picking winners and losers...\n");
printf(" Starting with %d voids\n", (int) voids.size());
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
voids[iVoid].accepted = 1;
}
// toss out voids that are obviously wrong
int iGood = 0;
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (isnan(voids[iVoid].vol) || isinf(voids[iVoid].vol)) {
numWrong++;
} else {
voids[iGood++] = voids[iVoid];
}
}
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) {
numTooSmall++;
} else {
voids[iGood++] = voids[iVoid];
}
}
voids.resize(iGood);
printf(" 2nd filter: rejected %d too small\n", numTooSmall);
// mark top-level voids
numAreParents = 0;
iGood = 0;
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].parentIndx != -1) {
numAreParents++;
voids[iVoid].isLeaf = true;
} else {
voids[iVoid].isLeaf = false;
}
}
// mark high-density voids
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].centralDen > args.maxCentralDen_arg) {
voids[iVoid].accepted = -1;
voids[iVoid].hasHighCentralDen = true;
numHighDen++;
} else {
voids[iVoid].hasHighCentralDen = false;
}
}
// count voids near survey edges
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].voidType == CENTRAL_VOID) {
numCentral++;
} else {
numEdge++;
}
}
printf(" Number kept: %d (out of %d)\n", (int) voids.size(), numVoids);
printf(" We have %d edge voids\n", numEdge);
printf(" We have %d central voids\n", numCentral);
printf(" We have %d too high central density\n", numHighDen);
printf(" We have %d that are not leaf nodes\n", numAreParents);
outputDir = string(args.outputDir_arg);
sampleName = (string(args.sampleName_arg)+".out");
dataPortions[0] = "central";
dataPortions[1] = "all";
printf(" Output fully trimmed catalog...\n");
prefix = "";
for (int i = 0; i < 2; i++) {
dataPortion = dataPortions[i];
outputVoids(outputDir, sampleName, prefix, dataPortion,
numPartTot,
voids,
args.isObservation_flag, boxLen, true, true);
}
printf(" Output fully untrimmed catalog...\n");
prefix = "untrimmed_";
for (int i = 0; i < 2; i++) {
dataPortion = dataPortions[i];
outputVoids(outputDir, sampleName, prefix, dataPortion,
numPartTot,
voids,
args.isObservation_flag, boxLen, false, false);
}
clock2 = clock();
printf(" Time: %f sec (for %d voids)\n",
(1.*clock2-clock1)/CLOCKS_PER_SEC, numVoids);
clock1 = clock();
printf("Done!\n");
} // end main
// ----------------------------------------------------------------------------
void openFiles(string outputDir, string sampleName,
string prefix, string dataPortion,
int numPartTot, int numKept,
FILE** fpZobov, FILE** fpCenters,
FILE** fpBarycenter, FILE** fpShapes,
FILE** fpExtra,
FILE** fpSkyPositions) {
*fpZobov = fopen((outputDir+"/"+prefix+"voidDesc_"+dataPortion+"_"+sampleName).c_str(), "w");
fprintf(*fpZobov, "%d particles, %d voids.\n", numPartTot, numKept);
fprintf(*fpZobov, "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
*fpBarycenter = fopen((outputDir+"/"+prefix+"macrocenters_"+dataPortion+"_"+sampleName).c_str(), "w");
*fpCenters = fopen((outputDir+"/"+prefix+"centers_"+dataPortion+"_"+sampleName).c_str(), "w");
fprintf(*fpCenters, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part, parent index, tree level, number of children, central density\n");
*fpSkyPositions = fopen((outputDir+"/"+prefix+"sky_positions_"+dataPortion+"_"+sampleName).c_str(), "w");
fprintf(*fpSkyPositions, "# RA, dec, redshift, radius (Mpc/h), void ID\n");
*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
// ----------------------------------------------------------------------------
void closeFiles(FILE* fpZobov, FILE* fpCenters,
FILE* fpBarycenter, FILE* fpShapes,
FILE* fpExtra,
FILE* fpSkyPositions) {
fclose(fpZobov);
fclose(fpCenters);
fclose(fpBarycenter);
//fclose(fpDistances);
fclose(fpShapes);
fclose(fpExtra);
fclose(fpSkyPositions);
} // end closeFile
// ----------------------------------------------------------------------------
void outputVoids(string outputDir, string sampleName, string prefix,
string dataPortion, int numPartTot,
vector<VOID> voids,
bool isObservation, double *boxLen, bool doTrim,
bool doCentralDenCut) {
int iVoid;
VOID outVoid;
FILE *fp, *fpZobov, *fpCenters, *fpCentersNoCut, *fpBarycenter,
*fpShapes, *fpExtra, *fpSkyPositions;
openFiles(outputDir, sampleName, prefix, dataPortion,
numPartTot, voids.size(),
&fpZobov, &fpCenters, &fpBarycenter,
&fpShapes, &fpExtra, &fpSkyPositions);
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
outVoid = voids[iVoid];
if (dataPortion == "central" && outVoid.voidType == EDGE_VOID) {
continue;
}
if (doTrim && outVoid.isLeaf) {
continue;
}
if (doCentralDenCut && outVoid.hasHighCentralDen) {
continue;
}
double outCenter[3];
outCenter[0] = outVoid.macrocenter[0];
outCenter[1] = outVoid.macrocenter[1];
outCenter[2] = outVoid.macrocenter[2];
//if (isObservation) {
// outCenter[0] = (outVoid.macrocenter[0]-boxLen[0]/2.)*100.;
// outCenter[1] = (outVoid.macrocenter[1]-boxLen[1]/2.)*100.;
// outCenter[2] = (outVoid.macrocenter[2]-boxLen[2]/2.)*100.;
//}
fprintf(fpZobov, "%d %d %d %f %f %d %d %f %d %f %f\n",
iVoid,
outVoid.voidID,
outVoid.coreParticle,
outVoid.coreDens,
outVoid.zoneVol,
outVoid.zoneNumPart,
outVoid.numZones,
outVoid.vol,
outVoid.numPart,
outVoid.densCon,
outVoid.voidProb);
fprintf(fpBarycenter, "%d %e %e %e\n",
outVoid.voidID,
outVoid.macrocenter[0],
outVoid.macrocenter[1],
outVoid.macrocenter[2]);
fprintf(fpCenters, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d %d %d %d %.2f\n",
outCenter[0],
outCenter[1],
outCenter[2],
outVoid.vol,
outVoid.radius,
outVoid.redshift,
4./3.*M_PI*pow(outVoid.radius, 3),
outVoid.voidID,
outVoid.densCon,
outVoid.numPart,
outVoid.parentIndx,
outVoid.level,
outVoid.numChildren,
outVoid.centralDen);
double phi = atan2(outVoid.macrocenter[1]-boxLen[1]/2.,
outVoid.macrocenter[0]-boxLen[0]/2.);
if (phi < 0) phi += 2.*M_PI;
double RA = phi * 180./M_PI;
double theta = acos((outVoid.macrocenter[2]-boxLen[2]/2.) /
outVoid.redshiftInMpc);
double dec = (M_PI/2. - theta) * 180./M_PI;
fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n",
RA,
dec,
outVoid.redshift,
outVoid.radius,
outVoid.voidID);
fprintf(fpShapes, "%d %.6f %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e\n",
outVoid.voidID,
outVoid.ellip,
gsl_vector_get(outVoid.eval, 0),
gsl_vector_get(outVoid.eval, 1),
gsl_vector_get(outVoid.eval, 2),
gsl_matrix_get(outVoid.evec, 0 ,0),
gsl_matrix_get(outVoid.evec, 1 ,0),
gsl_matrix_get(outVoid.evec, 2 ,0),
gsl_matrix_get(outVoid.evec, 0 ,1),
gsl_matrix_get(outVoid.evec, 1 ,1),
gsl_matrix_get(outVoid.evec, 2 ,1),
gsl_matrix_get(outVoid.evec, 0 ,2),
gsl_matrix_get(outVoid.evec, 1 ,2),
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, fpExtra, fpSkyPositions);
} // end outputVoids
Reference in a new issue View git blame Copy permalink