From 4df341f9ba1c56d926b7f42a627d73e4c347a2a8 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 28 Nov 2011 15:33:00 -0500 Subject: [PATCH] Added Healpix detection facilities --- mytools/FindHealpix.cmake | 42 +++++++ mytools/contour_pixels.cpp | 6 +- mytools/contour_pixels.hpp | 2 +- mytools/generateFromCatalog.cpp | 208 +++++++++++++++++++++++++++++++- vozutil.c | 2 +- 5 files changed, 252 insertions(+), 8 deletions(-) create mode 100644 mytools/FindHealpix.cmake diff --git a/mytools/FindHealpix.cmake b/mytools/FindHealpix.cmake new file mode 100644 index 0000000..b693ad0 --- /dev/null +++ b/mytools/FindHealpix.cmake @@ -0,0 +1,42 @@ +OPTION(ENABLE_OPENMP "Set to Yes if Healpix and/or you need openMP" OFF) + +# +# OpenMP handling +# + + +INCLUDE(FindOpenMP) + +IF(ENABLE_OPENMP) + + IF (NOT OPENMP_FOUND) + MESSAGE(ERROR "No known compiler option for enabling OpenMP") + ENDIF(NOT OPENMP_FOUND) + + SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKED_FLAGS} ${OpenMP_C_FLAGS}") + +ENDIF(ENABLE_OPENMP) + + + +SET(HEALPIX_BASE_PATH $ENV{HEALPIX}) +SET(HEALPIX_TARGET $ENV{HEALPIX_TARGET}) +SET(HEALPIX_PATH ${HEALPIX_BASE_PATH}/src/cxx/${HEALPIX_TARGET}) +SET(HEALPIX_HINT_INCLUDE ${HEALPIX_PATH}/include) +SET(HEALPIX_HINT_LIB ${HEALPIX_PATH}/lib) + +find_path(HEALPIX_INCLUDE_PATH NAMES healpix_map.h PATHS ${HEALPIX_HINT_INCLUDE}) + +find_library(HEALPIX_LIBRARY healpix_cxx PATHS ${HEALPIX_HINT_LIB}) +find_library(FFTPACK_LIBRARY fftpack PATHS ${HEALPIX_HINT_LIB}) +find_library(CFITSIO_LIBRARY cfitsio PATHS ${HEALPIX_HINT_LIB}) +find_library(CXXSUPPORT_LIBRARY cxxsupport PATHS ${HEALPIX_HINT_LIB}) +find_library(PSHT_LIBRARY psht PATHS ${HEALPIX_HINT_LIB}) +find_library(CUTILS_LIBRARY c_utils PATHS ${HEALPIX_HINT_LIB}) + +SET(HEALPIX_LIBRARIES + ${HEALPIX_LIBRARY} ${EXTRA_HEALPIX_LIBRARIES} + ${FFTPACK_LIBRARY} ${CXXSUPPORT_LIBRARY} ${CFITSIO_LIBRARY} +) \ No newline at end of file diff --git a/mytools/contour_pixels.cpp b/mytools/contour_pixels.cpp index 725b78c..725e6a5 100644 --- a/mytools/contour_pixels.cpp +++ b/mytools/contour_pixels.cpp @@ -7,7 +7,7 @@ using namespace std; static const bool DEBUG = true; -void computeContourPixels(Healpix_Map& m, vector contour) +void computeContourPixels(Healpix_Map& m, vector& contour) { for (int p = 0; p < m.Npix(); p++) { @@ -19,7 +19,7 @@ void computeContourPixels(Healpix_Map& m, vector contour) if (result[q] < 0) continue; - delta = (m[p]-0.5)*(m[result[q]]-0.5); + float delta = (m[p]-0.5)*(m[result[q]]-0.5); if (delta < 0) { contour.push_back(p); @@ -33,7 +33,7 @@ void computeContourPixels(Healpix_Map& m, vector contour) { Healpix_Map contour_map; - contour_map.SetNside(RING, m.Nside()); + contour_map.SetNside(m.Nside(), RING); contour_map.fill(0); for (int p = 0; p < contour.size(); p++) { diff --git a/mytools/contour_pixels.hpp b/mytools/contour_pixels.hpp index 341e409..bdf209c 100644 --- a/mytools/contour_pixels.hpp +++ b/mytools/contour_pixels.hpp @@ -4,6 +4,6 @@ #include #include -void computeContourPixels(Healpix_Map& map, std::vector contour); +void computeContourPixels(Healpix_Map& map, std::vector& contour); #endif diff --git a/mytools/generateFromCatalog.cpp b/mytools/generateFromCatalog.cpp index 7b41380..e3fbfb2 100644 --- a/mytools/generateFromCatalog.cpp +++ b/mytools/generateFromCatalog.cpp @@ -7,9 +7,12 @@ #include #include "generateFromCatalog_conf.h" #include "contour_pixels.hpp" +#include +#include using namespace std; using boost::format; +using namespace CosmoTool; struct NYU_Data { @@ -22,7 +25,23 @@ struct NYU_Data double phi_z; }; -void loadData(const string& fname, vector & data) +struct Position +{ + double xyz[3]; +}; + +struct ParticleData +{ + vector id_gal; + int id_mask; + vector pos; + double box[3][2]; + double Lmax; +}; + +typedef vector NYU_VData; + +void loadData(const string& fname, NYU_VData & data) { ifstream f(fname.c_str()); @@ -34,6 +53,173 @@ void loadData(const string& fname, vector & data) } } +void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data) +{ + double d2r = M_PI/180; + + output_data.pos.resize(data.size()); + output_data.id_gal.resize(data.size()); + + for (int j = 0; j < 3; j++) + { + output_data.box[j][0] = -INFINITY; + output_data.box[j][1] = INFINITY; + } + + for (int i = 0; i < data.size(); i++) + { + double ra = data[i].ra*d2r, dec = data[i].dec*d2r; + Position& p = output_data.pos[i]; + + 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); + output_data.id_gal[i] = data[i].index; + + 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]; + } + } + 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; + +} + +void generateSurfaceMask(generateFromCatalog_info& args , + Healpix_Map& mask, + vector& pixel_list, + NYU_VData& data, + ParticleData& output_data) +{ + // 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; + + 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; + + + 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); + cout << format("3d volume to fill: %g (Mpc/h)^3") % volume << endl; + + cout << format("Will insert %d particles") % numToInsert << endl; + + 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); + + output_data.pos.push_back(p); + output_data.id_gal.push_back(idx); + insertion++; + } + cout << format("Done. Inserted %d particles.") % insertion << endl; +} + +void saveData(ParticleData& pdata) +{ + NcFile f("particles.nc", NcFile::Replace); + + assert(f.is_valid()); + + NcDim *d = f.add_dim("space", 3); + NcDim *p = f.add_dim("Np", pdata.pos.size()); + NcVar *v = f.add_var("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->put_rec(d, x, j); + } + + v = f.add_var("id_gal", ncInt, p); + v->put(&pdata.id_gal[0], pdata.id_gal.size()); + + delete[] x; + +} + +void saveForZobov(ParticleData& pdata, const string& fname) +{ + UnformattedWrite f(fname); + static const char axis[] = { 'X', 'Y', 'Z' }; + double Lmax = pdata.Lmax; + + f.beginCheckpoint(); + f.writeInt32(pdata.pos.size()); + f.endCheckpoint(); + + for (int j = 0; j < 3; j++) + { + cout << format("Writing %c components...") % axis[j] << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < pdata.pos.size(); i++) + { + f.writeReal32((pdata.pos[i].xyz[j]+Lmax)/(2*Lmax)); + } + f.endCheckpoint(); + } +} int main(int argc, char **argv) { @@ -66,14 +252,30 @@ int main(int argc, char **argv) cout << "Loading NYU data..." << endl; vector data; - Healpix_Map mask; + Healpix_Map o_mask; vector pixel_list; + ParticleData output_data; + loadData(args_info.catalog_arg, data); + cout << "Loading mask..." << endl; - read_Healpix_map_from_fits(args_info.mask_arg, mask); + read_Healpix_map_from_fits(args_info.mask_arg, o_mask); + + Healpix_Map mask; + + mask.SetNside(128, RING); + mask.Import(o_mask); computeContourPixels(mask,pixel_list); + // We compute a cube holding all the galaxies + the survey surface mask + + generateGalaxiesInCube(data, output_data); + generateSurfaceMask(args_info, mask, pixel_list, data, output_data); + + saveForZobov(output_data, args_info.output_arg); + // saveData(output_data); + return 0; } diff --git a/vozutil.c b/vozutil.c index 3ab8fb8..ed60125 100644 --- a/vozutil.c +++ b/vozutil.c @@ -44,7 +44,7 @@ int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) } /* Delaunay triangulation*/ - sprintf (flags, "qhull s d"); + sprintf (flags, "qhull s d QJ"); exitcode= qh_new_qhull (dim, nvpall, points, ismalloc, flags, outfile, errfile);