From 616b1439ada18319248b0754247fd3a814d93af6 Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Tue, 10 Apr 2012 18:23:53 -0400
Subject: [PATCH 1/7] New tool (unfinished) to compare two void trees


From 324322271464392d894506ca8d97aa1086c1867b Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Tue, 10 Apr 2012 18:24:03 -0400
Subject: [PATCH 2/7] Added missing cmake commands


From f9bc27e88d4437e4b550adcb8c6ba955a3adba5b Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Thu, 19 Apr 2012 10:39:31 -0400
Subject: [PATCH 3/7] Made a copy of present jozov to jozov_persistent

---
 jozov_persistent.c | 550 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 550 insertions(+)
 create mode 100644 jozov_persistent.c

diff --git a/jozov_persistent.c b/jozov_persistent.c
new file mode 100644
index 0000000..f5fb4e7
--- /dev/null
+++ b/jozov_persistent.c
@@ -0,0 +1,550 @@
+#include <assert.h>
+/* jovoz.c by Mark Neyrinck */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#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 */
+} ZONE;
+
+typedef struct ZoneT {
+  int nadj; /* Number of zones on border */
+  int *adj; /* Each adjacent zone, with ... */
+  float *slv; /* Smallest Linking Volume */
+} 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<np;i++) {
+    fread(&p[i].nadj,1,sizeof(int),adj); 
+    /* The number of adjacencies per particle */
+    if (p[i].nadj > 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<np;i++) {
+    fread(&nin,1,sizeof(int),adj);
+    if (nin > 0)
+      for (k=0;k<nin;k++) {
+	fread(&j,1,sizeof(int),adj);
+	if (j < np) {
+	  /* Set both halves of the pair */
+	assert(i < j);
+          if (p[i].ncnt == p[i].nadj)
+            {
+		p[i].adj = (int *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(int));
+                 p[i].nadj++;
+            }
+          if (p[j].ncnt == p[j].nadj)
+            {
+              p[j].adj = (int *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(int));
+              p[j].nadj++;
+            }
+	  p[i].adj[p[i].ncnt] = j;
+	  p[j].adj[p[j].ncnt] = i;
+	  p[i].ncnt++; p[j].ncnt++;
+	} else {
+	  printf("%d: adj = %d\n",i,j);
+	}
+      }
+  }
+  fclose(adj);
+
+  /* Check that we got all the pairs */
+//  adj = fopen(adjfile, "r");
+//  fread(&np,1, sizeof(int),adj);
+  for (i=0;i<np;i++) {
+//    fread(&nin,1,sizeof(int),adj); /* actually nadj */
+    if (p[i].ncnt != p[i].nadj) {
+      p[i].nadj = p[i].ncnt;
+      printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
+      /*exit(0);*/
+    }
+  }
+//  fclose(adj);
+
+  /* Volumes */
+  vol = fopen(volfile, "r");
+  if (vol == NULL) {
+    printf("Unable to open volume file %s.\n\n",volfile);
+    exit(0);
+  }
+  fread(&np2,1, sizeof(int),adj);
+  if (np != np2) {
+    printf("Number of particles doesn't match! %d != %d\n",np,np2);
+    exit(0);
+  }
+  printf("%d particles\n", np);
+  FF;
+  for (i=0;i<np;i++) {
+    fread(&p[i].dens,1,sizeof(float),vol);
+    if ((p[i].dens < 1e-30) || (p[i].dens > 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<nzones;h++) {
+    zt[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
+    if (zt[h].adj == NULL) {
+      printf("Unable to allocate %d adj's of zone %d\n",zt[h].nadj,h);
+      exit(0);
+    }
+    zt[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
+    if (zt[h].slv == NULL) {
+      printf("Unable to allocate %d slv's of zone %d\n",zt[h].nadj,h);
+      exit(0);
+    }
+    zt[h].nadj = 0;
+  }
+
+  /* Find "weakest links" */
+  for (i = 0; i < np; i++) {
+    h = zonenum[jumped[i]];
+    for (j = 0; j < p[i].nadj; j++) {
+      testpart = p[i].adj[j];
+      if (h != zonenum[jumped[testpart]]) {
+	if (p[testpart].dens > 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) {
+	    zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
+	    zt[h].slv[zt[h].nadj] = p[i].dens;
+	    zt[h].nadj++;
+	  }
+	}
+      }
+    }
+  }
+  printf("Found zone adjacencies\n"); FF;
+
+  /* Free particle adjacencies */
+  for (i=0;i<np; i++) if (p[i].adj != 0) free(p[i].adj);
+
+  /* Use z instead of zt */
+  for (h=0;h<nzones;h++) {
+    /*printf("%d ",zt[h].nadj);*/
+    z[h].nadj = zt[h].nadj;
+    z[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
+    z[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
+    for (za = 0; za<zt[h].nadj; za++) {
+      z[h].adj[za] = zt[h].adj[za];
+      z[h].slv[za] = zt[h].slv[za];
+    }
+    free(zt[h].adj);
+    free(zt[h].slv);
+    z[h].np = numinh[z[h].core];
+  }
+  free(zt);
+  free(numinh);
+
+  m = (int **)malloc(nzones*sizeof(int *));
+  /* Not in the zone struct since it'll be freed up (contiguously, we hope)
+     soon */
+  nm = (int *)malloc(nzones*sizeof(int));
+  for (h=0; h<nzones; h++) {
+    m[h] = (int *)malloc(z[h].np*sizeof(int));
+    nm[h] = 0;
+    z[h].vol = 0.;
+  }
+  for (i=0; i<np; i++) {
+    h = zonenum[jumped[i]];
+    if (i == z[h].core) {
+      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;
+    }
+    z[h].vol += 1.0/(double)p[i].dens;
+    nm[h] ++;
+  }
+  free(nm);
+
+
+  zon = fopen(zonfile,"w");
+  if (zon == NULL) {
+    printf("Problem opening zonefile %s.\n\n",zonfile);
+    exit(0);
+  }
+  fwrite(&np,1,4,zon);
+  fwrite(&nzones,1,4,zon);
+  for (h=0; h<nzones; h++) {
+    fwrite(&(z[h].np),1,4,zon);
+    fwrite(m[h],z[h].np,4,zon);
+    free(m[h]);
+  }
+  free(m);
+  close(zon);
+
+  inyet = (char *)malloc(nzones*sizeof(char));
+  inyet2 = (char *)malloc(nzones*sizeof(char));
+  zonelist = (int *)malloc(nzones*sizeof(int));
+  zonelist2 = (int *)malloc(nzones*sizeof(int));
+  sorter = (double *)malloc((nzones+1)*sizeof(double));
+
+  for (h = 0; h< nzones; h++) {
+    inyet[h] = 0;
+    inyet2[h] = 0;
+  }
+
+  nhl = 0; 
+
+  maxvol = 0.;
+  minvol = BIGFLT;
+  maxdenscontrast = 0.;
+  for(i=0;i<np; i++){
+    if (p[i].dens > 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);
+
+  for (h = 0; h<nzones; h++) {
+    nhlcount = 0;
+    for (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 */
+      lowvol = BIGFLT; nl = 0; beaten = 0;
+      for (hl = 0; hl < nhl; hl++) {
+	h2 = zonelist[hl];
+	if (inyet[h2] == 1) { /* If it's not already identified as 
+				 an interior zone, with inyet=2 */
+	  interior = 1;
+	  for (za = 0; za < z[h2].nadj; za ++) {
+	    if (inyet[z[h2].adj[za]] == 0) {
+	      interior = 0;
+	      if (z[h2].slv[za] == lowvol) {
+		link[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];
+		link[0] = z[h2].adj[za];
+		nl = 1;
+	      }
+	    }
+	  }
+	  if (interior == 1) inyet[h2] = 2; /* No bordering exter. zones */
+	}
+      }
+
+      if (nl == 0) {
+	beaten = 1;
+	z[h].leak = maxvol;
+	continue;
+      }
+	
+      if (lowvol > voltol) {
+	beaten = 1;
+	z[h].leak = lowvol;
+	continue;
+      }
+
+      for (l=0; l < nl; l++)
+	if (p[z[link[l]].core].dens < p[z[h].core].dens)
+	  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;
+	    for (hl = 0; (hl < nhl2) && (beaten == 0); hl++) {
+	      h2 = zonelist2[hl];
+	      if (inyet2[h2] == 1) {
+		interior = 1; /* Guilty until proven innocent */
+		for (za = 0; za < z[h2].nadj; za ++) {
+		  link2 = z[h2].adj[za];
+		  if ((inyet[link2]+inyet2[link2]) == 0) {
+		    interior = 0;
+		    if (z[h2].slv[za] <= lowvol) {
+		      if (p[z[link2].core].dens < p[z[h].core].dens) {
+			beaten = 1;
+			break;
+		      }
+		      zonelist2[nhl2] = link2;
+		      inyet2[link2] = 1;
+		      nhl2++;
+		      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<nhl; q++) {
+      z[h].voljoin += z[zonelist[q]].vol;
+    }
+
+    z[h].nhl = nhl;
+
+    fwrite(&nhl,1,4,zon2);
+    fwrite(zonelist,nhl,4,zon2);
+  }
+  fclose(zon2);
+
+  printf("Maxdenscontrast = %f.\n",maxdenscontrast);
+
+  /* Assign sorter by probability (could use volume instead) */
+  for (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 = fopen(txtfile,"w");
+  fprintf(txt,"%d particles, %d voids.\n", np, nzones);
+  fprintf(txt,"Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
+  for (h=0; h<nzones; h++) {
+    i = iord[h];
+    prob = exp(-5.12*(z[i].denscontrast-1.) - 0.8*pow(z[i].denscontrast-1.,2.8));
+    fprintf(txt,"%d %d %d %e %e %d %d %e %d %f %6.2e\n",
+	    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);
+
+  } /* h+1 to start from 1, not zero */
+  fclose(txt);
+
+  return(0);
+}

From 3d81d499061c757ace9e92f2e8849b4db9a42a3b Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Thu, 19 Apr 2012 15:47:53 -0400
Subject: [PATCH 4/7] Added comments. Cleaned up the code

---
 jozov_persistent.c | 87 ++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 73 insertions(+), 14 deletions(-)

diff --git a/jozov_persistent.c b/jozov_persistent.c
index f5fb4e7..9a76e8a 100644
--- a/jozov_persistent.c
+++ b/jozov_persistent.c
@@ -27,12 +27,14 @@ 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 */
+  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);
@@ -114,11 +116,11 @@ int main(int argc,char **argv) {
 	fread(&j,1,sizeof(int),adj);
 	if (j < np) {
 	  /* Set both halves of the pair */
-	assert(i < j);
+	  assert(i < j);
           if (p[i].ncnt == p[i].nadj)
             {
-		p[i].adj = (int *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(int));
-                 p[i].nadj++;
+	      p[i].adj = (int *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(int));
+	      p[i].nadj++;
             }
           if (p[j].ncnt == p[j].nadj)
             {
@@ -136,17 +138,17 @@ int main(int argc,char **argv) {
   fclose(adj);
 
   /* Check that we got all the pairs */
-//  adj = fopen(adjfile, "r");
-//  fread(&np,1, sizeof(int),adj);
+  //  adj = fopen(adjfile, "r");
+  //  fread(&np,1, sizeof(int),adj);
   for (i=0;i<np;i++) {
-//    fread(&nin,1,sizeof(int),adj); /* actually nadj */
+    //    fread(&nin,1,sizeof(int),adj); /* actually nadj */
     if (p[i].ncnt != p[i].nadj) {
       p[i].nadj = p[i].ncnt;
       printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
       /*exit(0);*/
     }
   }
-//  fclose(adj);
+  //  fclose(adj);
 
   /* Volumes */
   vol = fopen(volfile, "r");
@@ -249,13 +251,17 @@ int main(int argc,char **argv) {
       exit(0);
     }
     zt[h].nadj = 0;
+    zt[h].smallest_saddle = BIGFLT;
   }
 
   /* Find "weakest links" */
   for (i = 0; i < np; i++) {
     h = zonenum[jumped[i]];
+
     for (j = 0; j < p[i].nadj; j++) {
+
       testpart = p[i].adj[j];
+
       if (h != zonenum[jumped[testpart]]) {
 	if (p[testpart].dens > p[i].dens) {
 	  /* there could be a weakest link through testpart */
@@ -282,8 +288,13 @@ int main(int argc,char **argv) {
 	      }
 	    }
 	  if (already == 0) {
-	    zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
-	    zt[h].slv[zt[h].nadj] = p[i].dens;
+	    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++;
 	  }
 	}
@@ -308,6 +319,7 @@ int main(int argc,char **argv) {
     free(zt[h].adj);
     free(zt[h].slv);
     z[h].np = numinh[z[h].core];
+    z[h].smallest_saddle = zt[h].smallest_saddle;
   }
   free(zt);
   free(numinh);
@@ -379,6 +391,12 @@ int main(int argc,char **argv) {
   }
   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<nzones; h++) {
     nhlcount = 0;
     for (hl = 0; hl < nhl; hl++)
@@ -391,15 +409,25 @@ int main(int argc,char **argv) {
     do {
       /* Find the lowest-volume (highest-density) adjacency */
       lowvol = BIGFLT; nl = 0; beaten = 0;
+
+      /* Look at all zones that have been investigated in all previous iterations. */
       for (hl = 0; hl < nhl; hl++) {
+
 	h2 = zonelist[hl];
+	/* Is this zone been processed and deemed interior in this round ? */
 	if (inyet[h2] == 1) { /* If it's not already identified as 
 				 an interior zone, with inyet=2 */
+	  /* Yes. Let's have a look at the adjacents... */
 	  interior = 1;
+
 	  for (za = 0; za < z[h2].nadj; za ++) {
+	    
+	    /* This zones has not been deemed as interior */
 	    if (inyet[z[h2].adj[za]] == 0) {
 	      interior = 0;
+	      /* Just at the threshold ? */
 	      if (z[h2].slv[za] == lowvol) {
+		/* Yes. Queue it for adding the zones in bunch */
 		link[nl] = z[h2].adj[za];
 		nl ++;
 		if (nl == NLINKS) {
@@ -407,14 +435,17 @@ int main(int argc,char **argv) {
 		  exit(0);
 		}
 	      }
+	      /* Just below the threshold ? */
 	      if (z[h2].slv[za] < lowvol) {
+		/* Set the new threshold and reset the list because we may have beater */
 		lowvol = z[h2].slv[za];
 		link[0] = z[h2].adj[za];
 		nl = 1;
 	      }
 	    }
 	  }
-	  if (interior == 1) inyet[h2] = 2; /* No bordering exter. zones */
+	  if (interior == 1)
+	    inyet[h2] = 2; /* No bordering exter. zones */
 	}
       }
 
@@ -430,47 +461,74 @@ int main(int argc,char **argv) {
 	continue;
       }
 
+
+      /* Check if among the found adjacent candidates there is a beater */
       for (l=0; l < nl; l++)
-	if (p[z[link[l]].core].dens < p[z[h].core].dens)
+	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) {
-		      if (p[z[link2].core].dens < p[z[h].core].dens) {
+		      /* 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;
+		if (interior == 1) 
+		  inyet2[h2] = 2;
 	      }
 	    }
 	  }
+
 	}
       }
       for (hl = 0; hl < nhl2; hl++)
@@ -494,7 +552,8 @@ int main(int argc,char **argv) {
     } 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.;
+    if (z[h].denscontrast < 1.) 
+      z[h].denscontrast = 1.;
     
     /* find biggest denscontrast */
     if (z[h].denscontrast > maxdenscontrast) {

From e9a6b3afa04ccd2bdcb18f7dab0e10af3f29b1f1 Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Thu, 19 Apr 2012 16:24:19 -0400
Subject: [PATCH 5/7] Rework the void tree to extract zone list in the other
 tree


From 9246e04a6457f145e8c9df1512af4ffc2fbc6a8e Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Thu, 19 Apr 2012 16:35:32 -0400
Subject: [PATCH 6/7] Do the loading of the segmented list


From bbee74fff05829ebb6c614cae2ee74e410fef80a Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Thu, 19 Apr 2012 16:35:46 -0400
Subject: [PATCH 7/7] Missing increment...