From 9274da1587f5c3fbf20dec5f087c87cda2f67fdd Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 15 Jan 2012 15:15:30 -0500 Subject: [PATCH 1/7] Add more mask --- mytools/contour_pixels.cpp | 9 +++++ mytools/contour_pixels.hpp | 2 ++ mytools/generateFromCatalog.cpp | 61 +++++++++++++++++++++++++++++++++ mytools/generateFromCatalog.ggo | 5 ++- voztie.c | 6 ++-- 5 files changed, 79 insertions(+), 4 deletions(-) diff --git a/mytools/contour_pixels.cpp b/mytools/contour_pixels.cpp index 725e6a5..6486ae7 100644 --- a/mytools/contour_pixels.cpp +++ b/mytools/contour_pixels.cpp @@ -7,8 +7,17 @@ using namespace std; static const bool DEBUG = true; +void computeFilledPixels(Healpix_Map& m, vector& filled) +{ + filled.clear(); + for (int p = 0; p < m.Npix(); p++) + if (m[p] > 0) + filled.push_back(p); +} + void computeContourPixels(Healpix_Map& m, vector& contour) { + contour.clear(); for (int p = 0; p < m.Npix(); p++) { fix_arr result; diff --git a/mytools/contour_pixels.hpp b/mytools/contour_pixels.hpp index bdf209c..d3a893b 100644 --- a/mytools/contour_pixels.hpp +++ b/mytools/contour_pixels.hpp @@ -4,6 +4,8 @@ #include #include + void computeContourPixels(Healpix_Map& map, std::vector& contour); +void computeFilledPixels(Healpix_Map& map, std::vector& contour); #endif diff --git a/mytools/generateFromCatalog.cpp b/mytools/generateFromCatalog.cpp index c41f150..62279ef 100644 --- a/mytools/generateFromCatalog.cpp +++ b/mytools/generateFromCatalog.cpp @@ -92,6 +92,65 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data) } +static double cube(double x) +{ + return x*x*x; +} + +void generateBoxMask(generateFromCatalog_info& args , + Healpix_Map& mask, + vector& pixel_list, + NYU_VData& data, + ParticleData& output_data) +{ + int idx = -1; + int insertion = 0; + double volume = pixel_list.size()*1.0/mask.Npix()*4*M_PI; + int numToInsert; + + idx = output_data.id_mask; + + cout << "Generate box mask..." << endl; + + double Rmax = output_data.Lmax, t = args.box_thickness_arg; + output_data.Lmax += args.box_thickness_arg; + + + volume *= Rmax*Rmax/1e6 * args.box_thickness_arg; + numToInsert = (int)floor(volume*args.box_density_arg); + + for (int i = 0; i < numToInsert; i++) + { + 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()*(cube(1+t/Rmax)-1) + 1,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++; + } +} + void generateSurfaceMask(generateFromCatalog_info& args , Healpix_Map& mask, vector& pixel_list, @@ -296,6 +355,8 @@ int main(int argc, char **argv) generateGalaxiesInCube(data, output_data); generateSurfaceMask(args_info, mask, pixel_list, data, output_data); + computeFilledPixels(mask,pixel_list); + generateBoxMask(args_info, mask, pixel_list, data, output_data); saveForZobov(output_data, args_info.output_arg, args_info.params_arg); // saveData(output_data); diff --git a/mytools/generateFromCatalog.ggo b/mytools/generateFromCatalog.ggo index bda567d..2a41f4f 100644 --- a/mytools/generateFromCatalog.ggo +++ b/mytools/generateFromCatalog.ggo @@ -8,4 +8,7 @@ option "mask" - "Healpix mask of unobserved data (in Equatorial coordinates)" st option "density_fake" - "Number density of boundary fake tracers (1 h^3/ Mpc^3)" double optional default="1" option "output" - "Filename of particle datafile" string required -option "params" - "Output parameters of the datacube" string required \ No newline at end of file +option "params" - "Output parameters of the datacube" string required + +option "box_thickness" - "Thickness of the limiting mask" double optional default="1" +option "box_density" - "Fake Galaxy density on the box" double optional default="1" diff --git a/voztie.c b/voztie.c index 242bf88..d779989 100644 --- a/voztie.c +++ b/voztie.c @@ -94,9 +94,9 @@ int main(int argc, char *argv[]) { fread(&volstemp,1,sizeof(float),part); if (vols[orig[p]] > -1.) if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) { - printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", - orig[p],vols[orig[p]],volstemp); - exit(0); +// printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", +// orig[p],vols[orig[p]],volstemp); +// exit(0); } vols[orig[p]] = volstemp; } From 8850969187b18ca759201631190ca9f001d4edc4 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 15 Jan 2012 16:39:32 -0500 Subject: [PATCH 2/7] Hack the voz1b1 mesh --- mytools/hack_mesh.cpp | 88 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100644 mytools/hack_mesh.cpp diff --git a/mytools/hack_mesh.cpp b/mytools/hack_mesh.cpp new file mode 100644 index 0000000..6cec2e1 --- /dev/null +++ b/mytools/hack_mesh.cpp @@ -0,0 +1,88 @@ +#include +#include +#include +#include +#include + +using namespace std; + +int main(int argc, char **argv) +{ + char *n = argv[1]; + char *n2 = argv[2]; + ifstream f(n); + ofstream f2(n2); + int np, nvp; + int *ids; + float *vols; + set mask_id; + + NcFile fp("params.nc"); + NcVar *v_id = fp.get_var("particle_ids"); + long int *e_N = v_id->edges(); + int *pid = new int[e_N[0]]; + + v_id->get(pid, e_N[0]); + + for (int i = 0; i < e_N[0]; i++) + mask_id.insert(pid[i]); + delete[] pid; + delete[] e_N; + + f.read((char*)&np, sizeof(int)); + f.read((char*)&nvp,sizeof(int)); + f2.write((char*)&np, sizeof(int)); + + ids = new int[nvp]; + vols = new float[nvp]; + f.read((char*)ids, sizeof(int)*nvp); + f.read((char *)vols, sizeof(float)*nvp); + + vector< vector > adjs; + vector adj0, new_ids; + vector new_vols; + for (int i = 0; i < nvp; i++) + { + int nadj; + bool masked_particle = false; + + + if (mask_id.find(ids[i]) != mask_id.end()) + { + masked_particle = true; + } + + f.read((char*)&nadj, sizeof(int)); + adj0.resize(nadj); + + if (nadj != 0) + { + f.read((char*)&adj0[0], sizeof(float)*nadj); + for (int j = 0; j < nadj; j++) + if (mask_id.find(adj0[j]) != mask_id.end()) + masked_particle = true; + } + + if (!masked_particle) + { + new_vols.push_back(vols[i]); + new_ids.push_back(ids[i]); + adjs.push_back(adj0); + } + } + + nvp = new_ids.size(); + f2.write((char*)&nvp, sizeof(int)); + f2.write((char*)&new_ids[0], sizeof(int)*nvp); + f2.write((char*)&new_vols[0], sizeof(float)*nvp); + + for (int i = 0; i < nvp; i++) + { + int nadj = adjs[i].size(); + f2.write((char*)&nadj, sizeof(int)); + if (nadj != 0) + f2.write((char*)&adjs[i], sizeof(int)*nadj); + } + + return 0; +} From ed7b2c3332d3f465d1577591f9a444be30cd9786 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 15 Jan 2012 16:39:37 -0500 Subject: [PATCH 3/7] Hack the voz1b1 mesh From 863ad3be410c9ee3b2e8472afd65728f902ebab4 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 15 Jan 2012 16:40:30 -0500 Subject: [PATCH 4/7] Revert "Add more mask" This reverts commit 9274da1587f5c3fbf20dec5f087c87cda2f67fdd. --- mytools/contour_pixels.cpp | 9 ----- mytools/contour_pixels.hpp | 2 -- mytools/generateFromCatalog.cpp | 61 --------------------------------- mytools/generateFromCatalog.ggo | 5 +-- voztie.c | 6 ++-- 5 files changed, 4 insertions(+), 79 deletions(-) diff --git a/mytools/contour_pixels.cpp b/mytools/contour_pixels.cpp index 6486ae7..725e6a5 100644 --- a/mytools/contour_pixels.cpp +++ b/mytools/contour_pixels.cpp @@ -7,17 +7,8 @@ using namespace std; static const bool DEBUG = true; -void computeFilledPixels(Healpix_Map& m, vector& filled) -{ - filled.clear(); - for (int p = 0; p < m.Npix(); p++) - if (m[p] > 0) - filled.push_back(p); -} - void computeContourPixels(Healpix_Map& m, vector& contour) { - contour.clear(); for (int p = 0; p < m.Npix(); p++) { fix_arr result; diff --git a/mytools/contour_pixels.hpp b/mytools/contour_pixels.hpp index d3a893b..bdf209c 100644 --- a/mytools/contour_pixels.hpp +++ b/mytools/contour_pixels.hpp @@ -4,8 +4,6 @@ #include #include - void computeContourPixels(Healpix_Map& map, std::vector& contour); -void computeFilledPixels(Healpix_Map& map, std::vector& contour); #endif diff --git a/mytools/generateFromCatalog.cpp b/mytools/generateFromCatalog.cpp index 62279ef..c41f150 100644 --- a/mytools/generateFromCatalog.cpp +++ b/mytools/generateFromCatalog.cpp @@ -92,65 +92,6 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data) } -static double cube(double x) -{ - return x*x*x; -} - -void generateBoxMask(generateFromCatalog_info& args , - Healpix_Map& mask, - vector& pixel_list, - NYU_VData& data, - ParticleData& output_data) -{ - int idx = -1; - int insertion = 0; - double volume = pixel_list.size()*1.0/mask.Npix()*4*M_PI; - int numToInsert; - - idx = output_data.id_mask; - - cout << "Generate box mask..." << endl; - - double Rmax = output_data.Lmax, t = args.box_thickness_arg; - output_data.Lmax += args.box_thickness_arg; - - - volume *= Rmax*Rmax/1e6 * args.box_thickness_arg; - numToInsert = (int)floor(volume*args.box_density_arg); - - for (int i = 0; i < numToInsert; i++) - { - 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()*(cube(1+t/Rmax)-1) + 1,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++; - } -} - void generateSurfaceMask(generateFromCatalog_info& args , Healpix_Map& mask, vector& pixel_list, @@ -355,8 +296,6 @@ int main(int argc, char **argv) generateGalaxiesInCube(data, output_data); generateSurfaceMask(args_info, mask, pixel_list, data, output_data); - computeFilledPixels(mask,pixel_list); - generateBoxMask(args_info, mask, pixel_list, data, output_data); saveForZobov(output_data, args_info.output_arg, args_info.params_arg); // saveData(output_data); diff --git a/mytools/generateFromCatalog.ggo b/mytools/generateFromCatalog.ggo index 2a41f4f..bda567d 100644 --- a/mytools/generateFromCatalog.ggo +++ b/mytools/generateFromCatalog.ggo @@ -8,7 +8,4 @@ option "mask" - "Healpix mask of unobserved data (in Equatorial coordinates)" st option "density_fake" - "Number density of boundary fake tracers (1 h^3/ Mpc^3)" double optional default="1" option "output" - "Filename of particle datafile" string required -option "params" - "Output parameters of the datacube" string required - -option "box_thickness" - "Thickness of the limiting mask" double optional default="1" -option "box_density" - "Fake Galaxy density on the box" double optional default="1" +option "params" - "Output parameters of the datacube" string required \ No newline at end of file diff --git a/voztie.c b/voztie.c index d779989..242bf88 100644 --- a/voztie.c +++ b/voztie.c @@ -94,9 +94,9 @@ int main(int argc, char *argv[]) { fread(&volstemp,1,sizeof(float),part); if (vols[orig[p]] > -1.) if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) { -// printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", -// orig[p],vols[orig[p]],volstemp); -// exit(0); + printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", + orig[p],vols[orig[p]],volstemp); + exit(0); } vols[orig[p]] = volstemp; } From 84b451d3ad97bb191cc2666dbf5e3d03e80270f7 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 15 Jan 2012 16:41:12 -0500 Subject: [PATCH 5/7] Wrong revert Revert "Revert "Add more mask"" This reverts commit 863ad3be410c9ee3b2e8472afd65728f902ebab4. --- mytools/contour_pixels.cpp | 9 +++++ mytools/contour_pixels.hpp | 2 ++ mytools/generateFromCatalog.cpp | 61 +++++++++++++++++++++++++++++++++ mytools/generateFromCatalog.ggo | 5 ++- voztie.c | 6 ++-- 5 files changed, 79 insertions(+), 4 deletions(-) diff --git a/mytools/contour_pixels.cpp b/mytools/contour_pixels.cpp index 725e6a5..6486ae7 100644 --- a/mytools/contour_pixels.cpp +++ b/mytools/contour_pixels.cpp @@ -7,8 +7,17 @@ using namespace std; static const bool DEBUG = true; +void computeFilledPixels(Healpix_Map& m, vector& filled) +{ + filled.clear(); + for (int p = 0; p < m.Npix(); p++) + if (m[p] > 0) + filled.push_back(p); +} + void computeContourPixels(Healpix_Map& m, vector& contour) { + contour.clear(); for (int p = 0; p < m.Npix(); p++) { fix_arr result; diff --git a/mytools/contour_pixels.hpp b/mytools/contour_pixels.hpp index bdf209c..d3a893b 100644 --- a/mytools/contour_pixels.hpp +++ b/mytools/contour_pixels.hpp @@ -4,6 +4,8 @@ #include #include + void computeContourPixels(Healpix_Map& map, std::vector& contour); +void computeFilledPixels(Healpix_Map& map, std::vector& contour); #endif diff --git a/mytools/generateFromCatalog.cpp b/mytools/generateFromCatalog.cpp index c41f150..62279ef 100644 --- a/mytools/generateFromCatalog.cpp +++ b/mytools/generateFromCatalog.cpp @@ -92,6 +92,65 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data) } +static double cube(double x) +{ + return x*x*x; +} + +void generateBoxMask(generateFromCatalog_info& args , + Healpix_Map& mask, + vector& pixel_list, + NYU_VData& data, + ParticleData& output_data) +{ + int idx = -1; + int insertion = 0; + double volume = pixel_list.size()*1.0/mask.Npix()*4*M_PI; + int numToInsert; + + idx = output_data.id_mask; + + cout << "Generate box mask..." << endl; + + double Rmax = output_data.Lmax, t = args.box_thickness_arg; + output_data.Lmax += args.box_thickness_arg; + + + volume *= Rmax*Rmax/1e6 * args.box_thickness_arg; + numToInsert = (int)floor(volume*args.box_density_arg); + + for (int i = 0; i < numToInsert; i++) + { + 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()*(cube(1+t/Rmax)-1) + 1,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++; + } +} + void generateSurfaceMask(generateFromCatalog_info& args , Healpix_Map& mask, vector& pixel_list, @@ -296,6 +355,8 @@ int main(int argc, char **argv) generateGalaxiesInCube(data, output_data); generateSurfaceMask(args_info, mask, pixel_list, data, output_data); + computeFilledPixels(mask,pixel_list); + generateBoxMask(args_info, mask, pixel_list, data, output_data); saveForZobov(output_data, args_info.output_arg, args_info.params_arg); // saveData(output_data); diff --git a/mytools/generateFromCatalog.ggo b/mytools/generateFromCatalog.ggo index bda567d..2a41f4f 100644 --- a/mytools/generateFromCatalog.ggo +++ b/mytools/generateFromCatalog.ggo @@ -8,4 +8,7 @@ option "mask" - "Healpix mask of unobserved data (in Equatorial coordinates)" st option "density_fake" - "Number density of boundary fake tracers (1 h^3/ Mpc^3)" double optional default="1" option "output" - "Filename of particle datafile" string required -option "params" - "Output parameters of the datacube" string required \ No newline at end of file +option "params" - "Output parameters of the datacube" string required + +option "box_thickness" - "Thickness of the limiting mask" double optional default="1" +option "box_density" - "Fake Galaxy density on the box" double optional default="1" diff --git a/voztie.c b/voztie.c index 242bf88..d779989 100644 --- a/voztie.c +++ b/voztie.c @@ -94,9 +94,9 @@ int main(int argc, char *argv[]) { fread(&volstemp,1,sizeof(float),part); if (vols[orig[p]] > -1.) if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) { - printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", - orig[p],vols[orig[p]],volstemp); - exit(0); +// printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", +// orig[p],vols[orig[p]],volstemp); +// exit(0); } vols[orig[p]] = volstemp; } From 63989236439435ce4474658ea5266398ed518ccf Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 15 Jan 2012 16:41:31 -0500 Subject: [PATCH 6/7] Wrong branch Revert "Hack the voz1b1 mesh" This reverts commit 8850969187b18ca759201631190ca9f001d4edc4. --- mytools/hack_mesh.cpp | 88 ------------------------------------------- 1 file changed, 88 deletions(-) delete mode 100644 mytools/hack_mesh.cpp diff --git a/mytools/hack_mesh.cpp b/mytools/hack_mesh.cpp deleted file mode 100644 index 6cec2e1..0000000 --- a/mytools/hack_mesh.cpp +++ /dev/null @@ -1,88 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; - -int main(int argc, char **argv) -{ - char *n = argv[1]; - char *n2 = argv[2]; - ifstream f(n); - ofstream f2(n2); - int np, nvp; - int *ids; - float *vols; - set mask_id; - - NcFile fp("params.nc"); - NcVar *v_id = fp.get_var("particle_ids"); - long int *e_N = v_id->edges(); - int *pid = new int[e_N[0]]; - - v_id->get(pid, e_N[0]); - - for (int i = 0; i < e_N[0]; i++) - mask_id.insert(pid[i]); - delete[] pid; - delete[] e_N; - - f.read((char*)&np, sizeof(int)); - f.read((char*)&nvp,sizeof(int)); - f2.write((char*)&np, sizeof(int)); - - ids = new int[nvp]; - vols = new float[nvp]; - f.read((char*)ids, sizeof(int)*nvp); - f.read((char *)vols, sizeof(float)*nvp); - - vector< vector > adjs; - vector adj0, new_ids; - vector new_vols; - for (int i = 0; i < nvp; i++) - { - int nadj; - bool masked_particle = false; - - - if (mask_id.find(ids[i]) != mask_id.end()) - { - masked_particle = true; - } - - f.read((char*)&nadj, sizeof(int)); - adj0.resize(nadj); - - if (nadj != 0) - { - f.read((char*)&adj0[0], sizeof(float)*nadj); - for (int j = 0; j < nadj; j++) - if (mask_id.find(adj0[j]) != mask_id.end()) - masked_particle = true; - } - - if (!masked_particle) - { - new_vols.push_back(vols[i]); - new_ids.push_back(ids[i]); - adjs.push_back(adj0); - } - } - - nvp = new_ids.size(); - f2.write((char*)&nvp, sizeof(int)); - f2.write((char*)&new_ids[0], sizeof(int)*nvp); - f2.write((char*)&new_vols[0], sizeof(float)*nvp); - - for (int i = 0; i < nvp; i++) - { - int nadj = adjs[i].size(); - f2.write((char*)&nadj, sizeof(int)); - if (nadj != 0) - f2.write((char*)&adjs[i], sizeof(int)*nadj); - } - - return 0; -} From ee18ace00020141464c469f8f081a837b4a43df4 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 15 Jan 2012 16:41:42 -0500 Subject: [PATCH 7/7] Wrong branch Revert "Hack the voz1b1 mesh" This reverts commit ed7b2c3332d3f465d1577591f9a444be30cd9786.