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