From 3dce2593d939cbf03b4e3f4ed4cbd418573dd13e Mon Sep 17 00:00:00 2001 From: "Paul M. Sutter" Date: Tue, 7 Jan 2025 20:04:29 +0800 Subject: [PATCH] Implemented (yet another) new boundary handling scheme, whereby we scan radially along survey edge while flagging nearest galaxies. The prepObservation routine was significantly cleaned up to accommodate this, but it was ultimately implemented in python (surveyTools.py) for ease of prototyping, with the intent to move it back into C later. Some general housekeeping, making sure some new parameters are passed around correctly, and removing the storage of some unused files. This update is considered HIGHLY UNSTABLE. It will almost certainly break somewhere for simulations. Still under active development. --- README.md | 26 +- c_source/prep/prepObservation.cpp | 427 +++++------------- c_source/prep/prepObservation.ggo | 4 + c_source/pruning/pruneVoids.cpp | 28 +- .../example_observation.py | 18 +- python_source/backend/classes.py | 12 +- python_source/backend/launchers.py | 117 +++-- python_source/backend/surveyTools.py | 147 ++++-- python_source/vide_pipeline/__main__.py | 23 +- 9 files changed, 348 insertions(+), 454 deletions(-) diff --git a/README.md b/README.md index c697861..64319b4 100644 --- a/README.md +++ b/README.md @@ -89,26 +89,24 @@ Usage: python3 -m vide_pipeline parameter_file.py The VIDE tools are all packaged in the `vide` package. +Running with observation +----------------------- + +An example parameter file and dataset is given in the examples/example_observation directory. The parameter file contains all the information VIDE needs to run: where to find inputs and place outputs, tolerances for managing boundary handling, and information about your particular datasets, like redshift boundaries. To see how this works, here is an example: + +cd examples/example_observation +python3 -m vide_pipeline example_observation.py + Running with simulation ----------------------- -Using simulation requires a preliminary step, consisting in using the script -`vide_prepare_simulation` which is installed during the installation procedure. -The script generates mock catalog and a default pipeline to handle simulations. -An example of the complete procedure is given here-below: -``` -mkdir $HOME/my_vide_test -cp python_tools/vide_pipeline/datasets/example_simulation.py $HOME/my_vide_test -mkdir $HOME/my_vide_test/examples -cp examples/example_simulation_z0.0.dat $HOME/my_vide_test/examples -cd $HOME/my_vide_test -vide_prepare_simulation --all --parm example_simulation.py +Working with simulations requires a preliminary step, consisting in using the script "vide_prepare_simulation" which is installed automatically. This script performs necessary processing on your simulation file, such as extracting slices, performing subsampling, placing particles on a lightcone, and so on. For a demonstration, see the "example_simulation.py" parameter file in the examples/example_simulation/ directory. Running this script creates a series of auxillary parameter files that can then be run individually for void finding. Here is an example of this procedure: + +cd examples/example_simulation +vide_prepare_simulation --all --parm example_simulation.py python3 -m vide_pipeline example_simulation/sim_ss1.0.py ``` -The example copies the required data in a separate directory. Then, we execute -the `vide_prepare_simulation` script to generate the auxiliary pipeline. The -`vide_pipeline` is finally executed on this generated script. Notes for CONDA --------------- diff --git a/c_source/prep/prepObservation.cpp b/c_source/prep/prepObservation.cpp index bf74b60..15f3e43 100644 --- a/c_source/prep/prepObservation.cpp +++ b/c_source/prep/prepObservation.cpp @@ -18,9 +18,9 @@ +*/ - #include #include +#include #include #include #include @@ -30,7 +30,6 @@ #include "contour_pixels.hpp" #include #include -#include #include #define LIGHT_SPEED 299792.458 @@ -65,10 +64,10 @@ struct ParticleData vector redshift; vector catalogID; vector uniqueID; - int id_mask; // PMS int mask_index; // END PMS + int edgeFlag = 0; vector pos; double box[3][2]; double Lmax; @@ -122,15 +121,13 @@ void loadData(const string& fname, NYU_VData & data) } void placeGalaxiesInCube(NYU_VData& data, ParticleData& output_data, - bool useComoving, double omegaM) -{ + bool useComoving, double omegaM) { double d2r = M_PI/180; gsl_function expanF; expanF.function = &expanFun; struct my_expan_params expanParams; - double maxZ = 2.0, z, result, error, *dL, *redshifts; - int numZ = 1000, iZ; + double result, error; size_t nEval; expanParams.Om = omegaM; @@ -138,20 +135,6 @@ void placeGalaxiesInCube(NYU_VData& data, ParticleData& output_data, expanParams.wa = 0.0; expanF.params = &expanParams; - dL = (double *) malloc(numZ * sizeof(double)); - redshifts = (double *) malloc(numZ * sizeof(double)); - - for (iZ = 0; iZ < numZ; iZ++) { - z = iZ * maxZ/numZ; - //gsl_integration_qng(&expanF, 0.0, z, 1.e-6, 1.e-6, &result, &error, &nEval); - dL[iZ] = result; - redshifts[iZ] = z; - } - - gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, numZ); - gsl_interp_init(interp, redshifts, dL, numZ); - gsl_interp_accel *acc = gsl_interp_accel_alloc(); - output_data.pos.resize(data.size()); output_data.id_gal.resize(data.size()); output_data.ra.resize(data.size()); @@ -168,43 +151,36 @@ void placeGalaxiesInCube(NYU_VData& data, ParticleData& output_data, for (int i = 0; i < data.size(); i++) { double ra = data[i].ra*d2r, dec = data[i].dec*d2r; + double Dc = data[i].cz; Position& p = output_data.pos[i]; if (useComoving) { - //double pos = gsl_interp_eval(interp, redshifts, dL, data[i].cz, acc); - // Maubert - Lower boundary in redshift set to 0 to be consistent with pruneVoids (was 1.e-6 before). gsl_integration_qng(&expanF, 0.0, data[i].cz/LIGHT_SPEED, - 1.e-6, - 1.e-6, &result, &error, &nEval); - double Dc = result*LIGHT_SPEED; - p.xyz[0] = Dc*cos(ra)*cos(dec); - p.xyz[1] = Dc*sin(ra)*cos(dec); - p.xyz[2] = Dc*sin(dec); - } else { - p.xyz[0] = data[i].cz*cos(ra)*cos(dec); - p.xyz[1] = data[i].cz*sin(ra)*cos(dec); - p.xyz[2] = data[i].cz*sin(dec); + 1.e-6, 1.e-6, &result, &error, &nEval); + Dc = result*LIGHT_SPEED; } -//printf("CREATE %e %e\n", data[i].cz, sqrt(p.xyz[0]*p.xyz[0] + p.xyz[1]*p.xyz[1] + p.xyz[2]*p.xyz[2])); + + p.xyz[0] = Dc*cos(ra)*cos(dec); + p.xyz[1] = Dc*sin(ra)*cos(dec); + p.xyz[2] = Dc*sin(dec); + output_data.id_gal[i] = data[i].index; 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++) - { - if (p.xyz[j] > output_data.box[j][0]) - output_data.box[j][0] = p.xyz[j]; - if (p.xyz[j] < output_data.box[j][1]) - output_data.box[j][1] = p.xyz[j]; - } -//printf("INSERT GAL %d %e %e %e\n", output_data.id_gal[i], p.xyz[0], p.xyz[1], p.xyz[2]); + for (int j = 0; j < 3; j++) { + if (p.xyz[j] > output_data.box[j][0]) + output_data.box[j][0] = p.xyz[j]; + if (p.xyz[j] < output_data.box[j][1]) + output_data.box[j][1] = p.xyz[j]; + } } - // normalize box - float left = 1.e99; - float right = -1.e99; + // normalize the box volume + float left = INFINITY; + float right = -INFINITY; for (int j = 0; j < 3; j++) { if (output_data.box[j][1] < left) left = output_data.box[j][1]; if (output_data.box[j][0] > right) right = output_data.box[j][0]; @@ -214,24 +190,27 @@ void placeGalaxiesInCube(NYU_VData& data, ParticleData& output_data, output_data.box[j][0] = right; } + double Rmax = -1; + for (int j = 0; j < 3; j++) { + Rmax = max(Rmax, max(output_data.box[j][0], -output_data.box[j][1])); + } + output_data.Lmax = Rmax; + cout << format("Galaxy position generated: %d galaxies") % output_data.pos.size() << endl; cout << format("box is %g < x < %g; %g < y < %g; %g < z < %g") % (1e-2*output_data.box[0][1]) % (1e-2*output_data.box[0][0]) % (1e-2*output_data.box[1][1]) % (1e-2*output_data.box[1][0]) % (1e-2*output_data.box[2][1]) % (1e-2*output_data.box[2][0]) << endl; - gsl_interp_free(interp); } -void generateSurfaceMask(prepObservation_info& args , +void flagEdgeGalaxies(prepObservation_info& args , Healpix_Map& mask, - vector& pixel_list, - vector& full_mask_list, + vector& contourPixels, NYU_VData& data, ParticleData& output_data, bool useComoving, - double omegaM) -{ + double omegaM) { //Maubert - Needed for comobile distance in mock_sphere gsl_function expanF; @@ -242,38 +221,14 @@ void generateSurfaceMask(prepObservation_info& args , expanParams.wa = 0.0; expanF.params = &expanParams; - double result, error ; + double result, error; size_t nEval; //End Maubert - Needed for comobile distance in mock_sphere - // Find the first free index - int idx = -1; - int insertion = 0; - double volume = pixel_list.size()*1.0/mask.Npix()*4*M_PI; - int numToInsert; - - for (int i = 0; i < output_data.id_gal.size(); i++) - { - if (idx < output_data.id_gal[i]) - idx = output_data.id_gal[i]+1; - } - - output_data.id_mask = idx; - -// PMS + // TODO - REMOVE THIS output_data.mask_index = output_data.id_gal.size(); -// END PMS - cout << "Generate surface mask..." << endl; - double Rmax = -1; - for (int j = 0; j < 3; j++) - { - Rmax = max(Rmax, max(output_data.box[j][0], -output_data.box[j][1])); - } - - output_data.Lmax = Rmax; - -// PMS - write a small text file with galaxy position (for diagnostic purposes) + // write a small text file with galaxy position (for diagnostic purposes) FILE *fp; fp = fopen("galaxies.txt", "w"); for (int i = 0; i < data.size(); i++) { @@ -284,208 +239,82 @@ void generateSurfaceMask(prepObservation_info& args , (p.xyz[2])); } fclose(fp); -// END PMS - cout << format("Rmax is %g, surface volume is %g") % (Rmax/100) % (volume/(4*M_PI)) << endl; - volume *= Rmax*Rmax*Rmax/3/1e6; - numToInsert = (int)floor(volume*args.density_fake_arg); - // TEST NOT USING MOCK PARTICLES - numToInsert = 0; - // END TEST - cout << format("3d volume to fill: %g (Mpc/h)^3") % volume << endl; +/* NOTE: temporarily moved to python for quick debugging. Will move back to + here once it's all sorted - cout << format("Will insert %d particles") % numToInsert << endl; - - fp = fopen("mock_galaxies.txt", "w"); - - double pct = 0; - for (int i = 0; i < numToInsert; i++) { - double new_pct = i*100./numToInsert; - - if (new_pct-pct > 5.) { - pct = new_pct; - cout << format(" .. %3.0f %%") % pct << endl; - } - - Position p; - bool stop_here; - - do { - int p0 = (int)floor(drand48()*pixel_list.size()); - vec3 v = mask.pix2vec(pixel_list[p0]); - double r = Rmax*pow(drand48(),1./3); - - p.xyz[0] = v.x * r; - p.xyz[1] = v.y * r; - p.xyz[2] = v.z * r; - - stop_here = true; - for (int j = 0; j < 3; j++) { - if (p.xyz[j] > output_data.box[j][0] || - p.xyz[j] < output_data.box[j][1]) - stop_here = false; - } - } - while (!stop_here); - -// PMS : write mock galaxies to a small file for diagnostic purposes - fprintf(fp, "%e %e %e\n", - (p.xyz[0]), - (p.xyz[1]), - (p.xyz[2])); -// END PMS - output_data.pos.push_back(p); - output_data.id_gal.push_back(idx); - 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++; - } - - fclose(fp); - - // PMS - // TEST - insert mock galaxies along box edge - this is for tesselation safety - fp = fopen("mock_boundary.txt", "w"); - double dx[3]; - dx[0] = output_data.box[0][1] - output_data.box[0][0]; - dx[1] = output_data.box[1][1] - output_data.box[1][0]; - dx[2] = output_data.box[2][1] - output_data.box[2][0]; - - 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++) { - - int iy = (iDir + 1) % 3; - int iz = (iDir + 2) % 3; - - for (int i = 0; i < nPart; i++) { - for (int j = 0; j < nPart; j++) { - Position p; - - p.xyz[iDir] = output_data.box[iDir][iFace]; - p.xyz[iy] = i * dx[iy]/nPart + output_data.box[iy][0]; - p.xyz[iz] = j * dx[iz]/nPart + output_data.box[iz][0]; - - output_data.pos.push_back(p); - output_data.id_gal.push_back(idx); - 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]), - (p.xyz[1]), - (p.xyz[2])); - } - } - } - } - fclose(fp); - // END PMS - - // PMS - // TEST - insert mock galaxies along spheres of survey redshift boundaries - fp = fopen("mock_sphere.txt", "w"); - //Maubert - insert mock galaxies according to useComoving specification - // Added & compute rmin & rmax out of loop - double rmin ; - double rmax ; + // convert redshift boundaries to covmoving if necessary + double rmin; + double rmax; if (useComoving) { - gsl_integration_qng(&expanF, 0.0, args.zMin_arg, - 1.e-6, - 1.e-6, &result, &error, &nEval); - rmin = result* LIGHT_SPEED; + gsl_integration_qng(&expanF, 0.0, args.zMin_arg, 1.e-6, 1.e-6, &result, + &error, &nEval); + rmin = result*LIGHT_SPEED; - gsl_integration_qng(&expanF, 0.0, args.zMax_arg, - 1.e-6, - 1.e-6, &result, &error, &nEval); - rmax = result* LIGHT_SPEED; + gsl_integration_qng(&expanF, 0.0, args.zMax_arg, 1.e-6, 1.e-6, &result, + &error, &nEval); + rmax = result*LIGHT_SPEED; } else { - rmin = args.zMin_arg * LIGHT_SPEED; - rmax = args.zMax_arg * LIGHT_SPEED; + rmin = args.zMin_arg * LIGHT_SPEED; + rmax = args.zMax_arg * LIGHT_SPEED; } - - // TEST NOT USING BOUNDARY PARTICLES - for (int q = 0; q < 0; q++) { - //for (int q = 0; q < full_mask_list.size(); q++) { - vec3 v = mask.pix2vec(full_mask_list[q]); - - Position p; - - - if (rmin > 0.) { - p.xyz[0] = v.x * rmin; - p.xyz[1] = v.y * rmin; - p.xyz[2] = v.z * rmin; - output_data.pos.push_back(p); - output_data.id_gal.push_back(idx); - 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]), - (p.xyz[1]), - (p.xyz[2])); - } - - - p.xyz[0] = v.x * rmax; - p.xyz[1] = v.y * rmax; - p.xyz[2] = v.z * rmax; - output_data.pos.push_back(p); - output_data.id_gal.push_back(idx); - 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]), - (p.xyz[1]), - (p.xyz[2])); - } - fclose(fp); - // END PMS - - cout << format("Done. Inserted %d particles.") % insertion << endl; -} - -void saveData(ParticleData& pdata) -{ - NcFile f("particles.nc", NcFile::replace); - - NcDim d = f.addDim("space", 3); - NcDim p = f.addDim("Np", pdata.pos.size()); - NcVar v = f.addVar("particles", ncDouble, {d, p}); - double *x = new double[pdata.pos.size()]; - - for (int j = 0; j < 3; j++) - { - - for (int i = 0; i < pdata.pos.size(); i++) - x[i] = pdata.pos[i].xyz[j]; - - v.putVar({size_t(j), 0}, {1, pdata.pos.size()}, x); - } - - v = f.addVar("id_gal", ncInt, std::vector({p})); - v.putVar(&pdata.id_gal[0]); - - delete[] x; -} + double dx = args.meanPartSep_arg; + int nSteps = floor( (rmax - rmin) / dx); + cout << "Assumed resolution element: " << dx << " " << nSteps << endl; -void saveForZobov(ParticleData& pdata, const string& fname, const string& paramname) + // flag galaxies near mask edges + // using the "ray marching" algorithm: follow rays along lines of sight + // of all mask edges, flagging nearest neighbor galaxies as we go + // TODO - replace this with faster kd-tree search + cout << "Flagging galaxies on edges of survey..." << endl; + + // explore rays along mask contours + for (int pixel : contourPixels) { + cout << "Working with pixel " << pixel << endl; + vec3 v = mask.pix2vec(pixel); + //cout << v*rmin << " " << v*rmax << endl; + + // march along single ray and find nearest neighbors + for (int n = 0; n <= nSteps; n++) { + double r = n*dx + rmin; + vec3 rayPos = v*r; + + double x = rayPos.x; + double y = rayPos.y; + double z = rayPos.z; + + //cout << "Step " << n << " " << rayPos << endl; + + // scan all galaxies + double minDist = INFINITY; + double dist = 0; + int closestGal = -1; + for (int i = 0; i < data.size(); i++) { + Position& galPos = output_data.pos[i]; + + dist = pow(galPos.xyz[0] - x, 2) + + pow(galPos.xyz[1] - y, 2) + + pow(galPos.xyz[2] - z, 2); + + if (dist < minDist) closestGal = i; + } // galaxy search + + + } // marching along one ray + + } // all contours + + + // flag galaxies near redshift boundaries + cout << "Flagging galaxies at redshift boundaries..." << endl; + +*/ + +} // end flagEdgeGalaxies + +void saveForZobov(ParticleData& pdata, const string& fname, + const string& paramname) { UnformattedWrite f(fname); static const char axis[] = { 'X', 'Y', 'Z' }; @@ -562,21 +391,6 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn v.putVar({0}, {size_t(nOutputPart)}, &pdata.id_gal[0]); //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) @@ -608,40 +422,40 @@ int main(int argc, char **argv) prepObservation_conf_print_version(); - cout << "Loading data " << args_info.catalog_arg << "..." << endl; + cout << "Loading galaxy data " << args_info.catalog_arg << "..." << endl; vector data; - Healpix_Map o_mask; - vector pixel_list; - vector full_mask_list; + vector contourPixels; ParticleData output_data; loadData(args_info.catalog_arg, data); - cout << "Loading mask..." << endl; + Healpix_Map mask; + Healpix_Map o_mask; + + int newNside = args_info.nsideForContour_arg; read_Healpix_map_from_fits(args_info.mask_arg, o_mask); - Healpix_Map mask; + if (newNside == -1) newNside = o_mask.Nside(); - mask.SetNside(128, RING); + mask.SetNside(newNside, RING); 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 + computeContourPixels(mask, contourPixels); cout << "Placing galaxies..." << endl; placeGalaxiesInCube(data, output_data, args_info.useComoving_flag, - args_info.omegaM_arg); - generateSurfaceMask(args_info, mask, pixel_list, full_mask_list, - data, output_data,args_info.useComoving_flag, - args_info.omegaM_arg); - + args_info.omegaM_arg); + + + //cout << "Flagging edge galaxies..." << endl; + flagEdgeGalaxies(args_info, mask, contourPixels, + data, output_data,args_info.useComoving_flag, + args_info.omegaM_arg); + saveForZobov(output_data, args_info.output_arg, args_info.params_arg); - // saveData(output_data); - // PMS + // PMS - TODO REMOVE THIS FILE *fp = fopen("mask_index.txt", "w"); fprintf(fp, "%d", output_data.mask_index); fclose(fp); @@ -650,6 +464,7 @@ int main(int argc, char **argv) fprintf(fp, "%d", output_data.pos.size()); fclose(fp); printf("Done!\n"); - // END PMS + // END PMS - TODO REMOVE THIS + return 0; } diff --git a/c_source/prep/prepObservation.ggo b/c_source/prep/prepObservation.ggo index e04cf1d..756a9b4 100644 --- a/c_source/prep/prepObservation.ggo +++ b/c_source/prep/prepObservation.ggo @@ -17,3 +17,7 @@ option "params" - "Output parameters of the datacube" string required option "useComoving" - "Convert to real space using LCDM cosmology" flag off option "omegaM" - "Omega Matter for fiducial cosmology" double optional default="0.27" + +option "nsideForContour" - "HEALPix NSIDE resolution for figuring out mask contours" int optional default="-1" + +option "meanPartSep" - "Estimated mean tracer seperation in h^3 / Mpc^3" double optional default="1" diff --git a/c_source/pruning/pruneVoids.cpp b/c_source/pruning/pruneVoids.cpp index f8de89a..f632d25 100644 --- a/c_source/pruning/pruneVoids.cpp +++ b/c_source/pruning/pruneVoids.cpp @@ -120,12 +120,12 @@ void openFiles(string outputDir, string sampleName, int mockIndex, int numKept, FILE** fpZobov, FILE** fpCenters, FILE** fpCentersNoCut, - FILE** fpBarycenter, FILE** fpDistances, FILE** fpShapes, + FILE** fpBarycenter, FILE** fpShapes, FILE** fpSkyPositions); void closeFiles(FILE* fpZobov, FILE* fpCenters, FILE* fpCentersNoCut, - FILE* fpBarycenter, FILE* fpDistances, FILE* fpShapes, + FILE* fpBarycenter, FILE* fpShapes, FILE* fpSkyPositions); void outputVoids(string outputDir, string sampleName, string prefix, @@ -374,7 +374,7 @@ int main(int argc, char **argv) { fclose(fp); // now the particles-zone - printf(" Loading particle-zone membership info...\n"); + printf(" Loading zone-particle membership info...\n"); fp = fopen(args.zone2Part_arg, "r"); fread(&dummy, 1, 4, fp); fread(&numZonesTot, 1, 4, fp); @@ -470,7 +470,7 @@ int main(int argc, char **argv) { // load voids *again* using Guilhem's code so we can get tree information clock3 = clock(); - printf(" Re-loading voids and building tree..\n"); + printf(" Re-loading voids and building tree...\n"); ZobovRep zobovCat; if (!loadZobov(args.voidDesc_arg, args.zone2Part_arg, args.void2Zone_arg, 0, zobovCat)) { @@ -691,12 +691,12 @@ int main(int argc, char **argv) { sqrt(pow(voids[iVoid].macrocenter[0] - boxLen[0]/2.,2) + pow(voids[iVoid].macrocenter[1] - boxLen[1]/2.,2) + pow(voids[iVoid].macrocenter[2] - boxLen[2]/2.,2)); - voids[iVoid].redshiftInMpc = voids[iVoid].redshiftInMpc; + //voids[iVoid].redshiftInMpc = voids[iVoid].redshiftInMpc; if (args.useComoving_flag) { redshift = gsl_interp_eval(interp, dL, redshifts, - voids[iVoid].redshiftInMpc, acc); + voids[iVoid].redshiftInMpc, acc); nearestEdge = fabs((redshift-args.zMax_arg)*LIGHT_SPEED/100.); voids[iVoid].redshift = redshift; } else { @@ -1002,7 +1002,7 @@ void openFiles(string outputDir, string sampleName, string prefix, string dataPortion, int mockIndex, int numKept, FILE** fpZobov, FILE** fpCenters, - FILE** fpBarycenter, FILE** fpDistances, FILE** fpShapes, + FILE** fpBarycenter, FILE** fpShapes, FILE** fpSkyPositions) { *fpZobov = fopen((outputDir+"/"+prefix+"voidDesc_"+dataPortion+"_"+sampleName).c_str(), "w"); @@ -1014,8 +1014,6 @@ void openFiles(string outputDir, string sampleName, *fpCenters = fopen((outputDir+"/"+prefix+"centers_"+dataPortion+"_"+sampleName).c_str(), "w"); fprintf(*fpCenters, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part, parent ID, tree level, number of children, central density\n"); - *fpDistances = fopen((outputDir+"/"+prefix+"boundaryDistances_"+dataPortion+"_"+sampleName).c_str(), "w"); - *fpSkyPositions = fopen((outputDir+"/"+prefix+"sky_positions_"+dataPortion+"_"+sampleName).c_str(), "w"); fprintf(*fpSkyPositions, "# RA, dec, redshift, radius (Mpc/h), void ID\n"); @@ -1027,13 +1025,13 @@ void openFiles(string outputDir, string sampleName, // ---------------------------------------------------------------------------- void closeFiles(FILE* fpZobov, FILE* fpCenters, - FILE* fpBarycenter, FILE* fpDistances, FILE* fpShapes, + FILE* fpBarycenter, FILE* fpShapes, FILE* fpSkyPositions) { fclose(fpZobov); fclose(fpCenters); fclose(fpBarycenter); - fclose(fpDistances); + //fclose(fpDistances); fclose(fpShapes); fclose(fpSkyPositions); @@ -1049,13 +1047,13 @@ void outputVoids(string outputDir, string sampleName, string prefix, int iVoid; VOID outVoid; FILE *fp, *fpZobov, *fpCenters, *fpCentersNoCut, *fpBarycenter, - *fpDistances, *fpShapes, *fpSkyPositions; + *fpShapes, *fpSkyPositions; openFiles(outputDir, sampleName, prefix, dataPortion, mockIndex, voids.size(), &fpZobov, &fpCenters, &fpBarycenter, - &fpDistances, &fpShapes, &fpSkyPositions); + &fpShapes, &fpSkyPositions); for (iVoid = 0; iVoid < voids.size(); iVoid++) { @@ -1104,6 +1102,7 @@ void outputVoids(string outputDir, string sampleName, string prefix, outVoid.macrocenter[1], outVoid.macrocenter[2]); + /* fprintf(fpDistances, "%d %e %e %e %e %e\n", outVoid.voidID, outVoid.nearestMock, @@ -1111,6 +1110,7 @@ void outputVoids(string outputDir, string sampleName, string prefix, outVoid.rescaledCoreDens, outVoid.nearestMockFromCore, outVoid.nearestGalFromCore); + */ fprintf(fpCenters, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d %d %d %d %.2f\n", outCenter[0], @@ -1164,6 +1164,6 @@ void outputVoids(string outputDir, string sampleName, string prefix, } // end iVoid closeFiles(fpZobov, fpCenters, fpBarycenter, - fpDistances, fpShapes, fpSkyPositions); + fpShapes, fpSkyPositions); } // end outputVoids diff --git a/examples/example_observation/example_observation.py b/examples/example_observation/example_observation.py index 67c8268..8064484 100644 --- a/examples/example_observation/example_observation.py +++ b/examples/example_observation/example_observation.py @@ -30,7 +30,7 @@ continueRun = False # 1 : extract redshift slices from data # 2 : void extraction using zobov # 3 : removal of small voids and voids near the edge -startCatalogStage = 1 +startCatalogStage = 2 endCatalogStage = 3 basePath = os.path.dirname(os.path.abspath(__file__)) @@ -50,7 +50,7 @@ figDir = os.path.join(workDir,"figs","example_observation") # optimization: maximum number of parallel threads to use numZobovThreads = 2 -# optimization: number of subdivisions of the box +# optimization: number of subdivisions of the volume numZobovDivisions = 2 # Maximum density for merging voids @@ -65,7 +65,7 @@ boundaryTolerance = 1.0 # don't change this dataSampleList = [] -# define your volume-limited samples +# define your data samples newSample = Sample( # path to galaxy file is inputDataDir+dataFile dataFile = "example_observation.dat", @@ -76,19 +76,22 @@ newSample = Sample( # a convenient nickname nickName = "exobs", - # don't change this + # don't change this or nothing will make sense dataType = "observation", # assume sample is volume-limited? volumeLimited = True, # HEALpix mask file - set to None to auto-compute + # NOTE: auto-computed masks are pretty terrible, so + # only do that if you have no other options #maskFile = "", maskFile = inputDataDir+"/example_observation_mask.fits", - # if maskFile blank, desired resolution for HEALpix - # mask mapping, otherwise pulled from maskFile - nsideForMask = 128, + # resolution for HEALpix mapping of survey edge contours + # Set to -1 to use nside from given fits file + # MUST be set if auto-computing mask + nsideForContour = 128, # radial selection function (if not volume limited) selFunFile = None, @@ -109,6 +112,7 @@ newSample = Sample( # density of mock particles in cubic Mpc/h # (make this as high as you can afford) + ### DEPRECATED fakeDensity = 0.05, # if true, convert to comoving space using LCDM cosmology diff --git a/python_source/backend/classes.py b/python_source/backend/classes.py index 6d179f9..347cd42 100644 --- a/python_source/backend/classes.py +++ b/python_source/backend/classes.py @@ -69,7 +69,7 @@ class Sample: nickName = "dim" outputDir = "" maskFile = "rast_window_512.fits" - nsideForMask = 128 + nsideForContour = 128 selFunFile = "czselfunc.all.dr72dim.dat" zBoundary = (0.0, 0.1) zBoundaryMpc = (0., 300) @@ -78,7 +78,8 @@ class Sample: zRange = (0.0, 0.1) omegaM = 0.27 minVoidRadius = -1 - fakeDensity = 0.01 + meanPartSep = 1 # calculated mean particle separation + fakeDensity = 0.01 # TODO - remove hasWeightedVolumes = False profileBinSize = 2 # Mpc autoNumInStack = -1 # set to >0 to automatically generate stacks of size N @@ -101,9 +102,9 @@ class Sample: stacks = [] def __init__(self, dataFile="", fullName="", dataUnit=1, - nickName="", maskFile="", nsideForMask=128, selFunFile="", + nickName="", maskFile="", nsideForContour=128, selFunFile="", zBoundary=(), zRange=(), zBoundaryMpc=(), boundaryWidth=0.1, - shiftSimZ=False, + shiftSimZ=False, meanPartSep = 1, minVoidRadius=-1, fakeDensity=0.01, volumeLimited=True, numAPSlices=1, hasWeightedVolumes=False, includeInHubble=True, partOfCombo=False, isCombo=False, @@ -118,7 +119,7 @@ class Sample: self.fullName = fullName self.nickName = nickName self.maskFile = maskFile - self.nsideForMask = nsideForMask + self.nsideForContour = nsideForContour self.selFunFile = selFunFile self.zBoundary = zBoundary self.zBoundaryMpc = zBoundaryMpc @@ -126,6 +127,7 @@ class Sample: self.shiftSimZ = shiftSimZ self.zRange = zRange self.minVoidRadius = minVoidRadius + self.meanPartSep = meanPartSep self.fakeDensity = fakeDensity self.hasWeightedVolumes = hasWeightedVolumes self.volumeLimited = volumeLimited diff --git a/python_source/backend/launchers.py b/python_source/backend/launchers.py index 0b6cf2b..4ab01e2 100644 --- a/python_source/backend/launchers.py +++ b/python_source/backend/launchers.py @@ -39,17 +39,20 @@ from backend.cosmologyTools import * from backend.surveyTools import * import pickle import scipy.interpolate as interpolate +import time NetCDFFile = Dataset ncFloat = 'f8' # Double precision -LIGHT_SPEED = 299792.458 +#LIGHT_SPEED = 299792.458 # ----------------------------------------------------------------------------- def launchPrep(sample, binPath, workDir=None, inputDataDir=None, outputDir=None, figDir=None, logFile=None, useComoving=False, continueRun=None, regenerate=False): + startTime = time.time() + if sample.dataType == "observation": sampleName = sample.fullName @@ -67,14 +70,29 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, datafile = inputDataDir+"/"+sample.dataFile if sample.maskFile == "": - sample.maskFile = outputDir + "/constructed_mask.fits" - figureOutMask(datafile, sample.nsideForMask, sample.maskFile) + if sample.nsideForContour == -1: + sample.nsideForContour = 128 + sample.maskFile = outputDir + "/constructed_mask.fits" + figureOutMask(datafile, sample.nsideForContour, sample.maskFile) + + # compute mean particle separation + (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0], + sample.zRange[1], sample.zRange[0], sample.zRange[1], "all", + sample.omegaM, useComoving=useComoving) + + numTracers = int(open(outputDir+"/mask_index.txt", "r").read()) + sample.meanPartSep = (1.*numTracers/boxVol/nbar)**(-1/3.) + + + # flag edge galaxies + galFile = outputDir + "galaxies.txt" edgeGalFile = outputDir + "/galaxy_edge_flags.txt" - edgeMaskFile = outputDir + "/mask_edge_map.fits" - findEdgeGalaxies(datafile, sample.maskFile, edgeGalFile, edgeMaskFile, + #edgeMaskFile = outputDir + "/mask_edge_map.fits" + contourFile = outputDir + "/contour_map.fits" + findEdgeGalaxies(galFile, sample.maskFile, edgeGalFile, contourFile, sample.zBoundary[0], sample.zBoundary[1], sample.omegaM, - useComoving, sample.boundaryWidth) + useComoving, sample.boundaryWidth, sample.meanPartSep) if useComoving: useComovingFlag = "useComoving" @@ -92,12 +110,15 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, %s %s omegaM %g + nsideForContour %g + meanPartSep %g """ % (datafile, sample.maskFile, outputFile, outputDir+"/zobov_slice_"+sampleName+".par", sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity, - useComovingFlag, inputParameterFlag, sample.omegaM) + useComovingFlag, inputParameterFlag, sample.omegaM, + sample.nsideForContour, sample.meanPartSep) - parmFile = os.getcwd()+"/generate_"+sample.fullName+".par" + parmFile = os.getcwd()+"/prep_"+sample.fullName+".par" if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")): with open(parmFile, mode="wt") as f: @@ -106,9 +127,11 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, with open(logFile, 'wt') as log: subprocess.call([binPath, arg1], stdout=log, stderr=log) if jobSuccessful(logFile, "Done!\n"): - print("done") + endTime = time.time() + walltime = endTime - startTime + print("done (%.2fs elapsed)" % walltime) else: - print("FAILED!") + print("FAILED! See log file for details.") exit(-1) else: @@ -118,7 +141,6 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, if os.access("contour_map.fits", os.F_OK): os.system("mv %s %s" % ("contour_map.fits", outputDir)) - os.system("mv %s %s" % ("mask_map.fits", outputDir)) if os.access("comoving_distance.txt", os.F_OK): os.system("mv %s %s" % ("comoving_distance.txt", outputDir)) @@ -129,15 +151,26 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, if os.access("galaxies.txt", os.F_OK): os.system("mv %s %s" % ("galaxies.txt", outputDir)) - os.system("mv %s %s" % ("mock_galaxies.txt", outputDir)) - os.system("mv %s %s" % ("mock_boundary.txt", outputDir)) - os.system("mv %s %s" % ("mock_sphere.txt", outputDir)) + #os.system("mv %s %s" % ("galaxy_edge_flags.txt", outputDir)) else: # simulation sampleName = sample.fullName datafile = inputDataDir+"/"+sample.dataFile + # compute mean particle separation + iX = float(sample.mySubvolume[0]) + iY = float(sample.mySubvolume[1]) + xMin = iX/sample.numSubvolumes * sample.boxLen + yMin = iY/sample.numSubvolumes * sample.boxLen + xMax = (iX+1)/sample.numSubvolumes * sample.boxLen + yMax = (iY+1)/sample.numSubvolumes * sample.boxLen + zMin = sample.zBoundaryMpc[0] + zMax = sample.zBoundaryMpc[1] + + boxVol = (xMax-xMin)*(yMax-yMin)*(zMax-zMin) + sample.meanPartSep = (1.*numTracers/boxVol)**(-1/3.) + # check if the final subsampling is done lastSample = sample.subsample.split(', ')[-1] doneLine = "Done! %5.2e\n" % float(lastSample) @@ -245,7 +278,7 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, cmd = "%s --configFile=%s" % (binPath,parmFile) log = open(logFile, 'a') arg1 = "--configFile=%s" % parmFile - + subprocess.call(cmd, stdout=log, stderr=log, shell=True) log.close() @@ -257,7 +290,7 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, doneLine = "Done! %5.2e\n" % keepFraction if not jobSuccessful(logFile, doneLine): - print("FAILED!") ### dies here for now + print("FAILED! See log file for details.") ### dies here for now exit(-1) prevSubSample = thisSubSample @@ -281,29 +314,6 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, os.system("mv %s %s" % ("total_particles.txt", outputDir)) #os.system("mv %s %s" % ("sample_info.txt", outputDir)) - # add to sample info file - if sample.dataType == "observation": - (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0], - sample.zRange[1], sample.zRange[0], sample.zRange[1], "all", - sample.omegaM, useComoving=useComoving) - else: - iX = float(sample.mySubvolume[0]) - iY = float(sample.mySubvolume[1]) - xMin = iX/sample.numSubvolumes * sample.boxLen - yMin = iY/sample.numSubvolumes * sample.boxLen - xMax = (iX+1)/sample.numSubvolumes * sample.boxLen - yMax = (iY+1)/sample.numSubvolumes * sample.boxLen - zMin = sample.zBoundaryMpc[0] - zMax = sample.zBoundaryMpc[1] - - boxVol = (xMax-xMin)*(yMax-yMin)*(zMax-zMin) - nbar = 1.0 - - numTracers = int(open(outputDir+"/mask_index.txt", "r").read()) - numTotal = int(open(outputDir+"/total_particles.txt", "r").read()) - - meanSep = (1.*numTracers/boxVol/nbar)**(-1/3.) - # save this sample's information with open(outputDir+"/sample_info.dat", mode='wb') as output: pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL) @@ -325,9 +335,8 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, fp.write("Number of simulation subvolumes: %s\n" % sample.numSubvolumes) fp.write("My subvolume index: %s\n" % sample.mySubvolume) fp.write("Estimated volume (cubic Mpc/h): %g\n" % boxVol) - fp.write("Number of real (non-boundary) tracers: %d\n" % numTracers) - fp.write("Total number of tracers: %d\n" % numTotal) - fp.write("Estimated mean tracer separation (Mpc/h): %g\n" % meanSep) + fp.write("Total number of tracers: %d\n" % numTracers) + fp.write("Estimated mean tracer separation (Mpc/h): %g\n" % sample.meanPartSep) fp.write("Minimum void size actually used (Mpc/h): %g\n" % sample.minVoidRadius) fp.close() @@ -336,6 +345,8 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None, numZobovDivisions=None, numZobovThreads=None, mergingThreshold=0.2): + startTime = time.time() + sampleName = sample.fullName datafile = outputDir+"zobov_slice_"+sampleName @@ -490,9 +501,11 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None, os.unlink(fileName) if jobSuccessful(logFile, "Done!\n"): - print("done") + endTime = time.time() + walltime = endTime - startTime + print("done (%.2fs elapsed)" % walltime) else: - print("FAILED!") + print("FAILED! See log file for details.") exit(-1) else: @@ -507,6 +520,8 @@ def launchPrune(sample, binPath, continueRun=None, useComoving=False, mergingThreshold=0.2, boundaryTolerance=1.0): + startTime = time.time() + sampleName = sample.fullName numVoids = sum(1 for line in \ @@ -580,9 +595,11 @@ def launchPrune(sample, binPath, if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or \ jobSuccessful(logFile, "Done!\n"): - print("done") + endTime = time.time() + walltime = endTime - startTime + print("done (%.2fs elapsed)" % walltime) else: - print("FAILED!") + print("FAILED! See log file for details.") #exit(-1) else: @@ -597,6 +614,8 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir, overlapFrac=0.25, matchMethod=None, strictMatch=False): + startTime = time.time() + sampleName1 = sample1.fullName sampleName2 = sample2.fullName @@ -663,7 +682,9 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir, log.close() #if jobSuccessful(logFile, "Done!\n"): - print("done") + endTime = time.time() + walltime = endTime - startTime + print("done (%.2fs elapsed)" % walltime) #else: # print "FAILED!" # exit(-1) @@ -707,7 +728,7 @@ def launchVelocityStack(sample, stack, binPath, if jobSuccessful(logFile, "Done!\n"): print("done") else: - print("FAILED!") + print("FAILED! See log file for details.") exit(-1) else: diff --git a/python_source/backend/surveyTools.py b/python_source/backend/surveyTools.py index b9f1423..8c63070 100644 --- a/python_source/backend/surveyTools.py +++ b/python_source/backend/surveyTools.py @@ -21,6 +21,7 @@ # distances, and expected void stretching import numpy as np +import scipy import healpy as healpy import os from backend import * @@ -127,63 +128,111 @@ def figureOutMask(galFile, nside, outMaskFile): return mask # ----------------------------------------------------------------------------- -# figures out which galaxies live on a mask edge, and also writes the edge -# map to an auxillary file -def findEdgeGalaxies(galFile, maskFile, edgeGalFile, edgeMaskFile, - zmin, zmax, omegaM, useComoving, boundaryWidth): +# figures out which galaxies live on a mask or redshift edge +def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile, + zmin, zmax, omegaM, useComoving, boundaryWidth, + meanPartSep): if useComoving: - zmin = comovingDistance(zmin, Om=omegaM) - zmax = comovingDistance(zmax, Om=omegaM) - #zmin = LIGHT_SPEED/100.*comovingDistance(zmin, Om=omegaM) - #zmax = LIGHT_SPEED/100.*comovingDistance(zmax, Om=omegaM) - #else: - # zmin *= LIGHT_SPEED/100. - # zmax *= LIGHT_SPEED/100. + zmin = comovingDistance(zmin, Om=omegaM)*LIGHT_SPEED + zmax = comovingDistance(zmax, Om=omegaM)*LIGHT_SPEED + else: + zmin *= LIGHT_SPEED + zmax *= LIGHT_SPEED - - mask = healpy.read_map(maskFile) - nside = healpy.get_nside(mask) + contourMap = healpy.read_map(contourFile) + nside = healpy.get_nside(contourMap) npix = healpy.nside2npix(nside) - edgeMask = np.zeros((npix)) - edgeFile = open(edgeGalFile, "w") + # load in galaxies + galPos = np.genfromtxt(galFile) + flagList = np.zeros(len(galPos[:,0])) + galTree = scipy.spatial.cKDTree(galPos) - for line in open(galFile): - line = line.split() - RA = float(line[3]) - Dec = float(line[4]) - z = float(line[5]) - - if useComoving: - z = comovingDistance(z/LIGHT_SPEED, Om=omegaM) - else: - z *= LIGHT_SPEED/100. - - phi, theta = convertAngle(RA, Dec) - - # check the mask edges - ipix = healpy.ang2pix(nside, theta, phi) - neighbors = healpy.get_all_neighbours(nside, ipix) - isOnMaskEdge = any(mask[p] == 0 for p in neighbors) + # flag galaxies near mask edges + # using the "ray marching" algorithm: follow rays along lines of sight + # of all mask edges, flagging nearest neighbor galaxies as we go - # check the redshift boundaries - zbuffer = (zmax-zmin)*boundaryWidth - isOnHighZEdge = (z >= zmax-zbuffer) - isOnLowZEdge = (z <= zmin+zbuffer) + raySteps = np.arange(zmin, zmax, meanPartSep) + + contourPixels = np.nonzero(contourMap)[0] + #print(contourPixels) + for pixel in contourPixels: + #print("Working with pixel %d" % pixel) + vec = healpy.pix2vec(nside,pixel) + x = raySteps * vec[0] + y = raySteps * vec[1] + z = raySteps * vec[2] + ray = np.array((x,y,z)).T + #print(ray) + + dist, nearest = galTree.query(ray) + flagList[nearest] = 1 + #print(nearest) - if isOnMaskEdge: - edgeFile.write("1\n") - edgeMask[ipix] = 1 - elif isOnHighZEdge: - edgeFile.write("2\n") - elif isOnLowZEdge: - edgeFile.write("3\n") - else: - edgeFile.write("0\n") + # flag galaxies near redsfhit boundaries + # TODO - save time by only covering portion of sphere with data + ds = np.sqrt(healpy.nside2pixarea(nside)) / 1000. + phi = np.arange(0, 2*np.pi, ds*2) + theta = np.arange(0, np.pi, ds) + vec = healpy.ang2vec(theta, phi) - edgeFile.close() - healpy.write_map(edgeMaskFile, edgeMask, overwrite=True, - dtype=np.dtype('float64')) + maxEdge = zmax * vec + dist, nearest = galTree.query(maxEdge) + #print(nearest) + #print(galPos[nearest]) + flagList[nearest] = 2 + + minEdge = zmin * vec + dist, nearest = galTree.query(minEdge) + #print(nearest) + #print(galPos[nearest]) + flagList[nearest] = 3 + + # output flag information + np.savetxt(edgeGalFile, flagList, fmt="%d") + + + +# # output galaxy edge flags +# edgeFile = open(edgeGalFile, "w") +# +# for line in open(galFile): +# line = line.split() +# RA = float(line[3]) +# Dec = float(line[4]) +# z = float(line[5]) +# +# if useComoving: +# z = comovingDistance(z/LIGHT_SPEED, Om=omegaM) +# else: +# z *= LIGHT_SPEED/100. +# +# phi, theta = convertAngle(RA, Dec) +# +# # check the mask edges +# ipix = healpy.ang2pix(nside, theta, phi) +# neighbors = healpy.get_all_neighbours(nside, ipix) +# isOnMaskEdge = any(mask[p] == 0 for p in neighbors) +# +# # check the redshift boundaries +# zbuffer = (zmax-zmin)*boundaryWidth +# isOnHighZEdge = (z >= zmax-zbuffer) +# isOnLowZEdge = (z <= zmin+zbuffer) +# +# if isOnMaskEdge: +# edgeFile.write("1\n") +# edgeMask[ipix] = 1 +# elif isOnHighZEdge: +# edgeFile.write("2\n") +# elif isOnLowZEdge: +# +#edgeFile.write("3\n") +# else: +# edgeFile.write("0\n") +# +# edgeFile.close() +# healpy.write_map(edgeMaskFile, edgeMask, overwrite=True, +# dtype=np.dtype('float64')) return diff --git a/python_source/vide_pipeline/__main__.py b/python_source/vide_pipeline/__main__.py index 716c1f1..840b4ca 100644 --- a/python_source/vide_pipeline/__main__.py +++ b/python_source/vide_pipeline/__main__.py @@ -36,6 +36,7 @@ if (len(sys.argv) == 1): exit(-1) if (len(sys.argv) > 1): + print("\n\n Welcome to VIDE!\n") filename = sys.argv[1] print(" Loading parameters from", filename) if not os.access(filename, os.F_OK): @@ -63,8 +64,8 @@ if not os.access(figDir, os.F_OK): if not continueRun: print(" Cleaning out log files...") - if startCatalogStage <= 1 and glob.glob(logDir+"/generate*") != []: - os.system("rm %s/generate*" % logDir) + if startCatalogStage <= 1 and glob.glob(logDir+"/prepare*") != []: + os.system("rm %s/prepare*" % logDir) if startCatalogStage <= 2 and glob.glob(logDir+"/zobov*") != []: os.system("rm %s/zobov*" % logDir) if startCatalogStage <= 3 and glob.glob(logDir+"/prune*") != []: @@ -83,24 +84,24 @@ for sample in dataSampleList: # --------------------------------------------------------------------------- if (startCatalogStage <= 1) and (endCatalogStage >= 1) and not sample.isCombo: - print(" Extracting tracers from catalog...", end=' ') + print(" Preparing tracers from catalog...", end='') sys.stdout.flush() - logFile = logDir+"/generate_"+sampleName+".out" + logFile = logDir+"/prepare_"+sampleName+".out" if sample.dataType == "observation": - GENERATE_PATH = CTOOLS_PATH+"/prepObservation" + PREPARE_PATH = CTOOLS_PATH+"/prepObservation" else: - GENERATE_PATH = CTOOLS_PATH+"/prepSimulation" + PREPARE_PATH = CTOOLS_PATH+"/prepSimulation" - launchPrep(sample, GENERATE_PATH, workDir=workDir, + launchPrep(sample, PREPARE_PATH, workDir=workDir, inputDataDir=inputDataDir, outputDir=outputDir, figDir=figDir, logFile=logFile, useComoving=sample.useComoving, continueRun=continueRun, regenerate=regenerateFlag) # -------------------------------------------------------------------------- if (startCatalogStage <= 2) and (endCatalogStage >= 2) and not sample.isCombo: - print(" Extracting voids with ZOBOV...", end=' ') + print(" Finding voids...", end='') sys.stdout.flush() launchZobov(sample, ZOBOV_PATH, outputDir=outputDir, logDir=logDir, @@ -111,7 +112,7 @@ for sample in dataSampleList: # ------------------------------------------------------------------------- if (startCatalogStage <= 3) and (endCatalogStage >= 3) and not sample.isCombo: - print(" Pruning void catalogs", "...", end=' ') + print(" Pruning void catalogs", "...", end='') sys.stdout.flush() logFile = logDir+"/pruneVoids_"+sampleName+".out" @@ -126,7 +127,7 @@ for sample in dataSampleList: # ------------------------------------------------------------------------- if (startCatalogStage <= 4) and (endCatalogStage >= 4): - print(" Plotting...", end=' ') + print(" Plotting...", end='') sys.stdout.flush() #for thisDataPortion in dataPortions: @@ -139,4 +140,4 @@ if (startCatalogStage <= 4) and (endCatalogStage >= 4): #plotVoidDistribution(workDir, dataSampleList, figDir, showPlot=False, # dataPortion=thisDataPortion, setName=setName) -print("\n Done!") +print("\n VIDE finished!")