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

Merge

Browse source
This commit is contained in:
Guilhem Lavaux 2012-11-18 21:57:36 -05:00
commit 9e79c784f3
29 changed files with 1679 additions and 2007 deletions
Download patch file Download diff file Expand all files Collapse all files

18
GetQhull.cmake
View file

@ -1,18 +0,0 @@
SET(QHULL_BASE_PATH CACHE PATH "Qhull base path")
find_path(QHULL_INCLUDE_PATH qhull_a.h HINTS ${QHULL_BASE_PATH}/src/libqhull)
find_path(QHULL_CPP_INCLUDE_PATH Qhull.h HINTS ${QHULL_BASE_PATH}/src/libqhullcpp)
find_library(QHULL_LIBRARY qhull_p HINTS ${QHULL_BASE_PATH}/lib)
find_library(QHULL_CPP_LIBRARY qhullcpp HINTS ${QHULL_BASE_PATH}/lib)
find_library(QHULL_P_LIBRARY qhullstatic_p HINTS ${QHULL_BASE_PATH}/lib)
if ((NOT QHULL_INCLUDE_PATH) OR (NOT QHULL_CPP_LIBRARY))
message(SEND_ERROR "Qhull library not found")
endif((NOT QHULL_INCLUDE_PATH) OR (NOT QHULL_CPP_LIBRARY))
set(QHULL_INCLUDES ${QHULL_INCLUDE_PATH} ${QHULL_INCLUDE_PATH}/.. ${QHULL_CPP_INCLUDE_PATH} ${QHULL_BASE_PATH}/src)
set(QHULL_LIBRARIES ${QHULL_CPP_LIBRARY} ${QHULL_P_LIBRARY})
add_definitions(-Dqh_QHpointer)
mark_as_advanced(QHULL_INCLUDE_PATH QHULL_CPP_INCLUDE_PATH QHULL_LIBRARY QHULL_CPP_LIBRARY QHULL_P_LIBRARY)

30
c_tools/libzobov/contour_pixels.cpp
View file

@ -54,3 +54,33 @@ void computeContourPixels(Healpix_Map<float>& m, vector<int>& contour)
write_Healpix_map_to_fits(h, contour_map, planckType<int>());
}
}
void computeMaskPixels(Healpix_Map<float>& m, vector<int>& contour)
{
for (int p = 0; p < m.Npix(); p++)
{
if (m[p]>0)
{
contour.push_back(p);
// This is boundary go to next pixel
}
}
if (DEBUG)
{
Healpix_Map<int> contour_map;
contour_map.SetNside(m.Nside(), RING);
contour_map.fill(0);
for (int p = 0; p < contour.size(); p++)
{
contour_map[contour[p]]=1;
}
fitshandle h;
h.create("!mask_map.fits");
write_Healpix_map_to_fits(h, contour_map, planckType<int>());
}
}

1
c_tools/libzobov/contour_pixels.hpp
View file

@ -7,5 +7,6 @@
void computeContourPixels(Healpix_Map<float>& map, std::vector<int>& contour);
void computeFilledPixels(Healpix_Map<float>& map, std::vector<int>& contour);
void computeMaskPixels(Healpix_Map<float>& map, std::vector<int>& contour);
#endif

4
c_tools/libzobov/voidTree.hpp
View file

@ -126,8 +126,8 @@ public:
// Go bigger, though I would say we should not to.
}
while (iter_candidate != candidateList.begin()) ;
if (vid_good_candidate < 0)
std::cout << "Failure to lookup parent (2)" << std::endl;
//if (vid_good_candidate < 0)
// std::cout << "Failure to lookup parent (2)" << std::endl;
return vid_good_candidate;
}

75
c_tools/mock/generateFromCatalog.cpp
View file

@ -27,6 +27,7 @@ struct NYU_Data
double cz;
double fgotten;
double phi_z;
double uniqueID;
};
struct Position
@ -40,6 +41,8 @@ struct ParticleData
vector<double> ra;
vector<double> dec;
vector<double> redshift;
vector<double> catalogID;
vector<double> uniqueID;
int id_mask;
// PMS
int mask_index;
@ -85,6 +88,7 @@ void loadData(const string& fname, NYU_VData & data)
{
NYU_Data d;
f >> d.index >> d.sector >> d.region >> d.ra >> d.dec >> d.cz >> d.fgotten >> d.phi_z;
d.uniqueID = d.index;
data.push_back(d);
}
}
@ -125,6 +129,7 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data,
output_data.ra.resize(data.size());
output_data.dec.resize(data.size());
output_data.redshift.resize(data.size());
output_data.uniqueID.resize(data.size());
for (int j = 0; j < 3; j++)
{
@ -143,7 +148,6 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data,
1.e-6,
1.e-6, &result, &error, &nEval);
double Dc = result*LIGHT_SPEED;
cout << "HELLO " << data[i].cz << " " << Dc << endl;
p.xyz[0] = Dc*cos(ra)*cos(dec);
p.xyz[1] = Dc*sin(ra)*cos(dec);
p.xyz[2] = Dc*sin(dec);
@ -157,6 +161,7 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data,
output_data.ra[i] = ra;
output_data.dec[i] = dec;
output_data.redshift[i] = data[i].cz;
output_data.uniqueID[i] = data[i].uniqueID;
for (int j = 0; j < 3; j++)
{
@ -273,6 +278,7 @@ void generateSurfaceMask(generateFromCatalog_info& args ,
output_data.ra.push_back(-1);
output_data.dec.push_back(-1);
output_data.redshift.push_back(-1);
output_data.uniqueID.push_back(-1);
//printf("INSERT MOCK %d %e %e %e\n", idx, p.xyz[0], p.xyz[1], p.xyz[2]);
insertion++;
}
@ -289,10 +295,10 @@ void generateSurfaceMask(generateFromCatalog_info& args ,
int nPart = 100;
// TEST
for (int iDir = 0; iDir < 0; iDir++) {
for (int iFace = 0; iFace < 0; iFace++) {
//for (int iDir = 0; iDir < 3; iDir++) {
//for (int iFace = 0; iFace < 2; iFace++) {
//for (int iDir = 0; iDir < 0; iDir++) {
//for (int iFace = 0; iFace < 0; iFace++) {
for (int iDir = 0; iDir < 3; iDir++) {
for (int iFace = 0; iFace < 2; iFace++) {
int iy = (iDir + 1) % 3;
int iz = (iDir + 2) % 3;
@ -310,6 +316,7 @@ void generateSurfaceMask(generateFromCatalog_info& args ,
output_data.ra.push_back(-1);
output_data.dec.push_back(-1);
output_data.redshift.push_back(-1);
output_data.uniqueID.push_back(-1);
insertion++;
fprintf(fp, "%e %e %e\n",
@ -341,6 +348,7 @@ void generateSurfaceMask(generateFromCatalog_info& args ,
output_data.ra.push_back(-1);
output_data.dec.push_back(-1);
output_data.redshift.push_back(-1);
output_data.uniqueID.push_back(-1);
insertion++;
fprintf(fp, "%e %e %e\n",
(p.xyz[0]),
@ -357,6 +365,7 @@ void generateSurfaceMask(generateFromCatalog_info& args ,
output_data.ra.push_back(-1);
output_data.dec.push_back(-1);
output_data.redshift.push_back(-1);
output_data.uniqueID.push_back(-1);
insertion++;
fprintf(fp, "%e %e %e\n",
(p.xyz[0]),
@ -401,6 +410,7 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn
UnformattedWrite f(fname);
static const char axis[] = { 'X', 'Y', 'Z' };
double Lmax = pdata.Lmax;
double r2d = 180./M_PI;
f.beginCheckpoint();
f.writeInt32(pdata.pos.size());
@ -418,6 +428,34 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn
f.endCheckpoint();
}
cout << format("Writing RA...") << endl;
f.beginCheckpoint();
for (uint32_t i = 0; i < pdata.pos.size(); i++) {
f.writeReal32(pdata.ra[i]*r2d);
}
f.endCheckpoint();
cout << format("Writing Dec...") << endl;
f.beginCheckpoint();
for (uint32_t i = 0; i < pdata.pos.size(); i++) {
f.writeReal32(pdata.dec[i]*r2d);
}
f.endCheckpoint();
cout << format("Writing Redshift...") << endl;
f.beginCheckpoint();
for (uint32_t i = 0; i < pdata.pos.size(); i++) {
f.writeReal32(pdata.redshift[i]);
}
f.endCheckpoint();
cout << format("Writing Unique ID...") << endl;
f.beginCheckpoint();
for (uint32_t i = 0; i < pdata.pos.size(); i++) {
f.writeReal32(pdata.uniqueID[i]);
}
f.endCheckpoint();
NcFile fp(paramname.c_str(), NcFile::Replace);
fp.add_att("range_x_min", -Lmax/100.);
@ -433,9 +471,6 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn
NcDim *NumPart_dim = fp.add_dim("numpart_dim", nOutputPart);
NcVar *v = fp.add_var("particle_ids", ncInt, NumPart_dim);
NcVar *vra = fp.add_var("RA", ncInt, NumPart_dim);
NcVar *vdec = fp.add_var("DEC", ncInt, NumPart_dim);
NcVar *vredshift = fp.add_var("z", ncInt, NumPart_dim);
//NcVar *v2 = fp.add_var("expansion", ncDouble, NumPart_dim);
//double *expansion_fac = new double[pdata.pos.size()];
@ -444,12 +479,23 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn
// expansion_fac[i] = 1.0;
v->put(&pdata.id_gal[0], nOutputPart);
vra->put(&pdata.ra[0], nOutputPart);
vdec->put(&pdata.dec[0], nOutputPart);
vredshift->put(&pdata.redshift[0], nOutputPart);
//v2->put(expansion_fac, pdata.pos.size());
//delete[] expansion_fac;
/*
FILE *infoFile = fopen("sample_info.txt", "w");
fprintf(infoFile, "x_min = %f\n", -Lmax/100.);
fprintf(infoFile, "x_max = %f\n", Lmax/100.);
fprintf(infoFile, "y_min = %f\n", -Lmax/100.);
fprintf(infoFile, "y_max = %f\n", Lmax/100.);
fprintf(infoFile, "z_min = %f\n", -Lmax/100.);
fprintf(infoFile, "z_max = %f\n", Lmax/100.);
fprintf(infoFile, "mask_index = %d\n", pdata.mask_index);
fprintf(infoFile, "total_particles = %d\n", pdata.pos.size());
fclose(infoFile);
*/
}
int main(int argc, char **argv)
@ -500,6 +546,7 @@ int main(int argc, char **argv)
mask.Import(o_mask);
computeContourPixels(mask,pixel_list);
computeMaskPixels(mask,full_mask_list);
// We compute a cube holding all the galaxies + the survey surface mask
@ -516,12 +563,6 @@ int main(int argc, char **argv)
fprintf(fp, "%d", output_data.mask_index);
fclose(fp);
fp = fopen("sample_info.txt", "w");
fprintf(fp, "Lmax = %f\n", output_data.Lmax);
fprintf(fp, "mask_index = %d\n", output_data.mask_index);
fprintf(fp, "total_particles = %d\n", output_data.pos.size());
fclose(fp);
fp = fopen("total_particles.txt", "w");
fprintf(fp, "%d", output_data.pos.size());
fclose(fp);

222
c_tools/mock/generateMock.cpp
View file

@ -18,9 +18,67 @@ using namespace CosmoTool;
#define LIGHT_SPEED 299792.458
static double gadgetUnit=1e-3;
SimuData *doLoadRamses(const char *basename, int baseid, int velAxis, bool goRedshift)
{
SimuData *d, *outd;
d = loadRamsesSimu(basename, baseid, -1, true, 0);
outd = new SimuData;
outd->NumPart = d->TotalNumPart;
outd->BoxSize = d->BoxSize;
outd->TotalNumPart = outd->NumPart;
outd->Hubble = d->Hubble;
outd->Omega_Lambda = d->Omega_Lambda;
outd->Omega_M = d->Omega_M;
outd->time = d->time;
for (int k = 0; k < 3; k++)
outd->Pos[k] = new float[outd->NumPart];
outd->Vel[2] = new float[outd->NumPart];
delete d;
int curCpu = 0;
cout << "loading cpu 0 " << endl;
while (d = loadRamsesSimu(basename, baseid, curCpu, true, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID))
{
for (int k = 0; k < 3; k++)
for (int i = 0; i < d->NumPart; i++)
{
assert(d->Id[i] >= 1);
assert(d->Id[i] <= outd->TotalNumPart);
outd->Pos[k][d->Id[i]-1] = d->Pos[k][i];
outd->Vel[2][d->Id[i]-1] = d->Vel[velAxis][i];
}
if (goRedshift)
for (int i = 0; i < d->NumPart; i++)
outd->Pos[velAxis][d->Id[i]-1] += d->Vel[velAxis][i]/100.;
delete d;
curCpu++;
cout << "loading cpu " << curCpu << endl;
}
return outd;
}
SimuData *myLoadGadget(const char *fname, int id, int flags)
{
return loadGadgetMulti(fname, id, flags);
SimuData *sim = loadGadgetMulti(fname, id, flags);
sim->BoxSize *= gadgetUnit*1000;
for (int j = 0; j < 3; j++)
{
if (sim->Pos[j] != 0) {
for (long i = 0; i < sim->NumPart; i++)
sim->Pos[j][i] *= gadgetUnit*1000;
}
}
return sim;
}
SimuData *doLoadSimulation(const char *gadgetname, int velAxis, bool goRedshift, SimuData *(*loadFunction)(const char *fname, int id, int flags))
@ -111,7 +169,7 @@ SimuData *doLoadMultidark(const char *multidarkname)
fscanf(fp, "%f\n", &outd->Omega_M);
fscanf(fp, "%f\n", &outd->Hubble);
fscanf(fp, "%f\n", &outd->time);
fscanf(fp, "%d\n", &outd->NumPart);
fscanf(fp, "%ld\n", &outd->NumPart);
outd->time = 1./(1.+outd->time); // convert to scale factor
outd->TotalNumPart = outd->NumPart;
@ -120,14 +178,17 @@ SimuData *doLoadMultidark(const char *multidarkname)
for (int k = 0; k < 3; k++)
outd->Pos[k] = new float[outd->NumPart];
outd->Vel[2] = new float[outd->NumPart];
outd->Id = new int[outd->NumPart];
cout << "loading multidark particles" << endl;
actualNumPart = 0;
for (int i = 0; i < outd->NumPart; i++) {
fscanf(fp, "%f %f %f %f\n", &outd->Pos[0][i], &outd->Pos[1][i],
fscanf(fp, "%d %d %f %f %f\n", &outd->Id[i],
&outd->Pos[0][i], &outd->Pos[1][i],
&outd->Pos[2][i], &outd->Vel[2][i]);
if (outd->Pos[0][i] == -99 && outd->Pos[1][i] == -99 &&
if (outd->Id[i] == -99 &&
outd->Pos[0][i] == -99 && outd->Pos[1][i] == -99 &&
outd->Pos[2][i] == -99 && outd->Vel[2][i] == -99) {
break;
} else {
@ -181,7 +242,7 @@ Interpolate make_cosmological_redshift(double OM, double OL, double z0, double z
return buildFromVector(pairs);
}
void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact)
void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact, bool cosmo_flag)
{
int x0, x1, x2;
@ -200,7 +261,11 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double
}
double z0 = 1/data->time - 1;
Interpolate z_vs_D = make_cosmological_redshift(data->Omega_M, data->Omega_Lambda, 0., z0+8*data->BoxSize*100/LIGHT_SPEED); // Redshift 2*z0 should be sufficient ?
Interpolate z_vs_D =
make_cosmological_redshift(data->Omega_M, data->Omega_Lambda,
0., z0+8*data->BoxSize*100/LIGHT_SPEED);
// Redshift 2*z0 should be sufficient ? This is fragile.
//A proper solver is needed here.
double z_base = reshift ? z0 : 0;
TotalExpansion e_computer;
double baseComovingDistance;
@ -229,13 +294,15 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double
// Distorted redshift
if (reduced_red == 0)
z = 0;
else {
z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.;
}
else if (cosmo_flag)
z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.;
else
z = reduced_red*LIGHT_SPEED/100.0;
expfact[i] = z / z_old;
// Add peculiar velocity
if (pecvel)
z += v/100;
// Add peculiar velocity
if (pecvel)
z += v/100;
}
catch(const InvalidRangeException& e) {
cout << "Trying to interpolate out of the tabulated range." << endl;
@ -296,6 +363,40 @@ void generateOutput(SimuData *data, int axis,
f.writeReal32(data->Pos[x2][i]);
}
f.endCheckpoint();
cout << "Writing RA..." << endl;
f.beginCheckpoint();
for (uint32_t i = 0; i < data->NumPart; i++)
{
f.writeReal32(data->Id[i]);
}
f.endCheckpoint();
cout << "Writing Dec..." << endl;
f.beginCheckpoint();
for (uint32_t i = 0; i < data->NumPart; i++)
{
f.writeReal32(data->Id[i]);
}
f.endCheckpoint();
cout << "Writing redshift..." << endl;
f.beginCheckpoint();
for (uint32_t i = 0; i < data->NumPart; i++)
{
f.writeReal32(data->Id[i]);
}
f.endCheckpoint();
cout << "Writing unique ID..." << endl;
f.beginCheckpoint();
for (uint32_t i = 0; i < data->NumPart; i++)
{
f.writeReal32(data->Id[i]);
}
f.endCheckpoint();
}
void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& args_info)
@ -391,6 +492,7 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info&
f.add_att("range_y_max", ranges[1][1]);
f.add_att("range_z_min", ranges[2][0]);
f.add_att("range_z_max", ranges[2][1]);
f.add_att("mask_index", -1);
NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim);
@ -401,6 +503,77 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info&
delete[] particle_id;
delete[] expansion_fac;
/*
FILE *fp = fopen("sample_info.txt", "w");
fprintf(fp, "x_min = %f\n", ranges[0][0]);
fprintf(fp, "x_max = %f\n", ranges[0][1]);
fprintf(fp, "y_min = %f\n", ranges[1][0]);
fprintf(fp, "y_max = %f\n", ranges[1][1]);
fprintf(fp, "z_min = %f\n", ranges[2][0]);
fprintf(fp, "z_max = %f\n", ranges[2][1]);
fprintf(fp, "mask_index = -1\n");
fprintf(fp, "total_particles = %d\n", boxed->NumPart);
fclose(fp);
*/
}
void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, generateMock_info& args_info)
{
NcFile f(args_info.inputParameter_arg);
NcVar *v;
int *particle_id;
double *expansion_fac;
boxed = new SimuData;
boxed->Hubble = simu->Hubble;
boxed->Omega_M = simu->Omega_M;
boxed->Omega_Lambda = simu->Omega_Lambda;
boxed->time = simu->time;
boxed->BoxSize = simu->BoxSize;
NcVar *v_id = f.get_var("particle_ids");
long *edges1;
double ranges[3][2];
double mul[3];
edges1 = v_id->edges();
assert(v_id->num_dims()==1);
boxed->NumPart = edges1[0];
delete[] edges1;
particle_id = new int[boxed->NumPart];
v_id->get(particle_id, boxed->NumPart);
ranges[0][0] = f.get_att("range_x_min")->as_double(0);
ranges[0][1] = f.get_att("range_x_max")->as_double(0);
ranges[1][0] = f.get_att("range_y_min")->as_double(0);
ranges[1][1] = f.get_att("range_y_max")->as_double(0);
ranges[2][0] = f.get_att("range_z_min")->as_double(0);
ranges[2][1] = f.get_att("range_z_max")->as_double(0);
for (int j = 0; j < 3; j++)
{
boxed->Pos[j] = new float[boxed->NumPart];
boxed->Vel[j] = 0;
mul[j] = 1.0/(ranges[j][1] - ranges[j][0]);
}
uint32_t k = 0;
for (uint32_t i = 0; i < boxed->NumPart; i++)
{
int id = particle_id[i];
for (int j = 0; j < 3; j++)
{
boxed->Pos[j][i] = (simu->Pos[j][id]-ranges[j][0])*mul[j];
}
}
delete[] particle_id;
}
int main(int argc, char **argv)
@ -433,13 +606,15 @@ int main(int argc, char **argv)
generateMock_conf_print_version();
gadgetUnit=args_info.gadgetUnit_arg;
if (args_info.ramsesBase_given || args_info.ramsesId_given)
{
if (args_info.ramsesBase_given && args_info.ramsesId_given) {
//simu = doLoadRamses(args_info.ramsesBase_arg,
// args_info.ramsesId_arg,
// args_info.axis_arg, false);
}
simu = doLoadRamses(args_info.ramsesBase_arg,
args_info.ramsesId_arg,
args_info.axis_arg, false);
}
else
{
cerr << "Both ramsesBase and ramsesId are required to be able to load snapshots" << endl;
@ -489,14 +664,15 @@ int main(int argc, char **argv)
double *expfact;
if (args_info.cosmo_flag){
metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag, args_info.peculiarVelocities_flag, expfact);
} else {
expfact = new double[simu->NumPart];
for (int j = 0; j < simu->NumPart; j++) expfact[j] = 1.0;
}
metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag,
args_info.peculiarVelocities_flag, expfact,
args_info.cosmo_flag);
if (args_info.inputParameter_given)
makeBoxFromParameter(simu, expfact, simuOut, args_info);
else
makeBox(simu, expfact, simuOut, args_info);
makeBox(simu, expfact, simuOut, args_info);
delete simu;
generateOutput(simuOut, args_info.axis_arg, args_info.output_arg);

6
c_tools/mock/generateMock.ggo
View file

@ -26,6 +26,10 @@ option "rangeZ_max" - "Maximum range in Z for making the box (after distorti
option "preReShift" - "Reshift the zero of the Z axis" flag off
option "peculiarVelocities" - "Added peculiar velocities distortion" flag off
option "cosmo" - "Apply cosmological redshift" flag on
option "cosmo" - "Apply cosmological redshift" flag off
option "subsample" - "Subsample the input simulation by the specified amount" double optional
option "inputParameter" - "Input geometry (optional, warning!)" string optional
option "gadgetUnit" - "Unit of length in gadget file in Mpc/h" double optional default="0.001"

148
c_tools/stacking/pruneVoids.cpp
View file

@ -39,7 +39,7 @@ typedef struct voidZoneStruct {
typedef struct voidStruct {
float vol, coreDens, zoneVol, densCon, voidProb, radius;
int voidID, numPart, numZones, coreParticle, zoneNumPart;
float maxRadius, nearestMock, centralDen, redshift;
float maxRadius, nearestMock, centralDen, redshift, redshiftInMpc;
float center[3], barycenter[3];
int accepted;
} VOID;
@ -89,15 +89,33 @@ int main(int argc, char **argv) {
double ranges[2][3], boxLen[3], mul;
double volNorm, radius;
int clock1, clock2;
int periodicX=0, periodicY=0, periodicZ=0;
numVoids = args_info.numVoids_arg;
mockIndex = args_info.mockIndex_arg;
tolerance = args_info.tolerance_arg;
clock1 = clock();
printf("Pruning parameters: %f %f %f\n", args_info.zMin_arg,
args_info.zMax_arg,
args_info.rMin_arg);
printf("Pruning parameters: %f %f %f %s\n", args_info.zMin_arg,
args_info.zMax_arg,
args_info.rMin_arg,
args_info.periodic_arg);
// check for periodic box
if (!args_info.isObservation_flag) {
if ( strchr(args_info.periodic_arg, 'x') != NULL) {
periodicX = 1;
printf("Will assume x-direction is periodic.\n");
}
if ( strchr(args_info.periodic_arg, 'y') != NULL) {
periodicY = 1;
printf("Will assume y-direction is periodic.\n");
}
if ( strchr(args_info.periodic_arg, 'z') != NULL) {
periodicZ = 1;
printf("Will assume z-direction is periodic.\n");
}
}
// load box size
printf("\n Getting info...\n");
@ -254,8 +272,9 @@ int main(int argc, char **argv) {
for (iVoid = 0; iVoid < numVoids; iVoid++) {
voidID = voids[iVoid].voidID;
//printf(" DOING %d (of %d) %d %d\n", iVoid, numVoids, voidID,
// voids[iVoid].numPart);
printf(" DOING %d (of %d) %d %d %f\n", iVoid, numVoids, voidID,
voids[iVoid].numPart,
voids[iVoid].radius);
voids[iVoid].center[0] = part[voids[iVoid].coreParticle].x;
voids[iVoid].center[1] = part[voids[iVoid].coreParticle].y;
@ -290,17 +309,14 @@ int main(int argc, char **argv) {
voids[iVoid].barycenter[1] = 0.;
voids[iVoid].barycenter[2] = 0.;
// TODO handle periodic boundaries?
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 (!args_info.isObservation_flag) {
// dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0]));
// dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1]));
// dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2]));
//}
if (periodicX) dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0]));
if (periodicY) dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1]));
if (periodicZ) dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2]));
voids[iVoid].barycenter[0] += voidPart[p].vol*(dist[0]);
voids[iVoid].barycenter[1] += voidPart[p].vol*(dist[1]);
@ -323,11 +339,9 @@ int main(int argc, char **argv) {
dist[1] = voidPart[p].y - voids[iVoid].barycenter[1];
dist[2] = voidPart[p].z - voids[iVoid].barycenter[2];
//if (!args_info.isObservation_flag) {
// dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0]));
// dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1]));
// dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2]));
//}
if (periodicX) dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0]));
if (periodicY) dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1]));
if (periodicZ) dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2]));
dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2);
if (dist2 < centralRad) centralDen += 1;
@ -350,13 +364,17 @@ int main(int argc, char **argv) {
}
voids[iVoid].maxRadius = sqrt(maxDist)/2.;
} else {
maxDist = 0.;
maxDist = 0.;
for (p = 0; p < voids[iVoid].numPart; p++) {
dist[0] = voidPart[p].x - voids[iVoid].barycenter[0];
dist[0] = voidPart[p].y - voids[iVoid].barycenter[1];
dist[0] = voidPart[p].z - voids[iVoid].barycenter[2];
if (periodicX) dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0]));
if (periodicY) dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1]));
if (periodicZ) dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2]));
dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2);
if (dist2 > maxDist) maxDist = dist2;
}
@ -381,28 +399,41 @@ int main(int argc, char **argv) {
}
if (args_info.isObservation_flag) {
voids[iVoid].redshift =
voids[iVoid].redshiftInMpc =
sqrt(pow(voids[iVoid].barycenter[0] - boxLen[0]/2.,2) +
pow(voids[iVoid].barycenter[1] - boxLen[1]/2.,2) +
pow(voids[iVoid].barycenter[2] - boxLen[2]/2.,2));
voids[iVoid].redshift = voids[iVoid].redshift;
redshift = voids[iVoid].redshift;
nearestEdge = fmin(fabs(redshift-args_info.zMin_arg*LIGHT_SPEED),
fabs(redshift-args_info.zMax_arg*LIGHT_SPEED));
voids[iVoid].redshiftInMpc = voids[iVoid].redshiftInMpc;
redshift = voids[iVoid].redshiftInMpc;
nearestEdge = fmin(fabs(redshift-args_info.zMin_arg*LIGHT_SPEED/100.),
fabs(redshift-args_info.zMax_arg*LIGHT_SPEED/100.));
voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.;
} else {
voids[iVoid].redshiftInMpc = voids[iVoid].barycenter[2];
voids[iVoid].redshift = voids[iVoid].barycenter[2]/LIGHT_SPEED*100.;
nearestEdge = fmin(
fabs(voids[iVoid].barycenter[0] - ranges[0][0]),
fabs(voids[iVoid].barycenter[0] - ranges[0][1]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[1] - ranges[1][0]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[1] - ranges[1][1]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[2] - ranges[2][0]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[2] - ranges[2][1]));
nearestEdge = 1.e99;
if (!periodicX) {
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[0] - ranges[0][0]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[0] - ranges[0][1]));
}
if (!periodicY) {
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[1] - ranges[1][0]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[1] - ranges[1][1]));
}
if (!periodicZ) {
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[2] - ranges[2][0]));
nearestEdge = fmin(nearestEdge,
fabs(voids[iVoid].barycenter[2] - ranges[2][1]));
}
}
if (nearestEdge < voids[iVoid].nearestMock) {
@ -411,30 +442,44 @@ int main(int argc, char **argv) {
} // iVoid
printf(" Picking winners and losers...\n");
numKept = numVoids;
for (iVoid = 0; iVoid < numVoids; iVoid++) {
voids[iVoid].accepted = 1;
}
for (iVoid = 0; iVoid < numVoids; iVoid++) {
if (voids[iVoid].densCon < 1.5) {
//voids[iVoid].accepted = 0;
}
// toss out voids that are obviously wrong
if (voids[iVoid].densCon > 1.e3) {
voids[iVoid].accepted = 0;
}
if (strcmp(args_info.dataPortion_arg, "edge") == 0 &&
tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) {
numKept--;
voids[iVoid].accepted = 0;
}
if (strcmp(args_info.dataPortion_arg, "central") == 0 &&
tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestMock) {
numKept--;
voids[iVoid].accepted = 0;
}
if (voids[iVoid].radius < args_info.rMin_arg) {
voids[iVoid].accepted = 0;
}
if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) {
numKept--;
voids[iVoid].accepted = -1;
}
if (voids[iVoid].radius < args_info.rMin_arg) {
numKept--;
voids[iVoid].accepted = 0;
}
}
numKept = 0;
for (iVoid = 0; iVoid < numVoids; iVoid++) {
if (voids[iVoid].accepted == 1) numKept++;
}
printf(" Number kept: %d (out of %d)\n", numKept, numVoids);
printf(" Output...\n");
@ -445,14 +490,14 @@ int main(int argc, char **argv) {
fpSkyPositions = fopen(args_info.outSkyPositions_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\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");
for (iVoid = 0; iVoid < numVoids; iVoid++) {
if (voids[iVoid].accepted != 1) continue;
fprintf(fp, "%d %d %d %f %f %d %d %f %d %f %f\n",
i,
iVoid,
voids[iVoid].voidID,
voids[iVoid].coreParticle,
voids[iVoid].coreDens,
@ -485,7 +530,7 @@ int main(int argc, char **argv) {
outCenter[2] = (voids[iVoid].barycenter[2]-boxLen[2]/2.)*100.;
}
fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d\n",
fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f\n",
outCenter[0],
outCenter[1],
outCenter[2],
@ -493,12 +538,12 @@ int main(int argc, char **argv) {
voids[iVoid].radius,
voids[iVoid].redshift,
4./3.*M_PI*pow(voids[iVoid].radius, 3),
voids[iVoid].voidID
);
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].redshift) * 180/M_PI,
asin((voids[iVoid].barycenter[2]-boxLen[2]/2.)/voids[iVoid].redshiftInMpc) * 180/M_PI,
voids[iVoid].redshift,
voids[iVoid].radius,
voids[iVoid].voidID);
@ -528,7 +573,7 @@ int main(int argc, char **argv) {
}
fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d\n",
fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f\n",
outCenter[0],
outCenter[1],
outCenter[2],
@ -536,7 +581,8 @@ int main(int argc, char **argv) {
voids[iVoid].radius,
voids[iVoid].redshift,
4./3.*M_PI*pow(voids[iVoid].radius, 3),
voids[iVoid].voidID);
voids[iVoid].voidID,
voids[iVoid].densCon);
}
fclose(fpInfo);

2
c_tools/stacking/pruneVoids.ggo
View file

@ -40,6 +40,8 @@ option "outSkyPositions" - "output sky positions 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 "tolerance" - "Fraction of void width to consider edge" double optional default="1.0"
option "centralRadFrac" - "Fraction of void radii to consider central region" double optional default="4"

5
external/cosmotool/src/loadSimu.hpp vendored
View file

@ -24,9 +24,10 @@ namespace CosmoTool
int *Id;
float *Pos[3];
float *Vel[3];
float *uniqueID;
int *type;
public:
SimuData() : Id(0),NumPart(0),type(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; }
SimuData() : Id(0),NumPart(0),type(0), uniqueID(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; uniqueID=0}
~SimuData()
{
for (int j = 0; j < 3; j++)
@ -40,6 +41,8 @@ namespace CosmoTool
delete[] type;
if (Id)
delete[] Id;
if (uniqueID)
delete[] uniqueID;
}
};

20
external/external_build.cmake vendored
View file

@ -22,28 +22,34 @@ SET(INTERNAL_QHULL ON)
IF(INTERNAL_GENGETOPT)
SET(GENGETOPT_URL "ftp://ftp.gnu.org/gnu/gengetopt/gengetopt-2.22.5.tar.gz" CACHE STRING "URL to download gengetopt from")
mark_as_advanced(GENGETOPT_URL)
ENDIF(INTERNAL_GENGETOPT)
IF(INTERNAL_HDF5)
SET(HDF5_URL "http://www.hdfgroup.org/ftp/HDF5/current/src/hdf5-1.8.9.tar.gz" CACHE STRING "URL to download HDF5 from")
mark_as_advanced(HDF5_URL)
ENDIF(INTERNAL_HDF5)
IF(INTERNAL_NETCDF)
SET(NETCDF_URL "http://www.unidata.ucar.edu/downloads/netcdf/ftp/netcdf-4.1.3.tar.gz" CACHE STRING "URL to download NetCDF from")
mark_as_advanced(NETCDF_URL)
ENDIF(INTERNAL_NETCDF)
IF(INTERNAL_BOOST)
SET(BOOST_URL "http://sourceforge.net/projects/boost/files/boost/1.49.0/boost_1_49_0.tar.gz/download" CACHE STRING "URL to download Boost from")
mark_as_advanced(BOOST_URL)
ELSE(INTERNAL_BOOST)
find_package(Boost 1.49.0 COMPONENTS format spirit phoenix python FATAL_ERROR)
ENDIF(INTERNAL_BOOST)
IF(INTERNAL_GSL)
SET(GSL_URL "ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz" CACHE STRING "URL to download GSL from ")
mark_as_advanced(GSL_URL)
ENDIF(INTERNAL_GSL)
IF(INTERNAL_QHULL)
SET(QHULL_URL "http://www.qhull.org/download/qhull-2012.1-src.tgz" CACHE STRING "URL to download QHull from")
mark_as_advanced(QHULL_URL)
ENDIF(INTERNAL_QHULL)
@ -72,10 +78,11 @@ if (INTERNAL_GENGETOPT)
BUILD_IN_SOURCE 1
INSTALL_COMMAND make install
)
SET(GENGETOPT ${GENGETOPT_BIN_DIR}/bin/gengetopt)
SET(GENGETOPT ${GENGETOPT_BIN_DIR}/bin/gengetopt CACHE FILEPATH "Path GenGetOpt binary")
else(INTERNAL_GENGETOPT)
find_program(GENGETOPT gengetopt)
endif(INTERNAL_GENGETOPT)
mark_as_advanced(GENGETOPT)
###############
# Build HDF5
@ -113,6 +120,7 @@ else(INTERNAL_HDF5)
find_library(HDF5HL_LIBRARY hdf5_hl)
endif (INTERNAL_HDF5)
SET(CONFIGURE_CPP_FLAGS "${CONFIGURE_CPP_FLAGS} -I${HDF5_INCLUDE_PATH}")
mark_as_advanced(HDF5_INCLUDE_PATH HDF5_LIBRARY HDF5_CPP_LIBRARY HDF5HL_LIBRARY HDF5HL_CPP_LIBRARY)
###############
# Build NetCDF
@ -154,6 +162,7 @@ ELSE(INTERNAL_NETCDF)
SET(CONFIGURE_CPP_FLAGS ${CONFIGURE_CPP_FLAGS}
-I${NETCDF_INCLUDE_PATH} -I${NETCDFCPP_INCLUDE_PATH})
endif (INTERNAL_NETCDF)
mark_as_advanced(NETCDF_LIBRARY NETCDFCPP_LIBRARY NETCDF_INCLUDE_PATH NETCDFCPP_INCLUDE_PATH)
##################
# Build BOOST
@ -171,8 +180,9 @@ if (INTERNAL_BOOST)
INSTALL_COMMAND echo "No install"
)
set(Boost_INCLUDE_DIRS ${BOOST_SOURCE_DIR} CACHE STRING "Boost path" FORCE)
set(Boost_LIBRARIES ${BOOST_SOURCE_DIR}/stage/lib/libboost_python.a)
set(Boost_LIBRARIES ${BOOST_SOURCE_DIR}/stage/lib/libboost_python.a CACHE STRING "Boost libraries" FORCE)
endif (INTERNAL_BOOST)
mark_as_advanced(Boost_INCLUDE_DIRS Boost_LIBRARIES)
##################
# Build GSl
@ -201,6 +211,7 @@ ELSE(INTERNAL_GSL)
find_library(GSLCBLAS_LIBRARY gslcblas)
find_path(GSL_INCLUDE_PATH NAMES gsl/gsl_blas.h)
ENDIF(INTERNAL_GSL)
mark_as_advanced(GSL_LIBRARY GSLCBLAS_LIBRARY GSL_INCLUDE_PATH)
##################
# Build CosmoTool
@ -242,7 +253,9 @@ ExternalProject_Add(cfitsio
BUILD_IN_SOURCE 1
INSTALL_COMMAND make install
)
SET(CFITSIO_LIBRARY ${CMAKE_BINARY_DIR}/ext_build/cfitsio/lib/libcfitsio.a)
SET(CFITSIO_PREFIX ${CMAKE_BINARY_DIR}/ext_build/cfitsio)
SET(CFITSIO_LIBRARY ${CFITSIO_PREFIX}/lib/libcfitsio.a)
SET(CFITSIO_INCLUDE_PATH ${CFITSIO_PREFIX}/include)
#################
# Build Healpix
@ -297,6 +310,7 @@ if (INTERNAL_QHULL)
add_definitions(-Dqh_QHpointer)
else(INTERNAL_QHULL)
message(FATAL_ERROR "Only packaged QHull is supported")
endif(INTERNAL_QHULL)
SET(QHULL_LIBRARIES ${QHULL_CPP_LIBRARY} ${QHULL_LIBRARY} )

71
external/external_python_build.cmake vendored
View file

@ -2,16 +2,36 @@ INCLUDE(FindPythonInterp)
SET(INTERNAL_NETCDF4_PYTHON ON)
SET(INTERNAL_CYTHON ON)
SET(INTERNAL_HEALPY ON)
IF (PYTHON_VERSION_STRING VERSION_LESS 2.7)
MESSAGE(STATUS "Python version is less than 2.7, argparse is needed.")
SET(INTERNAL_ARGPARSE ON)
ELSE (PYTHON_VERSION_STRING VERSION_LESS 2.7)
MESSAGE(STATUS "Python version is greater than 2.7, argparse is already bundled.")
ENDIF (PYTHON_VERSION_STRING VERSION_LESS 2.7)
IF(INTERNAL_CYTHON)
SET(CYTHON_URL "http://cython.org/release/Cython-0.17.1.tar.gz" CACHE STRING "URL to download Cython from")
mark_as_advanced(CYTHON_URL)
ENDIF(INTERNAL_CYTHON)
IF(INTERNAL_NETCDF4_PYTHON)
SET(NETCDF4_PYTHON_URL "http://netcdf4-python.googlecode.com/files/netCDF4-1.0.1.tar.gz" CACHE STRING "URL to download NetCDF4-python from")
mark_as_advanced(NETCDF4_PYTHON_URL)
ENDIF(INTERNAL_NETCDF4_PYTHON)
IF (INTERNAL_HEALPY)
SET(HEALPY_URL "http://github.com/healpy/healpy/archive/1.4.1.tar.gz" CACHE STRING "URL to download Healpy from")
mark_as_advanced(HEALPY_URL)
ENDIF(INTERNAL_HEALPY)
IF(INTERNAL_ARGPARSE)
SET(SETUPTOOLS_URL "http://pypi.python.org/packages/source/s/setuptools/setuptools-0.6c11.tar.gz" CACHE STRING "URL to download setuptools from")
SET(ARGPARSE_URL "http://argparse.googlecode.com/files/argparse-1.2.1.tar.gz" CACHE STRING "URL to download argparse from")
mark_as_advanced(ARGPARSE_URL SETUPTOOLS_URL)
ENDIF(INTERNAL_ARGPARSE)
execute_process(
COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_SOURCE_DIR}/external/detect_site.py ${CMAKE_BINARY_DIR}/ext_build/python
RESULT_VARIABLE RET_VALUE
@ -67,6 +87,57 @@ IF(INTERNAL_NETCDF4_PYTHON)
BUILD_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_build.cmake
INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake
)
SET(PREV_PYTHON_BUILD ${PREV_PYTHON_BUILD} netcdf4-python)
ENDIF(INTERNAL_NETCDF4_PYTHON)
IF(INTERNAL_HEALPY)
SET(BUILD_ENVIRONMENT
${CMAKE_COMMAND}
"-DPYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}"
"-DPYTHON_CPPFLAGS:STRING=${PYTHON_CPPFLAGS}"
"-DCFITSIO_EXT_LIB=${CFITSIO_LIBRARY}"
"-DCFITSIO_EXT_INC=${CFITSIO_INCLUDE_PATH}"
"-DCFITSIO_EXT_PREFIX=${CFITSIO_PREFIX}"
"-DNETCDF4_DIR=${NETCDF_BIN_DIR}"
"-DPYTHON_LDFLAGS:STRING=${PYTHON_LDFLAGS}"
"-DPYTHON_LOCAL_SITE_PACKAGE=${PYTHON_LOCAL_SITE_PACKAGE}"
"-DTARGET_PATH=${CMAKE_BINARY_DIR}/ext_build/python" "-P")
ExternalProject_Add(healpy
DEPENDS ${PREV_PYTHON_BUILD}
URL ${HEALPY_URL}
PREFIX ${BUILD_PREFIX}/healpy-prefix
CONFIGURE_COMMAND echo "No configure"
BUILD_IN_SOURCE 1
BUILD_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_build.cmake
INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake
)
ENDIF(INTERNAL_HEALPY)
IF(INTERNAL_ARGPARSE)
SET(BUILD_ENVIRONMENT
${CMAKE_COMMAND}
"-DPYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}"
"-DPYTHON_LOCAL_SITE_PACKAGE=${PYTHON_LOCAL_SITE_PACKAGE}"
"-DTARGET_PATH=${CMAKE_BINARY_DIR}/ext_build/python" "-P")
ExternalProject_Add(setuptools
URL ${SETUPTOOLS_URL}
PREFIX ${BUILD_PREFIX}/setuptools-prefix
CONFIGURE_COMMAND echo "No configure"
BUILD_IN_SOURCE 1
BUILD_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_build.cmake
INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake
)
ExternalProject_Add(argparse
DEPENDS setuptools
URL ${ARGPARSE_URL}
PREFIX ${BUILD_PREFIX}/argparse-prefix
CONFIGURE_COMMAND echo "No configure"
BUILD_IN_SOURCE 1
BUILD_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_build.cmake
INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake
)
SET(AUXILIARY_PYTHON_DEPEND ${AUXILIARY_PYTHON_DEPEND} argparse)
ENDIF(INTERNAL_ARGPARSE)

4
external/python_build.cmake vendored
View file

@ -4,6 +4,10 @@ SET(ENV{CPPFLAGS} ${PYTHON_CPPFLAGS})
SET(ENV{LDFLAGS} ${PYTHON_LDFLAGS})
SET(ENV{VOID_GSL} ${VOID_GSL})
SET(ENV{PYTHONPATH} ${PYTHON_LOCAL_SITE_PACKAGE}:$ENV{PYTHONPATH})
SET(ENV{CFITSIO_EXT_INC} ${CFITSIO_EXT_INC})
SET(ENV{CFITSIO_EXT_LIB} ${CFITSIO_EXT_LIB})
SET(ENV{CFITSIO_EXT_PREFIX} ${CFITSIO_EXT_PREFIX})
SET(PYTHON_BUILD_COMMAND ${PYTHON_EXECUTABLE} setup.py build)
MESSAGE(STATUS "Running ${PYTHON_BUILD_COMMAND}")
execute_process(

5
external/python_install.cmake vendored
View file

@ -3,8 +3,13 @@ SET(ENV{NETCDF4_DIR} ${NETCDF4_DIR})
SET(ENV{CPPFLAGS} ${PYTHON_CPPFLAGS})
SET(ENV{LDFLAGS} ${PYTHON_LDFLAGS})
SET(ENV{VOID_GSL} ${VOID_GSL})
SET(ENV{CFITSIO_EXT_INC} ${CFITSIO_EXT_INC})
SET(ENV{CFITSIO_EXT_PREFIX} ${CFITSIO_EXT_PREFIX})
SET(ENV{CFITSIO_EXT_LIB} ${CFITSIO_EXT_LIB})
SET(ENV{PYTHONPATH} ${PYTHON_LOCAL_SITE_PACKAGE}:$ENV{PYTHONPATH})
SET(PYTHON_INSTALL_COMMAND ${PYTHON_EXECUTABLE} setup.py install --prefix=${TARGET_PATH} --install-lib=${PYTHON_LOCAL_SITE_PACKAGE})
message(STATUS "Running ${PYTHON_INSTALL_COMMAND}")
execute_process(
COMMAND ${PYTHON_INSTALL_COMMAND}

321
pipeline/applyMaskToMock.py Executable file
View file

@ -0,0 +1,321 @@
#!/usr/bin/env python
# prepares input catalogs based on multidark simulations
# (borrows heavily from generateMock, but doesn't hold much in memory)
# also creates necessary analyzeVoids input files
import numpy as np
import os
import sys
import void_python_tools as vp
import argparse
import imp
import healpy as hp
# ------------------------------------------------------------------------------
def my_import(name):
mod = __import__(name)
components = name.split('.')
for comp in components[1:]:
mod = getattr(mod, comp)
return mod
# -----------------------------------------------------------------------------
LIGHT_SPEED = 299792.458
parser = argparse.ArgumentParser(description='options')
parser.add_argument('--scripts', dest='script', action='store_const',
const=True, default=False,
help='write scripts')
parser.add_argument('--parmFile', dest='parmFile',
default="",
help='path to parameter file')
args = parser.parse_args()
filename = args.parmFile
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
#------------------------------------------------------------------------------
def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1):
sampleName = setName
sampleName += "_z" + redshift
if iVol != -1: sampleName += "_d" + iVol
return sampleName
#------------------------------------------------------------------------------
# for given dataset parameters, outputs a script for use with analyzeVoids
def writeScript(setName, dataFileNameBase,
scriptDir, catalogDir, fileNums, redshifts, numSubvolumes,
numSlices, useVel, lbox, minRadius, omegaM, subsample=1.0,
suffix=".dat"):
if useVel: setName += "_pv"
scriptFileName = scriptDir + "/" + setName + ".py"
scriptFile = open(scriptFileName, 'w')
scriptFile.write("""#!/usr/bin/env/python
import os
from void_python_tools.backend.classes import *
continueRun = True # set to True to enable restarting aborted jobs
startCatalogStage = 1
endCatalogStage = 4
startAPStage = 1
endAPStage = 7
ZOBOV_PATH = os.getenv("PWD")+"/../zobov/"
CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/"
freshStack = True
errorBars = "CALCULATED"
numIncoherentRuns = 100
ranSeed = 101010
useLCDM = False
bias = 1.16
dataPortions = ["central"]
dataSampleList = []
""")
dataInfo = """
setName = "{setName}"
workDir = "{voidOutputDir}/{setName}/"
inputDataDir = "{inputDataDir}"
figDir = "{figDir}/{setName}/"
logDir = "{logDir}/{setName}/"
numZobovDivisions = {numZobovDivisions}
numZobovThreads = {numZobovThreads}
"""
scriptFile.write(dataInfo.format(setName=setName, figDir=figDir,
logDir=logDir, voidOutputDir=voidOutputDir,
inputDataDir=catalogDir,
numZobovDivisions=numZobovDivisions,
numZobovThreads=numZobovThreads))
sampleInfo = """
newSample = Sample(dataFile = "{dataFile}",
dataFormat = "{dataFormat}",
dataUnit = {dataUnit},
fullName = "{sampleName}",
nickName = "{sampleName}",
dataType = "simulation",
zBoundary = ({zMin}, {zMax}),
zRange = ({zMin}, {zMax}),
zBoundaryMpc = ({zMinMpc}, {zMaxMpc}),
omegaM = {omegaM},
minVoidRadius = {minRadius},
includeInHubble = True,
partOfCombo = False,
isCombo = False,
boxLen = {boxLen},
usePecVel = {usePecVel},
numSubvolumes = {numSubvolumes},
mySubvolume = "{mySubvolume}",
useLightCone = {useLightCone},
subsample = {subsample})
dataSampleList.append(newSample)
newSample.addStack({zMin}, {zMax}, {minRadius} , {minRadius}+2, True, False)
newSample.addStack({zMin}, {zMax}, {minRadius} , {minRadius}+4, True, False)
newSample.addStack({zMin}, {zMax}, {minRadius}+2, {minRadius}+6, True, False)
newSample.addStack({zMin}, {zMax}, {minRadius}+6, {minRadius}+10, True, False)
newSample.addStack({zMin}, {zMax}, {minRadius}+10, {minRadius}+18, True, False)
newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False)
"""
for (iFile, redshift) in enumerate(redshifts):
fileNum = fileNums[iFile]
zBox = float(redshift)
Om = float(omegaM)
zBoxMpc = LIGHT_SPEED/100.*vp.angularDiameter(zBox, Om=Om)
boxMaxMpc = zBoxMpc + lbox
# converter from redshift to comoving distance
zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000)
zVsDX = np.zeros(len(zVsDY))
for i in xrange(len(zVsDY)):
zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
if useLightCone:
boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \
LIGHT_SPEED*lbox, zVsDX, zVsDY)-zBox
dzSafe = 0.03
else:
boxWidthZ = lbox*100./LIGHT_SPEED
dzSafe = 0.0
for iSlice in xrange(numSlices):
sliceMin = zBox + dzSafe + iSlice*(boxWidthZ-dzSafe)/numSlices
sliceMax = zBox + dzSafe + (iSlice+1)*(boxWidthZ-dzSafe)/numSlices
sliceMinMpc = sliceMin*LIGHT_SPEED/100.
sliceMaxMpc = sliceMax*LIGHT_SPEED/100.
sliceMin = "%0.2f" % sliceMin
sliceMax = "%0.2f" % sliceMax
sliceMinMpc = "%0.1f" % sliceMinMpc
sliceMaxMpc = "%0.1f" % sliceMaxMpc
dataFileName = dataFileNameBase + fileNum + suffix
for iX in xrange(numSubvolumes):
for iY in xrange(numSubvolumes):
mySubvolume = "%d%d" % (iX, iY)
sampleName = getSampleName(setName, sliceMin, useVel,
iSlice=iSlice, iVol=mySubvolume)
scriptFile.write(sampleInfo.format(dataFile=dataFileName,
dataFormat=dataFormat,
dataUnit=dataUnit,
sampleName=sampleName,
zMin=sliceMin,
zMax=sliceMax,
zMinMpc=sliceMinMpc,
zMaxMpc=sliceMaxMpc,
omegaM=Om,
boxLen=lbox,
usePecVel=useVel,
minRadius=minRadius,
numSubvolumes=numSubvolumes,
mySubvolume=mySubvolume,
useLightCone=useLightCone,
subsample=subsample))
scriptFile.close()
return
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir)
if not os.access(catalogDir, os.F_OK): os.mkdir(catalogDir)
# -----------------------------------------------------------------------------
# now the SDSS HOD
parFileText = """
% cosmology
OMEGA_M {omegaM}
HUBBLE {hubble}
OMEGA_B 0.0469
SIGMA_8 0.82
SPECTRAL_INDX 0.95
ITRANS 5
REDSHIFT {redshift}
% halo definition
%DELTA_HALO 200
DELTA_HALO 740.74 % 200/Om_m
M_max 1.00E+16
% fit function types
pdfs 11
pdfc 2
EXCLUSION 4
% hod parameters
M_min {Mmin}
GALAXY_DENSITY 0.0111134 % computed automatically if M_min set, use for sanity
M1 {M1}
sigma_logM {sigma_logM}
alpha {alpha}
M_cut {Mcut}
% simulation info
real_space_xi 1
HOD 1
populate_sim 1
HaloFile {haloFile}
RESOLUTION {numPartPerSide}
BOX_SIZE {boxSize}
% output
root_filename hod
"""
print " Doing DR9 HOD"
# these parameters come from Manera et al 2012, eq. 26
parFileName = "./hod.par"
parFile = open(parFileName, 'w')
haloFile = catalogDir+haloFileBase+fileNums[iRedshift]
parFile.write(parFileText.format(omegaM=omegaM,
hubble=hubble,
redshift=redshift,
Mmin=1.23e13,
M1=1.e14,
sigma_logM=0.596,
alpha=1.0127,
Mcut=1.19399e13,
haloFile=haloFile,
numPartPerSide=numPart**(1/3.),
boxSize=lbox))
parFile.close()
os.system(hodPath+" "+parFileName+">& /dev/null")
# now place these particles on a lightcone, restrict redshift range, apply mask
mask = hp.read_map(maskFile)
nside = hp.get_nside(mask)
inFile = open('hod.mock', 'r')
outFile = open(catalogDir+"/mock.out"))
zBox = float(redshiftRange[0])
Om = float(omegaM)
# converter from redshift to comoving distance
zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000)
zVsDX = np.zeros(len(zVsDY))
for i in xrange(len(zVsDY)):
zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
for line in inFile:
line = line.split(',')
x = float(line[0]) - lbox/2.
y = float(line[1]) - lbox/2.
z = float(line[2]) - lbox/2.
vz = float(line[5])
zBoxInMpc = vp.angularDiameter(zBox, Om=Om)
redshift = np.sqrt(x*x + y*y + z*z)
redshift = np.interp(zBoxInMpc+100./LIGHT_SPEED*redshift, zVsDX, zVsDY)
if redshift < redshiftRange[0] or redshift > redshiftRange[1]: continue
RA = np.atan((y-lbox/2.)/(x-lbox/2.)) * 100/np.pi + 180.
Dec = np.asin((z-lboc/2.)/(redshift*LIGHT_SPEED/100.)) * 180/np.pi
phi = np.pi/180. * RA
theta = np.pi/2. - Dec*np.pi/180.
pos = np.zeros((3))
pix = hp.ang2pix(nside, theta, phi)
if mask[pix] <= 0: continue
print >> outFile, RA, Dec, redshift*LIGHT_SPEED, 0., x, y, z
inFile.close()
outFile.close()
os.system("rm ./hod.*")

62
pipeline/datasets/mock_dr9mid.py Normal file
View file

@ -0,0 +1,62 @@
import os
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# CONFIGURATION
# directory for the input simulation/observational particle files
catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/mock_dr9mid/"
# path to HOD code
hodPath = os.getenv("HOME")+"/projects/Voids/hod/HOD.x"
# path to mask
maskFile = os.getenv("HOME")+"/workspace/Voids/catalogs/boss/final_boss_mask.fits")
# where to put the final void catalog, figures, and output logs
voidOutputDir = os.getenv("HOME")+"/workspace/Voids/mock_dr9mid/"
figDir = os.getenv("PWD")+"/../figs/mock_dr9mid/"
logDir = os.getenv("PWD")+"/../logs/mock_dr9mid/"
# where to place the pipeline scripts
scriptDir = os.getenv("PWD")+"/mock_dr9mid/"
# simulation or observation?
dataType = "observation"
# available formats for simulation: gadget, multidark
dataFormat = "multidark"
dataUnit = 1 # as multiple of Mpc/h
# place particles on the lightcone?
useLightCone = True
# list of file numbers for the particle files
# to get particle file name, we take particleFileBase+fileNum
fileNums = (("0.53"))
# redshift range of the mock
redshiftRange = (0.53, 0.6)
# prefix to give all outputs
prefix = "mock_"
# common filename of halo files
haloFileBase = "mdr1_halos_z"
# adjust these two parameters given the memory contraints on your system:
# numZobovDivisions: how many sub-volumes per dimension will zobov process
# numZobovThreads: how many sub-volumes to process at once?
numZobovDivisions = 2
numZobovThreads = 2
# simulation information
numPart = 1024*1024*1024
lbox = 1000 # Mpc/h
omegaM = 0.27
hubble = 0.70
# END CONFIGURATION
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------

76
pipeline/datasets/multidark.py Normal file
View file

@ -0,0 +1,76 @@
import os
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# CONFIGURATION
# directory for the input simulation/observational particle files
catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/multidark/"
# path to HOD code
hodPath = os.getenv("HOME")+"/projects/Voids/hod/HOD.x"
# where to put the final void catalog, figures, and output logs
voidOutputDir = os.getenv("HOME")+"/workspace/Voids/multidark/"
figDir = os.getenv("PWD")+"/../figs/multidark/"
logDir = os.getenv("PWD")+"/../logs/multidark/"
# where to place the pipeline scripts
scriptDir = os.getenv("PWD")+"/multidark/"
# simulation or observation?
dataType = "simulation"
# available formats for simulation: gadget, multidark
dataFormat = "multidark"
dataUnit = 1 # as multiple of Mpc/h
# place particles on the lightcone?
useLightCone = True
# common filename of particle files
particleFileBase = "mdr1_particles_z"
# list of file numbers for the particle files
# to get particle file name, we take particleFileBase+fileNum
fileNums = (("0.0", "0.53", "1.0"))
# redshift of each file in the above list
redshifts = (("0.0", "0.53", "1.0"))
# how many independent slices along the z-axis?
numSlices = 4
# how many subdivisions along the x- and y- axis?
# ( = 2 will make 4 subvolumes for each slice, = 3 will make 9, etc.)
numSubvolumes = 1
# prefix to give all outputs
prefix = "md_"
# list of desired subsamples - these are in unts of h Mpc^-3!
#subSamples = [ 1.0 ]
subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002))
# common filename of halo files
haloFileBase = "mdr1_halos_z"
# minimum halo mass cuts to apply for the halo catalog
# use "none" to get all halos
minHaloMasses = (("none", 2e12, 1.23e13))
# adjust these two parameters given the memory contraints on your system:
# numZobovDivisions: how many sub-volumes per dimension will zobov process
# numZobovThreads: how many sub-volumes to process at once?
numZobovDivisions = 4
numZobovThreads = 2
# simulation information
numPart = 1024*1024*1024
lbox = 1000 # Mpc/h
omegaM = 0.27
hubble = 0.70
# END CONFIGURATION
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------

35
pipeline/generateCatalog.py
View file

@ -5,6 +5,7 @@
from void_python_tools.backend import *
from void_python_tools.plotting import *
import imp
import pickle
# ------------------------------------------------------------------------------
@ -23,9 +24,8 @@ if (len(sys.argv) > 1):
filename = sys.argv[1]
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file!"
print " Cannot find parameter file %s!" % filename
exit(-1)
#parms = __import__(filename[:-3], globals(), locals(), ['*'])
parms = imp.load_source("name", filename)
globals().update(vars(parms))
else:
@ -61,6 +61,27 @@ for sample in dataSampleList:
if not os.access(zobovDir, os.F_OK):
os.makedirs(zobovDir)
# save this sample's information
with open(zobovDir+"/sample_info.dat", 'w') as output:
pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL)
fp = open(zobovDir+"/sample_info.txt", 'w')
fp.write("Sample name: %s\n" % sample.fullName)
fp.write("Sample nickname: %s\n" % sample.nickName)
fp.write("Data type: %s\n" % sample.dataType)
fp.write("Redshift range: %f - %f\n" %(sample.zBoundary[0],sample.zBoundary[1]))
fp.write("Estimated mean particle separation: %g\n" % sample.minVoidRadius)
if (sample.dataType == "simulation"):
fp.write("Particles placed on lightcone: %g\n" % sample.useLightCone)
fp.write("Peculiar velocities included: %g\n" % sample.usePecVel)
fp.write("Additional subsampling fraction: %g\n" % sample.subsample)
fp.write("Simulation box length (Mpc/h): %g\n" % sample.boxLen)
fp.write("Simulation Omega_M: %g\n" % sample.omegaM)
fp.write("Number of simulation subvolumes: %g\n", sample.numSubvolumes)
fp.write("My subvolume index: %g\n", sample.mySubvolume)
fp.close()
# ---------------------------------------------------------------------------
if (startCatalogStage <= 1) and (endCatalogStage >= 1) and not sample.isCombo:
print " Extracting tracers from catalog...",
@ -84,7 +105,8 @@ for sample in dataSampleList:
sys.stdout.flush()
launchZobov(sample, ZOBOV_PATH, zobovDir=zobovDir, logDir=logDir,
continueRun=continueRun)
continueRun=continueRun, numZobovDivisions=numZobovDivisions,
numZobovThreads=numZobovThreads)
# -------------------------------------------------------------------------
if (startCatalogStage <= 3) and (endCatalogStage >= 3) and not sample.isCombo:
@ -108,6 +130,11 @@ if (startCatalogStage <= 4) and (endCatalogStage >= 4):
sys.stdout.flush()
for thisDataPortion in dataPortions:
plotNumberCounts(workDir, dataSampleList, figDir, showPlot=True, dataPortion=thisDataPortion, setName=setName)
plotRedshiftDistribution(workDir, dataSampleList, figDir, showPlot=False,
dataPortion=thisDataPortion, setName=setName)
plotSizeDistribution(workDir, dataSampleList, figDir, showPlot=False,
dataPortion=thisDataPortion, setName=setName)
plotNumberDistribution(workDir, dataSampleList, figDir, showPlot=False,
dataPortion=thisDataPortion, setName=setName)
print "\n Done!"

341
pipeline/prepareCatalogs.py
View file

@ -9,33 +9,19 @@ import os
import sys
import void_python_tools as vp
import argparse
import imp
catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/multidark/"
scriptDir = os.getenv("PWD")+"/multidark/"
hodPath = os.getenv("HOME")+"/projects/Voids/hod/HOD.x"
# ------------------------------------------------------------------------------
outputDir = os.getenv("HOME")+"/workspace/Voids/multidark/"
figDir = os.getenv("PWD")+"/../figs/multidark/"
logDir = os.getenv("PWD")+"/../logs/multidark/"
dataFormat = "multidark"
dataType = "simulation"
useLightCone = True
redshifts = (("0.0", "0.53", "1.0"))
redshiftFileBase = "mdr1_particles_z"
numSlices = 4
numSubvolumes = 1
prefix = "md_"
subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002))
numPart = 1024*1024*1024 # number of particles in base catalog
lbox = 1000
omegaM = 0.27
hubble = 0.70
def my_import(name):
mod = __import__(name)
components = name.split('.')
for comp in components[1:]:
mod = getattr(mod, comp)
return mod
# -----------------------------------------------------------------------------
LIGHT_SPEED = 299792.458
parser = argparse.ArgumentParser(description='options')
@ -50,54 +36,59 @@ parser.add_argument('--halos', dest='halos', action='store_const',
help='write halos')
parser.add_argument('--hod', dest='hod', action='store_const',
const=True, default=False,
help='write hos')
help='write hod')
parser.add_argument('--all', dest='all', action='store_const',
const=True, default=False,
help='write everything')
parser.add_argument('--parmFile', dest='parmFile',
default="",
help='path to parameter file')
args = parser.parse_args()
filename = args.parmFile
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
#------------------------------------------------------------------------------
def getSampleName(prefix, base, redshift, useVel, iSlice=-1, iVol=-1):
useVelString = ""
if useVel: useVelString = "_pv"
def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1):
sampleName = prefix + base + useVelString + "_z" + redshift
sampleName = setName
if iSlice != -1:
sampleName += "_s" + str(iSlice)
sampleName += "_z" + redshift
if iVol != -1:
sampleName += "_d" + iVol
if iVol != -1: sampleName += "_d" + iVol
return sampleName
#------------------------------------------------------------------------------
# for given dataset parameters, outputs a script for use with analyzeVoids
def writeScript(prefix, base, scriptDir, catalogDir, redshifts, numSubvolumes,
def writeScript(setName, dataFileNameBase,
scriptDir, catalogDir, fileNums, redshifts, numSubvolumes,
numSlices, useVel, lbox, minRadius, omegaM, subsample=1.0,
suffix=".dat"):
setName = prefix + base
dataFileNameBase = setName #+ ".dat"
if useVel: setName += "_pv"
scriptFileName = scriptDir + "/" + "md_" + base
if useVel: scriptFileName += "_pv"
scriptFileName += ".py"
scriptFileName = scriptDir + "/" + setName + ".py"
scriptFile = open(scriptFileName, 'w')
scriptFile.write("""#!/usr/bin/env/python
import os
from void_python_tools.backend.classes import *
continueRun = False
continueRun = True # set to True to enable restarting aborted jobs
startCatalogStage = 1
endCatalogStage = 3
endCatalogStage = 4
startAPStage = 1
endAPStage = 6
endAPStage = 7
ZOBOV_PATH = os.getenv("PWD")+"/../zobov/"
CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/"
@ -112,22 +103,28 @@ dataPortions = ["central"]
dataSampleList = []
""")
dataInfo = """
catalogName = "{sampleName}"
workDir = "{outputDir}/{sampleName}/"
inputDataDir = "{catalogDir}"
figDir = "{figDir}/{sampleName}/"
logDir = "{logDir}/{sampleName}/"
dataInfo = """
setName = "{setName}"
workDir = "{voidOutputDir}/{setName}/"
inputDataDir = "{inputDataDir}"
figDir = "{figDir}/{setName}/"
logDir = "{logDir}/{setName}/"
numZobovDivisions = {numZobovDivisions}
numZobovThreads = {numZobovThreads}
"""
scriptFile.write(dataInfo.format(sampleName=setName, figDir=figDir,
logDir=logDir, outputDir=outputDir,
catalogDir=catalogDir,
inputDataDir=catalogDir))
scriptFile.write(dataInfo.format(setName=setName, figDir=figDir,
logDir=logDir, voidOutputDir=voidOutputDir,
inputDataDir=catalogDir,
numZobovDivisions=numZobovDivisions,
numZobovThreads=numZobovThreads))
sampleInfo = """
newSample = Sample(dataFile = "{dataFile}",
dataFormat = "{dataFormat}",
dataUnit = {dataUnit},
fullName = "{sampleName}",
nickName = "{sampleName}",
dataType = "simulation",
@ -143,7 +140,6 @@ newSample = Sample(dataFile = "{dataFile}",
usePecVel = {usePecVel},
numSubvolumes = {numSubvolumes},
mySubvolume = "{mySubvolume}",
numSubDivisions = 4,
useLightCone = {useLightCone},
subsample = {subsample})
dataSampleList.append(newSample)
@ -154,8 +150,8 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+6, {minRadius}+10, True, False)
newSample.addStack({zMin}, {zMax}, {minRadius}+10, {minRadius}+18, True, False)
newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False)
"""
for redshift in redshifts:
for (iFile, redshift) in enumerate(redshifts):
fileNum = fileNums[iFile]
zBox = float(redshift)
Om = float(omegaM)
zBoxMpc = LIGHT_SPEED/100.*vp.angularDiameter(zBox, Om=Om)
@ -176,7 +172,6 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False)
dzSafe = 0.0
for iSlice in xrange(numSlices):
# trim off the first and last 10,000 km/s
sliceMin = zBox + dzSafe + iSlice*(boxWidthZ-dzSafe)/numSlices
sliceMax = zBox + dzSafe + (iSlice+1)*(boxWidthZ-dzSafe)/numSlices
@ -188,18 +183,19 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False)
sliceMinMpc = "%0.1f" % sliceMinMpc
sliceMaxMpc = "%0.1f" % sliceMaxMpc
dataFileName = dataFileNameBase + "_z" + redshift + suffix
dataFileName = dataFileNameBase + fileNum + suffix
for iX in xrange(numSubvolumes):
for iY in xrange(numSubvolumes):
mySubvolume = "%d%d" % (iX, iY)
sampleName = getSampleName(prefix, base, redshift, useVel,
sampleName = getSampleName(setName, sliceMin, useVel,
iSlice=iSlice, iVol=mySubvolume)
scriptFile.write(sampleInfo.format(dataFile=dataFileName,
dataFormat=dataFormat,
dataUnit=dataUnit,
sampleName=sampleName,
zMin=sliceMin,
zMax=sliceMax,
@ -221,12 +217,13 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False)
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir)
if not os.access(catalogDir, os.F_OK): os.mkdir(catalogDir)
#------------------------------------------------------------------------------
# first the directly downsampled runs
# Note: ss0.002 ~ SDSS DR7 dim2
# ss0.0004 ~ SDSS DR9 mid
baseResolution = numPart/lbox/lbox/lbox # particles/Mpc^3
baseResolution = float(numPart)/lbox/lbox/lbox # particles/Mpc^3
for thisSubSample in subSamples:
keepFraction = float(thisSubSample) / baseResolution
@ -235,110 +232,161 @@ for thisSubSample in subSamples:
if args.script or args.all:
print " Doing subsample", thisSubSample, " scripts"
setName = prefix+"ss"+str(thisSubSample)
if dataFormat == "multidark":
writeScript(prefix,
"ss"+str(thisSubSample), scriptDir, catalogDir, redshifts,
writeScript(setName, "md.ss"+str(thisSubSample)+"_z",
scriptDir, catalogDir, fileNums,
redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
subsample=1.0)
writeScript(prefix, "ss"+str(thisSubSample), scriptDir, catalogDir,
writeScript(setName, "md.ss"+str(thisSubSample)+"_z",
scriptDir, catalogDir,
fileNums,
redshifts, numSubvolumes, numSlices, True, lbox, minRadius,
omegaM, subsample=1.0)
elif dataFormat == "gadget":
writeScript("", redshiftFileBase, scriptDir, catalogDir, redshifts,
writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums,
redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
subsample=thisSubSample, suffix="")
writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums,
redshifts,
numSubvolumes, numSlices, True, lbox, minRadius, omegaM,
subsample=thisSubSample, suffix="")
elif dataFormat == "random":
writeScript(setName, "ran.ss"+str(thisSubSample)+"_z",
scriptDir, catalogDir, fileNums,
redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
subsample=1.0)
if args.subsample or args.all:
print " Doing subsample", thisSubSample
for redshift in redshifts:
for (iRedshift, redshift) in enumerate(redshifts):
print " redshift", redshift
dataFile = catalogDir+"/"+redshiftFileBase+redshift
inFile = open(dataFile, 'r')
sampleName = getSampleName("ss"+str(thisSubSample), redshift, False)
outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
if dataFormat == "multidark":
dataFile = catalogDir+"/"+particleFileBase+fileNums[iRedshift]
inFile = open(dataFile, 'r')
outFile.write("%f\n" %(lbox))
outFile.write("%s" %(omegaM))
outFile.write("%s" %(hubble))
outFile.write("%s\n" %(redshift))
outFile.write("%d\n" %(maxKeep))
sampleName = "md.ss"+str(thisSubSample)+"_z"+redshift
outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
numKept = 0
for line in inFile:
if np.random.uniform() > keepFraction: continue
numKept += 1
if numKept > maxKeep: break
line = line.split(',')
x = float(line[0])
y = float(line[1])
z = float(line[2])
vz = float(line[3])
outFile.write("%f\n" %(lbox))
outFile.write("%s\n" %(omegaM))
outFile.write("%s\n" %(hubble))
outFile.write("%s\n" %(redshift))
outFile.write("%d\n" %(maxKeep))
outFile.write("%e %e %e %e\n" %(x,y,z,vz))
numKept = 0
for (i,line) in enumerate(inFile):
if np.random.uniform() > keepFraction: continue
numKept += 1
if numKept > maxKeep: break
line = line.split(',')
x = float(line[0])
y = float(line[1])
z = float(line[2])
vz = float(line[3])
outFile.write("%d %e %e %e %e\n" %(i,x,y,z,vz))
outFile.write("-99 -99 -99 -99 -99\n")
inFile.close()
outFile.close()
elif dataFormat == "random":
sampleName = "ran.ss"+str(thisSubSample)+"_z"+redshift
outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
outFile.write("%f\n" %(lbox))
outFile.write("%s\n" %(omegaM))
outFile.write("%s\n" %(hubble))
outFile.write("%s\n" %(redshift))
outFile.write("%d\n" %(maxKeep))
for i in xrange(int(maxKeep)):
x = np.random.uniform()*lbox
y = np.random.uniform()*lbox
z = np.random.uniform()*lbox
outFile.write("%d %e %e %e %e\n" % (i, x,y,z, 0.))
outFile.write("-99 -99 -99 -99 -99\n")
outFile.close()
outFile.write("-99 -99 -99 -99\n")
print "KEEPING:", numKept, "...predicted:", maxKeep
inFile.close()
outFile.close()
# -----------------------------------------------------------------------------
# now halos
if args.script or args.all:
if (args.script or args.all) and dataFormat == "multidark":
print " Doing halo scripts"
dataFile = catalogDir+"/mdr1_halos_z"+redshifts[0]
inFile = open(dataFile, 'r')
numPart = 0
for line in inFile: numPart += 1
inFile.close()
for minHaloMass in minHaloMasses:
# estimate number of halos to get density
dataFile = catalogDir+haloFileBase+fileNums[0]
inFile = open(dataFile, 'r')
numPart = 0
for line in inFile:
line = line.split(',')
if minHaloMass == "none" or float(line[6]) > minHaloMass:
numPart += 1
inFile.close()
minRadius = int(np.ceil(lbox/numPart**(1./3.)))
minRadius = 2*int(np.ceil(lbox/numPart**(1./3.)))
if dataFormat == "multidark":
writeScript(prefix, "halos", scriptDir, catalogDir, redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM)
writeScript(prefix, "halos", scriptDir, catalogDir, redshifts, numSubvolumes,
numSlices, True, lbox, minRadius, omegaM)
if dataFormat == "multidark":
setName = prefix+"halos_min"+str(minHaloMass)
writeScript(setName, "md.halos_min"+str(minHaloMass)+"_z",
scriptDir, catalogDir, fileNums,
redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM)
writeScript(setName, "md.halos_min"+str(minHaloMass)+"_z",
scriptDir, catalogDir, fileNums,
redshifts,
numSubvolumes, numSlices, True, lbox, minRadius, omegaM)
if args.halos or args.all:
print " Doing halos"
for redshift in redshifts:
print " z = ", redshift
for minHaloMass in minHaloMasses:
print " min halo mass = ", minHaloMass
dataFile = catalogDir+"/mdr1_halos_z"+redshift
inFile = open(dataFile, 'r')
numPart = 0
for line in inFile: numPart += 1
inFile.close()
for (iRedshift, redshift) in enumerate(redshifts):
print " z = ", redshift
sampleName = getSampleName("halos", redshift, False)
dataFile = catalogDir+haloFileBase+fileNums[iRedshift]
inFile = open(dataFile, 'r')
numPart = 0
for line in inFile:
line = line.split(',')
if minHaloMass == "none" or float(line[6]) > minHaloMass:
numPart += 1
inFile.close()
outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
sampleName = "md.halos_z"+redshift
outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
outFile.write("%f\n" %(lbox))
outFile.write("%s" %(omegaM))
outFile.write("%s" %(hubble))
outFile.write("%s\n" %(redshift))
outFile.write("%d\n" %(numPart))
outFile.write("%f\n" %(lbox))
outFile.write("%s\n" %(omegaM))
outFile.write("%s\n" %(hubble))
outFile.write("%s\n" %(redshift))
outFile.write("%d\n" %(numPart))
inFile = open(dataFile, 'r')
numKept = 0
for line in inFile:
numKept += 1
line = line.split(',')
x = float(line[0])
y = float(line[1])
z = float(line[2])
vz = float(line[5])
inFile = open(dataFile, 'r')
for (iHalo,line) in enumerate(inFile):
line = line.split(',')
if minHaloMass == "none" or float(line[6]) > minHaloMass:
x = float(line[0])
y = float(line[1])
z = float(line[2])
vz = float(line[5])
# write to output file
outFile.write("%e %e %e %e\n" %(x,y,z,vz))
# write to output file
outFile.write("%d %e %e %e %e\n" %(iHalo,x,y,z,vz))
inFile.close()
outFile.close()
inFile.close()
outFile.close()
# -----------------------------------------------------------------------------
# now the SDSS HOD
@ -382,24 +430,25 @@ BOX_SIZE {boxSize}
root_filename hod
"""
if args.script or args.all:
if (args.script or args.all) and dataFormat == "multidark":
print " Doing DR7 HOD scripts"
if dataFormat == "multidark":
writeScript(prefix,
"hod_dr72dim2", scriptDir, catalogDir, redshifts,
setName = prefix+"hod_dr72dim2"
writeScript(setName, "md.hod_dr72dim2_z",
scriptDir, catalogDir, fileNums, redshifts,
numSubvolumes, numSlices, False, lbox, 5, omegaM)
writeScript(prefix,
"hod_dr72dim2", scriptDir, catalogDir, redshifts,
writeScript(setName, "md.hod_dr72dim2_z",
scriptDir, catalogDir, fileNums, redshifts,
numSubvolumes, numSlices, True, lbox, 5, omegaM)
if args.hod or args.all:
print " Doing DR7 HOD"
for redshift in redshifts:
for (iRedshift, redshift) in enumerate(redshifts):
print " z = ", redshift
parFileName = "./hod.par"
parFile = open(parFileName, 'w')
haloFile = catalogDir+"/mdr1_halos_z"+redshift
haloFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift]
parFile.write(parFileText.format(omegaM=omegaM,
hubble=hubble,
redshift=redshift,
@ -415,7 +464,7 @@ if args.hod or args.all:
os.system(hodPath+" "+parFileName+">& /dev/null")
sampleName = getSampleName("hod_dr72dim2", redshift, False)
sampleName = getSampleName("md.hod_dr72dim2", redshift, False)
outFileName = catalogDir+"/"+sampleName+".dat"
os.system("mv hod.mock" + " " + outFileName)
@ -423,24 +472,26 @@ if args.hod or args.all:
# -----------------------------------------------------------------------------
# now the BOSS HOD
if args.script or args.all:
if (args.script or args.all) and dataFormat == "multidark":
print " Doing DR9 HOD scripts"
if dataFormat == "multidark":
writeScript(prefix,
"hod_dr9mid", scriptDir, catalogDir, redshifts,
numSubvolumes, numSlices, False, lbox, 5, omegaM)
writeScript(prefix,
"hod_dr9mid", scriptDir, catalogDir, redshifts,
numSubvolumes, numSlices, True, lbox, 5, omegaM)
setName = prefix+"hod_dr9mid"
writeScript(setName, "md.hod_dr9mid_z",
scriptDir, catalogDir, fileNums, redshifts,
numSubvolumes, numSlices, False, lbox, 15, omegaM)
writeScript(setName, "md.hod_dr9mid_z",
scriptDir, catalogDir, fileNums, redshifts,
numSubvolumes, numSlices, True, lbox, 15, omegaM)
if args.hod or args.all:
print " Doing DR9 HOD"
for redshift in redshifts:
for (iRedshift, redshift) in enumerate(redshifts):
print " z = ", redshift
# these parameters come from Manera et al 2012, eq. 26
parFileName = "./hod.par"
parFile = open(parFileName, 'w')
haloFile = catalogDir+"/mdr1_halos_z"+redshift
haloFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift]
parFile.write(parFileText.format(omegaM=omegaM,
hubble=hubble,
redshift=redshift,
@ -456,7 +507,7 @@ if args.hod or args.all:
os.system(hodPath+" "+parFileName+">& /dev/null")
sampleName = getSampleName("hod_dr9mid", redshift, False)
sampleName = getSampleName("md.hod_dr9mid", redshift, False)
outFileName = catalogDir+"/"+sampleName+".dat"
os.system("mv hod.mock" + " " + outFileName)

237
pipeline/prepareGadgetCatalog.py
View file

@ -1,237 +0,0 @@
#!/usr/bin/env python
# script which prepares inputs for gadget-based void run
import numpy as np
import os
import sys
import void_python_tools as vp
import argparse
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# CONFIGURATION
# directory for the input simulation/observational particle files
catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/gadget/"
# where to put the final void catalog, figures, and output logs
voidOutputDir = os.getenv("HOME")+"/workspace/Voids/gadget/"
figDir = os.getenv("PWD")+"/../figs/gadget/"
logDir = os.getenv("PWD")+"/../logs/gadget/"
# where to place the pipeline scripts
scriptDir = os.getenv("PWD")+"/gadget/"
# simulation or observation?
dataType = "simulation"
# available formats for simulation: gadget, multidark
dataFormat = "gadget"
# place particles on the lightcone?
useLightCone = False
# common filename of particle files
redshiftFileBase = "mdr1_particles_z"
# list of redshifts for the particle files
# to get particle file name, we take redshiftFileBase+redshift
redshifts = (("0.0", "0.53", "1.0"))
# how many independent slices along the z-axis?
numSlices = 4
# how many subdivisions along the x- and y- axis?
# ( = 2 will make 4 subvolumes for each slice, = 3 will make 9, etc.)
numSubvolumes = 1
# prefix to give all outputs
prefix = "gadget_"
# list of desired subsamples
subSamples = [ 1.0, 0.1 ]
# simulation information
numPart = 1024*1024*1024
lbox = 1000 # Mpc/h
omegaM = 0.27
hubble = 0.70
# END CONFIGURATION
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
#------------------------------------------------------------------------------
LIGHT_SPEED = 299792.458
def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1):
sampleName = setName
if useVel: setName += "_pv"
if iSlice != -1: sampleName += "_s" + str(iSlice)
if iVol != -1: sampleName += "_d" + iVol
return sampleName
#------------------------------------------------------------------------------
# for given dataset parameters, outputs a script for use with analyzeVoids
def writeScript(setName, dataFileNameBase,
scriptDir, catalogDir, redshifts, numSubvolumes,
numSlices, useVel, lbox, minRadius, omegaM, subsample=1.0,
suffix=".dat"):
if useVel: setName += "_pv"
scriptFileName = scriptDir + "/" + setName + ".py"
scriptFile = open(scriptFileName, 'w')
scriptFile.write("""#!/usr/bin/env/python
import os
from void_python_tools.backend.classes import *
continueRun = False
startCatalogStage = 1
endCatalogStage = 3
startAPStage = 1
endAPStage = 6
ZOBOV_PATH = os.getenv("PWD")+"/../zobov/"
CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/"
freshStack = True
errorBars = "CALCULATED"
numIncoherentRuns = 100
ranSeed = 101010
useLCDM = False
bias = 1.16
dataPortions = ["central"]
dataSampleList = []
""")
dataInfo = """
setName = "{setName}"
workDir = "{voidOutputDir}/{setName}/"
inputDataDir = "{inputDataDir}"
figDir = "{figDir}/{setName}/"
logDir = "{logDir}/{setName}/"
"""
scriptFile.write(dataInfo.format(setName=setName, figDir=figDir,
logDir=logDir, voidOutputDir=voidOutputDir,
inputDataDir=catalogDir))
sampleInfo = """
newSample = Sample(dataFile = "{dataFile}",
dataFormat = "{dataFormat}",
fullName = "{sampleName}",
nickName = "{sampleName}",
dataType = "simulation",
zBoundary = ({zMin}, {zMax}),
zRange = ({zMin}, {zMax}),
zBoundaryMpc = ({zMinMpc}, {zMaxMpc}),
omegaM = {omegaM},
minVoidRadius = {minRadius},
includeInHubble = True,
partOfCombo = False,
isCombo = False,
boxLen = {boxLen},
usePecVel = {usePecVel},
numSubvolumes = {numSubvolumes},
mySubvolume = "{mySubvolume}",
numSubDivisions = 4,
useLightCone = {useLightCone},
subsample = {subsample})
dataSampleList.append(newSample)
"""
for redshift in redshifts:
zBox = float(redshift)
Om = float(omegaM)
zBoxMpc = LIGHT_SPEED/100.*vp.angularDiameter(zBox, Om=Om)
boxMaxMpc = zBoxMpc + lbox
# converter from redshift to comoving distance
zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000)
zVsDX = np.zeros(len(zVsDY))
for i in xrange(len(zVsDY)):
zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
if useLightCone:
boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \
LIGHT_SPEED*lbox, zVsDX, zVsDY)-zBox
dzSafe = 0.03
else:
boxWidthZ = lbox*100./LIGHT_SPEED
dzSafe = 0.0
for iSlice in xrange(numSlices):
sliceMin = zBox + dzSafe + iSlice*(boxWidthZ-dzSafe)/numSlices
sliceMax = zBox + dzSafe + (iSlice+1)*(boxWidthZ-dzSafe)/numSlices
sliceMinMpc = sliceMin*LIGHT_SPEED/100.
sliceMaxMpc = sliceMax*LIGHT_SPEED/100.
sliceMin = "%0.2f" % sliceMin
sliceMax = "%0.2f" % sliceMax
sliceMinMpc = "%0.1f" % sliceMinMpc
sliceMaxMpc = "%0.1f" % sliceMaxMpc
dataFileName = dataFileNameBase + redshift + suffix
for iX in xrange(numSubvolumes):
for iY in xrange(numSubvolumes):
mySubvolume = "%d%d" % (iX, iY)
sampleName = getSampleName(setName, redshift, useVel,
iSlice=iSlice, iVol=mySubvolume)
scriptFile.write(sampleInfo.format(dataFile=dataFileName,
dataFormat=dataFormat,
sampleName=sampleName,
zMin=sliceMin,
zMax=sliceMax,
zMinMpc=sliceMinMpc,
zMaxMpc=sliceMaxMpc,
omegaM=Om,
boxLen=lbox,
usePecVel=useVel,
minRadius=minRadius,
numSubvolumes=numSubvolumes,
mySubvolume=mySubvolume,
useLightCone=useLightCone,
subsample=subsample))
scriptFile.close()
return
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir)
#------------------------------------------------------------------------------
# first the directly downsampled runs
# Note: ss0.002 ~ SDSS DR7 dim2
# ss0.0004 ~ SDSS DR9 mid
baseResolution = numPart/lbox/lbox/lbox # particles/Mpc^3
for thisSubSample in subSamples:
keepFraction = float(thisSubSample) / baseResolution
maxKeep = keepFraction * numPart
minRadius = int(np.ceil(lbox/maxKeep**(1./3)))
print " Doing subsample", thisSubSample, " scripts"
setName = prefix+"ss"+str(thisSubSample)
writeScript(setName, redshiftFileBase, scriptDir, catalogDir, redshifts,
numSubvolumes, numSlices, True, lbox, minRadius, omegaM,
subsample=thisSubSample, suffix="")
writeScript(setName, redshiftFileBase, scriptDir, catalogDir, redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
subsample=thisSubSample, suffix="")

17
plotting/datasetsToPlot.py Executable file
View file

@ -0,0 +1,17 @@
#!/usr/bin/env python
workDir = "/home/psutter2/workspace/Voids/"
figDir = "./figs"
sampleDirList = [ "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/",
"multidark/md_ss01.0_pv/sample_md_ss1.0_pv_z0.56_d00/",
"multidark/md_halos_min1.23e13_pv/sample_md_halos_min1.23e13_pv_z0.56_d00/",
"random/ran_ss0.0004/sample_ran_ss0.0004_z0.56_d00/",
"random/ran_ss0.1/sample_ran_ss0.1_z0.56_d00/",
"multidark/md_hod_dr9mid_pv/sample_md_hod_dr9mid_pv_z0.56_d00/",
"multidark/md_ss0.0004_pv/sample_md_ss0.0004_pv_z0.56_d00/",
"sdss_dr9/sample_lss.dr9cmassmid.dat/" ]
dataPortion = "central"

93
plotting/plotCompareDensCon.py Executable file
View file

@ -0,0 +1,93 @@
#!/usr/bin/env python
# plots cumulative distributions of number counts versus density contrast
from void_python_tools.backend import *
from void_python_tools.plotting import *
import void_python_tools.apTools as vp
import imp
import pickle
import os
import matplotlib.pyplot as plt
import numpy as np
import argparse
# ------------------------------------------------------------------------------
from datasetsToPlot import *
plotNameBase = "compdenscon"
obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing?
parser = argparse.ArgumentParser(description='Plot.')
parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False,
help='display the plot (default: just write eps)')
args = parser.parse_args()
# ------------------------------------------------------------------------------
if not os.access(figDir, os.F_OK):
os.makedirs(figDir)
dataSampleList = []
for sampleDir in sampleDirList:
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
dataSampleList.append(pickle.load(input))
plt.clf()
plt.xlabel("Void Radius (Mpc/h)")
plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]")
plt.yscale('log')
plt.xlim(xmax=80.)
plotName = plotNameBase
for (iSample,sample) in enumerate(dataSampleList):
sampleName = sample.fullName
lineTitle = sampleName
if sample.dataType == "observation":
boxVol = vp.getSurveyProps(sample.maskFile,
sample.zBoundary[0], sample.zBoundary[1],
sample.zRange[0], sample.zRange[1], "all",
selectionFuncFile=sample.selFunFile)[0]
boxVol *= obsFudgeFactor
else:
boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] -
sample.zBoundaryMpc[0])
boxVol *= 1.e-9 # Mpc->Gpc
filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\
sampleName+".out"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
continue
data = np.loadtxt(filename, comments="#")
if data.ndim == 1:
print " Too few!"
continue
data = data[:,8]
indices = np.arange(0, len(data), 1)
sorted = np.sort(data)
plt.plot(sorted, indices[::-1]/boxVol, '-',
label=lineTitle, color=colorList[iSample],
linewidth=linewidth)
plt.legend(title = "Samples", loc = "upper right", prop={'size':8})
#plt.title(plotTitle)
plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
if args.showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")

93
plotting/plotCompareDist.py Executable file
View file

@ -0,0 +1,93 @@
#!/usr/bin/env python
# plots cumulative distributions of number counts
from void_python_tools.backend import *
from void_python_tools.plotting import *
import void_python_tools.apTools as vp
import imp
import pickle
import os
import matplotlib.pyplot as plt
import numpy as np
import argparse
# ------------------------------------------------------------------------------
from datasetsToPlot import *
plotNameBase = "compdist"
obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing?
parser = argparse.ArgumentParser(description='Plot.')
parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False,
help='display the plot (default: just write eps)')
args = parser.parse_args()
# ------------------------------------------------------------------------------
if not os.access(figDir, os.F_OK):
os.makedirs(figDir)
dataSampleList = []
for sampleDir in sampleDirList:
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
dataSampleList.append(pickle.load(input))
plt.clf()
plt.xlabel("Void Radius (Mpc/h)")
plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]")
plt.yscale('log')
plt.xlim(xmax=80.)
plotName = plotNameBase
for (iSample,sample) in enumerate(dataSampleList):
sampleName = sample.fullName
lineTitle = sampleName
if sample.dataType == "observation":
boxVol = vp.getSurveyProps(sample.maskFile,
sample.zBoundary[0], sample.zBoundary[1],
sample.zRange[0], sample.zRange[1], "all",
selectionFuncFile=sample.selFunFile)[0]
boxVol *= obsFudgeFactor
else:
boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] -
sample.zBoundaryMpc[0])
boxVol *= 1.e-9 # Mpc->Gpc
filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\
sampleName+".out"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
continue
data = np.loadtxt(filename, comments="#")
if data.ndim == 1:
print " Too few!"
continue
data = data[:,4]
indices = np.arange(0, len(data), 1)
sorted = np.sort(data)
plt.plot(sorted, indices[::-1]/boxVol, '-',
label=lineTitle, color=colorList[iSample],
linewidth=linewidth)
plt.legend(title = "Samples", loc = "upper right", prop={'size':8})
#plt.title(plotTitle)
plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
if args.showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")

1
python_tools/void_python_tools/backend/__init__.py
View file

@ -1,2 +1,3 @@
from classes import *
from launchers import *
from catalogPrep import *

1453
python_tools/void_python_tools/backend/backend.py
View file

File diff suppressed because it is too large Load diff

10
python_tools/void_python_tools/backend/classes.py
View file

@ -38,6 +38,7 @@ class Sample:
dataType = "observation"
dataFormat = "sdss"
dataFile = "lss.dr72dim.dat"
dataUNit = 1
fullName = "lss.dr72dim.dat"
nickName = "dim"
zobovDir = ""
@ -62,22 +63,21 @@ class Sample:
usePecVel = False
subsample = 1.0
useLightCone = True
numSubDivisions = 1
numSubvolumes = 1
mySubvolume = 1
stacks = []
def __init__(self, dataFile="", fullName="",
def __init__(self, dataFile="", fullName="", dataUnit=1,
nickName="", maskFile="", selFunFile="",
zBoundary=(), zRange=(), zBoundaryMpc=(),
minVoidRadius=0, fakeDensity=0.01, volumeLimited=True,
includeInHubble=True, partOfCombo=False, isCombo=False,
comboList=(), profileBinSize=2.0, skyFraction=0.19,
dataType="observation", numSubDivisions=2,
boxLen=1024, usePecVel=False, omegaM=0.27,
numSubvolumes=1, mySubvolume=1, dataFormat="sdss",
subsample="1.0", useLightCone=True):
dataType="observation",
subsample=1.0, useLightCone=True):
self.dataFile = dataFile
self.fullName = fullName
self.nickName = nickName
@ -97,7 +97,6 @@ class Sample:
self.profileBinSize = profileBinSize
self.skyFraction = skyFraction
self.dataType = dataType
self.numSubDivisions = numSubDivisions
self.boxLen = boxLen
self.usePecVel = usePecVel
self.omegaM = omegaM
@ -106,6 +105,7 @@ class Sample:
self.dataFormat = dataFormat
self.subsample = subsample
self.useLightCone = useLightCone
self.dataUnit = dataUnit
self.stacks = []

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

@ -84,7 +84,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
if os.access("mask_index.txt", os.F_OK):
os.system("mv %s %s" % ("mask_index.txt", zobovDir))
os.system("mv %s %s" % ("total_particles.txt", zobovDir))
os.system("mv %s %s" % ("sample_info.txt", zobovDir))
#os.system("mv %s %s" % ("sample_info.txt", zobovDir))
if os.access("galaxies.txt", os.F_OK):
os.system("mv %s %s" % ("galaxies.txt", zobovDir))
@ -102,7 +102,12 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
else:
includePecVelString = ""
if sample.dataFormat == "multidark":
if sample.useLightCone:
useLightConeString = "cosmo"
else:
useLightConeString = ""
if sample.dataFormat == "multidark" or sample.dataFormat == "random":
dataFileLine = "multidark " + datafile
elif sample.dataFormat == "gadget":
dataFileLine = "gadget " + datafile
@ -120,16 +125,20 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
output %s
outputParameter %s
%s
%s
gadgetUnit %g
rangeX_min %g
rangeX_max %g
rangeY_min %g
rangeY_max %g
rangeZ_min %g
rangeZ_max %g
subsample %g
subsample %s
""" % (dataFileLine, zobovDir+"/zobov_slice_"+sampleName,
zobovDir+"/zobov_slice_"+sampleName+".par",
includePecVelString,
useLightConeString,
sample.dataUnit,
xMin, xMax, yMin, yMax,
sample.zBoundaryMpc[0], sample.zBoundaryMpc[1],
sample.subsample)
@ -150,11 +159,16 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
else:
print "already done!"
if os.access("comoving_distance.txt", os.F_OK):
os.system("mv %s %s" % ("comoving_distance.txt", zobovDir))
#os.system("mv %s %s" % ("sample_info.txt", zobovDir))
if os.access(parmFile, os.F_OK):
os.unlink(parmFile)
# -----------------------------------------------------------------------------
def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None):
def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
numZobovDivisions=None, numZobovThreads=None):
sampleName = sample.fullName
@ -185,9 +199,9 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None):
if os.access(zobovDir+"/voidDesc_"+sampleName+".out", os.F_OK):
os.unlink(zobovDir+"/voidDesc_"+sampleName+".out")
cmd = "%s/vozinit %s 0.1 1.0 2 %s %g %s %s %s >& %s" % \
(binPath, datafile, \
"_"+sampleName, sample.numSubDivisions, \
cmd = "%s/vozinit %s 0.1 1.0 %g %s %g %s %s %s >& %s" % \
(binPath, datafile, numZobovDivisions, \
"_"+sampleName, numZobovThreads, \
binPath, zobovDir, maskIndex, logFile)
os.system(cmd)
@ -232,11 +246,18 @@ def launchPrune(sample, binPath, thisDataPortion=None,
totalPart = open(zobovDir+"/total_particles.txt", "r").read()
maxDen = 0.2*float(mockIndex)/float(totalPart)
observationLine = " --isObservation"
periodicLine = " --periodic=''"
else:
mockIndex = -1
maxDen = 0.2
observationLine = ""
periodicLine = " --periodic='"
if sample.numSubvolumes == 1: periodicLine += "xy"
if sample.zBoundaryMpc[0] == 0 and \
sample.zBoundaryMpc[1] == sample.boxLen : periodicLine += "z"
periodicLine += "' "
if not (continueRun and jobSuccessful(logFile, "NetCDF: Not a valid ID\n")):
cmd = binPath
cmd += " --partFile=" + zobovDir+"/zobov_slice_"+str(sampleName)
@ -255,6 +276,7 @@ def launchPrune(sample, binPath, thisDataPortion=None,
cmd += " --rMin=" + str(sample.minVoidRadius)
cmd += " --numVoids=" + str(numVoids)
cmd += observationLine
cmd += periodicLine
cmd += " --output=" + zobovDir+"/voidDesc_"+\
str(thisDataPortion)+"_"+\
str(sampleName)+".out"
@ -276,7 +298,8 @@ def launchPrune(sample, binPath, thisDataPortion=None,
cmd += " >& " + logFile
os.system(cmd)
if jobSuccessful(logFile, "NetCDF: Not a valid ID\n"):
if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or \
jobSuccessful(logFile, "Done!\n"):
print "done"
else:
print "FAILED!"
@ -356,8 +379,8 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
%s
ranSeed %d
dataPortion %s
barycenters %s
boundaryDistances %s
#barycenters %s
#boundaryDistances %s
%s
""" % \
(zobovDir+"/voidDesc_"+thisDataPortion+"_"+sampleName+".out",
@ -460,34 +483,37 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
return
# figure out box volume and average density
maskFile = sample.maskFile
sulFunFile = sample.selFunFile
if sample.dataType == "observation":
maskFile = sample.maskFile
sulFunFile = sample.selFunFile
if not os.access(sample.selFunFile, os.F_OK) and not volumeLimited:
print " Cannot find", selFunFile, "!"
exit(-1)
if not os.access(sample.selFunFile, os.F_OK) and not sample.volumeLimited:
print " Cannot find", selFunFile, "!"
exit(-1)
sys.stdout = open(logFile, 'a')
sys.stderr = open(logFile, 'a')
zMin = sample.zRange[0]
zMax = sample.zRange[1]
if not sample.volumeLimited:
props = vp.getSurveyProps(maskFile, stack.zMin,
stack.zMax, zMin, zMax, "all",
selectionFuncFile=sample.selFunFile)
sys.stdout = open(logFile, 'a')
sys.stderr = open(logFile, 'a')
zMin = sample.zRange[0]
zMax = sample.zRange[1]
if not sample.volumeLimited:
props = vp.getSurveyProps(maskFile, stack.zMin,
stack.zMax, zMin, zMax, "all",
selectionFuncFile=sample.selFunFile)
else:
zMinForVol = sample.zBoundary[0]
zMaxForVol = sample.zBoundary[1]
props = vp.getSurveyProps(maskFile, zMinForVol,
zMaxForVol, zMin, zMax, "all")
sys.stdout = sys.__stdout__
sys.stderr = sys.__stderr__
boxVol = props[0]
nbar = props[1]
if sample.volumeLimited:
nbar = 1.0
else:
zMinForVol = sample.zBoundary[0]
zMaxForVol = sample.zBoundary[1]
props = vp.getSurveyProps(maskFile, zMinForVol,
zMaxForVol, zMin, zMax, "all")
sys.stdout = sys.__stdout__
sys.stderr = sys.__stderr__
boxVol = props[0]
nbar = props[1]
if sample.volumeLimited:
nbar = 1.0
boxVol = sample.boxLen**3
summaryLine = runSuffix + " " + \
thisDataPortion + " " + \
@ -809,7 +835,7 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
maxtries = 5
while badChain:
Rexpect = (stack.rMin+stack.rMax)/2
Rtruncate = stack.rMax*3. + 1 # TEST
Rtruncate = stack.rMin*3. + 1 # TEST
ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect,
Niter=300000,
Nburn=100000,
@ -1173,7 +1199,11 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
voidDir = sample.zobovDir+"/stacks_" + runSuffix
centersFile = voidDir+"/centers.txt"
if os.access(centersFile, os.F_OK):
voidRedshifts = np.loadtxt(centersFile)[:,5]
voidRedshifts = np.loadtxt(centersFile)
if voidRedshifts.ndim > 1:
voidRedshifts = voidRedshifts[:,5]
else:
voidRedshifts = voidRedshifts[5]
#fp.write(str(len(voidRedshifts))+" ")
np.savetxt(fp, voidRedshifts[None])
else:

2
python_tools/void_python_tools/plotting/plotDefs.py
View file

@ -1,3 +1,5 @@
LIGHT_SPEED = 299792.458
colorList = ['r', 'b', 'g', 'y', 'c', 'm', 'y',
'brown', 'grey',
'darkred', 'orange', 'pink', 'darkblue',

218
python_tools/void_python_tools/plotting/plotTools.py
View file

@ -1,17 +1,20 @@
__all__=['plotNumberCounts']
__all__=['plotRedshiftDistribution', 'plotSizeDistribution', 'plot1dProfiles',
'plotMarg1d', 'plotNumberDistribution']
from void_python_tools.backend.classes import *
from plotDefs import *
import numpy as np
import os
import pylab as plt
import void_python_tools.apTools as vp
# -----------------------------------------------------------------------------
def plotNumberCounts(workDir=None, sampleList=None, figDir=None, plotNameBase="numbercount",
def plotRedshiftDistribution(workDir=None, sampleList=None, figDir=None,
plotNameBase="zdist",
showPlot=False, dataPortion=None, setName=None):
plt.clf()
plt.xlabel("Comoving Distance (Mpc/h)")
plt.xlabel("Redshift")
plt.ylabel("Number of Voids")
plotTitle = setName
@ -41,7 +44,7 @@ def plotNumberCounts(workDir=None, sampleList=None, figDir=None, plotNameBase="n
zMax = sample.zRange[1]
range = (zMin, zMax)
nbins = np.ceil((zMax-zMin)/0.1)
nbins = np.ceil((zMax-zMin)/0.02)
thisMax = np.max(data[:,5])
thisMin = np.min(data[:,5])
@ -67,3 +70,210 @@ def plotNumberCounts(workDir=None, sampleList=None, figDir=None, plotNameBase="n
os.system("display %s" % figDir+"/fig_"+plotName+".png")
# -----------------------------------------------------------------------------
def plotSizeDistribution(workDir=None, sampleList=None, figDir=None,
plotNameBase="sizedist",
showPlot=False, dataPortion=None, setName=None):
plt.clf()
plt.xlabel("Void Radius (Mpc/h)")
plt.ylabel("Number of Voids")
plotTitle = setName
plotName = plotNameBase
xMin = 1.e00
xMax = 0
for (iSample,sample) in enumerate(sampleList):
sampleName = sample.fullName
lineTitle = sampleName
filename = workDir+"/sample_"+sampleName+"/centers_"+dataPortion+"_"+\
sampleName+".out"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
continue
data = np.loadtxt(filename, comments="#")
if data.ndim == 1:
print " Too few!"
continue
xMin = 5
xMax = 140
range = (xMin, xMax)
nbins = np.ceil((xMax-xMin)/5)
#thisMax = np.max(data[:,5])
#thisMin = np.min(data[:,5])
#if thisMax > xMax: xMax = thisMax
#if thisMin < xMin: xMin = thisMin
plt.hist(data[:,4], bins=nbins,
label=lineTitle, color=colorList[iSample],
histtype = "step", range=range,
linewidth=linewidth)
plt.legend(title = "Samples", loc = "upper right")
plt.title(plotTitle)
plt.xlim(xMin, xMax)
#plt.xlim(xMin, xMax*1.4) # make room for legend
plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
if showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")
# -----------------------------------------------------------------------------
def plot1dProfiles(workDir=None, sampleList=None, figDir=None,
plotNameBase="1dprofile",
showPlot=False, dataPortion=None, setName=None):
plt.clf()
plt.xlabel(r"$R/R_{v,\mathrm{max}}$")
plt.ylabel(r"$n / \bar n$")
for (iSample,sample) in enumerate(sampleList):
sampleName = sample.fullName
for (iStack,stack) in enumerate(sample.stacks):
plotTitle = setName
plotName = plotNameBase
runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin,
stack.rMax, dataPortion)
plotTitle = sampleName + ", z = "+str(stack.zMin)+"-"+str(stack.zMax)+", R = "+str(stack.rMin)+"-"+str(stack.rMax)+ r" $h^{-1}$ Mpc"
filename = workDir+"/sample_"+sampleName+"/stacks_"+runSuffix+\
"profile_1d.txt"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
continue
data = np.loadtxt(filename, comments="#")
if data.ndim == 1:
print " Too few!"
continue
data[:,1] /= stack.rMax
plt.ylim(ymin=0.0, ymax=np.amax(data[:,2])+0.1)
plt.xlim(xmin=0.0, xmax=2.1)
plt.plot(data[:,1], data[:,2], label=lineTitle, color=colorList[0],
linewidth=linewidth)
plt.title(plotTitle)
plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
if showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")
# -----------------------------------------------------------------------------
def plotMarg1d(workDir=None, sampleList=None, figDir=None,
plotNameBase="marg1d",
showPlot=False, dataPortion=None, setName=None):
plotNames = ("Om", "w0", "wa")
plotTitles = ("$\Omega_M$", "$w_0$", "$w_a$")
files = ("Om", "w0", "wa")
for iPlot in range(len(plotNames)):
plt.clf()
plotName = plotNameBase+"_"+plotNames[iPlot]+"_"+dataPortion
plotTitle = plotTitles[iPlot]
dataFile = workDir + "/likelihoods_"+dataPortion+"_"+files[iPlot]+".dat"
plt.xlabel(plotTitle, fontsize="20")
plt.ylabel("Likelihood", fontsize="20")
plt.ylim(0.0, 1.0)
data = np.loadtxt(dataFile, comments="#")
plt.plot(data[:,0], data[:,1], color='k', linewidth=linewidth)
plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
if showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")
# -----------------------------------------------------------------------------
def plotNumberDistribution(workDir=None, sampleList=None, figDir=None,
plotNameBase="numberdist",
showPlot=False, dataPortion=None, setName=None):
plt.clf()
plt.xlabel("Void Radius (Mpc/h)")
plt.ylabel(r"N > R [$h^3$ Mpc$^{-3}$]")
plotTitle = setName
plotName = plotNameBase
plt.yscale('log')
for (iSample,sample) in enumerate(sampleList):
sampleName = sample.fullName
lineTitle = sampleName
if sample.dataType == "observation":
boxVol = vp.getSurveyProps(sample.maskFile,
sample.zBoundary[0], sample.zBoundary[1],
sample.zRange[0], sample.zRange[1], "all",
selectionFuncFile=sample.selFunFile)[0]
else:
boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] -
sample.zBoundaryMpc[0])
boxVol *= 1.e-9
filename = workDir+"/sample_"+sampleName+"/centers_"+dataPortion+"_"+\
sampleName+".out"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
continue
data = np.loadtxt(filename, comments="#")
if data.ndim == 1:
print " Too few!"
continue
data = data[:,4]
indices = np.arange(0, len(data), 1)
sorted = np.sort(data)
plt.plot(sorted, indices[::-1]/boxVol, '-',
label=lineTitle, color=colorList[iSample],
linewidth=linewidth)
plt.legend(title = "Samples", loc = "upper right")
plt.title(plotTitle)
plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
if showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")