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

cleaned up logic and output handling in prunevoids; added combined SDSS analysis script; added data preparation scripts

Browse source
This commit is contained in:
P.M. Sutter 2013-02-25 12:30:24 -06:00
parent 7f020f15e5
commit fa7264e4ee
3 changed files with 293 additions and 239 deletions
Download patch file Download diff file Expand all files Collapse all files

478
c_tools/stacking/pruneVoids.cpp
View file

@ -18,11 +18,15 @@
#include <stdio.h> #include <stdio.h>
#include <netcdfcpp.h> #include <netcdfcpp.h>
#include "pruneVoids_conf.h" #include "pruneVoids_conf.h"
#include <vector>
#define LIGHT_SPEED 299792.458 #define LIGHT_SPEED 299792.458
#define MPC2Z 100./LIGHT_SPEED #define MPC2Z 100./LIGHT_SPEED
#define Z2MPC LIGHT_SPEED/100. #define Z2MPC LIGHT_SPEED/100.
#define CENTRAL_VOID 1
#define EDGE_VOID 2
typedef struct partStruct { typedef struct partStruct {
float x, y, z, vol; float x, y, z, vol;
} PART; } PART;
@ -44,32 +48,38 @@ typedef struct voidStruct {
float nearestEdge; float nearestEdge;
float center[3], barycenter[3]; float center[3], barycenter[3];
int accepted; int accepted;
int voidType;
gsl_vector *eval; gsl_vector *eval;
gsl_matrix *evec; gsl_matrix *evec;
} VOID; } VOID;
void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters,
FILE* fpCenterNoCut,
FILE* fpSkyPositions, FILE* fpBarycenters, FILE* fpDistances,
FILE* fpShapes, bool isObservation, double *boxLen);
int main(int argc, char **argv) { int main(int argc, char **argv) {
// initialize arguments // initialize arguments
pruneVoids_info args_info; pruneVoids_info args;
pruneVoids_conf_params args_params; pruneVoids_conf_params args_params;
pruneVoids_conf_init(&args_info); pruneVoids_conf_init(&args);
pruneVoids_conf_params_init(&args_params); pruneVoids_conf_params_init(&args_params);
args_params.check_required = 0; args_params.check_required = 0;
if (pruneVoids_conf_ext (argc, argv, &args_info, &args_params)) if (pruneVoids_conf_ext (argc, argv, &args, &args_params))
return 1; return 1;
if (!args_info.configFile_given) { if (!args.configFile_given) {
if (pruneVoids_conf_required (&args_info, if (pruneVoids_conf_required (&args,
PRUNEVOIDS_CONF_PACKAGE)) PRUNEVOIDS_CONF_PACKAGE))
return 1; return 1;
} else { } else {
args_params.check_required = 1; args_params.check_required = 1;
args_params.initialize = 0; args_params.initialize = 0;
if (pruneVoids_conf_config_file (args_info.configFile_arg, if (pruneVoids_conf_config_file (args.configFile_arg,
&args_info, &args,
&args_params)) &args_params))
return 1; return 1;
} }
@ -77,12 +87,15 @@ int main(int argc, char **argv) {
int i, p, p2, numPartTot, numZonesTot, dummy, iVoid, iZ; int i, p, p2, numPartTot, numZonesTot, dummy, iVoid, iZ;
int numVoids, mockIndex, numKept; int numVoids, mockIndex, numKept;
double tolerance; double tolerance;
FILE *fp, *fpBarycenter, *fpDistances, *fpSkyPositions, *fpInfo; FILE *fp, *fpZobovCentral, *fpZobovAll, *fpCentersCentral, *fpCentersAll,
FILE *fpShapes; *fpCentersNoCutCentral, *fpCentersNoCutAll, *fpBarycenterCentral,
*fpBarycenterAll, *fpDistancesCentral, *fpDistancesAll,
*fpShapesCentral, *fpShapesAll, *fpSkyPositionsCentral,
*fpSkyPositionsAll;
PART *part, *voidPart; PART *part, *voidPart;
ZONE2PART *zones2Parts; ZONE2PART *zones2Parts;
VOID2ZONE *void2Zones; VOID2ZONE *void2Zones;
VOID *voids; std::vector<VOID> voids;
float *temp, junk, voidVol; float *temp, junk, voidVol;
int junkInt, voidID, numPart, numZones, zoneID, partID, maxNumPart; int junkInt, voidID, numPart, numZones, zoneID, partID, maxNumPart;
int coreParticle, zoneNumPart; int coreParticle, zoneNumPart;
@ -98,30 +111,30 @@ int main(int argc, char **argv) {
gsl_eigen_symmv_workspace *eigw = gsl_eigen_symmv_alloc(3); gsl_eigen_symmv_workspace *eigw = gsl_eigen_symmv_alloc(3);
numVoids = args_info.numVoids_arg; numVoids = args.numVoids_arg;
mockIndex = args_info.mockIndex_arg; mockIndex = args.mockIndex_arg;
tolerance = args_info.tolerance_arg; tolerance = args.tolerance_arg;
clock1 = clock(); clock1 = clock();
printf("Pruning parameters: %f %f %f %s\n", args_info.zMin_arg, printf("Pruning parameters: %f %f %f %s\n", args.zMin_arg,
args_info.zMax_arg, args.zMax_arg,
args_info.rMin_arg, args.rMin_arg,
args_info.periodic_arg); args.periodic_arg);
// check for periodic box // check for periodic box
periodicX = 0; periodicX = 0;
periodicY = 0; periodicY = 0;
periodicZ = 0; periodicZ = 0;
if (!args_info.isObservation_flag) { if (!args.isObservation_flag) {
if ( strchr(args_info.periodic_arg, 'x') != NULL) { if ( strchr(args.periodic_arg, 'x') != NULL) {
periodicX = 1; periodicX = 1;
printf("Will assume x-direction is periodic.\n"); printf("Will assume x-direction is periodic.\n");
} }
if ( strchr(args_info.periodic_arg, 'y') != NULL) { if ( strchr(args.periodic_arg, 'y') != NULL) {
periodicY = 1; periodicY = 1;
printf("Will assume y-direction is periodic.\n"); printf("Will assume y-direction is periodic.\n");
} }
if ( strchr(args_info.periodic_arg, 'z') != NULL) { if ( strchr(args.periodic_arg, 'z') != NULL) {
periodicZ = 1; periodicZ = 1;
printf("Will assume z-direction is periodic.\n"); printf("Will assume z-direction is periodic.\n");
} }
@ -129,7 +142,7 @@ int main(int argc, char **argv) {
// load box size // load box size
printf("\n Getting info...\n"); printf("\n Getting info...\n");
NcFile f_info(args_info.extraInfo_arg); NcFile f_info(args.extraInfo_arg);
ranges[0][0] = f_info.get_att("range_x_min")->as_double(0); ranges[0][0] = f_info.get_att("range_x_min")->as_double(0);
ranges[0][1] = f_info.get_att("range_x_max")->as_double(0); ranges[0][1] = f_info.get_att("range_x_max")->as_double(0);
ranges[1][0] = f_info.get_att("range_y_min")->as_double(0); ranges[1][0] = f_info.get_att("range_y_min")->as_double(0);
@ -143,7 +156,7 @@ int main(int argc, char **argv) {
// read in all particle positions // read in all particle positions
printf("\n Loading particles...\n"); printf("\n Loading particles...\n");
fp = fopen(args_info.partFile_arg, "r"); fp = fopen(args.partFile_arg, "r");
fread(&dummy, 1, 4, fp); fread(&dummy, 1, 4, fp);
fread(&numPartTot, 1, 4, fp); fread(&numPartTot, 1, 4, fp);
fread(&dummy, 1, 4, fp); fread(&dummy, 1, 4, fp);
@ -154,7 +167,7 @@ int main(int argc, char **argv) {
volNorm = numPartTot/(boxLen[0]*boxLen[1]*boxLen[2]); volNorm = numPartTot/(boxLen[0]*boxLen[1]*boxLen[2]);
printf(" VOL NORM = %f\n", volNorm); printf(" VOL NORM = %f\n", volNorm);
printf(" CENTRAL DEN = %f\n", args_info.maxCentralDen_arg); printf(" CENTRAL DEN = %f\n", args.maxCentralDen_arg);
fread(&dummy, 1, 4, fp); fread(&dummy, 1, 4, fp);
fread(temp, numPartTot, 4, fp); fread(temp, numPartTot, 4, fp);
@ -174,7 +187,7 @@ int main(int argc, char **argv) {
for (p = 0; p < numPartTot; p++) for (p = 0; p < numPartTot; p++)
part[p].z = mul*temp[p]; part[p].z = mul*temp[p];
if (!args_info.isObservation_flag) { if (!args.isObservation_flag) {
for (p = 0; p < numPartTot; p++) { for (p = 0; p < numPartTot; p++) {
part[p].x += ranges[0][0]; part[p].x += ranges[0][0];
part[p].y += ranges[1][0]; part[p].y += ranges[1][0];
@ -189,12 +202,12 @@ int main(int argc, char **argv) {
// read in desired voids // read in desired voids
printf(" Loading voids...\n"); printf(" Loading voids...\n");
fp = fopen(args_info.voidDesc_arg ,"r"); fp = fopen(args.voidDesc_arg ,"r");
fgets(line, sizeof(line), fp); fgets(line, sizeof(line), fp);
sscanf(line, "%d %s %d %s", &junkInt, junkStr, &junkInt, junkStr); sscanf(line, "%d %s %d %s", &junkInt, junkStr, &junkInt, junkStr);
fgets(line, sizeof(line), fp); fgets(line, sizeof(line), fp);
voids = (VOID *) malloc(numVoids * sizeof(VOID)); voids.resize(numVoids);
i = 0; i = 0;
while (fgets(line, sizeof(line), fp) != NULL) { while (fgets(line, sizeof(line), fp) != NULL) {
sscanf(line, "%d %d %d %f %f %d %d %f %d %f %f\n", &iVoid, &voidID, sscanf(line, "%d %d %d %f %f %d %d %f %d %f %f\n", &iVoid, &voidID,
@ -223,7 +236,7 @@ int main(int argc, char **argv) {
// load up the zone membership for each void // load up the zone membership for each void
printf(" Loading void-zone membership info...\n"); printf(" Loading void-zone membership info...\n");
fp = fopen(args_info.void2Zone_arg, "r"); fp = fopen(args.void2Zone_arg, "r");
fread(&numZonesTot, 1, 4, fp); fread(&numZonesTot, 1, 4, fp);
void2Zones = (VOID2ZONE *) malloc(numZonesTot * sizeof(VOID2ZONE)); void2Zones = (VOID2ZONE *) malloc(numZonesTot * sizeof(VOID2ZONE));
@ -241,7 +254,7 @@ int main(int argc, char **argv) {
// now the particles-zone // now the particles-zone
printf(" Loading particle-zone membership info...\n"); printf(" Loading particle-zone membership info...\n");
fp = fopen(args_info.zone2Part_arg, "r"); fp = fopen(args.zone2Part_arg, "r");
fread(&dummy, 1, 4, fp); fread(&dummy, 1, 4, fp);
fread(&numZonesTot, 1, 4, fp); fread(&numZonesTot, 1, 4, fp);
@ -259,7 +272,7 @@ int main(int argc, char **argv) {
// and finally volumes // and finally volumes
printf(" Loading particle volumes...\n"); printf(" Loading particle volumes...\n");
fp = fopen(args_info.partVol_arg, "r"); fp = fopen(args.partVol_arg, "r");
fread(&mask_index, 1, 4, fp); fread(&mask_index, 1, 4, fp);
if (mask_index != mockIndex) { if (mask_index != mockIndex) {
printf("NON-MATCHING MOCK INDICES!? %d %d\n", mask_index, mockIndex); printf("NON-MATCHING MOCK INDICES!? %d %d\n", mask_index, mockIndex);
@ -366,7 +379,7 @@ int main(int argc, char **argv) {
} }
// compute central density // compute central density
centralRad = voids[iVoid].radius/args_info.centralRadFrac_arg; centralRad = voids[iVoid].radius/args.centralRadFrac_arg;
centralDen = 0.; centralDen = 0.;
int numCentral = 0; int numCentral = 0;
for (p = 0; p < voids[iVoid].numPart; p++) { for (p = 0; p < voids[iVoid].numPart; p++) {
@ -386,7 +399,7 @@ int main(int argc, char **argv) {
// compute maximum extent // compute maximum extent
/* /*
if (args_info.isObservation_flag) { if (args.isObservation_flag) {
maxDist = 0.; maxDist = 0.;
for (p = 0; p < voids[iVoid].numPart; p++) { for (p = 0; p < voids[iVoid].numPart; p++) {
for (p2 = p; p2 < voids[iVoid].numPart; p2++) { for (p2 = p; p2 < voids[iVoid].numPart; p2++) {
@ -419,7 +432,7 @@ int main(int argc, char **argv) {
voids[iVoid].maxRadius = sqrt(maxDist); voids[iVoid].maxRadius = sqrt(maxDist);
// } // }
if (args_info.isObservation_flag) { if (args.isObservation_flag) {
// compute distance from center to nearest mock // compute distance from center to nearest mock
minDist = 1.e99; minDist = 1.e99;
for (p = mockIndex; p < numPartTot; p++) { for (p = mockIndex; p < numPartTot; p++) {
@ -436,16 +449,16 @@ int main(int argc, char **argv) {
voids[iVoid].nearestMock = 1.e99; voids[iVoid].nearestMock = 1.e99;
} }
if (args_info.isObservation_flag) { if (args.isObservation_flag) {
voids[iVoid].redshiftInMpc = voids[iVoid].redshiftInMpc =
sqrt(pow(voids[iVoid].barycenter[0] - boxLen[0]/2.,2) + sqrt(pow(voids[iVoid].barycenter[0] - boxLen[0]/2.,2) +
pow(voids[iVoid].barycenter[1] - boxLen[1]/2.,2) + pow(voids[iVoid].barycenter[1] - boxLen[1]/2.,2) +
pow(voids[iVoid].barycenter[2] - boxLen[2]/2.,2)); pow(voids[iVoid].barycenter[2] - boxLen[2]/2.,2));
voids[iVoid].redshiftInMpc = voids[iVoid].redshiftInMpc; voids[iVoid].redshiftInMpc = voids[iVoid].redshiftInMpc;
redshift = voids[iVoid].redshiftInMpc; redshift = voids[iVoid].redshiftInMpc;
nearestEdge = fabs(redshift-args_info.zMax_arg*LIGHT_SPEED/100.); nearestEdge = fabs(redshift-args.zMax_arg*LIGHT_SPEED/100.);
//nearestEdge = fmin(fabs(redshift-args_info.zMin_arg*LIGHT_SPEED/100.), //nearestEdge = fmin(fabs(redshift-args.zMin_arg*LIGHT_SPEED/100.),
// fabs(redshift-args_info.zMax_arg*LIGHT_SPEED/100.)); // fabs(redshift-args.zMax_arg*LIGHT_SPEED/100.));
voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.; voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.;
} else { } else {
@ -504,200 +517,271 @@ int main(int argc, char **argv) {
int numWrong = 0; int numWrong = 0;
int numHighDen = 0; int numHighDen = 0;
int numCentral = 0;
int numEdge = 0; int numEdge = 0;
int numNearZ = 0;
int numTooSmall = 0; int numTooSmall = 0;
printf(" Picking winners and losers...\n"); printf(" Picking winners and losers...\n");
for (iVoid = 0; iVoid < numVoids; iVoid++) { printf(" Starting with %d voids\n", voids.size());
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
voids[iVoid].accepted = 1; voids[iVoid].accepted = 1;
} }
for (iVoid = 0; iVoid < numVoids; iVoid++) { /*
int j = 0;
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].densCon < 1.5) { if (voids[iVoid].densCon < 1.5) {
// voids[iVoid].accepted = -4; // voids[iVoid].accepted = -4;
} }
}
*/
if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) { // toss out voids that are obviously wrong
int iGood = 0;
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].densCon > 1.e4) {
numWrong++;
} else {
voids[iGood++] = voids[iVoid];
}
}
voids.resize(iGood);
printf(" 1st filter: reiGoodected %d obviously bad\n", numWrong);
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: reiGoodected %d too small\n", numTooSmall);
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) {
numNearZ++;
} else {
voids[iGood++] = voids[iVoid];
}
}
voids.resize(iGood);
printf(" 3rd filter: reiGoodected %d too close to high redshift boundaries\n", numNearZ);
numNearZ = 0;
iGood = 0;
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
// assume the lower z-boundary is "soft" in observations
if (args.isObservation_flag &&
voids[iVoid].redshift < args.zMin_arg) {
numNearZ++;
} else {
voids[iGood++] = voids[iVoid];
}
}
voids.resize(iGood);
printf(" 4th filter: reiGoodected %d too close to low redshift boundaries\n", numNearZ);
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].centralDen > args.maxCentralDen_arg) {
voids[iVoid].accepted = -1; voids[iVoid].accepted = -1;
numHighDen++; numHighDen++;
} }
}
// toss out voids that are obviously wrong for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].densCon > 1.e4) { if (tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) {
voids[iVoid].accepted = -4; voids[iVoid].voidType = CENTRAL_VOID;
numWrong++; numCentral++;
} } else {
voids[iVoid].voidType = EDGE_VOID;
if (strcmp(args_info.dataPortion_arg, "edge") == 0 &&
tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) {
voids[iVoid].accepted = -3;
numEdge++; numEdge++;
} }
if (strcmp(args_info.dataPortion_arg, "central") == 0 &&
tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestMock) {
voids[iVoid].accepted = -3;
numEdge++;
}
if (voids[iVoid].radius < args_info.rMin_arg) {
voids[iVoid].accepted = -2;
numTooSmall++;
}
// *always* clean out near edges since there are no mocks there
if (tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestEdge) {
voids[iVoid].accepted = -3;
if (voids[iVoid].accepted == 1) numEdge++;
}
// assume the lower z-boundary is "soft" in observations
if (args_info.isObservation_flag &&
voids[iVoid].redshift < args_info.zMin_arg) {
voids[iVoid].accepted = -3;
if (voids[iVoid].accepted == 1) numEdge++;
}
} }
numKept = 0; printf(" Number kept: %d (out of %d)\n", voids.size(), numVoids);
for (iVoid = 0; iVoid < numVoids; iVoid++) { printf(" We have %d edge voids\n", numEdge);
if (voids[iVoid].accepted == 1) numKept++; printf(" We have %d central voids\n", numCentral);
} printf(" We have %d too high central density\n", numHighDen);
printf(" Number kept: %d (out of %d)\n", numKept, numVoids);
printf(" Rejected %d near the edge\n", numEdge);
printf(" Rejected %d too small\n", numTooSmall);
printf(" Rejected %d obviously bad\n", numWrong);
printf(" Rejected %d too high central density\n", numHighDen);
printf(" Output...\n"); printf(" Output...\n");
fp = fopen(args_info.output_arg, "w"); fpZobovCentral = fopen((std::string(args.outputDir_arg)+"/voidDesc_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fpBarycenter = fopen(args_info.outCenters_arg, "w"); fprintf(fpZobovCentral, "%d particles, %d voids.\n", mockIndex, numKept);
fpInfo = fopen(args_info.outInfo_arg, "w"); fprintf(fpZobovCentral, "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
fpDistances = fopen(args_info.outDistances_arg, "w");
fpSkyPositions = fopen(args_info.outSkyPositions_arg, "w");
fpShapes = fopen(args_info.outShapes_arg, "w");
fprintf(fp, "%d particles, %d voids.\n", mockIndex, numKept);
fprintf(fp, "see column in master void file\n");
fprintf(fpInfo, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast\n");
fprintf(fpSkyPositions, "# RA, dec, redshift, radius (Mpc/h), void ID\n");
fprintf(fpShapes, "# void ID, 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");
for (iVoid = 0; iVoid < numVoids; iVoid++) {
if (voids[iVoid].accepted != 1) continue; fpZobovAll = fopen((std::string(args.outputDir_arg)+"/voidDesc_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpZobovAll, "%d particles, %d voids.\n", mockIndex, numKept);
fprintf(fpZobovAll, "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
fprintf(fp, "%d %d %d %f %f %d %d %f %d %f %f\n", fpBarycenterCentral = fopen((std::string(args.outputDir_arg)+"/barycenters_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
iVoid, fpBarycenterAll = fopen((std::string(args.outputDir_arg)+"/barycenters_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
voids[iVoid].voidID,
voids[iVoid].coreParticle,
voids[iVoid].coreDens,
voids[iVoid].zoneVol,
voids[iVoid].zoneNumPart,
voids[iVoid].numZones,
voids[iVoid].vol,
voids[iVoid].numPart,
voids[iVoid].densCon,
voids[iVoid].voidProb);
fprintf(fpBarycenter, "%d %e %e %e\n", fpCentersCentral = fopen((std::string(args.outputDir_arg)+"/centers_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
voids[iVoid].voidID, fprintf(fpCentersCentral, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast\n");
voids[iVoid].barycenter[0],
voids[iVoid].barycenter[1],
voids[iVoid].barycenter[2]);
fprintf(fpDistances, "%d %e\n", fpCentersAll = fopen((std::string(args.outputDir_arg)+"/centers_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
voids[iVoid].voidID, fprintf(fpCentersAll, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast\n");
voids[iVoid].nearestMock);
double outCenter[3]; fpCentersNoCutCentral = fopen((std::string(args.outputDir_arg)+"/centers_nocut_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
outCenter[0] = voids[iVoid].barycenter[0]; fprintf(fpCentersNoCutCentral, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast\n");
outCenter[1] = voids[iVoid].barycenter[1];
outCenter[2] = voids[iVoid].barycenter[2];
if (args_info.isObservation_flag) { fpCentersNoCutAll = fopen((std::string(args.outputDir_arg)+"/centers_nocut_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
outCenter[0] = (voids[iVoid].barycenter[0]-boxLen[0]/2.)*100.; fprintf(fpCentersNoCutAll, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast\n");
outCenter[1] = (voids[iVoid].barycenter[1]-boxLen[1]/2.)*100.;
outCenter[2] = (voids[iVoid].barycenter[2]-boxLen[2]/2.)*100.;
}
fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f\n",
outCenter[0],
outCenter[1],
outCenter[2],
voids[iVoid].vol,
voids[iVoid].radius,
voids[iVoid].redshift,
4./3.*M_PI*pow(voids[iVoid].radius, 3),
voids[iVoid].voidID,
voids[iVoid].densCon);
fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n",
atan((voids[iVoid].barycenter[1]-boxLen[1]/2.) /
(voids[iVoid].barycenter[0]-boxLen[0]/2.)) * 180/M_PI + 180,
asin((voids[iVoid].barycenter[2]-boxLen[2]/2.) /
voids[iVoid].redshiftInMpc) * 180/M_PI,
voids[iVoid].redshift,
voids[iVoid].radius,
voids[iVoid].voidID);
fprintf(fpShapes, "%d %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",
voids[iVoid].voidID,
gsl_vector_get(voids[iVoid].eval, 0),
gsl_vector_get(voids[iVoid].eval, 1),
gsl_vector_get(voids[iVoid].eval, 2),
gsl_matrix_get(voids[iVoid].evec, 0 ,0),
gsl_matrix_get(voids[iVoid].evec, 0 ,1),
gsl_matrix_get(voids[iVoid].evec, 0 ,2),
gsl_matrix_get(voids[iVoid].evec, 1 ,0),
gsl_matrix_get(voids[iVoid].evec, 1 ,1),
gsl_matrix_get(voids[iVoid].evec, 1 ,2),
gsl_matrix_get(voids[iVoid].evec, 2 ,0),
gsl_matrix_get(voids[iVoid].evec, 2 ,1),
gsl_matrix_get(voids[iVoid].evec, 2 ,2)
);
}
fclose(fp);
fclose(fpInfo);
fclose(fpBarycenter);
fclose(fpDistances);
// print the centers catalog again but without central density cuts
fpInfo = fopen(args_info.outNoCutInfo_arg, "w");
fprintf(fpInfo, "# center x,y,z (km/s), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID\n");
for (iVoid = 0; iVoid < numVoids; iVoid++) {
if (voids[iVoid].accepted < -1) continue;
double outCenter[3];
outCenter[0] = voids[iVoid].barycenter[0];
outCenter[1] = voids[iVoid].barycenter[1];
outCenter[2] = voids[iVoid].barycenter[2];
if (args_info.isObservation_flag) {
outCenter[0] = (voids[iVoid].barycenter[0]-boxLen[0]/2.)*100.;
outCenter[1] = (voids[iVoid].barycenter[1]-boxLen[1]/2.)*100.;
outCenter[2] = (voids[iVoid].barycenter[2]-boxLen[2]/2.)*100.;
}
fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f\n", fpDistancesCentral = fopen((std::string(args.outputDir_arg)+"boundaryDistancesCentral_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
outCenter[0], fpDistancesAll = fopen((std::string(args.outputDir_arg)+"boundaryDistancesAll_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
outCenter[1],
outCenter[2], fpSkyPositionsCentral = fopen((std::string(args.outputDir_arg)+"/sky_positions_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
voids[iVoid].vol, fprintf(fpSkyPositionsCentral, "# RA, dec, redshift, radius (Mpc/h), void ID\n");
voids[iVoid].radius,
voids[iVoid].redshift, fpSkyPositionsAll = fopen((std::string(args.outputDir_arg)+"/sky_positions_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
4./3.*M_PI*pow(voids[iVoid].radius, 3), fprintf(fpSkyPositionsAll, "# RA, dec, redshift, radius (Mpc/h), void ID\n");
voids[iVoid].voidID,
voids[iVoid].densCon); fpShapesCentral = fopen((std::string(args.outputDir_arg)+"/shapes_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
} fprintf(fpShapesCentral, "# void ID, 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");
fclose(fpInfo);
fpShapesAll = fopen((std::string(args.outputDir_arg)+"/shapes_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpShapesAll, "# void ID, 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");
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].voidType == CENTRAL_VOID) {
outputVoid(iVoid, voids[iVoid], fpZobovCentral, fpCentersCentral,
fpCentersNoCutCentral, fpSkyPositionsCentral,
fpBarycenterCentral, fpDistancesCentral, fpShapesCentral,
args.isObservation_flag, boxLen);
}
if (voids[iVoid].voidType == EDGE_VOID ||
voids[iVoid].voidType == CENTRAL_VOID) {
outputVoid(iVoid, voids[iVoid], fpZobovAll, fpCentersAll,
fpCentersNoCutAll, fpSkyPositionsAll,
fpBarycenterAll, fpDistancesAll, fpShapesAll,
args.isObservation_flag, boxLen);
}
}
fclose(fpZobovCentral);
fclose(fpZobovAll);
fclose(fpCentersCentral);
fclose(fpCentersAll);
fclose(fpCentersNoCutCentral);
fclose(fpCentersNoCutAll);
fclose(fpBarycenterCentral);
fclose(fpBarycenterAll);
fclose(fpDistancesCentral);
fclose(fpDistancesAll);
fclose(fpShapesCentral);
fclose(fpShapesAll);
fclose(fpSkyPositionsCentral);
fclose(fpSkyPositionsAll);
clock2 = clock(); clock2 = clock();
printf(" Time: %f sec (for %d voids)\n", (1.*clock2-clock1)/CLOCKS_PER_SEC, numVoids); printf(" Time: %f sec (for %d voids)\n",
(1.*clock2-clock1)/CLOCKS_PER_SEC, numVoids);
clock1 = clock(); clock1 = clock();
printf("Done!\n"); printf("Done!\n");
return 0; return 0;
} // end main } // end main
// ----------------------------------------------------------------------------
void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters,
FILE* fpCenterNoCut, FILE* fpSkyPositions,
FILE* fpBarycenters, FILE* fpDistances, FILE* fpShapes,
bool isObservation, double *boxLen) {
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(fpBarycenters, "%d %e %e %e\n",
outVoid.voidID,
outVoid.barycenter[0],
outVoid.barycenter[1],
outVoid.barycenter[2]);
fprintf(fpDistances, "%d %e\n",
outVoid.voidID,
outVoid.nearestMock);
double outCenter[3];
outCenter[0] = outVoid.barycenter[0];
outCenter[1] = outVoid.barycenter[1];
outCenter[2] = outVoid.barycenter[2];
if (isObservation) {
outCenter[0] = (outVoid.barycenter[0]-boxLen[0]/2.)*100.;
outCenter[1] = (outVoid.barycenter[1]-boxLen[1]/2.)*100.;
outCenter[2] = (outVoid.barycenter[2]-boxLen[2]/2.)*100.;
}
if (outVoid.accepted == 1) {
fprintf(fpCenters, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f\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);
}
fprintf(fpCenterNoCut, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f\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);
fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n",
atan((outVoid.barycenter[1]-boxLen[1]/2.) /
(outVoid.barycenter[0]-boxLen[0]/2.)) * 180/M_PI + 180,
asin((outVoid.barycenter[2]-boxLen[2]/2.) /
outVoid.redshiftInMpc) * 180/M_PI,
outVoid.redshift,
outVoid.radius,
outVoid.voidID);
fprintf(fpShapes, "%d %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",
outVoid.voidID,
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, 0 ,1),
gsl_matrix_get(outVoid.evec, 0 ,2),
gsl_matrix_get(outVoid.evec, 1 ,0),
gsl_matrix_get(outVoid.evec, 1 ,1),
gsl_matrix_get(outVoid.evec, 1 ,2),
gsl_matrix_get(outVoid.evec, 2 ,0),
gsl_matrix_get(outVoid.evec, 2 ,1),
gsl_matrix_get(outVoid.evec, 2 ,2)
);
} // end outputVoid

17
c_tools/stacking/pruneVoids.ggo
View file

@ -26,21 +26,8 @@ option "zMax" - "Maximum redshift of sample" double optional default="10.0"
option "rMin" - "Minimum allowable void radius" double optional default="0.0" option "rMin" - "Minimum allowable void radius" double optional default="0.0"
option "outputDir" - "Directory to place outputs" string required
option "output" - "Output void file" string required option "sampleName" - "unique string to assign to outputs" string required
option "outDistances" - "output of distances from centers to nearest mock particle" string required
option "outCenters" - "output barycenters of voids" string required
option "outInfo" - "output info of voids" string required
option "outNoCutInfo" - "output info of voids" string required
option "outSkyPositions" - "output sky positions of voids" string required
option "outShapes" - "output shape information of voids" string required
option "dataPortion" - "all, central, or edge" string required
option "periodic" - "Set of edges which are periodic" string optional default="xy" option "periodic" - "Set of edges which are periodic" string optional default="xy"

37
python_tools/void_python_tools/backend/launchers.py
View file

@ -247,7 +247,7 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
os.unlink(vozScript) os.unlink(vozScript)
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
def launchPrune(sample, binPath, thisDataPortion=None, def launchPrune(sample, binPath,
summaryFile=None, logFile=None, zobovDir=None, summaryFile=None, logFile=None, zobovDir=None,
continueRun=None): continueRun=None):
@ -287,7 +287,6 @@ def launchPrune(sample, binPath, thisDataPortion=None,
cmd += " --extraInfo=" + zobovDir+"/zobov_slice_"+str(sampleName)+\ cmd += " --extraInfo=" + zobovDir+"/zobov_slice_"+str(sampleName)+\
".par" ".par"
cmd += " --tolerance=1.0" cmd += " --tolerance=1.0"
cmd += " --dataPortion=" + thisDataPortion
cmd += " --mockIndex=" + str(mockIndex) cmd += " --mockIndex=" + str(mockIndex)
cmd += " --maxCentralDen=" + str(maxDen) cmd += " --maxCentralDen=" + str(maxDen)
cmd += " --zMin=" + str(sample.zRange[0]) cmd += " --zMin=" + str(sample.zRange[0])
@ -296,27 +295,8 @@ def launchPrune(sample, binPath, thisDataPortion=None,
cmd += " --numVoids=" + str(numVoids) cmd += " --numVoids=" + str(numVoids)
cmd += observationLine cmd += observationLine
cmd += periodicLine cmd += periodicLine
cmd += " --output=" + zobovDir+"/voidDesc_"+\ cmd += " --outputDir=" + zobovDir
str(thisDataPortion)+"_"+\ cmd += " --sampleName=" + str(sampleName)
str(sampleName)+".out"
cmd += " --outCenters=" + zobovDir+"/barycenters_"+\
str(thisDataPortion)+"_"+\
str(sampleName)+".out"
cmd += " --outInfo=" + zobovDir+"/centers_"+\
str(thisDataPortion)+"_"+\
str(sampleName)+".out"
cmd += " --outNoCutInfo=" + zobovDir+"/centers_nocut_"+\
str(thisDataPortion)+"_"+\
str(sampleName)+".out"
cmd += " --outSkyPositions=" + zobovDir+"/sky_positions_"+\
str(thisDataPortion)+"_"+\
str(sampleName)+".out"
cmd += " --outShapes=" + zobovDir+"/shapes_"+\
str(thisDataPortion)+"_"+\
str(sampleName)+".out"
cmd += " --outDistances=" + zobovDir+"/boundaryDistances_"+\
str(thisDataPortion)+"_"+\
str(sampleName)+".out"
cmd += " &> " + logFile cmd += " &> " + logFile
f=file("run_prune.sh",mode="w") f=file("run_prune.sh",mode="w")
f.write(cmd) f.write(cmd)
@ -956,7 +936,7 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
return return
numVoids = int(open(voidDir+"/num_voids.txt", "r").readline()) numVoids = int(open(voidDir+"/num_voids.txt", "r").readline())
if numVoids < 15: if numVoids < 10:
print "not enough voids to fit; skipping!" print "not enough voids to fit; skipping!"
fp = open(voidDir+"/NOFIT", "w") fp = open(voidDir+"/NOFIT", "w")
fp.write("not enough voids: %d \n" % numVoids) fp.write("not enough voids: %d \n" % numVoids)
@ -1334,8 +1314,9 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
plotTitle = "all samples, incoherent "+\ plotTitle = "all samples, incoherent "+\
thisDataPortion+" voids (systematics corrected)" thisDataPortion+" voids (systematics corrected)"
else: else:
plotTitle = "all samples, "+thisDataPortion+\ #plotTitle = "all samples, "+thisDataPortion+\
" voids (systematics corrected)" # " voids (systematics corrected)"
plotTitle = setName + "(sysematics corrected)"
vp.do_all_obs(zbase, allExpList, aveDistList, vp.do_all_obs(zbase, allExpList, aveDistList,
rlist, plotTitle=plotTitle, sampleNames=shortSampleNames, rlist, plotTitle=plotTitle, sampleNames=shortSampleNames,
plotAve=True, mulfac = 1.16, plotAve=True, mulfac = 1.16,
@ -1429,7 +1410,9 @@ def launchLikelihood(dataPortions=None, logDir=None, workDir=None,
OmStart = 0.0, OmStart = 0.0,
OmEnd = 1.0, OmEnd = 1.0,
biasStart = 1.0, biasStart = 1.0,
biasEnd = 1.2, biasEnd = 1.32,
#biasStart = 1.15,
#biasEnd = 1.17,
outputBase = workDir+"/1dlikelihoods_"+thisDataPortion+"_", outputBase = workDir+"/1dlikelihoods_"+thisDataPortion+"_",
useBinAve = False) useBinAve = False)