diff --git a/c_tools/mock/loaders/sdf_loader.cpp b/c_tools/mock/loaders/sdf_loader.cpp index cab171d..f2b0c32 100644 --- a/c_tools/mock/loaders/sdf_loader.cpp +++ b/c_tools/mock/loaders/sdf_loader.cpp @@ -216,11 +216,14 @@ SimulationLoader *sdfLoader(const std::string& snapshot, int flags, hdr = new SimuData; SDFgetintOrDefault(sdfp, "version", &fileversion, 1); + double h0; if (fileversion == 1) { SDFgetfloatOrDie(sdfp, "Omega0", &hdr->Omega_M); SDFgetfloatOrDie(sdfp, "Lambda_prime", &hdr->Omega_Lambda); SDFgetfloatOrDie(sdfp, "H0", &hdr->Hubble); + hdr->Hubble *= 1000.*(one_kpc/one_Gyr); + h0 = hdr->Hubble / 100.; } else { @@ -232,12 +235,10 @@ SimulationLoader *sdfLoader(const std::string& snapshot, int flags, hdr->Omega_M += Or; hdr->Omega_Lambda += Of; - SDFgetfloatOrDie(sdfp, "H0", &hdr->Hubble); - + SDFgetfloatOrDie(sdfp, "h_100", &hdr->Hubble); + hdr->Hubble *= 100.; + h0 = hdr->Hubble / 100.; } - double h0; - hdr->Hubble *= 1000.0*(one_kpc/one_Gyr); - h0 = hdr->Hubble/100.; if (SDFhasname("R0", sdfp)) { diff --git a/c_tools/zobov2/findrtop.c b/c_tools/zobov2/findrtop.c new file mode 100644 index 0000000..8dbeeb3 --- /dev/null +++ b/c_tools/zobov2/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/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2.cpp new file mode 100644 index 0000000..718b8e9 --- /dev/null +++ b/c_tools/zobov2/jozov2.cpp @@ -0,0 +1,156 @@ +#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/zobov2/jozov2_io.cpp b/c_tools/zobov2/jozov2_io.cpp new file mode 100644 index 0000000..1fe97f4 --- /dev/null +++ b/c_tools/zobov2/jozov2_io.cpp @@ -0,0 +1,182 @@ +#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(); + } + + 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(); + } + + 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; + z[h].vol = 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/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp new file mode 100644 index 0000000..78f6694 --- /dev/null +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -0,0 +1,252 @@ +#ifdef OPENMP +#include +#endif + +#include +#include +#include +#include +#include "jozov2.hpp" +#include "zobov.hpp" + +using namespace std; +using boost::format; + +void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol) +{ + int *iord; + + double *sorter = new double[numZones+1]; + /* Assign sorter by probability (could use volume instead) */ + for (int h = 0; h < numZones; h++) + sorter[h] = (double)z[h].core; + + /* Text output file */ + + printf("about to sort (pre-watershed sort) ...\n");FF; + + iord = new int[numZones]; + + findrtop(sorter, numZones, iord, numZones); + delete[] sorter; + +#pragma omp parallel + { + char *inyet, *inyet2; + int *zonelist, *zonelist2; + int nhl; + int *links = new int[NLINKS]; + bool *done_zones; + + inyet = new char[numZones]; + inyet2 = new char[numZones]; + zonelist = new int[numZones]; + zonelist2 = new int[numZones]; + done_zones = new bool[numZones]; + + fill(inyet, inyet + numZones, 0); + fill(inyet2, inyet2 + numZones, 0); + fill(done_zones, done_zones + numZones, false); + + nhl = 0; +#pragma omp for schedule(dynamic,100) + for (int ii = 0; ii < numZones; ii++) + { + int nhlcount = 0; + int h = iord[ii]; + float lowvol; + bool beaten; + + for (int hl = 0; hl < nhl; hl++) + inyet[zonelist[hl]] = 0; + + zonelist[0] = h; + inyet[h] = 1; + nhl = 1; + z[h].npjoin = z[h].np; + do { + /* Find the lowest-volume (highest-density) adjacency */ + int nl = 0; + + beaten = false; + lowvol = BIGFLT; + + for (int hl = 0; hl < nhl; hl++) + { + int h2 = zonelist[hl]; + if (inyet[h2] == 1) { /* If it's not already identified as + an interior zone, with inyet=2 */ + bool interior = true; + for (int za = 0; za < z[h2].nadj; za++) + { + if (inyet[z[h2].adj[za]] == 0) + { + interior = false; + if (z[h2].slv[za] == lowvol) + { + links[nl] = z[h2].adj[za]; + nl ++; + if (nl == NLINKS) + { + printf("Too many links with the same rho_sl! Increase NLINKS from %d\n",nl); + exit(0); + } + } + if (z[h2].slv[za] < lowvol) + { + lowvol = z[h2].slv[za]; + links[0] = z[h2].adj[za]; + nl = 1; + } + } + } + if (interior) + inyet[h2] = 2; /* No bordering exter. zones */ + } + } + + if (nl == 0) + { + beaten = true; + z[h].leak = maxvol; + continue; + } + + if (lowvol > voltol) + { + beaten = true; + z[h].leak = lowvol; + continue; + } + + for (int l = 0; l < nl; l++) + { + //if (!done_zones[links[l]]) /* Equivalent to p[z[links[l]].core].dens < p[z[h].core].dens) as zones are sorted. */ + if (p[z[links[l]].core].dens < p[z[h].core].dens) + beaten = true; + } + + if (beaten) + { + z[h].leak = lowvol; + continue; + } + + /* Add everything linked to the link(s) */ + int nhl2 = 0; + for (int l = 0; l < nl; l++) + { + if (inyet2[links[l]] == 0) + { + zonelist2[nhl2] = links[l]; + inyet2[links[l]] = 1; + nhl2 ++; + 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 (!done_zones[link2]) { // Equivalent to p[z[link2].core].dens < p[z[h].core].dens) + if (p[z[link2].core].dens < p[z[h].core].dens) + beaten = true; + break; + } + zonelist2[nhl2] = link2; + inyet2[link2] = 1; + nhl2++; + added = true; + } + } + } + if (interior) + inyet2[h2] = 2; + } + } + } + } + } + for (int hl = 0; hl < nhl2; hl++) + inyet2[zonelist2[hl]] = 0; + + /* See if there's a beater */ + if (beaten) { + z[h].leak = lowvol; + } else { + for (int h2 = 0; h2 < nhl2; h2++) { + zonelist[nhl] = zonelist2[h2]; + inyet[zonelist2[h2]] = 1; + nhl++; + z[h].npjoin += z[zonelist2[h2]].np; + } + } + if (nhl/10000 > nhlcount) { + nhlcount = nhl/10000; + printf(" %d",nhl); FF; + } + } while((lowvol < BIGFLT) && (!beaten)); + + 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].voljoin = 0.; + z[h].zonelist = new int[nhl]; + z[h].numzones = nhl; + for (int q = 0; q < nhl; q++) { + z[h].voljoin += z[zonelist[q]].vol; + z[h].zonelist[q] = zonelist[q]; + } + + done_zones[h] = true; + z[h].nhl = nhl; + } + delete[] zonelist; + delete[] zonelist2; + delete[] links; + delete[] inyet; + delete[] inyet2; + delete[] done_zones; + + } + delete[] iord; + + 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/zobov2/jozov2_zones.cpp b/c_tools/zobov2/jozov2_zones.cpp new file mode 100644 index 0000000..ca59768 --- /dev/null +++ b/c_tools/zobov2/jozov2_zones.cpp @@ -0,0 +1,197 @@ +#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(); + + numZones = 0; + for (pid_t i = 0; i < np; i++) + if (numinh[i] > 0) + 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++; + } + + try + { + for (int h = 0; h < numZones; h++) { + zt[h].adj = new pid_t[zt[h].nadj]; + zt[h].slv = new float[zt[h].nadj]; + zt[h].nadj = 0; + } + } + catch(const std::bad_alloc& a) + { + cout << "Could not allocate memory for zone adjacencies." << 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; + 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; + 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; + zonenum[i] = h; + h++; + } else { + zonenum[i] = -1; + } + } + + buildZoneAdjacencies(p, np, z, zt, + nzones, 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/zobov2/zobov.hpp b/c_tools/zobov2/zobov.hpp new file mode 100644 index 0000000..832f071 --- /dev/null +++ b/c_tools/zobov2/zobov.hpp @@ -0,0 +1,34 @@ +#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/zobov/jozov_persistent.c b/zobov/jozov_persistent.c deleted file mode 100644 index 9a76e8a..0000000 --- a/zobov/jozov_persistent.c +++ /dev/null @@ -1,609 +0,0 @@ -#include -/* jovoz.c by Mark Neyrinck */ -#include -#include -#include - -#define BIGFLT 1e30 /* Biggest possible floating-point number */ -#define NLINKS 100 /* Number of possible links with the same rho_sl */ -#define FF fflush(stdout) - -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 */ - float smallest_saddle; /* Smallest saddle */ -} ZONE; - -typedef struct ZoneT { - int nadj; /* Number of zones on border */ - int *adj; /* Each adjacent zone, with ... */ - float *slv; /* Smallest Linking Volume */ - float smallest_saddle; /* Smallest saddle point value */ -} ZONET; - -void findrtop(double *a, int na, int *iord, int nb); - -int main(int argc,char **argv) { - - FILE *adj, *vol, *zon, *zon2, *txt; - PARTICLE *p; - ZONE *z; - ZONET *zt; - char *adjfile, *volfile, *zonfile, *zonfile2, *txtfile; - int i, j,k,l, h, h2,hl,n,np, np2, nzones, nhl, nhlcount, nhl2; - int *jumper, *jumped, *numinh; - int *zonenum, *zonelist, *zonelist2; - int link[NLINKS], link2, nl; - float lowvol, voltol, prob; - - int q,q2; - - int za, nin; - int testpart; - char already, interior, *inyet, *inyet2, added, beaten; - int *nm, **m; - - float maxvol, minvol; - double *sorter, e1,maxdenscontrast; - int *iord; - - e1 = exp(1.)-1.; - - if (argc != 7) { - 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\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 (voltol <= 0.) { - printf("Proceeding without a density threshold.\n"); - voltol = 1e30; - } - - adj = fopen(adjfile, "r"); - if (adj == NULL) { - printf("Unable to open %s\n",adjfile); - exit(0); - } - fread(&np,1, sizeof(int),adj); - - printf("adj: %d particles\n", np); - FF; - - p = (PARTICLE *)malloc(np*sizeof(PARTICLE)); - /* Adjacencies*/ - for (i=0;i 0) - p[i].adj = (int *)malloc(p[i].nadj*sizeof(int)); - else p[i].adj = 0; - p[i].ncnt = 0; /* Temporarily, it's an adj counter */ - } - for (i=0;i 0) - for (k=0;k 1e30)) { - printf("Whacked-out volume found, of particle %d: %f\n",i,p[i].dens); - p[i].dens = 1.; - } - p[i].dens = 1./p[i].dens; /* Get density from volume */ - } - fclose(vol); - - jumped = (int *)malloc(np*sizeof(int)); - jumper = (int *)malloc(np*sizeof(int)); - numinh = (int *)malloc(np*sizeof(int)); - - /* find jumper */ - for (i = 0; i < np; i++) { - minvol = p[i].dens; jumper[i] = -1; - for (j=0; j< p[i].nadj; j++) { - if (p[p[i].adj[j]].dens < minvol) { - jumper[i] = p[i].adj[j]; - minvol = p[jumper[i]].dens; - } - } - numinh[i] = 0; - } - - printf("About to jump ...\n"); FF; - - /* Jump */ - for (i = 0; i < np; i++) { - jumped[i] = i; - while (jumper[jumped[i]] > -1) - jumped[i] = jumper[jumped[i]]; - numinh[jumped[i]]++; - } - printf("Post-jump ...\n"); FF; - - nzones = 0; - for (i = 0; i < np; i++) - if (numinh[i] > 0) nzones++; - printf("%d initial zones found\n", nzones); - - z = (ZONE *)malloc(nzones*sizeof(ZONE)); - if (z == NULL) { - printf("Unable to allocate z\n"); - exit(0); - } - zt = (ZONET *)malloc(nzones*sizeof(ZONET)); - if (zt == NULL) { - printf("Unable to allocate zt\n"); - exit(0); - } - zonenum = (int *)malloc(np*sizeof(int)); - if (zonenum == NULL) { - printf("Unable to allocate zonenum\n"); - exit(0); - } - - h = 0; - for (i = 0; i < np; i++) - if (numinh[i] > 0) { - z[h].core = i; - zonenum[i] = h; - h++; - } else { - zonenum[i] = -1; - } - - /* Count border particles */ - for (i = 0; i < np; i++) - for (j = 0; j < p[i].nadj; j++) { - testpart = p[i].adj[j]; - if (jumped[i] != jumped[testpart]) - zt[zonenum[jumped[i]]].nadj++; - } - - for (h=0;h p[i].dens) { - /* there could be a weakest link through testpart */ - already = 0; - for (za = 0; za < zt[h].nadj; za++) - if (zt[h].adj[za] == zonenum[jumped[testpart]]) { - already = 1; - if (p[testpart].dens < zt[h].slv[za]) { - zt[h].slv[za] = p[testpart].dens; - } - } - if (already == 0) { - 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 = 0; - for (za = 0; za < zt[h].nadj; za++) - if (zt[h].adj[za] == zonenum[jumped[testpart]]) { - already = 1; - if (p[i].dens < zt[h].slv[za]) { - zt[h].slv[za] = p[i].dens; - } - } - if (already == 0) { - int q = zt[h].nadj; - - zt[h].adj[q] = zonenum[jumped[testpart]]; - zt[h].slv[q] = p[i].dens; - if (zt[h].smallest_saddle > zt[h].slv[q]) - zt[h].smallest_saddle = zt[h].slv[q]; - - zt[h].nadj++; - } - } - } - } - } - printf("Found zone adjacencies\n"); FF; - - /* Free particle adjacencies */ - for (i=0;i maxvol) maxvol = p[i].dens; - if (p[i].dens < minvol) minvol = p[i].dens; - } - printf("Densities range from %e to %e.\n",minvol,maxvol);FF; - - zon2 = fopen(zonfile2,"w"); - if (zon2 == NULL) { - printf("Problem opening zonefile %s.\n\n",zonfile2); - exit(0); - } - fwrite(&nzones,1,4,zon2); - - -#define CORE_DENSITY(a) p[z[a].core].dens -#define SADDLE_DENSITY(a) z[a].smallest_saddle -#define COMPARE_ZONES_BEATEN(a,b) (CORE_DENSITY(a) < CORE_DENSITY(b)) - - - for (h = 0; h voltol) { - beaten = 1; - z[h].leak = lowvol; - continue; - } - - - /* Check if among the found adjacent candidates there is a beater */ - for (l=0; l < nl; l++) - if (COMPARE_ZONES_BEATEN(link[l], h)) - beaten = 1; - - if (beaten == 1) { - z[h].leak = lowvol; - continue; - } - - /* Add everything linked to the link(s) */ - nhl2 = 0; - for (l=0; l < nl; l++) { - - if (inyet2[link[l]] == 0) { - zonelist2[nhl2] = link[l]; - inyet2[link[l]] = 1; - nhl2 ++; - added = 1; - - while ((added == 1) && (beaten == 0)) { - - added = 0; - - /* Loop over the newly found zones */ - for (hl = 0; (hl < nhl2) && (beaten == 0); hl++) { - h2 = zonelist2[hl]; - - /* This hole has already been added to zonelist2 ? */ - if (inyet2[h2] == 1) { - /* No, process it */ - - interior = 1; /* Guilty until proven innocent */ - - /* Take all adjacent zones of h2 */ - for (za = 0; za < z[h2].nadj; za ++) { - - link2 = z[h2].adj[za]; - - /* If this zone (link2) has not been already processed, process it */ - if ((inyet[link2]+inyet2[link2]) == 0) { - - interior = 0; - - /* Does it have a low saddle point ? */ - if (z[h2].slv[za] <= lowvol) { - /* Is this void beaten by the zone link2 ? */ - if (COMPARE_ZONES_BEATEN(link2, h)) { - /* Yes. Stop here. */ - beaten = 1; - break; - } - /* Insert it. */ - zonelist2[nhl2] = link2; - inyet2[link2] = 1; - nhl2++; - /* Note that something was added so please continue */ - added = 1; - } - } - } - if (interior == 1) - inyet2[h2] = 2; - } - } - } - - } - } - for (hl = 0; hl < nhl2; hl++) - inyet2[zonelist2[hl]] = 0; - - /* See if there's a beater */ - if (beaten == 1) { - z[h].leak = lowvol; - } else { - for (h2 = 0; h2 < nhl2; h2++) { - zonelist[nhl] = zonelist2[h2]; - inyet[zonelist2[h2]] = 1; - nhl++; - z[h].npjoin += z[zonelist2[h2]].np; - } - } - if (nhl/10000 > nhlcount) { - nhlcount = nhl/10000; - printf(" %d",nhl); FF; - } - } while((lowvol < BIGFLT) && (beaten == 0)); - - z[h].denscontrast = z[h].leak/p[z[h].core].dens; - if (z[h].denscontrast < 1.) - z[h].denscontrast = 1.; - - /* find biggest denscontrast */ - if (z[h].denscontrast > maxdenscontrast) { - maxdenscontrast = (double)z[h].denscontrast; - } - - /* 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].voljoin = 0.; - for (q = 0; q= QHULL_MAX_PARTICLES) { - printf("Too many particles to tesselate per division (Qhull will overflow). Increase divisions or reduce buffer size."); + printf("Too many particles to tesselate per division (Qhull will overflow). Increase divisions or reduce buffer size.\n"); fflush(stdout); exit(1); }