guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 23:31:12 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions
vide_public/c_source/pruning/pruneVoids.cpp

1039 lines
33 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 <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 parentID, 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,
int numPartTot, int numKept,
FILE** fpOutput);
void closeFiles(FILE* fpOutput);
void outputVoids(string outputDir, string sampleName, int numPartTot,
vector<VOID> voids,
bool isObservation, double *boxLen);
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 volNormZobov, volNormObs, 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));
volNormZobov = args.volNormZobov_arg;
volNormObs = args.volNormObs_arg;
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/volNormZobov*3./4./M_PI, 1./3.);
voids[i-1].accepted = 1;
voids[i-1].isLeaf = true;
voids[i-1].hasHighCentralDen = false;
voids[i-1].parentID = -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].parentID = checkID;
voids[iCheck].numChildren++;
}
}
} // end building tree
// find level in tree
for (iVoid = 0; iVoid < numVoids; iVoid++) {
int level = 0;
int parentIndx = iVoid;
int parentID = voids[iVoid].parentID;
while (parentID != -1) {
level++;
parentID = voids[parentIndx].parentID;
// get the index of parent void
for (int iCheck = 0; iCheck < numVoids; iCheck++){
if (voids[iCheck].voidID == parentID) {
parentIndx = iCheck;
break;
}
}
}
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 / (volNormObs*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].parentID != -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 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 untrimmed catalog...\n");
outputVoids(outputDir, sampleName,
numPartTot,
voids,
args.isObservation_flag, boxLen);
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,
int numPartTot, int numKept,
FILE** fpOutput) {
*fpOutput = fopen((outputDir+"/voidDatabase_"+sampleName).c_str(), "w");
//fprintf(*fpZobov, "%d particles, %d voids.\n", numPartTot, numKept);
fprintf(*fpOutput, "Void ID, void type, center (x,y,z) (Mpc/h), volume (normalized), volume(Mpc/h^3), radius (Mpc/h), redshift, RA, Dec, density contrast, max extent, nearest edge, num part, parent ID, tree level, num children, central density, core particle, core density, zone vol, zone num part, num zones, void probability, ellipticity, 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");
} // end openFiles
// ----------------------------------------------------------------------------
void closeFiles(FILE* fpOutput) {
fclose(fpOutput);
//fclose(fpCenters);
//fclose(fpBarycenter);
//fclose(fpDistances);
//fclose(fpShapes);
//fclose(fpExtra);
//fclose(fpSkyPositions);
} // end closeFile
// ----------------------------------------------------------------------------
void outputVoids(string outputDir, string sampleName, int numPartTot,
vector<VOID> voids,
bool isObservation, double *boxLen) {
int iVoid;
VOID outVoid;
FILE *fp, *fpOutput;
openFiles(outputDir, sampleName,
numPartTot, voids.size(),
&fpOutput);
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
outVoid = voids[iVoid];
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(fpOutput, "%d %d %e %e %e %e %e %e %e %e %e %e %e %e %d %d %d %d %e %d %e %e %d %d %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n",
outVoid.voidID,
outVoid.voidType,
outVoid.macrocenter[0],
outVoid.macrocenter[1],
outVoid.macrocenter[2],
outVoid.vol,
4./3.*M_PI*pow(outVoid.radius, 3),
outVoid.radius,
outVoid.redshift,
RA,
dec,
outVoid.densCon,
outVoid.maxRadius,
outVoid.nearestFlag,
outVoid.numPart,
outVoid.parentID,
outVoid.level,
outVoid.numChildren,
outVoid.centralDen,
outVoid.coreParticle,
outVoid.coreDens,
outVoid.zoneVol,
outVoid.zoneNumPart,
outVoid.numZones,
outVoid.voidProb,
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)
);
} // end iVoid
closeFiles(fpOutput);
} // end outputVoids
Reference in a new issue View git blame Copy permalink