From 91a222a3c3d073552ce84e8a54ca3759d6537a7d Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 1 Apr 2013 21:30:45 -0500 Subject: [PATCH 1/4] Re-Implemented jozov2 using priority_queues (far less complicated and maybe better performing) --- c_tools/zobov2/jozov2.cpp | 4 +- c_tools/zobov2/jozov2_watershed.cpp | 244 ++++++++++++++-------------- 2 files changed, 124 insertions(+), 124 deletions(-) diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2.cpp index 718b8e9..b33a906 100644 --- a/c_tools/zobov2/jozov2.cpp +++ b/c_tools/zobov2/jozov2.cpp @@ -139,10 +139,10 @@ int main(int argc,char **argv) if (z[i].np == 1) continue; - txt << format("%d %d %d %e %e %d %d %e %d %f %6.2e") + txt << format("%d %d %d %e %e %d %d %e %d %f %6.2e %f") % (h+1) % i % z[i].core % p[z[i].core].dens % z[i].vol % z[i].np % z[i].nhl % z[i].voljoin % z[i].npjoin - % z[i].denscontrast % prob << endl; + % z[i].denscontrast % prob % z[i].leak << endl; } /* h+1 to start from 1, not zero */ txt.close(); diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp index d6ef31d..3c46e97 100644 --- a/c_tools/zobov2/jozov2_watershed.cpp +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -2,6 +2,7 @@ #include #endif +#include #include #include #include @@ -13,6 +14,40 @@ using namespace std; using boost::format; +struct ZoneDensityPair +{ + int h; + double density; + + bool operator<(const ZoneDensityPair& p2) const + { + return density > p2.density; + } +}; + +typedef priority_queue ZoneQueue; + +static void build_process_queue(ZoneQueue& q, ZONE *z, char *inyet, int h) +{ + ZoneDensityPair zdp; + ZONE& z_h = z[h]; + bool interior = true; + + assert(inyet[h] == 1); + for (int za = 0; za < z_h.nadj; za++) + { + zdp.h = z_h.adj[za]; + zdp.density = z_h.slv[za]; + if (inyet[zdp.h] == 0) + { + q.push(zdp); + interior = false; + } + } + if (interior) + inyet[h] = 2; +} + void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol) { int *iord; @@ -28,36 +63,33 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo iord = new int[numZones]; - findrtop(sorter, numZones, iord, numZones); +// findrtop(sorter, numZones, iord, numZones); delete[] sorter; -#pragma omp parallel + //#pragma omp parallel { char *inyet, *inyet2; int *zonelist, *zonelist2; int nhl; - int *links = new int[NLINKS]; - bool *done_zones; int prev_ii = -1; inyet = new char[numZones]; inyet2 = new char[numZones]; zonelist = new int[numZones]; zonelist2 = new int[numZones]; - done_zones = new bool[numZones]; fill(inyet, inyet + numZones, 0); fill(inyet2, inyet2 + numZones, 0); - fill(done_zones, done_zones + numZones, false); nhl = 0; -#pragma omp for schedule(dynamic,1) + //#pragma omp for schedule(dynamic,1) for (int h = 0; h < numZones; h++) { int nhlcount = 0; - float lowvol; + float lowvol, z_cur_core_dens; bool beaten; - set to_process; + priority_queue to_process; + int link0; for (int hl = 0; hl < nhl; hl++) inyet[zonelist[hl]] = 0; @@ -65,144 +97,113 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo zonelist[0] = h; inyet[h] = 1; nhl = 1; - to_process.insert(h); z[h].npjoin = z[h].np; + z_cur_core_dens = p[z[h].core].dens; + + build_process_queue(to_process, z, inyet, h); + do { /* Find the lowest-volume (highest-density) adjacency */ int nl = 0; beaten = false; - lowvol = BIGFLT; - set::iterator iter = to_process.begin(); - while (iter != to_process.end()) - { - int h2 = *iter; - ZONE& z_cur = z[h2]; - bool interior = true, touched = false; - assert(inyet[h2] == 1); - - for (int za = 0; za < z_cur.nadj; za++) - { - if (inyet[z_cur.adj[za]] == 0) - { - interior = false; - if (z_cur.slv[za] == lowvol) - { - links[nl] = z_cur.adj[za]; - touched = true; - nl ++; - if (nl == NLINKS) - { - printf("Too many links with the same rho_sl! Increase NLINKS from %d\n",nl); - exit(0); - } - } - else - if (z_cur.slv[za] < lowvol) - { - lowvol = z_cur.slv[za]; - links[0] = z_cur.adj[za]; - nl = 1; - touched = true; - } - } - } - ++iter; - if (interior || !touched) - to_process.erase(h2); - if (interior) - inyet[h2] = 2; /* No bordering exter. zones */ - } - - if (nl == 0) + if (to_process.empty()) { beaten = true; z[h].leak = maxvol; continue; } + + do + { + lowvol = to_process.top().density; + link0 = to_process.top().h; + to_process.pop(); + } + while ((inyet[link0] != 0) && (!to_process.empty())); + + if (to_process.empty()) + { + beaten = true; + z[h].leak = maxvol; + continue; + } + + nl++; - if (lowvol > voltol) + if (lowvol > voltol) { beaten = true; z[h].leak = lowvol; continue; } - for (int l = 0; l < nl; l++) - { - //if (!done_zones[links[l]]) /* Equivalent to p[z[links[l]].core].dens < p[z[h].core].dens) as zones are sorted. */ - if (p[z[links[l]].core].dens < p[z[h].core].dens) - beaten = true; - } - - if (beaten) + if (p[z[link0].core].dens < z_cur_core_dens) { + beaten = true; z[h].leak = lowvol; continue; } /* Add everything linked to the link(s) */ int nhl2 = 0; - for (int l = 0; l < nl; l++) - { - if (inyet2[links[l]] == 0) - { - zonelist2[nhl2] = links[l]; - inyet2[links[l]] = 1; - nhl2 ++; - bool added = true; - while (added && !beaten) - { - added = false; - for (int hl = 0; (hl < nhl2) && (!beaten); hl++) - { - int h2 = zonelist2[hl]; - - if (inyet2[h2] == 1) { - bool interior = true; /* Guilty until proven innocent */ - for (int za = 0; za < z[h2].nadj; za ++) { - int link2 = z[h2].adj[za]; - if ((inyet[link2]+inyet2[link2]) == 0) { - interior = false; - if (z[h2].slv[za] <= lowvol) { - //if (!done_zones[link2]) { // Equivalent to p[z[link2].core].dens < p[z[h].core].dens) - if (p[z[link2].core].dens < p[z[h].core].dens) { - beaten = true; - break; - } - zonelist2[nhl2] = link2; - inyet2[link2] = 1; - nhl2++; - added = true; - } - } - } - if (interior) - inyet2[h2] = 2; - } - } - } - } - } - - /* See if there's a beater */ - if (beaten) { - z[h].leak = lowvol; - } else { - for (int h2 = 0; h2 < nhl2; h2++) { - int new_h = zonelist2[h2]; - + zonelist2[0] = link0; + inyet2[link0] = 1; + nhl2=1; + + 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 < z_cur_core_dens) { + beaten = true; + break; + } + zonelist2[nhl2] = link2; + inyet2[link2] = 1; + nhl2++; + added = true; + } + } + } + + if (interior) + inyet2[h2] = 2; + } + } + } + + /* See if there's a beater */ + if (beaten) { + z[h].leak = lowvol; + } else { + for (int h2 = 0; h2 < nhl2; h2++) { + int new_h = zonelist2[h2]; + zonelist[nhl] = new_h; assert(inyet[new_h] == 0); inyet[new_h] = 1; if (inyet2[new_h] != 2) - to_process.insert(new_h); + build_process_queue(to_process, z, inyet, new_h); nhl++; z[h].npjoin += z[new_h].np; } } - + for (int hl = 0; hl < nhl2; hl++) inyet2[zonelist2[hl]] = 0; if (nhl/10000 > nhlcount) { @@ -212,10 +213,12 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo (cout << format(" %d [%d]") % nhl % to_process.size()).flush(); nhlcount = nhl/10000; } - } while((lowvol < BIGFLT) && (!beaten)); - - z[h].denscontrast = z[h].leak/p[z[h].core].dens; - if (z[h].denscontrast < 1.) z[h].denscontrast = 1.; + } + while((lowvol < BIGFLT) && (!beaten)); + + z[h].denscontrast = z[h].leak/p[z[h].core].dens; + if (z[h].denscontrast < 1.) + z[h].denscontrast = 1.; /* Don't sort; want the core zone to be first */ @@ -233,19 +236,16 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo z[h].zonelist[q] = zonelist[q]; } - done_zones[h] = true; z[h].nhl = nhl; } delete[] zonelist; delete[] zonelist2; - delete[] links; delete[] inyet; delete[] inyet2; - delete[] done_zones; - + } delete[] iord; - + double maxdenscontrast = 0; #pragma omp parallel shared(maxdenscontrast) { From 5c40a6c2268d4fbc4f14e58c40130a64a4a45f5a Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 2 Apr 2013 09:06:23 -0500 Subject: [PATCH 2/4] Allow for a particle duplicate in voztie if and only if the adjacencies are strictly equal. --- zobov/voztie.c | 35 +++++++++++++++++++++++++++++------ 1 file changed, 29 insertions(+), 6 deletions(-) diff --git a/zobov/voztie.c b/zobov/voztie.c index 4aaaf57..5a15448 100644 --- a/zobov/voztie.c +++ b/zobov/voztie.c @@ -2,8 +2,21 @@ #include #include #include +#include #include "voz.h" +int compare_int(const void *a, const void *b) +{ + const int *ia = (const int *)a; + const int *ib = (const int *)b; + if ((*ia) < (*ib)) + return -1; + else if ((*ia) > (*ib)) + return 1; + else + return 0; +} + int main(int argc, char *argv[]) { FILE *part, *adj, *vol; @@ -135,15 +148,25 @@ int main(int argc, char *argv[]) { fread(&na,1,sizeof(int),part); pid_t pid = orig[p]; if (na > 0) { - assert(adjs[pid].nadj == 0);// || adjs[pid].nadj == na); - adjs[orig[p]].nadj = na; - if (adjs[pid].adj == 0) - adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t)); - if (adjs[orig[p]].adj == 0) { + pid_t *previous_adj = adjs[pid].adj; + assert(adjs[pid].nadj == 0 || adjs[pid].nadj == na); + adjs[pid].nadj = na; + adjs[pid].adj = (pid_t *)malloc(na*sizeof(pid_t)); + if (adjs[pid].adj == 0) { printf("Couldn't allocate adjs[orig[%d]].adj.\n",p); exit(0); } - fread(adjs[orig[p]].adj,na,sizeof(pid_t),part); + fread(adjs[pid].adj,na,sizeof(pid_t),part); + if (previous_adj != 0) + { + qsort(previous_adj, na, sizeof(pid_t), compare_int); + qsort(adjs[pid].adj, na, sizeof(pid_t), compare_int); + if (memcmp(previous_adj, adjs[pid].adj, sizeof(pid_t)*na) != 0) + { + printf("Inconsistent adjacencies between two divisions.\n"); + abort(); + } + } } else { printf("0"); fflush(stdout); } From 85546379cf7b7a3848fba79fcf4d9bc16c087255 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 2 Apr 2013 09:06:56 -0500 Subject: [PATCH 3/4] Cleanup some development cruft and re-enable openmp clause in jozov2 --- c_tools/zobov2/jozov2.cpp | 4 ++-- c_tools/zobov2/jozov2_watershed.cpp | 19 ++----------------- c_tools/zobov2/jozov2_zones.cpp | 7 +++++-- 3 files changed, 9 insertions(+), 21 deletions(-) diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2.cpp index b33a906..718b8e9 100644 --- a/c_tools/zobov2/jozov2.cpp +++ b/c_tools/zobov2/jozov2.cpp @@ -139,10 +139,10 @@ int main(int argc,char **argv) if (z[i].np == 1) continue; - txt << format("%d %d %d %e %e %d %d %e %d %f %6.2e %f") + txt << format("%d %d %d %e %e %d %d %e %d %f %6.2e") % (h+1) % i % z[i].core % p[z[i].core].dens % z[i].vol % z[i].np % z[i].nhl % z[i].voljoin % z[i].npjoin - % z[i].denscontrast % prob % z[i].leak << endl; + % z[i].denscontrast % prob << endl; } /* h+1 to start from 1, not zero */ txt.close(); diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp index 3c46e97..c66ba48 100644 --- a/c_tools/zobov2/jozov2_watershed.cpp +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -50,23 +50,9 @@ static void build_process_queue(ZoneQueue& q, ZONE *z, char *inyet, int h) void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol) { - int *iord; - - double *sorter = new double[numZones+1]; - /* Assign sorter by probability (could use volume instead) */ - for (int h = 0; h < numZones; h++) - sorter[h] = (double)z[h].core; - /* Text output file */ - printf("about to sort (pre-watershed sort) ...\n");FF; - - iord = new int[numZones]; - -// findrtop(sorter, numZones, iord, numZones); - delete[] sorter; - - //#pragma omp parallel +#pragma omp parallel { char *inyet, *inyet2; int *zonelist, *zonelist2; @@ -82,7 +68,7 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo fill(inyet2, inyet2 + numZones, 0); nhl = 0; - //#pragma omp for schedule(dynamic,1) +#pragma omp for schedule(dynamic,1) for (int h = 0; h < numZones; h++) { int nhlcount = 0; @@ -244,7 +230,6 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo delete[] inyet2; } - delete[] iord; double maxdenscontrast = 0; #pragma omp parallel shared(maxdenscontrast) diff --git a/c_tools/zobov2/jozov2_zones.cpp b/c_tools/zobov2/jozov2_zones.cpp index ca59768..ebfc974 100644 --- a/c_tools/zobov2/jozov2_zones.cpp +++ b/c_tools/zobov2/jozov2_zones.cpp @@ -39,10 +39,13 @@ void buildInitialZones(PARTICLE *p, pid_t np, pid_t* jumped, } (cout << "Post-jump ..." << endl).flush(); - numZones = 0; + pid_t loc_NumZones = 0; +#pragma omp parallel for schedule(static) reduction(+:loc_NumZones) for (pid_t i = 0; i < np; i++) if (numinh[i] > 0) - numZones++; + loc_NumZones++; + + numZones = loc_NumZones; cout << format("%d initial zones found") % numZones << endl; delete[] jumper; From 96b0ba31545a2fab6688134eba39912a84a347b2 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 4 Apr 2013 12:04:52 -0500 Subject: [PATCH 4/4] Fixed erroneous handling of zones with exact same leakage value in jozov2. --- c_tools/zobov2/jozov2_watershed.cpp | 52 ++++++++++++++++++++--------- 1 file changed, 36 insertions(+), 16 deletions(-) diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp index c66ba48..365e836 100644 --- a/c_tools/zobov2/jozov2_watershed.cpp +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -18,16 +18,17 @@ struct ZoneDensityPair { int h; double density; + double core; bool operator<(const ZoneDensityPair& p2) const { - return density > p2.density; + return (density > p2.density) || (density==p2.density && core > p2.core); } }; typedef priority_queue ZoneQueue; -static void build_process_queue(ZoneQueue& q, ZONE *z, char *inyet, int h) +static void build_process_queue(ZoneQueue& q, ZONE *z, PARTICLE *p, char *inyet, int h) { ZoneDensityPair zdp; ZONE& z_h = z[h]; @@ -38,6 +39,7 @@ static void build_process_queue(ZoneQueue& q, ZONE *z, char *inyet, int h) { zdp.h = z_h.adj[za]; zdp.density = z_h.slv[za]; + zdp.core = p[z[zdp.h].core].dens; if (inyet[zdp.h] == 0) { q.push(zdp); @@ -72,32 +74,34 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo for (int h = 0; h < numZones; h++) { int nhlcount = 0; - float lowvol, z_cur_core_dens; + float previous_lowvol = BIGFLT, lowvol, z_cur_core_dens; bool beaten; priority_queue to_process; int link0; + int save_nhl; + pid_t save_npjoin; 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; + save_nhl = nhl = 1; + save_npjoin = z[h].npjoin = z[h].np; z_cur_core_dens = p[z[h].core].dens; - build_process_queue(to_process, z, inyet, h); + build_process_queue(to_process, z, p, inyet, h); do { /* Find the lowest-volume (highest-density) adjacency */ - int nl = 0; - beaten = false; if (to_process.empty()) { beaten = true; z[h].leak = maxvol; + save_npjoin = z[h].npjoin; + save_nhl = nhl; continue; } @@ -111,13 +115,21 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo if (to_process.empty()) { + save_npjoin = z[h].npjoin; + save_nhl = nhl; beaten = true; z[h].leak = maxvol; continue; } - nl++; - + /* See if there's a beater */ + if (previous_lowvol != lowvol) + { + save_npjoin = z[h].npjoin; + save_nhl = nhl; + previous_lowvol = lowvol; + } + if (lowvol > voltol) { beaten = true; @@ -177,16 +189,17 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo if (beaten) { z[h].leak = lowvol; } else { + for (int h2 = 0; h2 < nhl2; h2++) { int new_h = zonelist2[h2]; zonelist[nhl] = new_h; assert(inyet[new_h] == 0); + z[h].npjoin += z[new_h].np; inyet[new_h] = 1; if (inyet2[new_h] != 2) - build_process_queue(to_process, z, inyet, new_h); + build_process_queue(to_process, z, p, inyet, new_h); nhl++; - z[h].npjoin += z[new_h].np; } } @@ -201,6 +214,12 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo } } while((lowvol < BIGFLT) && (!beaten)); + + if (!beaten) + { + save_npjoin = z[h].npjoin; + save_nhl = nhl; + } z[h].denscontrast = z[h].leak/p[z[h].core].dens; if (z[h].denscontrast < 1.) @@ -214,15 +233,16 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo FF; } /* Calculate volume */ + z[h].npjoin = save_npjoin; z[h].voljoin = 0.; - z[h].zonelist = new int[nhl]; - z[h].numzones = nhl; - for (int q = 0; q < nhl; q++) { + z[h].zonelist = new int[save_nhl]; + z[h].numzones = save_nhl; + for (int q = 0; q < save_nhl; q++) { z[h].voljoin += z[zonelist[q]].vol; z[h].zonelist[q] = zonelist[q]; } - z[h].nhl = nhl; + z[h].nhl = save_nhl; } delete[] zonelist; delete[] zonelist2;