guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 15:21:11 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Splitted jozov2 into several files for better readability

Browse source
This commit is contained in:
Guilhem Lavaux 2013-04-01 09:32:09 -04:00
parent 60e692a0cb
commit 9b0e7e115e
5 changed files with 651 additions and 583 deletions
Download patch file Download diff file Expand all files Collapse all files

586
c_tools/zobov2/jozov2.cpp
View file

@ -7,594 +7,13 @@
#include <cmath>
#include <assert.h>
/* 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) */

47
c_tools/zobov2/jozov2.hpp Normal file
View file

@ -0,0 +1,47 @@
#ifndef __JOZOV2_HPP
#define __JOZOV2_HPP
#include <string>
#include <exception>
#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

182
c_tools/zobov2/jozov2_io.cpp Normal file
View file

@ -0,0 +1,182 @@
#include <iostream>
#include <fstream>
#include <string>
#include <boost/format.hpp>
#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));
}
}

222
c_tools/zobov2/jozov2_watershed.cpp Normal file
View file

@ -0,0 +1,222 @@
#include <iostream>
#include <fstream>
#include <boost/format.hpp>
#include <string>
#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;
}

197
c_tools/zobov2/jozov2_zones.cpp Normal file
View file

@ -0,0 +1,197 @@
#include <iostream>
#include <fstream>
#include "jozov2.hpp"
#include "zobov.hpp"
#include <boost/format.hpp>
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;
}
}