From 9c329086aff70226dcfb13a658cacde8b493c491 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 30 Mar 2013 15:26:37 -0500 Subject: [PATCH 01/11] Fixed computation of the maximum number of particles accepted by Qhull --- zobov/vozinit.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/zobov/vozinit.c b/zobov/vozinit.c index c93f3f6..c1e7865 100644 --- a/zobov/vozinit.c +++ b/zobov/vozinit.c @@ -5,7 +5,7 @@ #define DL for (d=0;d<3;d++) #define BF 1e30 -#define QHULL_MAX_PARTICLES (1L<<24-1) +#define QHULL_MAX_PARTICLES ((1L<<24)-1) int posread(char *posfile, float ***p, float fact); @@ -149,9 +149,10 @@ int main(int argc, char *argv[]) { printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax); numGuards = 6*(NGUARD+1)*(NGUARD+1); + printf("Total max particles: %d\n" , nvpbufmax+numGuards); if (nvpbufmax+numGuards >= 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); } From 9b7a931db9bfcb23f1534462cb22d8367ae8253f Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 31 Mar 2013 11:08:57 -0400 Subject: [PATCH 02/11] New C++ jozov. Preliminary work before rewrite of the watershed transform algorithm. --- c_tools/zobov2/jozov2.cpp | 701 ++++++++++++++++++++++++++++++++++++++ c_tools/zobov2/zobov.hpp | 31 ++ 2 files changed, 732 insertions(+) create mode 100644 c_tools/zobov2/jozov2.cpp create mode 100644 c_tools/zobov2/zobov.hpp diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2.cpp new file mode 100644 index 0000000..2a04a41 --- /dev/null +++ b/c_tools/zobov2/jozov2.cpp @@ -0,0 +1,701 @@ +#include +#include +#include +#include +#include +#include +#include +#include +/* jovoz.c by Mark Neyrinck */ + +#include "zobov.hpp" + +using namespace std; +using boost::format; + +#define BIGFLT 1e30 /* Biggest possible floating-point number */ +#define NLINKS 1000 /* Number of possible links with the same rho_sl */ +#define FF cout.flush() + +extern "C" void findrtop(double *a, int na, int *iord, int nb); + +class FileError: virtual std::exception +{ +}; + +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 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[i]; + z[h].vol += 1.0/(double)p[i].dens; + } +} + + +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 doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol) +{ + char *inyet, *inyet2; + int *zonelist, *zonelist2; + int nhl; + int *links = new int[NLINKS]; + float maxdenscontrast = 0; + + ofstream zon2(zonfile2.c_str()); + if (!zon2) + { + cout << format("Problem opening zonefile %s.)") % zonfile2 << endl; + throw FileError(); + } + zon2.write((char *)&numZones,sizeof(int)); + + 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; + for (int h = 0; h < numZones; h++) + { + int nhlcount = 0; + 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 (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 (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.; + + /* 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 (int q = 0; q < nhl; q++) { + z[h].voljoin += z[zonelist[q]].vol; + } + + z[h].nhl = nhl; + + zon2.write((char *)&nhl, sizeof(int)); + zon2.write((char *)zonelist, nhl*sizeof(int)); + } + delete[] links; + + cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl; +} + +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); + + sorter = new double[nzones+1]; + + 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(zonfile2, p, np, z, nzones, maxvol, voltol); + + 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 = (int *)malloc(nzones*sizeof(int)); + + findrtop(sorter, nzones, iord, nzones); + + 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 Date: Sun, 31 Mar 2013 11:09:34 -0400 Subject: [PATCH 03/11] Free sorter --- c_tools/zobov2/jozov2.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2.cpp index 2a04a41..e2e0000 100644 --- a/c_tools/zobov2/jozov2.cpp +++ b/c_tools/zobov2/jozov2.cpp @@ -650,8 +650,6 @@ int main(int argc,char **argv) buildZones(p, np, jumped, z, nzones, zonenum); writeZoneFile(zonfile, p, np, z, nzones, zonenum, jumped); - sorter = new double[nzones+1]; - maxvol = 0.; minvol = BIGFLT; for(pid_t i = 0; i < np; i++){ @@ -677,6 +675,7 @@ int main(int argc,char **argv) iord = (int *)malloc(nzones*sizeof(int)); findrtop(sorter, nzones, iord, nzones); + delete[] sorter; txt.open(txtfile.c_str()); txt << format("%d particles, %d voids.") % np % nzones << endl; From af425f345de3f23bfaac32dccc177a1441bed241 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 31 Mar 2013 17:01:44 -0400 Subject: [PATCH 04/11] Better hubble constant decoding in SDF --- c_tools/mock/loaders/sdf_loader.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/c_tools/mock/loaders/sdf_loader.cpp b/c_tools/mock/loaders/sdf_loader.cpp index b716536..bb1f5f3 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)) { From 6a15e595216413b68ba5be16ca0654939f5389c8 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 31 Mar 2013 19:35:03 -0400 Subject: [PATCH 05/11] Sort and use the sorted zone property to change some if statements --- c_tools/zobov2/jozov2.cpp | 41 +++++++++++++++++++++++++++++++++------ c_tools/zobov2/zobov.hpp | 3 +++ 2 files changed, 38 insertions(+), 6 deletions(-) diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2.cpp index e2e0000..14cbf30 100644 --- a/c_tools/zobov2/jozov2.cpp +++ b/c_tools/zobov2/jozov2.cpp @@ -314,6 +314,8 @@ void buildZones(PARTICLE *p, pid_t np, pid_t *&jumped, for (pid_t i=0; i < np; i++) { int h = zonenum[i]; z[h].vol += 1.0/(double)p[i].dens; + z[h].numzones = 0; + z[h].zonelist = 0; } } @@ -373,7 +375,9 @@ void doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, in int *zonelist, *zonelist2; int nhl; int *links = new int[NLINKS]; + int *iord; float maxdenscontrast = 0; + bool *done_zones; ofstream zon2(zonfile2.c_str()); if (!zon2) @@ -387,20 +391,37 @@ void doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, in inyet2 = new char[numZones]; zonelist = new int[numZones]; zonelist2 = new int[numZones]; + done_zones = new bool[nzones]; fill(inyet, inyet + numZones, 0); fill(inyet2, inyet2 + numZones, 0); + fill(done_zones, done_zones + numZones, false); + + 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].core; + + /* Text output file */ + + printf("about to sort (pre-watershed sort) ...\n");FF; + + iord = new int[nzones]; + + findrtop(sorter, numZones, iord, numZones); + delete[] sorter; nhl = 0; - for (int h = 0; h < numZones; h++) + for (int ii = 0; ii < numZones; h++) { 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; @@ -462,7 +483,7 @@ void doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, in for (int l = 0; l < nl; l++) { - if (p[z[links[l]].core].dens < p[z[h].core].dens) + if (!done_zones[links[l]]) /* Equivalent to p[z[links[l]].core].dens < p[z[h].core].dens) as zones are sorted. */ beaten = true; } @@ -496,7 +517,7 @@ void doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, in if ((inyet[link2]+inyet2[link2]) == 0) { interior = false; if (z[h2].slv[za] <= lowvol) { - if (p[z[link2].core].dens < p[z[h].core].dens) { + if (!done_zones[link2]) { // Equivalent to p[z[link2].core].dens < p[z[h].core].dens) beaten = true; break; } @@ -550,10 +571,14 @@ void doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, in } /* 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; zon2.write((char *)&nhl, sizeof(int)); @@ -672,7 +697,7 @@ int main(int argc,char **argv) printf("about to sort ...\n");FF; - iord = (int *)malloc(nzones*sizeof(int)); + iord = new int[nzones]; findrtop(sorter, nzones, iord, nzones); delete[] sorter; @@ -694,7 +719,11 @@ int main(int argc,char **argv) } /* h+1 to start from 1, not zero */ txt.close(); - + + delete[] iord; + delete[] z; + delete[] p; + cout << "Done!" << endl; return(0); } diff --git a/c_tools/zobov2/zobov.hpp b/c_tools/zobov2/zobov.hpp index ab9a9ec..832f071 100644 --- a/c_tools/zobov2/zobov.hpp +++ b/c_tools/zobov2/zobov.hpp @@ -20,6 +20,9 @@ typedef struct Zone { 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 { From 60e692a0cb3c7e9c7180d330a623abfc4c4ba6d3 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 1 Apr 2013 09:02:59 -0400 Subject: [PATCH 06/11] Removed old unused code (jozov_persistent). Activated jozov2 cleaned implementation of jozov in C++. --- c_tools/zobov2/findrtop.c | 81 +++++ c_tools/zobov2/jozov2.cpp | 23 +- zobov/jozov_persistent.c | 609 -------------------------------------- 3 files changed, 96 insertions(+), 617 deletions(-) create mode 100644 c_tools/zobov2/findrtop.c delete mode 100644 zobov/jozov_persistent.c 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 index 14cbf30..212fd31 100644 --- a/c_tools/zobov2/jozov2.cpp +++ b/c_tools/zobov2/jozov2.cpp @@ -391,28 +391,28 @@ void doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, in inyet2 = new char[numZones]; zonelist = new int[numZones]; zonelist2 = new int[numZones]; - done_zones = new bool[nzones]; + done_zones = new bool[numZones]; fill(inyet, inyet + numZones, 0); fill(inyet2, inyet2 + numZones, 0); fill(done_zones, done_zones + numZones, false); - sorter = new double[nzones+1]; + double *sorter = new double[numZones+1]; /* Assign sorter by probability (could use volume instead) */ - for (int h = 0; h < nzones; h++) + 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[nzones]; + iord = new int[numZones]; findrtop(sorter, numZones, iord, numZones); delete[] sorter; nhl = 0; - for (int ii = 0; ii < numZones; h++) + for (int ii = 0; ii < numZones; ii++) { int nhlcount = 0; int h = iord[ii]; @@ -580,11 +580,18 @@ void doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, in done_zones[h] = true; z[h].nhl = nhl; - - zon2.write((char *)&nhl, sizeof(int)); - zon2.write((char *)zonelist, nhl*sizeof(int)); } + delete[] zonelist; + delete[] zonelist2; delete[] links; + delete[] iord; + + 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)); + } cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl; } 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 Date: Mon, 1 Apr 2013 09:32:09 -0400 Subject: [PATCH 07/11] Splitted jozov2 into several files for better readability --- c_tools/zobov2/jozov2.cpp | 586 +--------------------------- c_tools/zobov2/jozov2.hpp | 47 +++ c_tools/zobov2/jozov2_io.cpp | 182 +++++++++ c_tools/zobov2/jozov2_watershed.cpp | 222 +++++++++++ c_tools/zobov2/jozov2_zones.cpp | 197 ++++++++++ 5 files changed, 651 insertions(+), 583 deletions(-) create mode 100644 c_tools/zobov2/jozov2.hpp create mode 100644 c_tools/zobov2/jozov2_io.cpp create mode 100644 c_tools/zobov2/jozov2_watershed.cpp create mode 100644 c_tools/zobov2/jozov2_zones.cpp diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2.cpp index 212fd31..718b8e9 100644 --- a/c_tools/zobov2/jozov2.cpp +++ b/c_tools/zobov2/jozov2.cpp @@ -7,594 +7,13 @@ #include #include /* jovoz.c by Mark Neyrinck */ - +#include "jozov2.hpp" #include "zobov.hpp" using namespace std; using boost::format; -#define BIGFLT 1e30 /* Biggest possible floating-point number */ -#define NLINKS 1000 /* Number of possible links with the same rho_sl */ -#define FF cout.flush() -extern "C" void findrtop(double *a, int na, int *iord, int nb); - -class FileError: virtual std::exception -{ -}; - -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 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[i]; - z[h].vol += 1.0/(double)p[i].dens; - z[h].numzones = 0; - z[h].zonelist = 0; - } -} - - -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 doWatershed(const std::string& zonfile2, PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol) -{ - char *inyet, *inyet2; - int *zonelist, *zonelist2; - int nhl; - int *links = new int[NLINKS]; - int *iord; - float maxdenscontrast = 0; - bool *done_zones; - - ofstream zon2(zonfile2.c_str()); - if (!zon2) - { - cout << format("Problem opening zonefile %s.)") % zonfile2 << endl; - throw FileError(); - } - zon2.write((char *)&numZones,sizeof(int)); - - 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); - - 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; - - nhl = 0; - 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. */ - 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) - 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.; - - /* 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.; - 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[] iord; - - 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)); - } - - cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl; -} int main(int argc,char **argv) { @@ -693,7 +112,8 @@ int main(int argc,char **argv) cout << format("Densities range from %e to %e.") % minvol % maxvol << endl; FF; - doWatershed(zonfile2, p, np, z, nzones, maxvol, voltol); + doWatershed(p, np, z, nzones, maxvol, voltol); + writeVoidFile(zonfile2, z, nzones); sorter = new double[nzones+1]; /* Assign sorter by probability (could use volume instead) */ diff --git a/c_tools/zobov2/jozov2.hpp b/c_tools/zobov2/jozov2.hpp new file mode 100644 index 0000000..f23d3c9 --- /dev/null +++ b/c_tools/zobov2/jozov2.hpp @@ -0,0 +1,47 @@ +#ifndef __JOZOV2_HPP +#define __JOZOV2_HPP + +#include +#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..c94de42 --- /dev/null +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -0,0 +1,222 @@ +#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) +{ + char *inyet, *inyet2; + int *zonelist, *zonelist2; + int nhl; + int *links = new int[NLINKS]; + int *iord; + float maxdenscontrast = 0; + 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); + + 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; + + nhl = 0; + 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. */ + 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) + 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.; + + /* 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.; + 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[] iord; + + + 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..32bd739 --- /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[i]; + z[h].vol += 1.0/(double)p[i].dens; + z[h].numzones = 0; + z[h].zonelist = 0; + } +} From b555533fa36c82e41480b4229df479c78d8b0d0d Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 1 Apr 2013 09:39:58 -0400 Subject: [PATCH 08/11] Added openmp-ization of jozov2_watershed --- c_tools/zobov2/jozov2_watershed.cpp | 355 +++++++++++++++------------- 1 file changed, 190 insertions(+), 165 deletions(-) diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp index c94de42..33a7dba 100644 --- a/c_tools/zobov2/jozov2_watershed.cpp +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -1,3 +1,7 @@ +#ifdef OPENMP +#include +#endif + #include #include #include @@ -18,16 +22,6 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo float maxdenscontrast = 0; 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); - double *sorter = new double[numZones+1]; /* Assign sorter by probability (could use volume instead) */ for (int h = 0; h < numZones; h++) @@ -42,181 +36,212 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo findrtop(sorter, numZones, iord, numZones); delete[] sorter; - nhl = 0; - 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; +#pragma omp parallel + { + inyet = new char[numZones]; + inyet2 = new char[numZones]; + zonelist = new int[numZones]; + zonelist2 = new int[numZones]; + done_zones = new bool[numZones]; - 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; + 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++) - { - 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++) + 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. */ + 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) { - if (inyet[z[h2].adj[za]] == 0) + zonelist2[nhl2] = links[l]; + inyet2[links[l]] = 1; + nhl2 ++; + bool added = true; + while (added && !beaten) { - interior = false; - if (z[h2].slv[za] == lowvol) + added = false; + for (int hl = 0; (hl < nhl2) && (!beaten); hl++) { - 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); + 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) + beaten = true; + break; + } + zonelist2[nhl2] = link2; + inyet2[link2] = 1; + nhl2++; + added = true; + } + } } - } - if (z[h2].slv[za] < lowvol) - { - lowvol = z[h2].slv[za]; - links[0] = z[h2].adj[za]; - nl = 1; + if (interior) + inyet2[h2] = 2; + } } } } - if (interior) - inyet[h2] = 2; /* No bordering exter. zones */ + } + 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)); - if (nl == 0) - { - beaten = true; - z[h].leak = maxvol; - continue; - } - - if (lowvol > voltol) - { - beaten = true; - z[h].leak = lowvol; - continue; - } + z[h].denscontrast = z[h].leak/p[z[h].core].dens; + if (z[h].denscontrast < 1.) z[h].denscontrast = 1.; - 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. */ - 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) - 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; - } + /* 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; } - if (nhl/10000 > nhlcount) { - nhlcount = nhl/10000; - printf(" %d",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]; } - } while((lowvol < BIGFLT) && (!beaten)); - - 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; + + done_zones[h] = true; + z[h].nhl = nhl; } - - /* 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; + delete[] zonelist; + delete[] zonelist2; + delete[] links; + delete[] iord; + + } + + 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; + } } - /* 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; + +#pragma omp critical + { + if (maxdenscontrast_local > maxdenscontrast) + maxdenscontrast = maxdenscontrast_local; } - delete[] zonelist; - delete[] zonelist2; - delete[] links; - delete[] iord; - + } cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl; } From 33b32ea43de839488acb15a277e512b918861de1 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 1 Apr 2013 09:43:05 -0400 Subject: [PATCH 09/11] Correct memory cleanup. --- c_tools/zobov2/jozov2_watershed.cpp | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp index 33a7dba..dcbe087 100644 --- a/c_tools/zobov2/jozov2_watershed.cpp +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -14,13 +14,7 @@ using boost::format; void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol) { - char *inyet, *inyet2; - int *zonelist, *zonelist2; - int nhl; - int *links = new int[NLINKS]; int *iord; - float maxdenscontrast = 0; - bool *done_zones; double *sorter = new double[numZones+1]; /* Assign sorter by probability (could use volume instead) */ @@ -38,6 +32,12 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo #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]; @@ -218,11 +218,14 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo delete[] zonelist; delete[] zonelist2; delete[] links; - delete[] iord; - - } + delete[] inyet; + delete[] inyet2; + delete[] done_zones; - maxdenscontrast = 0; + } + delete[] iord; + + double maxdenscontrast = 0; #pragma omp parallel shared(maxdenscontrast) { double maxdenscontrast_local = 0; From d0901e6d95fadb01baf55786c317bcb93daa9329 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 1 Apr 2013 10:15:11 -0400 Subject: [PATCH 10/11] Fixed zone lookup in jozov2 --- c_tools/zobov2/jozov2_zones.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/c_tools/zobov2/jozov2_zones.cpp b/c_tools/zobov2/jozov2_zones.cpp index 32bd739..ca59768 100644 --- a/c_tools/zobov2/jozov2_zones.cpp +++ b/c_tools/zobov2/jozov2_zones.cpp @@ -189,7 +189,7 @@ void buildZones(PARTICLE *p, pid_t np, pid_t *&jumped, delete[] numinh; for (pid_t i=0; i < np; i++) { - int h = zonenum[i]; + int h = zonenum[jumped[i]]; z[h].vol += 1.0/(double)p[i].dens; z[h].numzones = 0; z[h].zonelist = 0; From 45c9d4d62a80bbbac986da2bac1c0090bb14342d Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 1 Apr 2013 10:25:15 -0400 Subject: [PATCH 11/11] Temporarily disable new done_zones array --- c_tools/zobov2/jozov2_watershed.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp index dcbe087..78f6694 100644 --- a/c_tools/zobov2/jozov2_watershed.cpp +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -121,7 +121,8 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo 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 (!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; } @@ -155,7 +156,8 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo 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 (!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; }