From b1962100a8b80a04da581026e4ea9da76e81dbe7 Mon Sep 17 00:00:00 2001 From: "Paul M. Sutter" Date: Wed, 22 May 2024 15:13:55 -0400 Subject: [PATCH] re-added some essential routines; updated example sceipts to reflect new layout --- c_tools/jozov2/CMakeLists.txt | 26 + c_tools/jozov2/findrtop.c | 81 ++ c_tools/jozov2/jozov2.cpp | 174 +++ c_tools/jozov2/jozov2.hpp | 65 + c_tools/jozov2/jozov2_io.cpp | 199 +++ c_tools/jozov2/jozov2_watershed.cpp | 294 +++++ c_tools/jozov2/jozov2_zones.cpp | 243 ++++ c_tools/jozov2/zobov.hpp | 52 + c_tools/pruning/CMakeLists.txt | 8 + c_tools/pruning/pruneVoids.cpp | 1213 +++++++++++++++++++ c_tools/pruning/pruneVoids.ggo | 44 + python_tools/vide/backend/cosmologyTools.py | 112 ++ 12 files changed, 2511 insertions(+) create mode 100644 c_tools/jozov2/CMakeLists.txt create mode 100644 c_tools/jozov2/findrtop.c create mode 100644 c_tools/jozov2/jozov2.cpp create mode 100644 c_tools/jozov2/jozov2.hpp create mode 100644 c_tools/jozov2/jozov2_io.cpp create mode 100644 c_tools/jozov2/jozov2_watershed.cpp create mode 100644 c_tools/jozov2/jozov2_zones.cpp create mode 100644 c_tools/jozov2/zobov.hpp create mode 100644 c_tools/pruning/CMakeLists.txt create mode 100644 c_tools/pruning/pruneVoids.cpp create mode 100644 c_tools/pruning/pruneVoids.ggo create mode 100644 python_tools/vide/backend/cosmologyTools.py diff --git a/c_tools/jozov2/CMakeLists.txt b/c_tools/jozov2/CMakeLists.txt new file mode 100644 index 0000000..33deb8a --- /dev/null +++ b/c_tools/jozov2/CMakeLists.txt @@ -0,0 +1,26 @@ +function(omp_add_flags TARGET LANG) + if(OPENMP_FOUND) + if(NOT LANG) + get_target_property(LANG ${TARGET} LINKER_LANGUAGE) + endif() + if(LANG MATCHES "C(XX)?") + set_property(TARGET ${TARGET} APPEND + PROPERTY COMPILE_FLAGS ${OpenMP_${LANG}_FLAGS}) + set_property(TARGET ${TARGET} APPEND + PROPERTY LINK_FLAGS ${OpenMP_${LANG}_FLAGS}) + else() + message(WARNING "omp_add_flags: target '${TARGET}'" + " link language '${LANG}' is not supported") + endif() + endif() +endfunction() + +add_executable(jozov2 jozov2.cpp jozov2_io.cpp jozov2_zones.cpp jozov2_watershed.cpp findrtop.c) + +if (ENABLE_OPENMP) + Find_Package(OpenMP) + omp_add_flags(jozov2 CXX) + add_definitions(-DOPENMP) +ENDIF(ENABLE_OPENMP) + +install(TARGETS jozov2 DESTINATION ${VIDE_BIN}) diff --git a/c_tools/jozov2/findrtop.c b/c_tools/jozov2/findrtop.c new file mode 100644 index 0000000..8dbeeb3 --- /dev/null +++ b/c_tools/jozov2/findrtop.c @@ -0,0 +1,81 @@ +#include +#include + +/*------------------------------------------------------------------------------ + Find nb elements of Real array a having the largest value. + Returns index iord of these elements, ordered so iord[0] corresponds + to element a[iord[0]] having the largest value. + If nb > na, then the last nb-na elements of iord are undefined. + Elements of a that are equal are left in their original order. + + Courtesy Andrew Hamilton. +*/ +void findrtop(double *a, int na, int *iord, int nb) +{ +#undef ORDER +#define ORDER(ia, ja) a[ia] > a[ja] || (a[ia] == a[ja] && ia < ja) + + int i, ia, ib, it, n; + + n = (na <= nb)? na : nb; + if (n <= 0) return; + + /* heap first n elements, so smallest element is at top of heap */ + for (ib = (n >> 1); ib < n; ib++) { + iord[ib] = ib; + } + for (ia = (n >> 1) - 1; ia >= 0; ia--) { + i = ia; + for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) { + if (ib+1 < n) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(ia, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = ia; + } + + /* now compare rest of elements of array to heap */ + for (ia = n; ia < na; ia++) { + /* if new element is greater than smallest, sift it into heap */ + i = 0; + if (ORDER(ia, iord[i])) { + for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) { + if (ib+1 < n) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(ia, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = ia; + } + } + + /* unheap iord so largest element is at top */ + for (ia = n - 1; ia > 0; ia--) { + it = iord[ia]; + i = 0; + iord[ia] = iord[i]; + for (ib = (i << 1) + 1; ib < ia; ib = (i << 1) + 1) { + if (ib+1 < ia) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(it, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = it; + } +} diff --git a/c_tools/jozov2/jozov2.cpp b/c_tools/jozov2/jozov2.cpp new file mode 100644 index 0000000..f4cb142 --- /dev/null +++ b/c_tools/jozov2/jozov2.cpp @@ -0,0 +1,174 @@ +/*+ + VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2.cpp + Copyright (C) 2010-2014 Guilhem Lavaux + Copyright (C) 2011-2014 P. M. Sutter + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; version 2 of the License. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. ++*/ +#include +#include +#include +#include +#include +#include +#include +#include +/* jovoz.c by Mark Neyrinck */ +#include "jozov2.hpp" +#include "zobov.hpp" + +using namespace std; +using boost::format; + + + +int main(int argc,char **argv) +{ + ofstream txt; + PARTICLE *p; + pid_t np; + + ZONE *z; + int nzones; + + string adjfile, volfile, zonfile, zonfile2, txtfile; + int *jumped; + int *zonenum; + float voltol, prob; + + float maxvol, minvol; + double *sorter; + int *iord; + + int mockIndex; + + if (argc != 8) { + printf("Wrong number of arguments.\n"); + printf("arg1: adjacency file\n"); + printf("arg2: volume file\n"); + printf("arg3: output file containing particles in each zone\n"); + printf("arg4: output file containing zones in each void\n"); + printf("arg5: output text file\n"); + printf("arg6: Density threshold (0 for no threshold)\n"); + printf("arg7: Beginning index of mock galaxies\n\n"); + exit(0); + } + adjfile = argv[1]; + volfile = argv[2]; + zonfile = argv[3]; + zonfile2 = argv[4]; + txtfile = argv[5]; + if (sscanf(argv[6],"%f",&voltol) == 0) { + printf("Bad density threshold.\n"); + exit(0); + } + if (sscanf(argv[7],"%d",&mockIndex) == 0) { + printf("Bad mock galaxy index.\n"); + exit(0); + } + printf("TOLERANCE: %f\n", voltol); + if (voltol <= 0.) { + printf("Proceeding without a density threshold.\n"); + voltol = 1e30; + } + + try + { + readAdjacencyFile(adjfile, p, np); + } + catch(const FileError& e) + { + return 1; + } + if (mockIndex < 0) + mockIndex = np; + + /* Check that we got all the pairs */ + for (int i = 0; i < np; i++) + { + if (p[i].ncnt != p[i].nadj && i < mockIndex) { + cout + << format("We didn't get all of %d's adj's; %d != %d.") + % i % p[i].ncnt % p[i].nadj + << endl; + p[i].nadj = p[i].ncnt; + } + } + + /* Volumes */ + try + { + readVolumeFile(volfile, p, np, mockIndex); + } + catch (const FileError& e) + { + return 2; + } + + buildZones(p, np, jumped, z, nzones, zonenum); + writeZoneFile(zonfile, p, np, z, nzones, zonenum, jumped); + + maxvol = 0.; + minvol = BIGFLT; + for(pid_t i = 0; i < np; i++){ + if (p[i].dens > maxvol) + maxvol = p[i].dens; + if (p[i].dens < minvol) + minvol = p[i].dens; + } + cout << format("Densities range from %e to %e.") % minvol % maxvol << endl; + FF; + + doWatershed(p, np, z, nzones, maxvol, voltol); + writeVoidFile(zonfile2, z, nzones); + + sorter = new double[nzones+1]; + /* Assign sorter by probability (could use volume instead) */ + for (int h = 0; h < nzones; h++) + sorter[h] = (double)z[h].denscontrast; + + /* Text output file */ + + printf("about to sort ...\n");FF; + + iord = new int[nzones]; + + findrtop(sorter, nzones, iord, nzones); + delete[] sorter; + + txt.open(txtfile.c_str()); + txt << format("%d particles, %d voids.") % np % nzones << endl; + txt << "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb" << endl; + for (int h=0; h +#include +#include "zobov.hpp" + +#define BIGFLT 1e30 /* Biggest possible floating-point number */ +#define NLINKS 1000 /* Number of possible links with the same rho_sl */ +#define FF cout.flush() + +class FileError: virtual std::exception +{ +}; + +void readAdjacencyFile(const std::string& adjfile, PARTICLE*& p, pid_t& np) + throw(FileError); + +void readVolumeFile(const std::string& volfile, PARTICLE *p, pid_t np, + pid_t mockIndex) + throw(FileError); + +void buildInitialZones(PARTICLE *p, pid_t np, pid_t* jumped, + pid_t *numinh, pid_t& numZones); + +void buildZoneAdjacencies(PARTICLE *p, pid_t np, + ZONE *z, ZONET *zt, + int numZones, + pid_t *jumped, + int *zonenum, + int *numinh); + +void buildZones(PARTICLE *p, pid_t np, pid_t *&jumped, + ZONE*& z, int& nzones, + int*& zonenum); + +void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol); + +void writeZoneFile(const std::string& zonfile, PARTICLE* p, pid_t np, + ZONE *z, int numZones, int* zonenum, int *jumped); + +void writeVoidFile(const std::string& zonfile2, ZONE *z, int numZones); + + +extern "C" void findrtop(double *a, int na, int *iord, int nb); + +#endif diff --git a/c_tools/jozov2/jozov2_io.cpp b/c_tools/jozov2/jozov2_io.cpp new file mode 100644 index 0000000..cc6abf7 --- /dev/null +++ b/c_tools/jozov2/jozov2_io.cpp @@ -0,0 +1,199 @@ +/*+ + VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_io.cpp + Copyright (C) 2010-2014 Guilhem Lavaux + Copyright (C) 2011-2014 P. M. Sutter + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; version 2 of the License. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. ++*/ +#include +#include +#include +#include +#include "jozov2.hpp" +#include "zobov.hpp" + +using namespace std; +using boost::format; + +void readAdjacencyFile(const string& adjfile, PARTICLE*& p, pid_t& np) + throw(FileError) +{ + ifstream adj(adjfile.c_str()); + + if (!adj) { + cout << format("Unable to open %s") % adjfile << endl; + throw FileError(); + } + adj.read((char *)&np, sizeof(pid_t)); + + cout << format("adj: %d particles") % np << endl; + FF; + + p = new PARTICLE[np]; + + /* Adjacencies*/ + for (pid_t i=0; i < np; i++) { + adj.read((char*)&p[i].nadj, sizeof(pid_t)); + if (!adj) + throw FileError(); + + /* The number of adjacencies per particle */ + if (p[i].nadj > 0) + p[i].adj = new pid_t[p[i].nadj]; + else + p[i].adj = 0; + p[i].ncnt = 0; /* Temporarily, it's an adj counter */ + } + + for (pid_t i = 0; i < np; i++) + { + pid_t nin; + adj.read((char*)&nin, sizeof(pid_t)); + if (!adj) + throw FileError(); + + if (nin > 0) + { + for (int k=0; k < nin; k++) { + pid_t j; + + adj.read((char *)&j,sizeof(pid_t)); + if (j < np) + { + /* Set both halves of the pair */ + assert(i < j); + if (p[i].ncnt == p[i].nadj) + { + cout << format("OVERFLOW for particle %d (pending %d). List of accepted:") % i % j << endl; + for (int q = 0; q < p[i].nadj; q++) + cout << format(" %d") % p[i].adj[q] << endl; +// throw FileError(); + } + else + if (p[j].ncnt == p[j].nadj) + { + cout << format("OVERFLOW for particle %d (pending %d). List of accepted:") % j % i << endl; + for (int q = 0; q < p[j].nadj; q++) + cout << format(" %d\n") % p[j].adj[q] << endl; +// throw FileError(); + } + else{ + p[i].adj[p[i].ncnt] = j; + p[j].adj[p[j].ncnt] = i; + p[i].ncnt++; + p[j].ncnt++; } + } + else + { + cout << format("%d: adj = %d") % i % j << endl; + } + } + } + } +} + +void readVolumeFile(const std::string& volfile, PARTICLE *p, pid_t np, + pid_t mockIndex) + throw(FileError) +{ + ifstream vol(volfile.c_str()); + pid_t np2; + + if (!vol) + { + cout << "Unable to open volume file " << volfile << endl; + throw FileError(); + } + vol.read((char *)&np2, sizeof(pid_t)); + if (np != np2) { + cout << format("Number of particles doesn't match! %d != %d") % np %np2 << endl; + throw FileError(); + } + cout << format("%d particles") % np << endl; + FF; + + for (pid_t i = 0; i < np; i++) { + vol.read((char*)&p[i].dens, sizeof(float)); + if (((p[i].dens < 1e-30) || (p[i].dens > 1e30)) && (i < mockIndex)) { + cout << format("Whacked-out volume found, of particle %d: %f") % i % p[i].dens << endl; + p[i].dens = 1.; + } + p[i].dens = 1./p[i].dens; /* Get density from volume */ + } +} + +void writeZoneFile(const std::string& zonfile, PARTICLE* p, pid_t np, + ZONE *z, int numZones, int* zonenum, int *jumped) +{ + pid_t **m = new pid_t *[numZones]; + pid_t *nm = new pid_t[numZones]; + + /* Not in the zone struct since it'll be freed up (contiguously, we hope) + soon */ + for (int h=0; h < numZones; h++) + { + m[h] = new pid_t[z[h].np]; + nm[h] = 0; + } + + for (pid_t i = 0; i < np; i++) + { + int h = zonenum[jumped[i]]; + if (z[h].core == i) + { + m[h][nm[h]] = m[h][0]; + m[h][0] = i; /* Put the core particle at the top of the list */ + } + else + { + m[h][nm[h]] = i; + } + nm[h] ++; + } + delete[] nm; + + ofstream zon(zonfile.c_str()); + if (!zon) + { + cout << format("Problem opening zonefile %s.") % zonfile << endl; + throw FileError(); + } + + zon.write((char *)&np, sizeof(pid_t)); + zon.write((char *)&numZones, sizeof(int)); + for (int h = 0; h < numZones; h++) + { + zon.write((char *)&(z[h].np), sizeof(pid_t)); + zon.write((char *)m[h],z[h].np * sizeof(pid_t)); + delete[] m[h]; + } + delete[] m; +} + +void writeVoidFile(const string& zonfile2, ZONE *z, int numZones) +{ + ofstream zon2(zonfile2.c_str()); + if (!zon2) + { + cout << format("Problem opening zonefile %s.)") % zonfile2 << endl; + throw FileError(); + } + zon2.write((char *)&numZones,sizeof(int)); + cout << "Writing void/zone relations..." << endl; + for (int h = 0; h < numZones; h++) + { + zon2.write((char *)&z[h].nhl, sizeof(int)); + zon2.write((char *)z[h].zonelist, z[h].nhl*sizeof(int)); + } +} diff --git a/c_tools/jozov2/jozov2_watershed.cpp b/c_tools/jozov2/jozov2_watershed.cpp new file mode 100644 index 0000000..7382c93 --- /dev/null +++ b/c_tools/jozov2/jozov2_watershed.cpp @@ -0,0 +1,294 @@ +/*+ + VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_watershed.cpp + Copyright (C) 2010-2014 Guilhem Lavaux + Copyright (C) 2011-2014 P. M. Sutter + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; version 2 of the License. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. ++*/ +#ifdef OPENMP +#include +#endif + +#include +#include +#include +#include +#include +#include +#include "jozov2.hpp" +#include "zobov.hpp" + +using namespace std; +using boost::format; + +struct ZoneDensityPair +{ + int h; + double density; + double core; + + bool operator<(const ZoneDensityPair& p2) const + { + return (density > p2.density) || (density==p2.density && core > p2.core); + } +}; + +typedef priority_queue ZoneQueue; + +static void build_process_queue(ZoneQueue& q, ZONE *z, PARTICLE *p, char *inyet, int h) +{ + ZoneDensityPair zdp; + ZONE& z_h = z[h]; + bool interior = true; + + assert(inyet[h] == 1); + for (int za = 0; za < z_h.nadj; za++) + { + zdp.h = z_h.adj[za]; + zdp.density = z_h.slv[za]; + zdp.core = p[z[zdp.h].core].dens; + if (inyet[zdp.h] == 0) + { + q.push(zdp); + interior = false; + } + } + if (interior) + inyet[h] = 2; +} + +void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol) +{ + /* Text output file */ + +#pragma omp parallel + { + char *inyet, *inyet2; + int *zonelist, *zonelist2; + int nhl; + int prev_ii = -1; + + inyet = new char[numZones]; + inyet2 = new char[numZones]; + zonelist = new int[numZones]; + zonelist2 = new int[numZones]; + + fill(inyet, inyet + numZones, 0); + fill(inyet2, inyet2 + numZones, 0); + + nhl = 0; +#pragma omp for schedule(dynamic,1) + for (int h = 0; h < numZones; h++) + { + int nhlcount = 0; + float previous_lowvol = BIGFLT, lowvol, z_cur_core_dens; + bool beaten; + priority_queue to_process; + int link0; + int save_nhl; + pid_t save_npjoin; + + for (int hl = 0; hl < nhl; hl++) + inyet[zonelist[hl]] = 0; + + zonelist[0] = h; + inyet[h] = 1; + save_nhl = nhl = 1; + save_npjoin = z[h].npjoin = z[h].np; + z_cur_core_dens = p[z[h].core].dens; + + build_process_queue(to_process, z, p, inyet, h); + + do { + /* Find the lowest-volume (highest-density) adjacency */ + beaten = false; + + if (to_process.empty()) + { + beaten = true; + z[h].leak = maxvol; + save_npjoin = z[h].npjoin; + save_nhl = nhl; + continue; + } + + do + { + lowvol = to_process.top().density; + link0 = to_process.top().h; + to_process.pop(); + } + while ((inyet[link0] != 0) && (!to_process.empty())); + + if (to_process.empty()) + { + save_npjoin = z[h].npjoin; + save_nhl = nhl; + beaten = true; + z[h].leak = maxvol; + continue; + } + + /* See if there's a beater */ + if (previous_lowvol != lowvol) + { + save_npjoin = z[h].npjoin; + save_nhl = nhl; + previous_lowvol = lowvol; + } + + if (lowvol > voltol) + { + beaten = true; + z[h].leak = lowvol; + continue; + } + + if (p[z[link0].core].dens < z_cur_core_dens) + { + beaten = true; + z[h].leak = lowvol; + continue; + } + + /* Add everything linked to the link(s) */ + int nhl2 = 0; + zonelist2[0] = link0; + inyet2[link0] = 1; + nhl2=1; + + bool added = true; + while (added && !beaten) + { + added = false; + for (int hl = 0; (hl < nhl2) && (!beaten); hl++) + { + int h2 = zonelist2[hl]; + + if (inyet2[h2] == 1) { + bool interior = true; /* Guilty until proven innocent */ + + for (int za = 0; za < z[h2].nadj; za ++) { + int link2 = z[h2].adj[za]; + + if ((inyet[link2]+inyet2[link2]) == 0) { + interior = false; + if (z[h2].slv[za] <= lowvol) { + if (p[z[link2].core].dens < z_cur_core_dens) { + beaten = true; + break; + } + zonelist2[nhl2] = link2; + inyet2[link2] = 1; + nhl2++; + added = true; + } + } + } + + if (interior) + inyet2[h2] = 2; + } + } + } + + /* See if there's a beater */ + if (beaten) { + z[h].leak = lowvol; + } else { + + for (int h2 = 0; h2 < nhl2; h2++) { + int new_h = zonelist2[h2]; + + zonelist[nhl] = new_h; + assert(inyet[new_h] == 0); + z[h].npjoin += z[new_h].np; + inyet[new_h] = 1; + if (inyet2[new_h] != 2) + build_process_queue(to_process, z, p, inyet, new_h); + nhl++; + } + } + + for (int hl = 0; hl < nhl2; hl++) + inyet2[zonelist2[hl]] = 0; + if (nhl/10000 > nhlcount) { + if (nhlcount == 0) + (cout << format("Zone %d: %d") % h % nhl).flush(); + else + (cout << format(" %d [%d]") % nhl % to_process.size()).flush(); + nhlcount = nhl/10000; + } + } + while((lowvol < BIGFLT) && (!beaten)); + + if (!beaten) + { + save_npjoin = z[h].npjoin; + save_nhl = nhl; + } + + z[h].denscontrast = z[h].leak/p[z[h].core].dens; + if (z[h].denscontrast < 1.) + z[h].denscontrast = 1.; + + + /* Don't sort; want the core zone to be first */ + + if (nhlcount > 0) { /* Outputs the number of zones in large voids */ + printf(" h%d:%d\n",h,nhl); + FF; + } + /* Calculate volume */ + z[h].npjoin = save_npjoin; + z[h].voljoin = 0.; + z[h].zonelist = new int[save_nhl]; + z[h].numzones = save_nhl; + for (int q = 0; q < save_nhl; q++) { + z[h].voljoin += z[zonelist[q]].vol; + z[h].zonelist[q] = zonelist[q]; + } + + z[h].nhl = save_nhl; + } + delete[] zonelist; + delete[] zonelist2; + delete[] inyet; + delete[] inyet2; + + } + + double maxdenscontrast = 0; +#pragma omp parallel shared(maxdenscontrast) + { + double maxdenscontrast_local = 0; + +#pragma omp for schedule(static) + for (int h = 0; h < numZones; h++) + { + /* find biggest denscontrast */ + if (z[h].denscontrast > maxdenscontrast_local) { + maxdenscontrast_local = (double)z[h].denscontrast; + } + } + +#pragma omp critical + { + if (maxdenscontrast_local > maxdenscontrast) + maxdenscontrast = maxdenscontrast_local; + } + } + + cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl; +} diff --git a/c_tools/jozov2/jozov2_zones.cpp b/c_tools/jozov2/jozov2_zones.cpp new file mode 100644 index 0000000..2e2f4fe --- /dev/null +++ b/c_tools/jozov2/jozov2_zones.cpp @@ -0,0 +1,243 @@ +/*+ + VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_zones.cpp + Copyright (C) 2010-2014 Guilhem Lavaux + Copyright (C) 2011-2014 P. M. Sutter + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; version 2 of the License. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. ++*/ +#include +#include +#include "jozov2.hpp" +#include "zobov.hpp" +#include + +using namespace std; +using boost::format; + +void buildInitialZones(PARTICLE *p, pid_t np, pid_t* jumped, + pid_t *numinh, pid_t& numZones) +{ + pid_t *jumper = new pid_t[np]; + float minvol; + + /* find jumper */ + for (pid_t i = 0; i < np; i++) { + PARTICLE& cur_p = p[i]; + + minvol = cur_p.dens; + jumper[i] = -1; + for (int j = 0; j < cur_p.nadj; j++) { + if (p[cur_p.adj[j]].dens < minvol) { + jumper[i] = cur_p.adj[j]; + minvol = p[jumper[i]].dens; + } + } + numinh[i] = 0; + } + + (cout << "About to jump ... " << endl).flush(); + + /* Jump */ + for (pid_t i = 0; i < np; i++) { + jumped[i] = i; + while (jumper[jumped[i]] > -1) + jumped[i] = jumper[jumped[i]]; + numinh[jumped[i]]++; + } + (cout << "Post-jump ..." << endl).flush(); + + pid_t loc_NumZones = 0; +#pragma omp parallel for schedule(static) reduction(+:loc_NumZones) + for (pid_t i = 0; i < np; i++) + if (numinh[i] > 0) + loc_NumZones++; + + numZones = loc_NumZones; + cout << format("%d initial zones found") % numZones << endl; + + delete[] jumper; +} + +void buildZoneAdjacencies(PARTICLE *p, pid_t np, + ZONE *z, ZONET *zt, + int numZones, + pid_t *jumped, + int *zonenum, + int *numinh) +{ + /* Count border particles */ + for (pid_t i = 0; i < np; i++) + for (int j = 0; j < p[i].nadj; j++) { + pid_t testpart = p[i].adj[j]; + if (jumped[i] != jumped[testpart]) + zt[zonenum[jumped[i]]].nadj++; + } + + size_t nadjAlloced = 0; + cout << "Total numZones = " << numZones << endl; + try + { + for (int h = 0; h < numZones; h++) { + if (zt[h].nadj > 0) { + zt[h].adj = new pid_t[zt[h].nadj]; + zt[h].slv = new float[zt[h].nadj]; + nadjAlloced += zt[h].nadj; + zt[h].nadj = 0; + } else { + zt[h].adj = 0; + zt[h].slv = 0; + } + } + } + catch(const std::bad_alloc& a) + { + cout << "Could not allocate memory for zone adjacencies (nadj so far: " << nadjAlloced << ", memory needed: " << (nadjAlloced*(sizeof(pid_t)+sizeof(float))) << ")" << endl; + throw a; + } + + /* Find "weakest links" */ + for (pid_t i = 0; i < np; i++) + { + int h = zonenum[jumped[i]]; + ZONET& zt_h = zt[h]; + + for (int j = 0; j < p[i].nadj; j++) + { + pid_t testpart = p[i].adj[j]; + bool already; + + if (h != zonenum[jumped[testpart]]) { + if (p[testpart].dens > p[i].dens) { + /* there could be a weakest link through testpart */ + already = false; + for (int za = 0; za < zt_h.nadj; za++) + if (zt_h.adj[za] == zonenum[jumped[testpart]]) { + already = true; + if (p[testpart].dens < zt_h.slv[za]) { + zt_h.slv[za] = p[testpart].dens; + } + } + if (!already) { + zt_h.adj[zt_h.nadj] = zonenum[jumped[testpart]]; + zt_h.slv[zt_h.nadj] = p[testpart].dens; + zt_h.nadj++; + } + } else { /* There could be a weakest link through i */ + already = false; + for (int za = 0; za < zt_h.nadj; za++) + if (zt_h.adj[za] == zonenum[jumped[testpart]]) { + already = true; + if (p[i].dens < zt_h.slv[za]) { + zt_h.slv[za] = p[i].dens; + } + } + if (!already) { + zt_h.adj[zt_h.nadj] = zonenum[jumped[testpart]]; + zt_h.slv[zt_h.nadj] = p[i].dens; + zt_h.nadj++; + } + } + } + } + } + (cout <<"Found zone adjacencies" << endl).flush(); + + /* Free particle adjacencies */ + for (pid_t i = 0; i < np; i++) + { + if (p[i].adj != 0) + delete[] p[i].adj; + } + + /* Use z instead of zt */ + for (int h = 0; h < numZones; h++) { + z[h].nadj = zt[h].nadj; + if (zt[h].nadj > 0) { + z[h].adj = new pid_t[zt[h].nadj]; + z[h].slv = new float[zt[h].nadj]; + for (int za = 0; za < zt[h].nadj; za++) { + z[h].adj[za] = zt[h].adj[za]; + z[h].slv[za] = zt[h].slv[za]; + } + delete[] zt[h].adj; + delete[] zt[h].slv; + } else { + z[h].adj = 0; + z[h].slv = 0; + } + z[h].np = numinh[z[h].core]; + } +} + +void buildZones(PARTICLE *p, pid_t np, pid_t *&jumped, + ZONE*& z, int& nzones, + int*& zonenum) +{ + pid_t *numinh; + ZONET *zt; + + jumped = new pid_t[np]; + numinh = new pid_t[np]; + + buildInitialZones(p, np, jumped, numinh, nzones); + + try + { + z = new ZONE[nzones]; + zt = new ZONET[nzones]; + zonenum = new int[np]; + } + catch (const std::bad_alloc& e) + { + cout << "Unable to do zone allocations" << endl; + throw e; + } + + int h = 0; + for (pid_t i = 0; i < np; i++) + { + if (numinh[i] > 0) { + z[h].core = i; + z[h].vol = 0; + z[h].np = z[h].npjoin = z[h].nadj = z[h].nhl = 0; + z[h].leak = 0; + z[h].adj = 0; + z[h].slv = 0; + z[h].denscontrast = 0; + z[h].vol = z[h].voljoin = 0; + zonenum[i] = h; + h++; + } else { + zonenum[i] = -1; + } + } + for (size_t z = 0; z < nzones; z++) { + zt[z].nadj = 0; + zt[z].adj = 0; + zt[z].slv = 0; + } + + buildZoneAdjacencies(p, np, z, zt, + h, jumped, zonenum, numinh); + + delete[] zt; + delete[] numinh; + + for (pid_t i=0; i < np; i++) { + int h = zonenum[jumped[i]]; + z[h].vol += 1.0/(double)p[i].dens; + z[h].numzones = 0; + z[h].zonelist = 0; + } +} diff --git a/c_tools/jozov2/zobov.hpp b/c_tools/jozov2/zobov.hpp new file mode 100644 index 0000000..e083ffa --- /dev/null +++ b/c_tools/jozov2/zobov.hpp @@ -0,0 +1,52 @@ +/*+ + VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/zobov.hpp + Copyright (C) 2010-2014 Guilhem Lavaux + Copyright (C) 2011-2014 P. M. Sutter + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; version 2 of the License. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. ++*/ +#ifndef __ZOBOV_HPP +#define __ZOBOV_HPP + +typedef struct Particle { + float dens; + int nadj; + int ncnt; + int *adj; +} PARTICLE; + +typedef struct Zone { + int core; /* Identity of peak particle */ + int np; /* Number of particles in zone */ + int npjoin; /* Number of particles in the joined void */ + int nadj; /* Number of adjacent zones */ + int nhl; /* Number of zones in final joined void */ + float leak; /* Volume of last leak zone*/ + int *adj; /* Each adjacent zone, with ... */ + float *slv; /* Smallest Linking Volume */ + float denscontrast; /* density contrast */ + double vol; /* Total volume of all particles in the zone */ + double voljoin; /* Total volume of all particles in the joined void */ + + int *zonelist; /* Zones bound to the void. */ + int numzones; /* Number of zones bound. */ +} ZONE; + +typedef struct ZoneT { + int nadj; /* Number of zones on border */ + int *adj; /* Each adjacent zone, with ... */ + float *slv; /* Smallest Linking Volume */ +} ZONET; + +#endif diff --git a/c_tools/pruning/CMakeLists.txt b/c_tools/pruning/CMakeLists.txt new file mode 100644 index 0000000..fea19e3 --- /dev/null +++ b/c_tools/pruning/CMakeLists.txt @@ -0,0 +1,8 @@ +include_directories(${CMAKE_CURRENT_BINARY_DIR}) + +SET(pruneVoids_SRCS pruneVoids.cpp) +add_genopt(pruneVoids_SRCS pruneVoids.ggo pruneVoids_conf STRUCTNAME pruneVoids_info) +add_executable(pruneVoids ${pruneVoids_SRCS}) +target_link_libraries(pruneVoids ${ZOB_LIBS}) + +install(TARGETS pruneVoids DESTINATION ${VIDE_BIN}) diff --git a/c_tools/pruning/pruneVoids.cpp b/c_tools/pruning/pruneVoids.cpp new file mode 100644 index 0000000..a82ea72 --- /dev/null +++ b/c_tools/pruning/pruneVoids.cpp @@ -0,0 +1,1213 @@ +/*+ + VIDE -- Void IDentification and Examination -- ./c_tools/stacking/pruneVoids.cpp + Copyright (C) 2010-2014 Guilhem Lavaux + Copyright (C) 2011-2014 P. M. Sutter + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; version 2 of the License. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. ++*/ + + + +// Reads in the void catalog and removes any void that could potentially +// be affected by a mock particle. It does this by computing the longest +// particle distance within each void and comparing it to the distance +// of the nearest mock particle. If the void could potentially by rotated +// to include this particle, we throw out the void. + +// This is placed here instead of using the edgeAvoidance option in +// stackVoidsZero so that we can optionally filter the entire +// catalog at once before the stacking phase. This is useful +// for producing a "clean" void catalog for other purposes. + +#include "gsl/gsl_math.h" +#include "gsl/gsl_eigen.h" +#include "string.h" +#include "ctype.h" +#include "stdlib.h" +#include +#include +#include +#include "pruneVoids_conf.h" +#include +#include "assert.h" +#include "voidTree.hpp" +#include "loadZobov.hpp" +#include +#include + +#define LIGHT_SPEED 299792.458 +#define MPC2Z 100./LIGHT_SPEED +#define Z2MPC LIGHT_SPEED/100. + +#define CENTRAL_VOID 1 +#define EDGE_VOID 2 + +using namespace std; + +typedef struct partStruct { + float x, y, z, vol; + int nadj, ncnt; + int *adj; +} PART; + +typedef struct zoneStruct { + int numPart; + int *partIDs; +} ZONE2PART; + +typedef struct voidZoneStruct { + int numZones; + int *zoneIDs; +} VOID2ZONE; + +typedef struct voidStruct { + float vol, coreDens, zoneVol, densCon, voidProb, radius; + float rescaledCoreDens; + int voidID, numPart, numZones, coreParticle, zoneNumPart; + float maxRadius, nearestMock, centralDen, redshift, redshiftInMpc; + float nearestMockFromCore, nearestGalFromCore; + float nearestEdge; + float center[3], macrocenter[3]; + int accepted; + int voidType; + int parentID, numChildren, level; + bool isLeaf, hasHighCentralDen; + gsl_vector *eval; + gsl_matrix *evec; + float ellip; +} VOID; + +// this defines the expansion function that we will integrate +// Laveaux & Wandelt (2012) Eq. 24 +struct my_expan_params { double Om; double w0; double wa; }; +double expanFun (double z, void * p) { + struct my_expan_params * params = (struct my_expan_params *)p; + double Om = (params->Om); + double w0 = (params->w0); + double wa = (params->wa); + + //const double h0 = 1.0; + const double h0 = 0.71; + double ez; + + double wz = w0 + wa*z/(1+z); + + ez = Om*pow(1+z,3) + (1.-Om); + //ez = Om*pow(1+z,3) + pow(h0,2) * (1.-Om)*pow(1+z,3+3*wz); + + ez = sqrt(ez); + //ez = sqrt(ez)/h0; + + ez = 1./ez; + + return ez; +} + +void openFiles(string outputDir, string sampleName, + string prefix, string dataPortion, + int mockIndex, int numKept, + FILE** fpZobov, FILE** fpCenters, + FILE** fpCentersNoCut, + FILE** fpBarycenter, FILE** fpDistances, FILE** fpShapes, + FILE** fpSkyPositions); + +void closeFiles(FILE* fpZobov, FILE* fpCenters, + FILE* fpCentersNoCut, + FILE* fpBarycenter, FILE* fpDistances, FILE* fpShapes, + FILE* fpSkyPositions); + +void outputVoids(string outputDir, string sampleName, string prefix, + string dataPortion, int mockIndex, + vector voids, + bool isObservation, double *boxLen, + bool doTrim, bool doCentralDenCut); + +int main(int argc, char **argv) { + + // initialize arguments + pruneVoids_info args; + pruneVoids_conf_params args_params; + + pruneVoids_conf_init(&args); + pruneVoids_conf_params_init(&args_params); + + args_params.check_required = 0; + if (pruneVoids_conf_ext (argc, argv, &args, &args_params)) + return 1; + + if (!args.configFile_given) { + if (pruneVoids_conf_required (&args, + PRUNEVOIDS_CONF_PACKAGE)) + return 1; + } else { + args_params.check_required = 1; + args_params.initialize = 0; + if (pruneVoids_conf_config_file (args.configFile_arg, + &args, + &args_params)) + return 1; + } + + // initialize cosmology integrator and interpolator + gsl_function expanF; + expanF.function = &expanFun; + struct my_expan_params expanParams; + expanParams.Om = args.omegaM_arg; + expanParams.w0 = -1.0; + expanParams.wa = 0.0; + expanF.params = &expanParams; + double result, error; + size_t nEval; + + int iZ, numZ = 8000; + double maxZ = 10.0, z, *dL, *redshifts; + 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*LIGHT_SPEED/100.; + //printf("HERE %e %e\n", z, dL[iZ]); + redshifts[iZ] = z; + } + gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, numZ); + gsl_interp_init(interp, dL, redshifts, numZ); + gsl_interp_accel *acc = gsl_interp_accel_alloc(); + + + int i, p, p2, numPartTot, numZonesTot, dummy, iVoid; + int numVoids, mockIndex; + double tolerance; + FILE *fp; + PART *part, *voidPart; + ZONE2PART *zones2Parts; + VOID2ZONE *void2Zones; + vector voids; + float *temp, junk, voidVol; + int junkInt, voidID, numPart, numZones, zoneID, partID, maxNumPart; + int coreParticle, zoneNumPart; + float coreDens, zoneVol, densCon, voidProb, dist[3], dist2, minDist, maxDist; + float centralRad, centralDen; + double nearestEdge, redshift; + char line[500], junkStr[10]; + string outputDir, sampleName, dataPortion, prefix; + int mask_index; + double ranges[3][2], boxLen[3], mul; + double volNorm, radius; + int clock1, clock2, clock3, clock4; + double interval; + int periodicX=0, periodicY=0, periodicZ=0; + string dataPortions[2]; + + gsl_eigen_symmv_workspace *eigw = gsl_eigen_symmv_alloc(3); + + numVoids = args.numVoids_arg; + mockIndex = args.mockIndex_arg; + tolerance = args.tolerance_arg; + + clock1 = clock(); + printf("Pruning parameters: %f %f %f %s\n", args.zMin_arg, + args.zMax_arg, + args.rMin_arg, + args.periodic_arg); + + // check for periodic box + periodicX = 0; + periodicY = 0; + periodicZ = 0; + if (!args.isObservation_flag) { + if ( strchr(args.periodic_arg, 'x') != NULL) { + periodicX = 1; + printf("Will assume x-direction is periodic.\n"); + } + if ( strchr(args.periodic_arg, 'y') != NULL) { + periodicY = 1; + printf("Will assume y-direction is periodic.\n"); + } + if ( strchr(args.periodic_arg, 'z') != NULL) { + periodicZ = 1; + printf("Will assume z-direction is periodic.\n"); + } + } + + // load box size + printf("\n Getting info...\n"); + netCDF::NcFile f_info(args.extraInfo_arg, netCDF::NcFile::read); + f_info.getAtt("range_x_min").getValues(&ranges[0][0]); + f_info.getAtt("range_x_max").getValues(&ranges[0][1]); + f_info.getAtt("range_y_min").getValues(&ranges[1][0]); + f_info.getAtt("range_y_max").getValues(&ranges[1][1]); + f_info.getAtt("range_z_min").getValues(&ranges[2][0]); + f_info.getAtt("range_z_max").getValues(&ranges[2][1]); + + printf(" Range xmin %e\n", ranges[0][0]); + printf(" Range xmax %e\n", ranges[0][1]); + printf(" Range ymin %e\n", ranges[1][0]); + printf(" Range ymax %e\n", ranges[1][1]); + printf(" Range zmin %e\n", ranges[2][0]); + printf(" Range zmax %e\n", ranges[2][1]); + + boxLen[0] = ranges[0][1] - ranges[0][0]; + boxLen[1] = ranges[1][1] - ranges[1][0]; + boxLen[2] = ranges[2][1] - ranges[2][0]; + + // read in all particle positions + clock3 = clock(); + printf("\n Loading particles...\n"); + fp = fopen(args.partFile_arg, "r"); + fread(&dummy, 1, 4, fp); + fread(&numPartTot, 1, 4, fp); + fread(&dummy, 1, 4, fp); + + part = (PART *) malloc(numPartTot * sizeof(PART)); + temp = (float *) malloc(numPartTot * sizeof(float)); + + volNorm = numPartTot/(boxLen[0]*boxLen[1]*boxLen[2]); + printf(" VOL NORM = %f\n", volNorm); + + printf(" CENTRAL DEN = %f\n", args.maxCentralDen_arg); + + fread(&dummy, 1, 4, fp); + fread(temp, numPartTot, 4, fp); + mul = ranges[0][1] - ranges[0][0]; + for (p = 0; p < numPartTot; p++) + part[p].x = mul*temp[p]; + fread(&dummy, 1, 4, fp); + fread(&dummy, 1, 4, fp); + fread(temp, numPartTot, 4, fp); + mul = ranges[1][1] - ranges[1][0]; + for (p = 0; p < numPartTot; p++) + part[p].y = mul*temp[p]; + fread(&dummy, 1, 4, fp); + fread(&dummy, 1, 4, fp); + fread(temp, numPartTot, 4, fp); + mul = ranges[2][1] - ranges[2][0]; + for (p = 0; p < numPartTot; p++) + part[p].z = mul*temp[p]; + + if (!args.isObservation_flag) { + for (p = 0; p < numPartTot; p++) { + part[p].x += ranges[0][0]; + part[p].y += ranges[1][0]; + part[p].z += ranges[2][0]; + } + } + fclose(fp); + + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + printf(" Read %d particles (%.2f sec)...\n", numPartTot, interval); + + + if (mockIndex == -1) mockIndex = numPartTot; + + // read in desired voids + clock3 = clock(); + printf(" Loading voids...\n"); + fp = fopen(args.voidDesc_arg ,"r"); + fgets(line, sizeof(line), fp); + sscanf(line, "%d %s %d %s", &junkInt, junkStr, &junkInt, junkStr); + fgets(line, sizeof(line), fp); + + voids.resize(numVoids); + i = 0; + while (fgets(line, sizeof(line), fp) != NULL) { + sscanf(line, "%d %d %d %f %f %d %d %f %d %f %f\n", &iVoid, &voidID, + &coreParticle, &coreDens, &zoneVol, &zoneNumPart, &numZones, + &voidVol, &numPart, &densCon, &voidProb); + + i++; + voids[i-1].coreParticle = coreParticle; + voids[i-1].zoneNumPart = zoneNumPart; + voids[i-1].coreDens = coreDens; + voids[i-1].zoneVol = zoneVol; + voids[i-1].voidID = voidID; + voids[i-1].vol = voidVol; + voids[i-1].numPart = numPart; + voids[i-1].numZones = numZones; + voids[i-1].densCon = densCon; + voids[i-1].voidProb = voidProb; + + voids[i-1].radius = pow(voidVol/volNorm*3./4./M_PI, 1./3.); + voids[i-1].accepted = 1; + + voids[i-1].isLeaf = true; + voids[i-1].hasHighCentralDen = false; + voids[i-1].numChildren = 0; + voids[i-1].parentID = -1; + voids[i-1].level = 0; + + voids[i-1].eval = gsl_vector_alloc(3); + voids[i-1].evec = gsl_matrix_alloc(3,3); + voids[i-1].ellip = 0; + } + fclose(fp); + + // load up the zone membership for each void + printf(" Loading void-zone membership info...\n"); + fp = fopen(args.void2Zone_arg, "r"); + fread(&numZonesTot, 1, 4, fp); + void2Zones = (VOID2ZONE *) malloc(numZonesTot * sizeof(VOID2ZONE)); + + for (iZ = 0; iZ < numZonesTot; iZ++) { + fread(&numZones, 1, 4, fp); + + void2Zones[iZ].numZones = numZones; + void2Zones[iZ].zoneIDs = (int *) malloc(numZones * sizeof(int)); + + for (p = 0; p < numZones; p++) { + fread(&void2Zones[iZ].zoneIDs[p], 1, 4, fp); + } + } + fclose(fp); + + // now the particles-zone + printf(" Loading particle-zone membership info...\n"); + fp = fopen(args.zone2Part_arg, "r"); + fread(&dummy, 1, 4, fp); + fread(&numZonesTot, 1, 4, fp); + + zones2Parts = (ZONE2PART *) malloc(numZonesTot * sizeof(ZONE2PART)); + for (iZ = 0; iZ < numZonesTot; iZ++) { + fread(&numPart, 1, 4, fp); + + zones2Parts[iZ].numPart = numPart; + zones2Parts[iZ].partIDs = (int *) malloc(numPart * sizeof(int)); + + for (p = 0; p < numPart; p++) { + fread(&zones2Parts[iZ].partIDs[p], 1, 4, fp); + } + } + + // and finally volumes + printf(" Loading particle volumes...\n"); + fp = fopen(args.partVol_arg, "r"); + fread(&mask_index, 1, 4, fp); + if (mask_index != mockIndex) { + printf("NON-MATCHING MOCK INDICES!? %d %d\n", mask_index, mockIndex); + exit(-1); + } + for (p = 0; p < mask_index; p++) { + fread(&temp[0], 1, 4, fp); + part[p].vol = temp[0]; + } + fclose(fp); + +/* + // and finally finally adjacencies + printf(" Loading particle adjacencies...\n"); + fp = fopen(args.partAdj_arg, "r"); + fread(&mask_index, 1, 4, fp); + if (mask_index != mockIndex) { + printf("NON-MATCHING MOCK INDICES!? %d %d\n", mask_index, mockIndex); + exit(-1); + } + int tempInt; + for (p = 0; p < mask_index; p++) { + fread(&tempInt, 1, 4, fp); + part[p].nadj = tempInt; + part[p].ncnt = 0; + if (part[p].nadj > 0) + part[p].adj = (int *) malloc(part[p].nadj * sizeof(int)); + } + for (p = 0; p < mask_index; p++) { + fread(&tempInt, 1, 4, fp); + int nin = tempInt; + if (nin > 0) { + for (int nAdj = 0; nAdj < nin; nAdj++) { + int tempAdj; + fread(&tempAdj, 1, 4, fp); + + // this bit has been readjusted just in case we are + // accidentally still linking to mock particles + + //if (tempAdj < mask_index) { + assert(p < tempAdj); + //if (part[p].ncnt == part[p].nadj) { + // printf("OVERFLOW %d\n", p); + //} else if (part[tempAdj].ncnt == part[tempAdj].nadj) { + // printf("OVERFLOW %d\n", tempAdj); + //} else { + part[p].adj[part[p].ncnt] = tempAdj; + part[p].ncnt++; + if (tempAdj < mask_index) { + part[tempAdj].adj[part[tempAdj].ncnt] = p; + part[tempAdj].ncnt++; + } + //} + //} + } +//printf("ADJ %d %d %d %d %d\n", p, nin, part[p].nadj, nAdj, tempInt); + } + } + fclose(fp); +*/ + + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + printf(" Read voids (%.2f sec)...\n", interval); + + // load voids *again* using Guilhem's code so we can get tree + clock3 = clock(); + //if (!args.isObservation_flag) { + printf(" Re-loading voids and building tree..\n"); + ZobovRep zobovCat; + if (!loadZobov(args.voidDesc_arg, args.zone2Part_arg, + args.void2Zone_arg, + 0, zobovCat)) { + printf("Error loading catalog!\n"); + return -1; + } + VoidTree *tree; + tree = new VoidTree(zobovCat); + zobovCat.allZones.erase(zobovCat.allZones.begin(), zobovCat.allZones.end()); + + // copy tree information to our own data structures + for (iVoid = 0; iVoid < numVoids; iVoid++) { + voidID = voids[iVoid].voidID; + voids[iVoid].parentID = tree->getParent(voidID); + voids[iVoid].numChildren = tree->getChildren(voidID).size(); + + // compute level in tree + int level = 0; + int parentID = tree->getParent(voidID); + while (parentID != -1) { + level++; + parentID = tree->getParent(parentID); + } + voids[iVoid].level = level; + } + //} // end re-load + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + printf(" Re-read voids (%.2f sec)...\n", interval); + + // check boundaries + printf(" Computing void properties...\n"); + + maxNumPart = 0; + for (iVoid = 0; iVoid < numVoids; iVoid++) { + if (voids[iVoid].numPart > maxNumPart) maxNumPart = voids[iVoid].numPart; + } + + voidPart = (PART *) malloc(maxNumPart * sizeof(PART)); + + for (iVoid = 0; iVoid < numVoids; iVoid++) { + + voidID = voids[iVoid].voidID; + printf(" DOING %d (of %d) %d %d %f\n", iVoid, numVoids, voidID, + voids[iVoid].numPart, + voids[iVoid].radius); + + voids[iVoid].center[0] = part[voids[iVoid].coreParticle].x; + voids[iVoid].center[1] = part[voids[iVoid].coreParticle].y; + voids[iVoid].center[2] = part[voids[iVoid].coreParticle].z; + + // first load up particles into a buffer + clock3 = clock(); + i = 0; + for (iZ = 0; iZ < void2Zones[voidID].numZones; iZ++) { + zoneID = void2Zones[voidID].zoneIDs[iZ]; + + for (p = 0; p < zones2Parts[zoneID].numPart; p++) { + partID = zones2Parts[zoneID].partIDs[p]; + + if (partID > mask_index || + (part[partID].vol < 1.e-27 && part[partID].vol > 0.)) { + printf("BAD PART!? %d %d %e", partID, mask_index, part[partID].vol); + exit(-1); + } + + voidPart[i].x = part[partID].x; + voidPart[i].y = part[partID].y; + voidPart[i].z = part[partID].z; + voidPart[i].vol = part[partID].vol; + +/* + // testing for edge contamination + if (part[partID].vol < 1.e-27) { + printf("CONTAMINATED!! %d %d\n", iVoid, partID); + } else { + //printf("NORMAL!! %d %d %e\n", iVoid, partID, part[partID].vol); + } + for (int iAdj = 0; iAdj < part[partID].ncnt; iAdj++) { + if (part[partID].adj[iAdj] > mockIndex) { + printf("CONTAMINATED!! %d %d %d\n", iVoid, partID, iAdj); + } + } +*/ + i++; + } + } + + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + //printf(" %.2f for buffer\n", interval); + + // compute macrocenters + clock3 = clock(); + double weight = 0.; + voids[iVoid].macrocenter[0] = 0.; + voids[iVoid].macrocenter[1] = 0.; + voids[iVoid].macrocenter[2] = 0.; + + for (p = 0; p < voids[iVoid].numPart; p++) { + dist[0] = voidPart[p].x - voids[iVoid].center[0]; + dist[1] = voidPart[p].y - voids[iVoid].center[1]; + dist[2] = voidPart[p].z - voids[iVoid].center[2]; + + if (periodicX && fabs(dist[0]) > boxLen[0]/2.) + dist[0] = dist[0] - copysign(boxLen[0], dist[0]); + if (periodicY && fabs(dist[1]) > boxLen[1]/2.) + dist[1] = dist[1] - copysign(boxLen[1], dist[1]); + if (periodicZ && fabs(dist[2]) > boxLen[2]/2.) + dist[2] = dist[2] - copysign(boxLen[2], dist[2]); + + voids[iVoid].macrocenter[0] += voidPart[p].vol*(dist[0]); + voids[iVoid].macrocenter[1] += voidPart[p].vol*(dist[1]); + voids[iVoid].macrocenter[2] += voidPart[p].vol*(dist[2]); + weight += voidPart[p].vol; + } + voids[iVoid].macrocenter[0] /= weight; + voids[iVoid].macrocenter[1] /= weight; + voids[iVoid].macrocenter[2] /= weight; + voids[iVoid].macrocenter[0] += voids[iVoid].center[0]; + voids[iVoid].macrocenter[1] += voids[iVoid].center[1]; + voids[iVoid].macrocenter[2] += voids[iVoid].center[2]; + + if (periodicX) { + if (voids[iVoid].macrocenter[0] > ranges[0][1]) + voids[iVoid].macrocenter[0] = voids[iVoid].macrocenter[0] - boxLen[0]; + if (voids[iVoid].macrocenter[0] < ranges[0][0]) + voids[iVoid].macrocenter[0] = boxLen[0] + voids[iVoid].macrocenter[0]; + } + if (periodicY) { + if (voids[iVoid].macrocenter[1] > ranges[1][1]) + voids[iVoid].macrocenter[1] = voids[iVoid].macrocenter[1] - boxLen[1]; + if (voids[iVoid].macrocenter[1] < ranges[1][0]) + voids[iVoid].macrocenter[1] = boxLen[1] + voids[iVoid].macrocenter[1]; + } + if (periodicZ) { + if (voids[iVoid].macrocenter[2] > ranges[2][1]) + voids[iVoid].macrocenter[2] = voids[iVoid].macrocenter[2] - boxLen[2]; + if (voids[iVoid].macrocenter[2] < ranges[2][0]) + voids[iVoid].macrocenter[2] = boxLen[2] + voids[iVoid].macrocenter[2]; + } + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + //printf(" %.2f for macrocenter\n", interval); + + // compute central density + clock3 = clock(); + centralRad = voids[iVoid].radius/args.centralRadFrac_arg; + centralDen = 0.; + int numCentral = 0; + for (p = 0; p < voids[iVoid].numPart; p++) { + dist[0] = fabs(voidPart[p].x - voids[iVoid].macrocenter[0]); + dist[1] = fabs(voidPart[p].y - voids[iVoid].macrocenter[1]); + dist[2] = fabs(voidPart[p].z - voids[iVoid].macrocenter[2]); + + if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]); + if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]); + if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]); + + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + if (sqrt(dist2) < centralRad) numCentral += 1; + } + voids[iVoid].centralDen = numCentral / (volNorm*4./3. * M_PI * + pow(centralRad, 3.)); + + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + //printf(" %.2f for central density\n", interval); + + //coreParticle = voids[iVoid].coreParticle; + //voids[iVoid].rescaledCoreDens = voids[iVoid].coreDens*(pow(1.*mockIndex/numPartTot,3)); + // // compute distance from core to nearest mock + // minDist = 1.e99; + // for (p = mockIndex; p < numPartTot; p++) { + // dist[0] = part[coreParticle].x - part[p].x; + // dist[1] = part[coreParticle].y - part[p].y; + // dist[2] = part[coreParticle].z - part[p].z; + // + // dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + // if (dist2 < minDist) minDist = dist2; + // } + // voids[iVoid].nearestMockFromCore = sqrt(minDist); + // + // // compute distance from core to nearest mock + // minDist = 1.e99; + // for (p = 0; p < mockIndex; p++) { + // dist[0] = part[coreParticle].x - part[p].x; + // dist[1] = part[coreParticle].y - part[p].y; + // dist[2] = part[coreParticle].z - part[p].z; + // + // dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + // if (dist2 < minDist && dist2 > 1.e-10) minDist = dist2; + // } + // voids[iVoid].nearestGalFromCore = sqrt(minDist); + + // compute maximum extent +/* + if (args.isObservation_flag) { + maxDist = 0.; + for (p = 0; p < voids[iVoid].numPart; p++) { + for (p2 = p; p2 < voids[iVoid].numPart; p2++) { + + dist[0] = voidPart[p].x - voidPart[p2].x; + dist[1] = voidPart[p].y - voidPart[p2].y; + dist[2] = voidPart[p].z - voidPart[p2].z; + + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + if (dist2 > maxDist) maxDist = dist2; + } + } + voids[iVoid].maxRadius = sqrt(maxDist)/2.; + } else { +*/ + + clock3 = clock(); + maxDist = 0.; + for (p = 0; p < voids[iVoid].numPart; p++) { + + dist[0] = fabs(voidPart[p].x - voids[iVoid].macrocenter[0]); + dist[1] = fabs(voidPart[p].y - voids[iVoid].macrocenter[1]); + dist[2] = fabs(voidPart[p].z - voids[iVoid].macrocenter[2]); + + if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]); + if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]); + if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]); + + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + if (dist2 > maxDist) maxDist = dist2; + } + voids[iVoid].maxRadius = sqrt(maxDist); + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + //printf(" %.2f for maximum extent\n", interval); +// } + + clock3 = clock(); + if (args.isObservation_flag) { + // compute distance from center to nearest mock + minDist = 1.e99; + for (p = mockIndex; p < numPartTot; p++) { + dist[0] = voids[iVoid].macrocenter[0] - part[p].x; + dist[1] = voids[iVoid].macrocenter[1] - part[p].y; + dist[2] = voids[iVoid].macrocenter[2] - part[p].z; + + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + if (dist2 < minDist) minDist = dist2; + } + voids[iVoid].nearestMock = sqrt(minDist); + + } else { + voids[iVoid].nearestMock = 1.e99; + } + + if (args.isObservation_flag) { + voids[iVoid].redshiftInMpc = + 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; + + + if (args.useComoving_flag) { + redshift = gsl_interp_eval(interp, dL, redshifts, + voids[iVoid].redshiftInMpc, acc); + //printf("HELLO %e %e\n", redshift, args.zMax_arg); + nearestEdge = fabs((redshift-args.zMax_arg)*LIGHT_SPEED/100.); + voids[iVoid].redshift = redshift; + } else { + redshift = voids[iVoid].redshiftInMpc; + nearestEdge = fabs(redshift-args.zMax_arg*LIGHT_SPEED/100.); + voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.; + } + //nearestEdge = fmin(fabs(redshift-args.zMin_arg*LIGHT_SPEED/100.), + // fabs(redshift-args.zMax_arg*LIGHT_SPEED/100.)); + + } else { + + voids[iVoid].redshiftInMpc = voids[iVoid].macrocenter[2]; + if (args.useComoving_flag) { + voids[iVoid].redshift = gsl_interp_eval(interp, dL, redshifts, + voids[iVoid].redshiftInMpc, acc); + } else { + voids[iVoid].redshift = voids[iVoid].macrocenter[2]/LIGHT_SPEED*100.; + } + + nearestEdge = 1.e99; + + if (!periodicX) { + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].macrocenter[0] - ranges[0][0])); + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].macrocenter[0] - ranges[0][1])); + } + if (!periodicY) { + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].macrocenter[1] - ranges[1][0])); + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].macrocenter[1] - ranges[1][1])); + } + if (!periodicZ) { + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].macrocenter[2] - ranges[2][0])); + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].macrocenter[2] - ranges[2][1])); + } + } + + voids[iVoid].nearestEdge = nearestEdge; + + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + //printf(" %.2f for nearest edge\n", interval); + + // compute eigenvalues and vectors for orientation and shape + clock3 = clock(); + double inertia[9]; + for (int i = 0; i < 9; i++) inertia[i] = 0.; + + for (int p = 0; p < voids[iVoid].numPart; p++) { + dist[0] = voidPart[p].x - voids[iVoid].macrocenter[0]; + dist[1] = voidPart[p].y - voids[iVoid].macrocenter[1]; + dist[2] = voidPart[p].z - voids[iVoid].macrocenter[2]; + + if (periodicX && fabs(dist[0]) > boxLen[0]/2.) + dist[0] = dist[0] - copysign(boxLen[0], dist[0]); + if (periodicY && fabs(dist[1]) > boxLen[1]/2.) + dist[1] = dist[1] - copysign(boxLen[1], dist[1]); + if (periodicZ && fabs(dist[2]) > boxLen[2]/2.) + dist[2] = dist[2] - copysign(boxLen[2], dist[2]); + + for (int i = 0; i < 3; i++) { + for (int j = i; j < 3; j++) { + if (i == j) inertia[i*3+j] += dist[0]*dist[0] + dist[1]*dist[1] + + dist[2]*dist[2]; + inertia[i*3+j] -= dist[i]*dist[j]; + } + } + } + inertia[1*3+0] = inertia[0*3+1]; + inertia[2*3+0] = inertia[0*3+2]; + inertia[2*3+1] = inertia[1*3+2]; + + gsl_matrix_view m = gsl_matrix_view_array(inertia, 3, 3); + gsl_eigen_symmv(&m.matrix, voids[iVoid].eval, voids[iVoid].evec, eigw); + + float a = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,1) + + gsl_vector_get(voids[iVoid].eval,2) - + gsl_vector_get(voids[iVoid].eval,0))); + float b = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,2) + + gsl_vector_get(voids[iVoid].eval,0) - + gsl_vector_get(voids[iVoid].eval,1))); + float c = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,0) + + gsl_vector_get(voids[iVoid].eval,1) - + gsl_vector_get(voids[iVoid].eval,2))); + float ca; + float cb = c/b; + + float smallest = 1.e99; + float largest = 0.; + for (int i = 0; i < 3; i++) { + if (gsl_vector_get(voids[iVoid].eval,i) < smallest) + smallest = gsl_vector_get(voids[iVoid].eval,i); + if (gsl_vector_get(voids[iVoid].eval,i) > largest) + largest = gsl_vector_get(voids[iVoid].eval,i); + } + // TEST + voids[iVoid].ellip = 1.0 - sqrt(sqrt(fabs(smallest/largest))); + + //if (a < c) ca = a/c; + //if (a >= c) ca = c/a; + //voids[iVoid].ellip = fabs(1.0 - ca); + + //if (a < c) ca = a*a/(c*c); + //if (a >= c) ca = (c*c)/(a*a); + //voids[iVoid].ellip = sqrt(fabs(1.0 - ca)); + + clock4 = clock(); + interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC; + //printf(" %.2f for ellipticity\n", interval); + } // iVoid + + gsl_eigen_symmv_free(eigw); + + int numWrong = 0; + int numHighDen = 0; + int numCentral = 0; + int numEdge = 0; + int numNearZ = 0; + int numAreParents = 0; + int numTooSmall = 0; + + printf(" Picking winners and losers...\n"); + printf(" Starting with %d voids\n", (int) voids.size()); + + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + voids[iVoid].accepted = 1; + } + +/* + int j = 0; + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + if (voids[iVoid].densCon < 1.5) { +// voids[iVoid].accepted = -4; + } + } +*/ + + // toss out voids that are obviously wrong + int iGood = 0; + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + if (voids[iVoid].densCon > 1.e4 || isnan(voids[iVoid].vol) || + isinf(voids[iVoid].vol)) { + numWrong++; + } else { + voids[iGood++] = voids[iVoid]; + } + } + voids.resize(iGood); + printf(" 1st filter: rejected %d obviously bad\n", numWrong); + + iGood = 0; + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + if (voids[iVoid].radius < args.rMin_arg) { + numTooSmall++; + } else { + voids[iGood++] = voids[iVoid]; + } + } + voids.resize(iGood); + printf(" 2nd filter: rejected %d too small\n", numTooSmall); + + + iGood = 0; + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + // *always* clean out near edges since there are no mocks there + if (tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestEdge || + tolerance*voids[iVoid].radius > voids[iVoid].nearestEdge) { + numNearZ++; + } else { + voids[iGood++] = voids[iVoid]; + } + } + voids.resize(iGood); + printf(" 3rd filter: rejected %d too close to high redshift boundaries\n", numNearZ); + + numNearZ = 0; + iGood = 0; + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + // assume the lower z-boundary is "soft" in observations + if (args.isObservation_flag && + voids[iVoid].redshift < args.zMin_arg) { + numNearZ++; + } else { + voids[iGood++] = voids[iVoid]; + } + } + voids.resize(iGood); + + //Maubert - Uncommented this part : to be sure that voids do not cross maximum redshift asked for in zrange + iGood = 0; + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + // just in case + if (args.isObservation_flag && + voids[iVoid].redshift > args.zMax_arg) { + numNearZ++; + } else { + voids[iGood++] = voids[iVoid]; + } + } + voids.resize(iGood); + // Maubert - End of Uncommented part + + + printf(" 4th filter: rejected %d outside redshift boundaries\n", numNearZ); + + // take only top-level voids + numAreParents = 0; + iGood = 0; + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + if (voids[iVoid].parentID != -1) { + numAreParents++; + voids[iVoid].isLeaf = true; + } else { + voids[iVoid].isLeaf = false; + } + } + + + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + if (voids[iVoid].centralDen > args.maxCentralDen_arg) { + voids[iVoid].accepted = -1; + voids[iVoid].hasHighCentralDen = true; + numHighDen++; + } else { + voids[iVoid].hasHighCentralDen = false; + } + } + + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + if (tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) { + voids[iVoid].voidType = CENTRAL_VOID; + numCentral++; + } else { + voids[iVoid].voidType = EDGE_VOID; + numEdge++; + } + } + + printf(" Number kept: %d (out of %d)\n", (int) voids.size(), numVoids); + printf(" We have %d edge voids\n", numEdge); + printf(" We have %d central voids\n", numCentral); + printf(" We have %d too high central density\n", numHighDen); + printf(" We have %d that are not leaf nodes\n", numAreParents); + + + outputDir = string(args.outputDir_arg); + sampleName = (string(args.sampleName_arg)+".out"); + + dataPortions[0] = "central"; + dataPortions[1] = "all"; + + printf(" Output fully trimmed catalog...\n"); + prefix = ""; + for (int i = 0; i < 2; i++) { + dataPortion = dataPortions[i]; + + outputVoids(outputDir, sampleName, prefix, dataPortion, + mockIndex, + voids, + args.isObservation_flag, boxLen, true, true); + } + + printf(" Output fully untrimmed catalog...\n"); + prefix = "untrimmed_"; + for (int i = 0; i < 2; i++) { + dataPortion = dataPortions[i]; + + outputVoids(outputDir, sampleName, prefix, dataPortion, + mockIndex, + voids, + args.isObservation_flag, boxLen, false, false); + } + + printf(" Output partially-trimmed catalogs...\n"); + prefix = "untrimmed_dencut_"; + for (int i = 0; i < 2; i++) { + dataPortion = dataPortions[i]; + + outputVoids(outputDir, sampleName, prefix, dataPortion, + mockIndex, + voids, + args.isObservation_flag, boxLen, false, true); + } + + prefix = "trimmed_nodencut_"; + for (int i = 0; i < 2; i++) { + dataPortion = dataPortions[i]; + + outputVoids(outputDir, sampleName, prefix, dataPortion, + mockIndex, + voids, + args.isObservation_flag, boxLen, true, false); + } + + + + + + clock2 = clock(); + printf(" Time: %f sec (for %d voids)\n", + (1.*clock2-clock1)/CLOCKS_PER_SEC, numVoids); + clock1 = clock(); + + + printf("Done!\n"); +} // end main + + +// ---------------------------------------------------------------------------- +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** fpSkyPositions) { + + *fpZobov = fopen((outputDir+"/"+prefix+"voidDesc_"+dataPortion+"_"+sampleName).c_str(), "w"); + fprintf(*fpZobov, "%d particles, %d voids.\n", mockIndex, numKept); + fprintf(*fpZobov, "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n"); + + *fpBarycenter = fopen((outputDir+"/"+prefix+"macrocenters_"+dataPortion+"_"+sampleName).c_str(), "w"); + + *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"); + + *fpShapes = fopen((outputDir+"/"+prefix+"shapes_"+dataPortion+"_"+sampleName).c_str(), "w"); + fprintf(*fpShapes, "# void ID, ellip, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n"); + +} // end openFiles + + +// ---------------------------------------------------------------------------- +void closeFiles(FILE* fpZobov, FILE* fpCenters, + FILE* fpBarycenter, FILE* fpDistances, FILE* fpShapes, + FILE* fpSkyPositions) { + + fclose(fpZobov); + fclose(fpCenters); + fclose(fpBarycenter); + fclose(fpDistances); + fclose(fpShapes); + fclose(fpSkyPositions); + +} // end closeFile + +// ---------------------------------------------------------------------------- +void outputVoids(string outputDir, string sampleName, string prefix, + string dataPortion, int mockIndex, + vector voids, + bool isObservation, double *boxLen, bool doTrim, + bool doCentralDenCut) { + + int iVoid; + VOID outVoid; + FILE *fp, *fpZobov, *fpCenters, *fpCentersNoCut, *fpBarycenter, + *fpDistances, *fpShapes, *fpSkyPositions; + + + openFiles(outputDir, sampleName, prefix, dataPortion, + mockIndex, voids.size(), + &fpZobov, &fpCenters, &fpBarycenter, + &fpDistances, &fpShapes, &fpSkyPositions); + + + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + outVoid = voids[iVoid]; + + + if (dataPortion == "central" && outVoid.voidType == EDGE_VOID) { + continue; + } + + if (doTrim && outVoid.isLeaf) { + continue; + } + + if (doCentralDenCut && outVoid.hasHighCentralDen) { + continue; + } + + double outCenter[3]; + outCenter[0] = outVoid.macrocenter[0]; + outCenter[1] = outVoid.macrocenter[1]; + outCenter[2] = outVoid.macrocenter[2]; + + //if (isObservation) { + // outCenter[0] = (outVoid.macrocenter[0]-boxLen[0]/2.)*100.; + // outCenter[1] = (outVoid.macrocenter[1]-boxLen[1]/2.)*100.; + // outCenter[2] = (outVoid.macrocenter[2]-boxLen[2]/2.)*100.; + //} + + fprintf(fpZobov, "%d %d %d %f %f %d %d %f %d %f %f\n", + iVoid, + outVoid.voidID, + outVoid.coreParticle, + outVoid.coreDens, + outVoid.zoneVol, + outVoid.zoneNumPart, + outVoid.numZones, + outVoid.vol, + outVoid.numPart, + outVoid.densCon, + outVoid.voidProb); + + fprintf(fpBarycenter, "%d %e %e %e\n", + outVoid.voidID, + outVoid.macrocenter[0], + outVoid.macrocenter[1], + outVoid.macrocenter[2]); + + fprintf(fpDistances, "%d %e %e %e %e %e\n", + outVoid.voidID, + outVoid.nearestMock, + outVoid.radius, + 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], + outCenter[1], + outCenter[2], + outVoid.vol, + outVoid.radius, + outVoid.redshift, + 4./3.*M_PI*pow(outVoid.radius, 3), + outVoid.voidID, + outVoid.densCon, + outVoid.numPart, + outVoid.parentID, + outVoid.level, + outVoid.numChildren, + outVoid.centralDen); + + double phi = atan2(outVoid.macrocenter[1]-boxLen[1]/2., + outVoid.macrocenter[0]-boxLen[0]/2.); + if (phi < 0) phi += 2.*M_PI; + double RA = phi * 180./M_PI; + + double theta = acos((outVoid.macrocenter[2]-boxLen[2]/2.) / + outVoid.redshiftInMpc); + double dec = (M_PI/2. - theta) * 180./M_PI; + + fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n", + RA, + dec, + outVoid.redshift, + outVoid.radius, + outVoid.voidID); + + fprintf(fpShapes, "%d %.6f %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e\n", + outVoid.voidID, + outVoid.ellip, + gsl_vector_get(outVoid.eval, 0), + gsl_vector_get(outVoid.eval, 1), + gsl_vector_get(outVoid.eval, 2), + gsl_matrix_get(outVoid.evec, 0 ,0), + gsl_matrix_get(outVoid.evec, 1 ,0), + gsl_matrix_get(outVoid.evec, 2 ,0), + gsl_matrix_get(outVoid.evec, 0 ,1), + gsl_matrix_get(outVoid.evec, 1 ,1), + gsl_matrix_get(outVoid.evec, 2 ,1), + gsl_matrix_get(outVoid.evec, 0 ,2), + gsl_matrix_get(outVoid.evec, 1 ,2), + gsl_matrix_get(outVoid.evec, 2 ,2) + ); + + } // end iVoid + + closeFiles(fpZobov, fpCenters, fpBarycenter, + fpDistances, fpShapes, fpSkyPositions); + +} // end outputVoids diff --git a/c_tools/pruning/pruneVoids.ggo b/c_tools/pruning/pruneVoids.ggo new file mode 100644 index 0000000..77f4521 --- /dev/null +++ b/c_tools/pruning/pruneVoids.ggo @@ -0,0 +1,44 @@ +package "removeEdgeVoids" +version "alpha" + + +option "configFile" - "Configuration filename" string optional + +option "partFile" - "Particle file from generateFromCatalog" string required + +option "extraInfo" - "Extra particle file from generateFromCatalog" string required + +option "voidDesc" - "Void list from ZOBOV" string required + +option "void2Zone" - "Zone file from ZOBOV" string required + +option "partVol" - "Particle volume file from ZOBOV" string required + +option "partAdj" - "Adjacency file from ZOBOV" string required + +option "zone2Part" - "Particle file from ZOBOV" string required +option "mockIndex" - "Beginning index of mock particles" int required + +option "numVoids" - "Number of voids" int required + +option "isObservation" - "We are working with observational data" flag off + +option "useComoving" - "Void positions are in comoving coordinates" flag off + +option "omegaM" - "Omega_M for redshift convertion" double optional default="0.27" + +option "zMin" - "Minimum redshift of sample" double optional default="0.0" +option "zMax" - "Maximum redshift of sample" double optional default="10.0" + +option "rMin" - "Minimum allowable void radius" double optional default="0.0" + +option "outputDir" - "Directory to place outputs" string required +option "sampleName" - "unique string to assign to outputs" string required + +option "periodic" - "Set of edges which are periodic" string optional default="xy" + +option "tolerance" - "Fraction of void width to consider edge" double optional default="1.0" + +option "centralRadFrac" - "Fraction of void radii to consider central region" double optional default="4" + +option "maxCentralDen" - "Maximum central density to accept as a void" double optional default="0.2" diff --git a/python_tools/vide/backend/cosmologyTools.py b/python_tools/vide/backend/cosmologyTools.py new file mode 100644 index 0000000..c6edd24 --- /dev/null +++ b/python_tools/vide/backend/cosmologyTools.py @@ -0,0 +1,112 @@ +#+ +# VIDE -- Void IDentification and Examination -- ./python_tools/vide/apTools/chi2/cosmologyTools.py +# Copyright (C) 2010-2014 Guilhem Lavaux +# Copyright (C) 2011-2014 P. M. Sutter +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; version 2 of the License. +# +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +#+ +# a suite of functions to compute expansion rates, angular diameter +# distances, and expected void stretching + +import numpy as np +import scipy.integrate as integrate +import healpy as healpy +import os +from vide.backend import * + +__all__=['expansion', 'angularDiameter', 'aveExpansion', 'getSurveyProps'] + +# returns 1/E(z) for the given cosmology +def expansion(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): + ez = Om * (1+z)**3 + (Ot-Om)# * (1+z)**(3.+3*wz) + #ez = Om * (1+z)**3 + (Ot-Om)# * integrade.quad(eosDE, 0.0, z, args=(w0,wa))[0] + ez = 1./np.sqrt(ez) + return ez + +# returns DE value at redshift z +def eosDE(z, w0 = -1.0, wa = 0.0): + return w0 + wa*z/(1+z) + +# returns D_A(z) for the given cosmology +def angularDiameter(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): + da = integrate.quad(expansion, 0.0, z, args=(Om, Ot, w0, wa))[0] + return da + +# ----------------------------------------------------------------------------- +# returns average expected expansion for a given redshift range +def aveExpansion(zStart, zEnd, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): + if zStart == 0.0: zStart = 1.e-6 + ave = integrate.quad(expansion, zStart, zEnd, args=(Om, Ot, w0, wa))[0] + ave = (zEnd-zStart)/ave + return ave + +# ----------------------------------------------------------------------------- +# returns the volume and galaxy density for a given redshit slice +def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectionFuncFile=None, useComoving=False): + + LIGHT_SPEED = 299792.458 + + mask = healpy.read_map(maskFile) + area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi + + if useComoving: + zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=0.27) + zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=0.27) + selFunMin = LIGHT_SPEED/100.*angularDiameter(selFunMin, Om=0.27) + selFunMax = LIGHT_SPEED/100.*angularDiameter(selFunMax, Om=0.27) + else: + zmin = zmin * 3000 + zmax = zmax * 3000 + selFunMin *= 3000 + selFunMax *= 3000 + volume = area * (zmax**3 - zmin**3) / 3 + + if selectionFuncFile == None: + nbar = 1.0 + + elif not os.access(selectionFuncFile, os.F_OK): + print(" Warning, selection function file %s not found, using default of uniform selection." % selectionFuncFile) + nbar = 1.0 + + else: + selfunc = np.genfromtxt(selectionFuncFile) + selfunc = np.array(selfunc) + selfunc[:,0] = selfunc[:,0]/100. + selfuncUnity = selfunc + selfuncUnity[:,1] = 1.0 + selfuncMin = selfunc[0,0] + selfuncMax = selfunc[-1,0] + selfuncDx = selfunc[1,0] - selfunc[0,0] + selfuncN = np.size(selfunc[:,0]) + + selFunMin = max(selFunMin, selfuncMin) + selFunMax = min(selFunMax, selfuncMax) + + + def f(z): return selfunc[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 + def fTotal(z): return selfuncUnity[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 + + zrange = np.linspace(selFunMin, selFunMax) + + nbar = scipy.integrate.quad(f, selFunMin, selFunMax) + nbar = nbar[0] + + ntotal = scipy.integrate.quad(fTotal, 0.0, max(selfuncUnity[:,0])) + ntotal = ntotal[0] + + nbar = ntotal / area / nbar + + + return (volume, nbar)