Merge bergeron:ns/Software/void_identification

.gitignore

@@ -1,8 +1,4 @@
 *~
 *.o
 *.prog
-voz1b1
-jozov
-vozinit
-voztie
 .mydepend

c_tools/analysis/voidOverlap.cpp

@@ -177,24 +177,70 @@ int main(int argc, char **argv) {
       newMatch.commonVolProg = 0;
       newMatch.dist = rdist; 
 
-      if (rdist/catalog1.voids[iVoid1].radius > 1.5) continue;
+      //if (rdist > 1.5) continue;
+      if (rdist/catalog1.voids[iVoid1].radius > 2.0) continue;
 
-      if (catalog1.voids[iVoid1].matches.size() < MAX_MATCHES) {
+      // see if center is contained in void
+      match = false;
+      voidID1 = catalog1.voids[iVoid1].voidID;
+      for (iZ1 = 0; iZ1 < catalog1.void2Zones[voidID1].numZones; iZ1++) {
+        zoneID1 = catalog1.void2Zones[voidID1].zoneIDs[iZ1];
+        for (p1 = 0; p1 < catalog1.zones2Parts[zoneID1].numPart; p1++) {
+          partID1 = catalog1.zones2Parts[zoneID1].partIDs[p1];     
+
+          dist[0] = fabs(catalog1.part[partID1].x - 
+                         catalog2.voids[iVoid2].barycenter[0]);
+          dist[1] = fabs(catalog1.part[partID1].y - 
+                         catalog2.voids[iVoid2].barycenter[1]);
+          dist[2] = fabs(catalog1.part[partID1].z - 
+                         catalog2.voids[iVoid2].barycenter[2]);
+
+          if (periodicX) dist[0] = fmin(dist[0], 1.0 - dist[0]);
+          if (periodicY) dist[1] = fmin(dist[1], 1.0 - dist[1]);
+          if (periodicZ) dist[2] = fmin(dist[2], 1.0 - dist[2]);
+
+          rdist = sqrt(dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2]);
+
+          r1 = pow(3./4./M_PI*catalog1.part[partID1].volume / 
+                   catalog1.numPartTot, 1./3.);
+          if (rdist <= matchDist*4*r1) {
+            match = true;
+            break;
+          }
+        }
+      }
+       
+      if (!match) continue; 
+      //if (rdist/catalog1.voids[iVoid1].radius > 1.) continue;
+
+      //if (catalog1.voids[iVoid1].matches.size() < MAX_MATCHES) {
         catalog1.voids[iVoid1].matches.push_back(newMatch);
-      } else {
+      //} else {
-        // find the farthest match
+      //  // find the farthest match
-        float farthestMatchDist = 0;
+      //  float farthestMatchDist = 0;
-        int   farthestMatchID = 0;
+      //  int   farthestMatchID = 0;
-        for (iMatch = 0; iMatch < MAX_MATCHES; iMatch++) {
+      //  for (iMatch = 0; iMatch < MAX_MATCHES; iMatch++) {
-          if (catalog1.voids[iVoid1].matches[iMatch].dist > farthestMatchDist){
+      //    if (catalog1.voids[iVoid1].matches[iMatch].dist > farthestMatchDist){
-            farthestMatchDist = catalog1.voids[iVoid1].matches[iMatch].dist;
+      //      farthestMatchDist = catalog1.voids[iVoid1].matches[iMatch].dist;
-            farthestMatchID   = iMatch;
+       //     farthestMatchID   = iMatch;
+       //   }           
+       // } 
+       // if (rdist < farthestMatchDist)
+       //   catalog1.voids[iVoid1].matches[farthestMatchID] = newMatch;
+      //}
     }
   } 
-        if (rdist < farthestMatchDist)
+
-          catalog1.voids[iVoid1].matches[farthestMatchID] = newMatch;
+  // pick the closest matches to speed up computation
-      }
+  for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
+    int numPotential = catalog1.voids[iVoid1].matches.size();
+    catalog1.voids[iVoid1].numMatches = numPotential;
+    
+    if (numPotential > MAX_MATCHES) {
+      sortMatches(catalog1.voids[iVoid1].matches, catalog2, 1);
     }
+
+    catalog1.voids[iVoid1].matches.resize(MAX_MATCHES); 
   }
 
   printf(" Determining overlap...\n");
@@ -267,9 +313,13 @@ int main(int argc, char **argv) {
       int matchID = catalog1.voids[iVoid1].matches[iMatch].matchID;
       catalog1.voids[iVoid1].matches[iMatch].merit = 
         catalog1.voids[iVoid1].matches[iMatch].commonVolOrig * 
-        catalog1.voids[iVoid1].matches[iMatch].commonVolProg / 
+        catalog1.voids[iVoid1].matches[iMatch].commonVolOrig / 
         catalog1.voids[iVoid1].vol /
-        catalog2.voids[matchID].vol;
+        catalog2.voids[iVoid1].vol;
+        //catalog1.voids[iVoid1].matches[iMatch].commonVolOrig * 
+        //catalog1.voids[iVoid1].matches[iMatch].commonVolProg / 
+        //catalog1.voids[iVoid1].vol /
+        //catalog2.voids[matchID].vol;
         //pow(catalog1.voids[iVoid1].matches[iMatch].commonVol,2) / 
         //catalog1.voids[iVoid1].vol /
         //catalog2.voids[matchID].vol;
@@ -296,7 +346,6 @@ int main(int argc, char **argv) {
                      catalog1.voids[iVoid1].vol;
       if (commonVolRatio > 0.2) catalog1.voids[iVoid1].numBigMatches++;
     } 
-    catalog1.voids[iVoid1].numMatches = catalog1.voids[iVoid1].matches.size();
   }
 
   // output summary
@@ -305,7 +354,7 @@ int main(int argc, char **argv) {
   filename = string(args.outfile_arg);
   filename = filename.append("summary.out");
   fp = fopen(filename.c_str(), "w");
-  fprintf(fp, "# void ID, radius, radius ratio, common volume ratio (to original), common volume ratio (to match), relative dist, num matches, num significant matches, match ID, merit, ellipticity ratio\n");
+  fprintf(fp, "# void ID, radius, radius ratio, common volume ratio (to original), common volume ratio (to match), relative dist, num matches, num significant matches, match ID, merit, ellipticity ratio, density contrast\n");
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
     int voidID = catalog1.voids[iVoid1].voidID;
     if (catalog1.voids[iVoid1].numMatches > 0) {
@@ -322,7 +371,7 @@ int main(int argc, char **argv) {
       rdist = catalog1.voids[iVoid1].matches[0].dist;
       rdist /= catalog1.voids[iVoid1].radius;
 
-      fprintf(fp, "%d %.4f %.4f %e %e %.4f %d %d %d %e %e\n", voidID, 
+      fprintf(fp, "%d %.4f %.4f %e %e %.4f %d %d %d %e %e %e\n", voidID, 
       //fprintf(fp, "%d %.4f %.4f %.4f %.4f %.4f %d %d %d\n", voidID, 
                                    catalog1.voids[iVoid1].radiusMpc,
                                    rRatio, 
@@ -333,10 +382,11 @@ int main(int argc, char **argv) {
                                    catalog1.voids[iVoid1].numBigMatches,
                                    catalog2.voids[iVoid2].voidID,
                                    catalog1.voids[iVoid1].matches[0].merit,
-                                   ellipRatio);
+                                   ellipRatio,
+                                   catalog1.voids[iVoid1].densCon);
 
     } else {
-      fprintf(fp, "%d %.4f 0.0 0.0 0.0 0.0 0.0 0 0 0 0.0\n", voidID, 
+      fprintf(fp, "%d %.4f 0.0 0.0 0.0 0.0 0.0 0 0 0 0.0 0.0\n", voidID, 
                                            catalog1.voids[iVoid1].radiusMpc);
     }
   } // end printing
@@ -347,25 +397,26 @@ int main(int argc, char **argv) {
   filename = string(args.outfile_arg);
   filename = filename.append("detail.out");
   fp = fopen(filename.c_str(), "w");
-  int MAX_OUT = 10;
+  fprintf(fp, "# match ID, merit, relative dist, relative radius\n");
-  fprintf(fp, "# void ID, match common vol\n");
   for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
     int voidID = catalog1.voids[iVoid1].voidID;
-    fprintf(fp,"%d: ", voidID);
+    fprintf(fp,"#%d (%.4f Mpc/h):\n", voidID, catalog1.voids[iVoid1].radiusMpc);
-    for (iMatch = 0; iMatch < MAX_OUT; iMatch++) {
+    for (iMatch = 0; iMatch < catalog1.voids[iVoid1].matches.size(); iMatch++) {
-      if (iMatch < catalog1.voids[iVoid1].matches.size()) {
       commonVolRatio = catalog1.voids[iVoid1].matches[iMatch].commonVolOrig / 
                        //catalog1.voids[iVoid1].numPart;
                        catalog1.voids[iVoid1].vol;
 
-        fprintf(fp, "%.3f %.2f ", catalog1.voids[iVoid1].matches[iMatch].dist, commonVolRatio);
+      int matchID = catalog1.voids[iVoid1].matches[iMatch].matchID;
-        //fprintf(fp, "%.2f ", commonVolRatio);
+      fprintf(fp, "%d %e %.4f %.4f\n", matchID, 
-      } else {
+                                   catalog1.voids[iVoid1].matches[iMatch].merit,
-        fprintf(fp, "0.00 ");
+                                   catalog1.voids[iVoid1].matches[iMatch].dist/
-      }
+                                     catalog1.voids[iVoid1].radius,
+                                   catalog2.voids[matchID].radius/
+                                     catalog1.voids[iVoid1].radius);
     }
 
-    fprintf(fp, "\n");
+    fprintf(fp, "\n ");
+    fprintf(fp, "\n ");
   } // end printing detail
   fclose(fp);
 

c_tools/stacking/pruneVoids.cpp

@@ -578,6 +578,8 @@ int main(int argc, char **argv) {
 
     // compute eigenvalues and vectors for orientation and shape
     double inertia[9];
+    for (int i = 0; i < 9; i++) inertia[i] = 0.;
+
     for (int p = 0; p < voids[iVoid].numPart; p++) {
       dist[0] = voidPart[p].x - voids[iVoid].barycenter[0];
       dist[1] = voidPart[p].y - voids[iVoid].barycenter[1];
@@ -614,9 +616,11 @@ int main(int argc, char **argv) {
     float c = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,0) +
                         gsl_vector_get(voids[iVoid].eval,1) - 
                         gsl_vector_get(voids[iVoid].eval,2)));
-    float ca = c/a;
+    float ca;
     float cb = c/b;
-    voids[iVoid].ellip = 1.0 - c/a;
+    if (a < c)  ca = a/c;
+    if (a >= c) ca = c/a;
+    voids[iVoid].ellip = fabs(1.0 - ca);
 
   } // iVoid
 
@@ -938,7 +942,7 @@ void outputVoids(string outputDir, string sampleName, string prefix,
              outVoid.radius,
              outVoid.voidID);
 
-     fprintf(fpShapes, "%d %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n", 
+     fprintf(fpShapes, "%d %.6f %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e\n", 
              outVoid.voidID,
              outVoid.ellip,
              gsl_vector_get(outVoid.eval, 0),

c_tools/zobov2/voz1b1/voz.h

@@ -0,0 +1,29 @@
+#ifndef __VOZ_H
+#define __VOZ_H
+
+#define MAXVERVER 100000
+#define NGUARD 84 /*Actually, the number of SPACES between guard points
+##define NGUARD 42 /*Actually, the number of SPACES between guard points
+		    in each dim */
+
+typedef int pid_t;
+
+typedef struct Partadj {
+  pid_t nadj;
+  pid_t *adj;
+} PARTADJ;
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs);
+int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol);
+int posread(char *posfile, float ***p, float fact);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif

c_tools/zobov2/voz1b1/voz1b1.cpp

@@ -0,0 +1,395 @@
+#include <boost/format.hpp>
+#include <iostream>
+#include <fstream>
+#include <CosmoTool/miniargs.hpp>
+#include <cmath>
+#include "libqhull/qhull_a.h"
+#include "voz.h"
+#include "voz_io.hpp"
+
+using CosmoTool::MINIARG_STRING;
+using CosmoTool::MINIARG_DOUBLE;
+using CosmoTool::MINIARG_INT;
+using CosmoTool::MINIARG_NULL;
+using boost::format;
+using namespace std;
+
+#define DL for (d=0;d<3;d++)
+
+bool checkParameters(int *numdiv, int *b)
+{
+  for (int i = 0; i < 3; i++)
+    {
+      if (numdiv[i] < 0 || b[i] < 0 || b[i] >= numdiv[i])
+        return false;
+    }
+  return true;
+}
+
+
+struct BoxData
+{
+  float width[3], totwidth[3];
+  float width2[3], totwidth2[3];
+  float bf, s[3], g[3], c[3];
+  coordT *parts;
+  pid_t *orig;
+  pid_t nvpall, nvp, nvpbuf;
+  bool guard_added;
+  double xyz_min[3], xyz_max[3];
+
+  void prepareBox(const PositionData& pdata, int *b_id, int *numdivs);
+
+  void checkParticle(double *xyz, bool& in_main, bool& in_buf);
+  void addGuardPoints();
+  void findBoundary();
+};
+
+void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf)
+{
+  in_main = in_buf = true;
+  for (int d = 0; d < 3; d++)
+    {
+      xyz[d] -= (double)c[d];
+      if (xyz[d] > width2[d])
+        xyz[d] -= width[d];
+      if (xyz[d] < -width2[d])
+        xyz[d] += width[d];
+
+      in_buf = in_buf && (fabs(xyz[d]) < totwidth2[d]);
+      in_main = in_main && (fabs(xyz[d]) <= width2[d]);
+    }
+}
+
+void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs)
+{
+  guard_added = false;
+
+  for (int i = 0; i < 3; i++)
+    {
+      width[i] = (pdata.xyz_max[i] - pdata.xyz_min[i])/numdivs[i];
+      width2[i] = 0.5*width[i];
+
+      totwidth[i] = width[i]+2.*bf;
+      totwidth2[i] = width2[i] + bf;
+
+      s[i] = width[i]/(float)NGUARD;
+      if ((bf*bf - 2.*s[i]*s[i]) < 0.)
+        {
+          printf("bf = %f, s = %f.\n",bf,s[i]);
+          printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.",
+                 sqrt(2.)*width[i]/bf);
+          exit(0);
+        }
+      g[i] = (bf / 2.)*(1. + sqrt(1 - 2.*s[i]*s[i]/(bf*bf)));
+      cout << format("s[%d] = %f, bf = %f, g[%d] = %f.") % s[i] % i % bf % g[i] % i << endl;
+
+      c[i] = b_id[i] * width[i];
+    }
+  cout.flush();
+
+  cout << format("c: %f,%f,%f") % c[0] % c[1] % c[2] << endl;
+
+  /* Assign temporary array*/
+  nvpbuf = 0; /* Number of particles to tesselate, including buffer */
+  nvp = 0; /* Without the buffer */
+
+  for (pid_t i = 0; i < pdata.np; i++)
+    {
+      double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] };
+      bool is_it_in_buf, is_it_in_main;
+
+      checkParticle(xyz, is_it_in_buf, is_it_in_main);
+
+      if (is_it_in_buf)
+        nvpbuf++;
+      if (is_it_in_main)
+        nvp++;
+    }
+
+  nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */
+
+  parts = new coordT[3*nvpbuf];
+  orig = new pid_t[nvpbuf];
+
+  if (parts == 0)
+    {
+      cout << "Unable to allocate parts" << endl;
+      exit(0);
+    }
+  if (orig == 0)
+    {
+      cout << "Unable to allocate orig" << endl;
+      exit(0);
+    }
+
+  nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */
+
+  for (int j = 0; j < 3; j++)
+    {
+      xyz_min[j] = -std::numeric_limits<double>::max();
+      xyz_max[j] = std::numeric_limits<double>::max();
+    }
+  for (pid_t i = 0; i < pdata.np; i++)
+    {
+      bool is_it_in_main, is_it_in_buf;
+      double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] };
+
+      checkParticle(xyz, is_it_in_main, is_it_in_buf);
+
+      if (is_it_in_main) {
+        for (int j = 0; j < 3; j++)
+          {
+            parts[3*nvp+j] = xyz[j];
+            xyz_min[j] = min(xyz_min[j], xyz[j]);
+            xyz_max[j] = max(xyz_max[j], xyz[j]);
+          }
+        orig[nvp] = i;
+        nvp++;
+      }
+  }
+  cout << format("nvp = %d") %nvp << endl;
+  cout << format("x: %f,%f; y: %f,%f; z:%f,%f") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl;
+  nvpbuf = nvp;
+  for (pid_t i = 0; i < pdata.np; i++)
+    {
+      bool is_it_in_main, is_it_in_buf;
+      double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] };
+
+      checkParticle(xyz, is_it_in_main, is_it_in_buf);
+
+      if (is_it_in_buf && !is_it_in_main)
+        {
+          for (int j = 0; j < 3; j++)
+            {
+              parts[3*nvpbuf+j] = xyz[j];
+              xyz_min[j] = min(xyz_min[j], xyz[j]);
+              xyz_max[j] = max(xyz_max[j], xyz[j]);
+            }
+          orig[nvpbuf] = i;
+
+          nvpbuf++;
+        }
+    }
+  cout << format("nvpbuf = %d") % nvpbuf << endl;
+  cout << format("x: %f,%f; y: %f,%f; z:%f,%f\n") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl;
+  nvpall = nvpbuf;
+
+  double predict = 1;
+  for (int j = 0; j < 3; j++)
+    predict *= (width[j]+2*bf);
+  predict *= pdata.np;
+  cout << format("There should be ~ %g points; there are %d\n") % predict % nvpbuf << endl;
+}
+
+void BoxData::addGuardPoints()
+{
+  int rot_g[3][3]  = { {0, 0, 1}, {0,1,0}, {1,0,0} };
+  int rot_si[3][3] = { {1, 0, 0}, {1,0,0}, {0,1,0} };
+  int rot_sj[3][3] = { {0, 1, 0}, {0,0,1}, {0,0,1} };
+
+  for (int k = 0; k < 3; k++)
+    {
+
+      /* Add guard points */
+      for (int i=0; i <= NGUARD; i++)
+        {
+          for (int j=0; j <= NGUARD; j++)
+            {
+              /* Bottom */
+              for (int a = 0; a < 3; a++)
+                parts[3*nvpall+a] = -width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] - rot_g[k][a] * g[a];
+
+              nvpall++;
+
+              /* Top */
+              for (int a = 0; a < 3; a++)
+                parts[3*nvpall+a] = (2*rot_g[k][a]-1)*width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] + rot_g[k][a] * g[a];
+
+              nvpall++;
+            }
+        }
+    }
+}
+
+void BoxData::findBoundary()
+{
+  double BF = std::numeric_limits<double>::max();
+
+  for (int j = 0; j < 3; j++)
+    {
+      xyz_min[j] = -BF;
+      xyz_max[j] = BF;
+    }
+
+  for (pid_t i = nvpbuf; i < nvpall; i++) {
+    for (int j = 0; j < 3; j++)
+      {
+        xyz_min[j] = std::min(xyz_min[j], parts[3*i+j]);
+        xyz_max[j] = std::min(xyz_max[j], parts[3*i+j]);
+      }
+  }
+}
+
+void saveTesselation(const string& outfile, PositionData& pdata, BoxData& boxdata, PARTADJ *adjs, float *vols)
+{
+  ofstream out(outfile.c_str());
+  if (!out)
+    {
+      cout << format("Unable to open %s") % outfile << endl;
+      exit(0);
+    }
+  out.write((char *)&pdata.np, sizeof(pid_t));
+  out.write((char *)&boxdata.nvp, sizeof(pid_t));
+  cout << format("nvp = %d") % boxdata.nvp << endl;
+
+  /* Tell us where the original particles were */
+  out.write((char *)boxdata.orig, sizeof(pid_t)*boxdata.nvp);
+  /* Volumes*/
+  out.write((char *)vols,sizeof(float)*boxdata.nvp);
+  /* Adjacencies */
+  for (pid_t i = 0; i < boxdata.nvp; i++)
+    {
+      out.write((char*)&adjs[i].nadj, sizeof(pid_t));
+      if (adjs[i].nadj > 0)
+        out.write((char *)adjs[i].adj, adjs[i].nadj*sizeof(pid_t));
+      else
+        (cout << "0").flush();
+    }
+  out.close();
+}
+
+int main(int argc, char *argv[]) {
+  PositionData pdata;
+  BoxData boxdata;
+
+  coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
+  pointT intpoints[3*MAXVERVER];
+  string outfile;
+  ofstream out;
+
+  char *suffix, *outDir, *posfile;
+  PARTADJ *adjs;
+  float *vols;
+  pid_t *orig;
+  double border, boxsize[3];
+  int b[3], numdiv[3];
+  double totalvol;
+
+  CosmoTool::MiniArgDesc args[] = {
+    { "POSITION FILE", &posfile, MINIARG_STRING },
+    { "BORDER SIZE", &border, MINIARG_DOUBLE },
+    { "BOX_X", &boxsize[0], MINIARG_DOUBLE },
+    { "BOX_Y", &boxsize[1], MINIARG_DOUBLE },
+    { "BOX_Z", &boxsize[2], MINIARG_DOUBLE },
+    { "SUFFIX", &suffix, MINIARG_STRING },
+    { "NUM_DIVISION_X", &numdiv[0], MINIARG_INT },
+    { "NUM_DIVISION_Y", &numdiv[1], MINIARG_INT },
+    { "NUM_DIVISION_Z", &numdiv[2], MINIARG_INT },
+    { "B0", &b[0], MINIARG_INT },
+    { "B1", &b[1], MINIARG_INT },
+    { "B2", &b[2], MINIARG_INT },
+    { "OUTPUT DIRECTORY", &outDir, MINIARG_STRING },
+    { 0, 0, MINIARG_NULL }
+  };
+
+  if (!CosmoTool::parseMiniArgs(argc, argv, args))
+    return 1;
+
+  if (!checkParameters(numdiv, b))
+    return 2;
+
+  /* Boxsize should be the range in r, yielding a range 0-1 */
+  if (!pdata.readFrom(posfile))
+    return 3;
+
+  (cout << pdata.np << " particles" << endl).flush();
+
+  pdata.findExtrema();
+
+  (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f")
+    % pdata.np
+    % pdata.xyz_min[0] % pdata.xyz_max[0]
+    % pdata.xyz_min[1] % pdata.xyz_max[1]
+    % pdata.xyz_min[2] % pdata.xyz_max[2]).flush();
+
+
+  if (border > 0.)
+    boxdata.bf = border;
+  else
+    boxdata.bf = 0.1;
+
+  boxdata.prepareBox(pdata, b, numdiv);
+  pdata.destroy();
+  boxdata.addGuardPoints();
+
+  adjs = new PARTADJ[boxdata.nvpall];
+  if (adjs == 0)
+    {
+      cout << "Unable to allocate adjs" << endl;
+      return 0;
+    }
+
+  boxdata.findBoundary();
+
+  cout << format("Added guard points to total %d points (should be %d)")
+            % boxdata.nvpall % (boxdata.nvpbuf + 6*(NGUARD+1)*(NGUARD+1)) << endl;
+  cout << format("x: %f,%f; y: %f,%f; z:%f,%f") % boxdata.xyz_min[0] % boxdata.xyz_max[0] % boxdata.xyz_min[1] % boxdata.xyz_max[1] % boxdata.xyz_min[2] % boxdata.xyz_max[2] << endl;
+
+  /* Do tesselation*/
+  printf("File read.  Tessellating ...\n"); fflush(stdout);
+  int exitcode = delaunadj(boxdata.parts, boxdata.nvp, boxdata.nvpbuf, boxdata.nvpall, &adjs);
+  if (exitcode != 0)
+   {
+     printf("Error while tesselating. Stopping here."); fflush(stdout);
+     return 4;
+   }
+
+  /* Now calculate volumes*/
+  printf("Now finding volumes ...\n"); fflush(stdout);
+  vols = new float[boxdata.nvp];
+
+  for (pid_t i = 0; i < boxdata.nvp; i++)
+    { /* Just the original particles
+       * Assign adjacency coordinate array*/
+      /* Volumes */
+      for (int j = 0; j < adjs[i].nadj; j++)
+        {
+          for (int d = 0; d < 3; d++)
+            {
+              deladjs[3*j + d] = boxdata.parts[3*adjs[i].adj[j]+d] - boxdata.parts[3*i+d];
+
+              if (deladjs[3*j+d] < -boxdata.width2[d])
+                deladjs[3*j+d] += boxdata.width[d];
+              if (deladjs[3*j+d] > boxdata.width2[d])
+                deladjs[3*j+d] -= boxdata.width[d];
+            }
+        }
+
+      exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i]));
+      vols[i] *= pdata.np/pdata.V0;
+      if ((i % 1000) == 0)
+        cout << format("%d: %d, %f") % i % adjs[i].nadj % vols[i] << endl;
+    }
+
+  /* Get the adjacencies back to their original values */
+  for (pid_t i=0; i<boxdata.nvp; i++)
+    for (int j=0; j < adjs[i].nadj; j++)
+      adjs[i].adj[j] = orig[adjs[i].adj[j]];
+
+  totalvol = 0.;
+  for (pid_t i=0;i<boxdata.nvp; i++)
+    totalvol += vols[i];
+
+  cout << format("Average volume = %g") % (totalvol/boxdata.nvp) << endl;
+
+  /* Now the output!
+     First number of particles */
+  outfile = str(format("%s/part.%s.%02d.%02d.%02d") % outDir % suffix % b[0] % b[1] % b[2]);
+
+  cout << format("Output to %s") %outfile << endl << endl;
+  saveTesselation(outfile, pdata, boxdata, adjs, vols);
+  delete[] adjs;
+
+  return(0);
+}

c_tools/zobov2/voz1b1/voz_io.cpp

@@ -0,0 +1,67 @@
+#include <limits>
+#include <iostream>
+#include <fstream>
+#include "voz_io.hpp"
+#include "CosmoTool/fortran.hpp"
+#include <cmath>
+
+using namespace CosmoTool;
+using namespace std;
+
+bool PositionData::readFrom(const string& fname)
+{
+  try
+    {
+      UnformattedRead f(fname.c_str());
+
+      f.beginCheckpoint();
+      np = f.readInt32();
+      f.endCheckpoint();
+
+      for (int j = 0; j < 3; j++)
+        {
+          xyz[j] = new float[np];
+
+          f.beginCheckpoint();
+          for (int p = 0; p < np; p++)
+            xyz[j][p] = f.readReal32();
+          f.endCheckpoint();
+        }
+    }
+  catch (const NoSuchFileException& e)
+    {
+      return false;
+    }
+  catch (const InvalidUnformattedAccess& e)
+    {
+      return false;
+    }
+  catch (const EndOfFileException& e)
+    {
+      return false;
+    }
+
+  return true;
+}
+
+void PositionData::findExtrema()
+{
+  const float BF = std::numeric_limits<float>::max();
+
+  for (int i = 0; i < 3; i++)
+    {
+      xyz_min[i] = BF;
+      xyz_max[i] = -BF;
+    }
+
+  V0=1;
+  for (int i = 0; i < 3; i++)
+    {
+      for (pid_t p = 0; p < np; p++)
+        {
+          xyz_min[i] = min(xyz_min[i], xyz[i][p]);
+          xyz_max[i] = max(xyz_max[i], xyz[i][p]);
+        }
+      V0 *= (xyz_max[i]-xyz_min[i]);
+    }
+}

c_tools/zobov2/voz1b1/voz_io.hpp

@@ -0,0 +1,24 @@
+#ifndef __VOZ_IO_HPP
+#define __VOZ_IO_HPP
+
+#include <string>
+
+struct PositionData
+{
+  float *xyz[3];
+  pid_t np;
+  float xyz_min[3], xyz_max[3];
+  float V0;
+
+  void destroy()
+  {
+    for (int j = 0; j < 3; j++)
+      delete[] xyz[j];
+  }
+
+  void findExtrema();
+
+  bool readFrom(const std::string& fname);
+};
+
+#endif

c_tools/zobov2/voz1b1/vozutil.c

@@ -0,0 +1,245 @@
+#include "libqhull/qhull_a.h"
+#include "voz.h"
+
+#define FOREACHvertex2_(vertices) FOREACHsetelement_(vertexT, vertices2,vertex2)
+
+/*char qh_version[] = "user_eg 3.1 2001/10/04";  */
+
+int compar(const void * n1, const void * n2) {
+  int i1,i2;
+
+  i1 = *(int *)n1;
+  i2 = *(int *)n2;
+  return 2*(i1 > i2) - 1 + (i1 == i2);
+}
+
+/* Finds Delaunay adjacencies of a set of points */
+int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) {
+  
+  int dim= 3;	            /* dimension of points */
+  boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
+  char flags[250];          /* option flags for qhull, see qh_opt.htm */
+  FILE *outfile= stdout;    /* output from qh_produce_output()
+			       use NULL to skip qh_produce_output() */
+  FILE *errfile= stderr;    /* error messages from qhull code */
+  int exitcode;             /* 0 if no error from qhull */
+  int curlong, totlong;	    /* memory remaining after qh_memfreeshort */
+  int i, ver, count;
+  int numfacets, numsimplicial, numridges, totneighbors, numneighbors,
+    numcoplanars, numtricoplanars;
+  setT *vertices, *vertices2, *vertex_points, *coplanar_points;
+  vertexT *vertex, **vertexp;
+  vertexT *vertex2, **vertex2p;
+  int vertex_i, vertex_n;
+  facetT *facet, **facetp, *neighbor, **neighborp;
+  pointT *point, **pointp;
+  int numdiv;
+  
+  PARTADJ adjst;
+  
+  int errorReported = 0;
+
+  adjst.adj = (int *)malloc(MAXVERVER*sizeof(int));
+  if (adjst.adj == NULL) {
+    printf("Unable to allocate adjst.adj\n");
+    exit(0);
+  }
+
+  /* Delaunay triangulation*/
+  sprintf (flags, "qhull s d Qt");
+  exitcode= qh_new_qhull (dim, nvpall, points, ismalloc,
+                      flags, outfile, errfile); 
+  
+
+  if (!exitcode) {                  /* if no error */
+    /* 'qh facet_list' contains the convex hull */
+    
+    /* From qh_printvneighbors */
+    qh_countfacets(qh facet_list, NULL, 0, &numfacets, &numsimplicial, 
+		   &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
+    qh_vertexneighbors();
+    vertices= qh_facetvertices (qh facet_list, NULL, 0);
+    vertex_points= qh_settemp (nvpall);
+    coplanar_points= qh_settemp (nvpall);
+    qh_setzero (vertex_points, 0, nvpall);
+    qh_setzero (coplanar_points, 0, nvpall);
+    FOREACHvertex_(vertices)
+      qh_point_add (vertex_points, vertex->point, vertex);
+    FORALLfacet_(qh facet_list) {
+      FOREACHpoint_(facet->coplanarset)
+	qh_point_add (coplanar_points, point, facet);
+    }
+    ver = 0;
+    FOREACHvertex_i_(vertex_points) {
+      (*adjs)[ver].nadj = 0;
+      if (vertex) { 
+	/* Count the neighboring vertices, check that all are real
+	   neighbors */
+	adjst.nadj = 0;
+	FOREACHneighbor_(vertex) {
+	  if ((*adjs)[ver].nadj > -1) {
+	    if (neighbor->visitid) {
+	      vertices2 = neighbor->vertices;
+	      FOREACHvertex2_(vertices2) {
+		if (ver != qh_pointid(vertex2->point)) {
+		  adjst.adj[adjst.nadj] = (int)qh_pointid(vertex2->point);
+		  adjst.nadj ++;
+		  if (adjst.nadj > MAXVERVER) {
+		    printf("Increase MAXVERVER to at least %d!\n",adjst.nadj);
+		    exit(0);
+		  }
+		}
+	      }
+	    } else {
+	      printf(" %d",ver);
+	      (*adjs)[ver].nadj = -1; /* There are unreal vertices here */
+	    }
+	  }
+	}
+      } else (*adjs)[ver].nadj = -2;
+
+      /* Enumerate the unique adjacencies*/
+      if (adjst.nadj >= 4) {
+	qsort((void *)adjst.adj, adjst.nadj, sizeof(int), &compar);
+	count = 1;
+
+	for (i=1; i<adjst.nadj; i++)
+	  if (adjst.adj[i] != adjst.adj[i-1]) {
+	    if (adjst.adj[i] >= nvpbuf && !errorReported) {
+        errorReported = 1;
+	      printf("Guard point encountered.  Increase border and/or nguard.\n");
+	      printf("P:(%f,%f,%f), G: (%f,%f,%f)\n",points[3*ver],points[3*ver+1],points[3*ver+2],
+		     points[3*adjst.adj[i]],points[3*adjst.adj[i]+1],points[3*adjst.adj[i]+2]);
+	    }
+	    count++;
+	  }
+	(*adjs)[ver].adj = (int *)malloc(count*sizeof(int));
+	if ((*adjs)[ver].adj == NULL) {
+	  printf("Unable to allocate (*adjs)[ver].adj\n");
+	  exit(0);
+	}
+	(*adjs)[ver].adj[0] = adjst.adj[0];
+	count = 1;
+	for (i=1; i<adjst.nadj; i++)
+	  if (adjst.adj[i] != adjst.adj[i-1]) {
+	    (*adjs)[ver].adj[count] = adjst.adj[i];
+	    count++;
+	  }
+	(*adjs)[ver].nadj = count;
+      } else {
+	printf("Number of adjacencies %d < 4, particle %d -> %d\n",adjst.nadj,ver,ver);
+	exit(0);
+      }
+      ver++;
+      if (ver == nvp) break;
+    }
+    qh_settempfree (&coplanar_points);
+    qh_settempfree (&vertex_points);
+    qh_settempfree (&vertices);
+  }
+  qh_freeqhull(!qh_ALL);                 /* free long memory */
+  qh_memfreeshort (&curlong, &totlong);  /* free short memory and memory allocator */
+  if (curlong || totlong) 
+    fprintf (errfile, "qhull internal warning (delaunadj): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
+  free(adjst.adj);
+  return exitcode;
+}
+  
+/* Calculates the Voronoi volume from a set of Delaunay adjacencies */
+int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol) {
+  int dim= 3;	            /* dimension of points */
+  boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
+  char flags[250];          /* option flags for qhull, see qh_opt.htm */
+  FILE *outfile= NULL;    /* output from qh_produce_output()
+			       use NULL to skip qh_produce_output() */
+  FILE *errfile= stderr;    /* error messages from qhull code */
+  int exitcode;             /* 0 if no error from qhull */
+  facetT *facet;	    /* set by FORALLfacets */
+  int curlong, totlong;	    /* memory remaining after qh_memfreeshort */
+
+  coordT *point, *normp, *coordp, *feasiblep, *deladj;
+  int i, j, k;
+  boolT zerodiv;
+  float runsum;
+  char region;
+  /*coordT *points;
+    pointT *intpoints;*/
+  
+  /* make point array from adjacency coordinates (add offset)*/
+  /*points = (coordT *)malloc(4*numpoints*sizeof(coordT));
+  if (points == NULL) {
+    printf("Unable to allocate points\n");
+    exit(0);
+    }*/
+  for (i=0; i<numpoints; i++) {
+    runsum = 0.;
+    deladj = deladjs + 3*i;
+    point = points + 4*i;
+    for (j=0;j<3;j++) {
+      runsum += deladj[j]*deladj[j];
+      point[j] = deladj[j];
+    }
+    point[3] = -0.5*runsum;
+  }
+  sprintf (flags, "qhull H0");
+
+  exitcode= qh_new_qhull (4, numpoints, points, ismalloc,
+			  flags, outfile, errfile);
+  
+  numpoints = 0;
+  if (!exitcode) {                  /* if no error */
+    FORALLfacets {
+      numpoints++;
+    }
+    /* Now we know how many points */
+    /*intpoints = (pointT *)malloc(dim*numpoints*sizeof(pointT));
+    if (intpoints == NULL) {
+    printf("Unable to allocate intpoints\n");
+    exit(0);
+    }*/
+
+    j = 0;
+    FORALLfacets {
+      if (!qh feasible_point) {
+	fprintf (stdout, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n");
+	qh_errexit( qh_ERRinput, NULL, NULL);
+      }
+      point= coordp= intpoints + j*3;
+      j++;
+      normp= facet->normal;
+      feasiblep= qh feasible_point;
+      if (facet->offset < -qh MINdenom) {
+	for (k= qh hull_dim; k--; )
+	  *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++);
+      }else {
+	for (k= qh hull_dim; k--; ) {
+	  *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1,
+				   &zerodiv) + *(feasiblep++);
+	  if (zerodiv) {
+	    qh_memfree (point, qh normal_size);
+	    printf("LABELprintinfinite\n");
+	    exit(0);
+	  }
+	}
+      }
+    }
+  }
+  qh_freeqhull (!qh_ALL);
+  qh_memfreeshort (&curlong, &totlong);
+
+  /* Now we calculate the volume */
+  sprintf (flags, "qhull FA");
+  exitcode= qh_new_qhull (dim, numpoints, intpoints, ismalloc,
+			  flags, outfile, errfile);
+  
+  qh_getarea(qh facet_list);
+  *vol = qh totvol;
+  
+  qh_freeqhull (!qh_ALL);
+  qh_memfreeshort (&curlong, &totlong);
+  if (curlong || totlong) 
+    fprintf (errfile, "qhull internal warning (vorvol): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
+  /*free(points); free(intpoints);*/
+  
+  return exitcode;
+}

python_tools/pipeline_source/defaults.py

@@ -23,6 +23,8 @@ import os
 # -----------------------------------------------------------------------------
 # DEFAULT CONFIGURATION
 
+datasetName = ""
+
 startCatalogStage = 1
 endCatalogStage = 3
 

python_tools/pipeline_source/prepareCatalogs.in.py

@@ -29,6 +29,7 @@ import sys
 import void_python_tools as vp
 import argparse
 import imp
+import subprocess
 
 # -----------------------------------------------------------------------------
 
@@ -81,31 +82,39 @@ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1):
   return sampleName
 
 #------------------------------------------------------------------------------
-def getNickName(sampleName):
+def getNickName(setName, sampleName):
 
   splitName = sampleName.split('_')
 
+  nickName = datasetName  
+
   if "ss" in splitName[1]:
-    nickName = "Subsample = " + splitName[1].replace("ss","")
+    nickName += " SS " + splitName[1].replace("ss","")
+    #nickName = "Subsample = " + splitName[1].replace("ss","")
     if "pv" in splitName[2]: 
-      nickName += ", z = " + splitName[3].replace("z","") + ", Pev.Vel."
+      nickName += ", z = " + splitName[3].replace("z","") + " (w/ PV)"
     else:
       nickName += ", z = " + splitName[2].replace("z","")
   elif "hod" in splitName[1]:
-    nickName = "HOD = " + splitName[2]
+    nickName += " HOD " + splitName[2]
     if "pv" in splitName[3]: 
-      nickName += ", z = " + splitName[4].replace("z","") + ", Pec.Vel."
+      nickName += ", z = " + splitName[4].replace("z","") + " (w/ PV)"
     else:
       nickName += ", z = " + splitName[3].replace("z","")
   elif "halos" in splitName[1]:
-    nickName = "Halos, min mass = " + splitName[2].replace("min","")
+    if "none" in splitName[2]:
+      nickName += " All Halos"
+    else:
+      nickName += " Halos > " + splitName[2].replace("min","")
     if "pv" in splitName[3]: 
-      nickName += ", z = " + splitName[4].replace("z","") + ", Pec.Vel."
+      nickName += ", z = " + splitName[4].replace("z","") + " (w/ PV)"
     else:
       nickName += ", z = " + splitName[3].replace("z","")
   else:
     nickName = sampleName
 
+  if "ran" in setName: nickName = "Random" + nickName
+
   return nickName
 
 
@@ -274,7 +283,7 @@ newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, Fal
 
           sampleName = getSampleName(setName, sliceMin, useVel,
                                      iSlice=iSlice, iVol=mySubvolume)
-          nickName = getNickName(sampleName)
+          nickName = getNickName(setName, sampleName)
 
 
           scriptFile.write(sampleInfo.format(dataFile=dataFileName,
@@ -428,7 +437,7 @@ for iSubSample in xrange(len(subSamples)):
           rescale_position = hubble/1000./scale
           shift = lbox/2.
           rescale_velocity = 3.08567802e16/3.1558149984e16
-          command = "%s %s x y z vz vy vx mass | awk '{print $1*%g+%g, $2*%g+%g, $3*%g+%g, $4*%g, $5*%g, $6*%g, $7}' > %s" % (SDFcvt_PATH, dataFile,
+          command = "%s -a 200000 %s x y z vz vy vx mass | awk '{print $1*%g+%g, $2*%g+%g, $3*%g+%g, $4*%g, $5*%g, $6*%g, $7}' > %s" % (SDFcvt_PATH, dataFile,
                                      rescale_position,
                                      shift,
                                      rescale_position,
@@ -521,38 +530,51 @@ if (args.script or args.all) and haloFileBase != "":
     sys.stdout.flush()
 
     # estimate number of halos to get density
-    #if haloFileDummy == '':
+    if haloFileDummy == '':
-    #  dataFile = catalogDir+haloFileBase+fileNums[0]
+      dataFile = catalogDir+haloFileBase+fileNums[0]
-    #else:
+    else:
-    #  dataFile = catalogDir+haloFileBase.replace(haloFileDummy, fileNums[0])
+      dataFile = catalogDir+haloFileBase.replace(haloFileDummy, 
-    #
+                                                 fileNums[0])
-    #inFile = open(dataFile, 'r')
+    numPart = 0
-    #numPart = 0
+    if dataFormat == "sdf":
-    #for (iLine, line) in enumerate(inFile): 
+      SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64"
-    #  if iLine < haloFileNumComLines: continue
+      if minHaloMass == "none":
-    #  line = line.split(haloFileColSep)
+        command = "%s -a 200000 %s mass | wc" % (SDFcvt_PATH, dataFile)
-    #  if minHaloMass == "none" or float(line[haloFileMCol]) > minHaloMass:
+      else:
-    #    numPart += 1
+        command = "%s -a 200000 %s mass | awk '{if ($1>%g) print $1}' | wc" % (SDFcvt_PATH, dataFile, minHaloMass)
-    #inFile.close()
+      numPart = subprocess.check_output(command, shell=True)
+      numPart = int(numPart.split()[0])
+    else:
+      inFile = open(dataFile, 'r')
+      for (iHalo,line) in enumerate(inFile):
+        if iHalo < haloFileNumComLines: continue
+        line = line.split(haloFileColSep)
+        if minHaloMass == "none" or float(line[haloFileMCol]) > minHaloMass:
+          numPart += 1
+        inFile.close()
 
-    #minRadius = 2*int(np.ceil(lbox/numPart**(1./3.)))
+    minRadius = int(np.ceil(lbox/numPart**(1./3.)))
-    minRadies = 10
   
-    setName = prefix+"halos_min"+str(minHaloMass)
+    if minHaloMass != "none": 
+      strMinHaloMass = "%.2e" % minHaloMass
+    else:
+      strMinHaloMass = "none"
+ 
+    setName = prefix+"halos_min"+strMinHaloMass
     fileList = []
     for (iRedshift, redshift) in enumerate(redshifts):
       sampleName = getSampleName(setName, redshift, False)
       outFileName = sampleName+".dat"
       fileList.append(outFileName) 
 
-    writeScript(setName, prefix+"halos_min"+str(minHaloMass)+"_z", "multidark",
+    writeScript(setName, prefix+"halos_min"+strMinHaloMass+"_z", "multidark",
                 scriptDir, catalogDir, fileNums, 
                 redshifts, 
                 numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
                 dataFileNameList = fileList)
    
     if doPecVel:
-      writeScript(setName, prefix+"halos_min"+str(minHaloMass)+"_z", 
+      writeScript(setName, prefix+"halos_min"+strMinHaloMass+"_z", 
                   "multidark",
                   scriptDir, catalogDir, fileNums, 
                   redshifts, 
@@ -583,6 +605,9 @@ if (args.halos or args.all) and haloFileBase != "":
           if "nhalos" in line:
             numPart = int(line.split()[3].strip(';'))
             break
+          if "npart" in line:
+            numPart = int(line.split()[3].strip(';'))
+            break
         inFile.close()
       else:
         for (iLine, line) in enumerate(inFile):
@@ -592,7 +617,12 @@ if (args.halos or args.all) and haloFileBase != "":
             numPart += 1
         inFile.close()
 
-      sampleName = prefix+"halos_min"+str(minHaloMass)+"_z"+redshifts[iRedshift]
+      if minHaloMass != "none": 
+        strMinHaloMass = "%.2e" % minHaloMass
+      else:
+        strMinHaloMass = "none"
+
+      sampleName = prefix+"halos_min"+strMinHaloMass+"_z"+redshifts[iRedshift]
       outFileName = catalogDir+"/"+sampleName+".dat"
       outFile = open(outFileName, 'w')
       outFile.write("%f\n" %(lbox))
@@ -604,8 +634,10 @@ if (args.halos or args.all) and haloFileBase != "":
 
       if dataFormat == "sdf":
         SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64"
-        if minHaloMass == "none": minHaloMass = 0.0
+        if minHaloMass == "none":
-        command = "%s %s mass id x y z vz vy vx | awk '{if ($1>%g) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, minHaloMass, outFileName )
+          command = "%s -a 200000 %s mass ident x y z vz vy vx | awk '{print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, outFileName )
+        else:
+          command = "%s -a 200000 %s mass ident x y z vz vy vx | awk '{if ($1>%g) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, minHaloMass, outFileName )
         os.system(command)
         outFile = open(outFileName, 'a')
         outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n")
@@ -688,15 +720,38 @@ if (args.script or args.all) and haloFileBase != "":
       fileList.append(outFileName) 
 
     print "  ", thisHod['name']
+
+    # estimate number of halos to get density
+    if haloFileDummy == '':
+      dataFile = catalogDir+haloFileBase+fileNums[0]
+    else:
+      dataFile = catalogDir+haloFileBase.replace(haloFileDummy,
+                                                 fileNums[0])
+    numPart = 0
+    if dataFormat == "sdf":
+      SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64"
+      command = "%s -a 200000 %s mass | awk '{if ($1>%g) print $1}' | wc" % (SDFcvt_PATH, dataFile, thisHod['Mcut'])
+      numPart = subprocess.check_output(command, shell=True)
+      numPart = int(numPart.split()[0])
+    else:
+      inFile = open(dataFile, 'r')
+      for (iHalo,line) in enumerate(inFile):
+        if iHalo < haloFileNumComLines: continue
+        line = line.split(haloFileColSep)
+        if float(line[haloFileMCol]) > thisHod['Mcut']: numPart += 1
+        inFile.close()
+
+    minRadius = int(np.ceil(lbox/numPart**(1./3.)))
+
     setName = prefix+"hod_"+thisHod['name']
     writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark",
                 scriptDir, catalogDir, fileNums, redshifts, 
-                numSubvolumes, numSlices, False, lbox, 15, omegaM,
+                numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
                 dataFileNameList = fileList)
     if doPecVel:
       writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark",
                   scriptDir, catalogDir, fileNums, redshifts, 
-                  numSubvolumes, numSlices, True, lbox, 15, omegaM,
+                  numSubvolumes, numSlices, True, lbox, minRadius, omegaM,
                   dataFileNameList = fileList)
 
 if (args.hod or args.all) and haloFileBase != "":
@@ -716,7 +771,7 @@ if (args.hod or args.all) and haloFileBase != "":
       inFile = haloFile
       outFile = haloFile+"_temp"
       SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64"
-      command = "%s %s mass x y z vx vy vz>>%s" % (SDFcvt_PATH, inFile, outFile)
+      command = "%s -a 200000 %s mass x y z vx vy vz>>%s" % (SDFcvt_PATH, inFile, outFile)
       os.system(command)
       haloFile = outFile
 
@@ -730,13 +785,11 @@ if (args.hod or args.all) and haloFileBase != "":
                                      hubble=hubble,
                                      redshift=redshift,
                                      Mmin=thisHod['Mmin'],
-                                     #Mmin=1.23e13,
                                      M1=thisHod['M1'],
                                      sigma_logM=thisHod['sigma_logM'],
                                      alpha=thisHod['alpha'],
                                      Mcut=thisHod['Mcut'],
                                      galden=thisHod['galDens'],
-                                     #galden=0.000225,
                                      haloFile=haloFile,
                                      haloFileFormat=dataFormat,
                                      numPartPerSide=numPart**(1/3.),
@@ -744,7 +797,13 @@ if (args.hod or args.all) and haloFileBase != "":
                                      workDir=catalogDir))
       parFile.close()
 
-      os.system(hodPath+" "+parFileName+">& /dev/null")
+      tempFile = "./hod.out"
+      os.system(hodPath+" "+parFileName+">& " + tempFile)
+      for line in open(tempFile):
+        if "MLO" in line:
+          print "     (minimum halo mass = ", line.split()[1], ")"
+          break
+      os.unlink(tempFile)
 
       sampleName = getSampleName(prefix+"hod_"+thisHod['name'], redshift, False)
       outFileName = catalogDir+"/"+sampleName+".dat"

python_tools/void_python_tools/plotting/plotDefs.py

@@ -20,9 +20,8 @@
 LIGHT_SPEED = 299792.458 
 
 colorList = ['r', 'b', 'g', 'y', 'c', 'm', 
-             'brown', 'grey',
+             'darkred', 'grey',
-             'darkred', 'orange', 'pink', 'darkblue',
+             'orange', 'darkblue',
-             'lightblue', 'chocolate',
              'indigo', 'lightseagreen', 'maroon', 'olive',
              'royalblue', 'palevioletred', 'seagreen', 'tomato',
              'aquamarine', 'darkslateblue',