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!")