From 9b0e7e115ea2576048fadbb14ce8e4394ec2a96c Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 1 Apr 2013 09:32:09 -0400 Subject: [PATCH] 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; + } +}