re-added some essential routines; updated example sceipts to reflect new layout

c_tools/jozov2/CMakeLists.txt

@@ -0,0 +1,26 @@
+function(omp_add_flags TARGET LANG)
+   if(OPENMP_FOUND)
+     if(NOT LANG)
+       get_target_property(LANG ${TARGET} LINKER_LANGUAGE)
+     endif()
+     if(LANG MATCHES "C(XX)?")
+       set_property(TARGET ${TARGET} APPEND
+                    PROPERTY COMPILE_FLAGS ${OpenMP_${LANG}_FLAGS})
+       set_property(TARGET ${TARGET} APPEND
+                    PROPERTY LINK_FLAGS ${OpenMP_${LANG}_FLAGS})
+     else()
+       message(WARNING "omp_add_flags: target '${TARGET}'"
+                      " link language '${LANG}' is not supported")
+    endif()
+  endif()
+endfunction()
+
+add_executable(jozov2 jozov2.cpp jozov2_io.cpp jozov2_zones.cpp jozov2_watershed.cpp findrtop.c)
+
+if (ENABLE_OPENMP)
+  Find_Package(OpenMP)
+  omp_add_flags(jozov2 CXX)
+  add_definitions(-DOPENMP)
+ENDIF(ENABLE_OPENMP)
+
+install(TARGETS jozov2 DESTINATION ${VIDE_BIN})

c_tools/jozov2/findrtop.c

@@ -0,0 +1,81 @@
+#include <math.h>
+#include <stdlib.h>
+
+/*------------------------------------------------------------------------------
+  Find nb elements of Real array a having the largest value.
+  Returns index iord of these elements, ordered so iord[0] corresponds
+  to element a[iord[0]] having the largest value.
+  If nb > na, then the last nb-na elements of iord are undefined.
+  Elements of a that are equal are left in their original order.
+
+  Courtesy Andrew Hamilton.
+*/
+void findrtop(double *a, int na, int *iord, int nb)
+{
+#undef	ORDER
+#define	ORDER(ia, ja)	a[ia] > a[ja] || (a[ia] == a[ja] && ia < ja)
+
+  int i, ia, ib, it, n;
+  
+  n = (na <= nb)? na : nb;
+  if (n <= 0) return;
+  
+  /* heap first n elements, so smallest element is at top of heap */
+  for (ib = (n >> 1); ib < n; ib++) {
+    iord[ib] = ib;
+  }
+  for (ia = (n >> 1) - 1; ia >= 0; ia--) {
+    i = ia;
+    for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) {
+      if (ib+1 < n) {
+	if (ORDER(iord[ib], iord[ib+1])) ib++;
+      }
+      if (ORDER(ia, iord[ib])) {
+	iord[i] = iord[ib];
+	i = ib;
+      } else {
+	break;
+      }
+    }
+    iord[i] = ia;
+  }
+  
+  /* now compare rest of elements of array to heap */
+  for (ia = n; ia < na; ia++) {
+    /* if new element is greater than smallest, sift it into heap */
+    i = 0;
+    if (ORDER(ia, iord[i])) {
+      for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) {
+	if (ib+1 < n) {
+	  if (ORDER(iord[ib], iord[ib+1])) ib++;
+	}
+	if (ORDER(ia, iord[ib])) {
+	  iord[i] = iord[ib];
+	  i = ib;
+	} else {
+	  break;
+	}
+      }
+      iord[i] = ia;
+    }
+  }
+  
+  /* unheap iord so largest element is at top */
+  for (ia = n - 1; ia > 0; ia--) {
+    it = iord[ia];
+    i = 0;
+    iord[ia] = iord[i];
+    for (ib = (i << 1) + 1; ib < ia; ib = (i << 1) + 1) {
+      if (ib+1 < ia) {
+	if (ORDER(iord[ib], iord[ib+1])) ib++;
+      }
+      if (ORDER(it, iord[ib])) {
+	iord[i] = iord[ib];
+	i = ib;
+      } else {
+	break;
+      }
+    }
+    iord[i] = it;
+  }
+}

c_tools/jozov2/jozov2.cpp

@@ -0,0 +1,174 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#include <boost/format.hpp>
+#include <exception>
+#include <iostream>
+#include <fstream>
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <assert.h>
+/* jovoz.c by Mark Neyrinck */
+#include "jozov2.hpp"
+#include "zobov.hpp"
+
+using namespace std;
+using boost::format;
+
+
+
+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);
+
+  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(p, np, z, nzones, maxvol, voltol);
+  writeVoidFile(zonfile2, z, nzones);
+
+  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 = new int[nzones];
+
+  findrtop(sorter, nzones, iord, nzones);
+  delete[] sorter;
+
+  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<nzones; h++)
+    {
+      int i = iord[h];
+      prob = exp(-5.12*(z[i].denscontrast-1.) - 0.8*pow(z[i].denscontrast-1.,2.8));
+      if (z[i].np == 1) 
+	continue;
+      
+      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 << endl;
+
+    } /* h+1 to start from 1, not zero */
+  txt.close();
+
+  delete[] iord;
+  delete[] z;
+  delete[] p;
+
+  cout << "Done!" << endl;
+  return(0);
+}

c_tools/jozov2/jozov2.hpp

@@ -0,0 +1,65 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#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

c_tools/jozov2/jozov2_io.cpp

@@ -0,0 +1,199 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_io.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#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();
+                  }
+                else               
+                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();
+                  }
+                else{
+                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;
+    }
+
+  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));
+    }
+}

c_tools/jozov2/jozov2_watershed.cpp

@@ -0,0 +1,294 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_watershed.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#ifdef OPENMP
+#include <omp.h>
+#endif
+
+#include <queue>
+#include <set>
+#include <iostream>
+#include <fstream>
+#include <boost/format.hpp>
+#include <string>
+#include "jozov2.hpp"
+#include "zobov.hpp"
+
+using namespace std;
+using boost::format;
+
+struct ZoneDensityPair
+{
+  int h;
+  double density;
+  double core;
+
+  bool operator<(const ZoneDensityPair& p2) const
+  {
+    return (density > p2.density) || (density==p2.density && core > p2.core);
+  }
+};
+
+typedef priority_queue<ZoneDensityPair> ZoneQueue;
+
+static void build_process_queue(ZoneQueue& q, ZONE *z, PARTICLE *p, 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];
+     zdp.core = p[z[zdp.h].core].dens;
+     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)
+{
+  /* Text output file */
+
+#pragma omp parallel
+  {
+    char *inyet, *inyet2;
+    int *zonelist, *zonelist2;
+    int nhl;
+    int prev_ii = -1;
+
+    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;
+#pragma omp for schedule(dynamic,1)
+    for (int h = 0; h < numZones; h++)
+      {
+        int nhlcount = 0;
+        float previous_lowvol = BIGFLT, lowvol, z_cur_core_dens;
+        bool beaten;
+        priority_queue<ZoneDensityPair> 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;
+        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, p, inyet, h);
+	
+        do {
+          /* Find the lowest-volume (highest-density) adjacency */
+          beaten = false;
+          
+          if (to_process.empty())
+            {
+              beaten = true;
+              z[h].leak = maxvol;
+              save_npjoin = z[h].npjoin;
+              save_nhl = nhl;
+              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())
+	    {
+              save_npjoin = z[h].npjoin;
+              save_nhl = nhl;
+	      beaten = true;
+	      z[h].leak = maxvol;
+	      continue;
+	    }
+
+          /* 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;
+              z[h].leak = lowvol;
+              continue;
+            }
+          
+          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;
+	  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);
+              z[h].npjoin += z[new_h].np;
+              inyet[new_h] = 1;
+              if (inyet2[new_h] != 2)
+                build_process_queue(to_process, z, p, inyet, new_h);
+              nhl++;
+            }
+          }
+	  
+          for (int hl = 0; hl < nhl2; hl++)
+            inyet2[zonelist2[hl]] = 0;
+          if (nhl/10000 > nhlcount) {
+            if (nhlcount == 0)
+              (cout << format("Zone %d: %d") % h % nhl).flush(); 
+            else
+              (cout << format(" %d [%d]") % nhl % to_process.size()).flush();
+            nhlcount = nhl/10000;
+          }
+	}
+	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.) 
+	  z[h].denscontrast = 1.;
+        
+        
+        /* 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].npjoin = save_npjoin;
+        z[h].voljoin = 0.;
+        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 = save_nhl;
+      }
+    delete[] zonelist;
+    delete[] zonelist2;
+    delete[] inyet;
+    delete[] inyet2;
+    
+  }
+  
+  double maxdenscontrast = 0;
+#pragma omp parallel shared(maxdenscontrast)
+  {
+    double maxdenscontrast_local = 0;
+
+#pragma omp for schedule(static)
+    for (int h = 0; h < numZones; h++)
+      {
+        /* find biggest denscontrast */
+        if (z[h].denscontrast > maxdenscontrast_local) {
+          maxdenscontrast_local = (double)z[h].denscontrast;
+        }
+      }
+
+#pragma omp critical
+    {
+      if (maxdenscontrast_local > maxdenscontrast)
+        maxdenscontrast = maxdenscontrast_local;
+    }
+  }
+
+  cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl;
+}

c_tools/jozov2/jozov2_zones.cpp

@@ -0,0 +1,243 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_zones.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#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();
+  
+  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) 
+      loc_NumZones++;
+
+  numZones = loc_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++;
+    }
+
+  size_t nadjAlloced = 0;
+  cout << "Total numZones = " << numZones << endl;
+  try
+    {
+      for (int h = 0; h < numZones; h++) {
+        if (zt[h].nadj > 0) {
+          zt[h].adj = new pid_t[zt[h].nadj];
+          zt[h].slv = new float[zt[h].nadj];
+          nadjAlloced += zt[h].nadj;
+          zt[h].nadj = 0;
+        } else {
+          zt[h].adj = 0;
+          zt[h].slv = 0;
+        }
+      }
+    }
+  catch(const std::bad_alloc& a)
+    {
+      cout << "Could not allocate memory for zone adjacencies (nadj so far: " << nadjAlloced << ", memory needed: " << (nadjAlloced*(sizeof(pid_t)+sizeof(float))) << ")" << 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;
+    if (zt[h].nadj > 0) {
+      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;
+    } else {
+      z[h].adj = 0;
+      z[h].slv = 0;
+    }
+    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;
+        z[h].vol = 0;
+        z[h].np = z[h].npjoin = z[h].nadj = z[h].nhl = 0;
+        z[h].leak = 0;
+        z[h].adj = 0;
+        z[h].slv = 0;
+        z[h].denscontrast = 0;
+        z[h].vol = z[h].voljoin = 0;
+        zonenum[i] = h;
+        h++;
+      } else {
+        zonenum[i] = -1;
+      }
+    }
+  for (size_t z = 0; z < nzones; z++) {
+    zt[z].nadj = 0;
+    zt[z].adj = 0;
+    zt[z].slv = 0;
+  }
+
+  buildZoneAdjacencies(p, np, z, zt,
+                       h, jumped, zonenum, numinh);
+
+  delete[] zt;
+  delete[] numinh;
+
+  for (pid_t i=0; i < np; i++) {
+    int h = zonenum[jumped[i]];
+    z[h].vol += 1.0/(double)p[i].dens;
+    z[h].numzones = 0;
+    z[h].zonelist = 0;
+  }
+}

c_tools/jozov2/zobov.hpp

@@ -0,0 +1,52 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/zobov.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#ifndef __ZOBOV_HPP
+#define __ZOBOV_HPP
+
+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 */
+
+  int *zonelist; /* Zones bound to the void. */
+  int numzones; /* Number of zones bound. */
+} ZONE;
+
+typedef struct ZoneT {
+  int nadj; /* Number of zones on border */
+  int *adj; /* Each adjacent zone, with ... */
+  float *slv; /* Smallest Linking Volume */
+} ZONET;
+
+#endif

c_tools/pruning/CMakeLists.txt

@@ -0,0 +1,8 @@
+include_directories(${CMAKE_CURRENT_BINARY_DIR})
+
+SET(pruneVoids_SRCS pruneVoids.cpp)
+add_genopt(pruneVoids_SRCS pruneVoids.ggo pruneVoids_conf STRUCTNAME pruneVoids_info)
+add_executable(pruneVoids ${pruneVoids_SRCS})
+target_link_libraries(pruneVoids ${ZOB_LIBS})
+
+install(TARGETS pruneVoids DESTINATION ${VIDE_BIN})

c_tools/pruning/pruneVoids.ggo

@@ -0,0 +1,44 @@
+package "removeEdgeVoids"
+version "alpha"
+
+
+option "configFile" - "Configuration filename" string optional
+
+option "partFile" - "Particle file from generateFromCatalog" string required
+
+option "extraInfo" - "Extra particle file from generateFromCatalog" string required
+
+option "voidDesc" - "Void list from ZOBOV" string required
+
+option "void2Zone" - "Zone file from ZOBOV" string required
+
+option "partVol" - "Particle volume file from ZOBOV" string required
+
+option "partAdj" - "Adjacency file from ZOBOV" string required
+
+option "zone2Part" - "Particle file from ZOBOV" string required
+option "mockIndex" - "Beginning index of mock particles" int required
+
+option "numVoids" - "Number of voids" int required
+
+option "isObservation" - "We are working with observational data" flag off
+
+option "useComoving" - "Void positions are in comoving coordinates" flag off
+
+option "omegaM" - "Omega_M for redshift convertion" double optional default="0.27"
+
+option "zMin" - "Minimum redshift of sample" double optional default="0.0"
+option "zMax" - "Maximum redshift of sample" double optional default="10.0"
+
+option "rMin" - "Minimum allowable void radius" double optional default="0.0"
+
+option "outputDir" - "Directory to place outputs" string required
+option "sampleName" - "unique string to assign to outputs" string required
+
+option "periodic" - "Set of edges which are periodic" string optional default="xy"
+
+option "tolerance" - "Fraction of void width to consider edge" double optional default="1.0"
+
+option "centralRadFrac" - "Fraction of void radii to consider central region" double optional default="4"
+
+option "maxCentralDen" - "Maximum central density to accept as a void" double optional default="0.2" 

python_tools/vide/backend/cosmologyTools.py

@@ -0,0 +1,112 @@
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/apTools/chi2/cosmologyTools.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+# 
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+# a suite of functions to compute expansion rates, angular diameter 
+# distances, and expected void stretching
+
+import numpy as np
+import scipy.integrate as integrate
+import healpy as healpy
+import os
+from vide.backend import *
+
+__all__=['expansion', 'angularDiameter', 'aveExpansion', 'getSurveyProps']
+
+# returns 1/E(z) for the given cosmology
+def expansion(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
+  ez = Om * (1+z)**3 + (Ot-Om)# * (1+z)**(3.+3*wz)
+  #ez = Om * (1+z)**3 + (Ot-Om)# * integrade.quad(eosDE, 0.0, z, args=(w0,wa))[0]
+  ez = 1./np.sqrt(ez)
+  return ez
+
+# returns DE value at redshift z
+def eosDE(z, w0 = -1.0, wa = 0.0):
+  return w0 + wa*z/(1+z)
+  
+# returns D_A(z) for the given cosmology
+def angularDiameter(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
+  da = integrate.quad(expansion, 0.0, z, args=(Om, Ot, w0, wa))[0]
+  return da
+  
+# -----------------------------------------------------------------------------
+# returns average expected expansion for a given redshift range
+def aveExpansion(zStart, zEnd, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
+  if zStart == 0.0: zStart = 1.e-6
+  ave = integrate.quad(expansion, zStart, zEnd, args=(Om, Ot, w0, wa))[0]
+  ave = (zEnd-zStart)/ave
+  return ave
+
+# -----------------------------------------------------------------------------
+# returns the volume and galaxy density for a given redshit slice
+def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectionFuncFile=None, useComoving=False):
+
+  LIGHT_SPEED = 299792.458
+
+  mask = healpy.read_map(maskFile)
+  area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi
+
+  if useComoving:
+    zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=0.27)
+    zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=0.27)
+    selFunMin = LIGHT_SPEED/100.*angularDiameter(selFunMin, Om=0.27)
+    selFunMax = LIGHT_SPEED/100.*angularDiameter(selFunMax, Om=0.27)
+  else:
+    zmin = zmin * 3000
+    zmax = zmax * 3000
+    selFunMin *= 3000
+    selFunMax *= 3000
+  volume = area * (zmax**3  - zmin**3) / 3
+
+  if selectionFuncFile == None:
+    nbar = 1.0
+
+  elif not os.access(selectionFuncFile, os.F_OK):
+    print("   Warning, selection function file %s not found, using default of uniform selection." % selectionFuncFile)
+    nbar = 1.0
+
+  else:
+    selfunc = np.genfromtxt(selectionFuncFile)
+    selfunc = np.array(selfunc)
+    selfunc[:,0] = selfunc[:,0]/100.
+    selfuncUnity = selfunc
+    selfuncUnity[:,1] = 1.0
+    selfuncMin = selfunc[0,0]
+    selfuncMax = selfunc[-1,0]
+    selfuncDx = selfunc[1,0] - selfunc[0,0]
+    selfuncN = np.size(selfunc[:,0])
+  
+    selFunMin = max(selFunMin, selfuncMin)
+    selFunMax = min(selFunMax, selfuncMax)
+  
+  
+    def f(z): return selfunc[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 
+    def fTotal(z): return selfuncUnity[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 
+  
+    zrange = np.linspace(selFunMin, selFunMax)
+  
+    nbar = scipy.integrate.quad(f, selFunMin, selFunMax)
+    nbar = nbar[0]
+  
+    ntotal = scipy.integrate.quad(fTotal, 0.0, max(selfuncUnity[:,0]))
+    ntotal = ntotal[0]
+  
+    nbar = ntotal / area / nbar
+
+
+  return (volume, nbar)