merged guilhem's and my updates

c_tools/libzobov/loadZobov.cpp

@@ -1,3 +1,4 @@
+#include <cassert>
 #include <string>
 #include <fstream>
 #include <iostream>
@@ -152,6 +153,7 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName,
       for (int j = 0; j < z.allVoids[volId].zId.size(); j++)
         {
           int zzid = z.allVoids[volId].zId[j];
+          assert(zzid < z.allZones.size());
 	  actualNumber += z.allZones[zzid].pId.size();
         }
       

c_tools/libzobov/particleInfo.cpp

@@ -1,3 +1,4 @@
+#include <cstdlib>
 #include <netcdfcpp.h>
 #include <CosmoTool/fortran.hpp>
 #include "particleInfo.hpp"
@@ -5,6 +6,34 @@
 using namespace std;
 using namespace CosmoTool;
 
+template<bool failure>
+double retrieve_attr_safe_double(NcFile& f, const char *name, double defval)
+{
+  NcAtt *a = f.get_att(name);
+  if (a == 0)
+    {
+      if (failure)
+        abort();
+      return defval;
+    }
+  return a->as_double(0);
+}
+
+template<bool failure>
+int retrieve_attr_safe_int(NcFile& f, const char *name, int defval)
+{
+  NcAtt *a = f.get_att(name);
+  if (a == 0)
+    {
+      if (failure)
+        abort();
+      return defval;
+    }
+  return a->as_int(0);
+}
+
+  
+
 bool loadParticleInfo(ParticleInfo& info,
 		      const std::string& particles, 
 		      const std::string& extra_info)
@@ -13,18 +42,19 @@ bool loadParticleInfo(ParticleInfo& info,
       int isObservation;
 
   NcFile f_info(extra_info.c_str());
+  NcError nerr(NcError::verbose_nonfatal);
  
   if (!f_info.is_valid())
     return false;
 
-  info.ranges[0][0] = f_info.get_att("range_x_min")->as_double(0);
+  info.ranges[0][0] = retrieve_attr_safe_double<true>(f_info, "range_x_min", 0);
-  info.ranges[0][1] = f_info.get_att("range_x_max")->as_double(0);
+  info.ranges[0][1] = retrieve_attr_safe_double<true>(f_info, "range_x_max", 0);
-  info.ranges[1][0] = f_info.get_att("range_y_min")->as_double(0);
+  info.ranges[1][0] = retrieve_attr_safe_double<true>(f_info, "range_y_min", 0);
-  info.ranges[1][1] = f_info.get_att("range_y_max")->as_double(0);
+  info.ranges[1][1] = retrieve_attr_safe_double<true>(f_info, "range_y_max", 0);
-  info.ranges[2][0] = f_info.get_att("range_z_min")->as_double(0);
+  info.ranges[2][0] = retrieve_attr_safe_double<true>(f_info, "range_z_min", 0);
-  info.ranges[2][1] = f_info.get_att("range_z_max")->as_double(0);
+  info.ranges[2][1] = retrieve_attr_safe_double<true>(f_info, "range_z_max", 0);
-  info.mask_index = f_info.get_att("mask_index")->as_int(0); //PMS
+  info.mask_index = retrieve_attr_safe_int<true>(f_info, "mask_index", 0);
-  isObservation = f_info.get_att("is_observation")->as_int(0); //PMS
+  isObservation = retrieve_attr_safe_int<false>(f_info, "is_observation", 0);
 
   for (int i = 0; i < 3; i++)
     info.length[i] = info.ranges[i][1] - info.ranges[i][0];

c_tools/mock/generateMock.cpp

@@ -317,9 +317,9 @@ void buildBox(SimuData *simu, long num_targets, long loaded,
       
       for (int j = 0; j < 3; j++)
 	{
-	  boxed->Pos[j][loaded] = (simu->Pos[j][pid]-ranges[j*2])*mul[j];
+	  boxed->Pos[j][loaded] = max(min((simu->Pos[j][pid]-ranges[j*2])*mul[j], double(1)), double(0));
-	  assert(boxed->Pos[j][loaded] > 0);
+	  assert(boxed->Pos[j][loaded] >= 0);
-	  assert(boxed->Pos[j][loaded] < 1);
+	  assert(boxed->Pos[j][loaded] <= 1);
 	}
       uniqueID[loaded] = (simu_uniqueID != 0) ? simu_uniqueID[pid] : 0;
       expansion_fac[loaded] = efac[pid];
@@ -400,6 +400,7 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a
   mul = new double[3];
   ranges = new double[6];
   snapshot_split = new long[*num_snapshots];
+  expansion_fac = new double[boxed->NumPart];
 
 
   boxed->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
@@ -408,6 +409,7 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a
   boxed->new_attribute("particle_id", particle_id, delete_adaptor<long>);
   boxed->new_attribute("num_snapshots", num_snapshots, delete_adaptor<int>);
   boxed->new_attribute("snapshot_split", snapshot_split, delete_adaptor<long>);
+  boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor<double>);
 
   v_id->get(particle_id, boxed->NumPart);
   v_snap->get(snapshot_split, *num_snapshots);

c_tools/mock/loaders/gadget_loader.cpp

@@ -77,7 +77,7 @@ public:
 
 SimulationLoader *gadgetLoader(const std::string& snapshot, double Mpc_unitLength, int flags)
 {
-  bool singleFile;
+  bool singleFile = false;
   int num_files;
   SimuData *d;
 

c_tools/stacking/pruneVoids.cpp

@@ -91,7 +91,7 @@ int main(int argc, char **argv) {
   double nearestEdge, redshift;
   char line[500], junkStr[10];
   int mask_index;
-  double ranges[2][3], boxLen[3], mul; 
+  double ranges[3][2], boxLen[3], mul; 
   double volNorm, radius;
   int clock1, clock2;
   int periodicX=0, periodicY=0, periodicZ=0;
@@ -180,7 +180,7 @@ int main(int argc, char **argv) {
       part[p].y += ranges[1][0];
       part[p].z += ranges[2][0];
     }
-  } 
+  }
   fclose(fp); 
 
   printf(" Read %d particles...\n", numPartTot);
@@ -350,19 +350,19 @@ int main(int argc, char **argv) {
       if (voids[iVoid].barycenter[0] > boxLen[0])
         voids[iVoid].barycenter[0] = voids[iVoid].barycenter[0] - boxLen[0];
       if (voids[iVoid].barycenter[0] < 0)
-        voids[iVoid].barycenter[0] = boxLen[0] - voids[iVoid].barycenter[0];
+        voids[iVoid].barycenter[0] = boxLen[0] + voids[iVoid].barycenter[0];
     }
     if (periodicY) {
       if (voids[iVoid].barycenter[1] > boxLen[1])
         voids[iVoid].barycenter[1] = voids[iVoid].barycenter[1] - boxLen[1];
-      if (voids[iVoid].barycenter[1] < 1)
+      if (voids[iVoid].barycenter[1] < 0)
-        voids[iVoid].barycenter[1] = boxLen[1] - voids[iVoid].barycenter[1];
+        voids[iVoid].barycenter[1] = boxLen[1] + voids[iVoid].barycenter[1];
     }
     if (periodicZ) {
       if (voids[iVoid].barycenter[2] > boxLen[2])
         voids[iVoid].barycenter[2] = voids[iVoid].barycenter[2] - boxLen[2];
-      if (voids[iVoid].barycenter[2] < 2)
+      if (voids[iVoid].barycenter[2] < 0)
-        voids[iVoid].barycenter[2] = boxLen[2] - voids[iVoid].barycenter[2];
+        voids[iVoid].barycenter[2] = boxLen[2] + voids[iVoid].barycenter[2];
     }
 
     // compute central density
@@ -406,8 +406,8 @@ int main(int argc, char **argv) {
       for (p = 0; p < voids[iVoid].numPart; p++) {
   
         dist[0] = fabs(voidPart[p].x - voids[iVoid].barycenter[0]);
-        dist[0] = fabs(voidPart[p].y - voids[iVoid].barycenter[1]);
+        dist[1] = fabs(voidPart[p].y - voids[iVoid].barycenter[1]);
-        dist[0] = fabs(voidPart[p].z - voids[iVoid].barycenter[2]);
+        dist[2] = fabs(voidPart[p].z - voids[iVoid].barycenter[2]);
 
         if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]);
         if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]);
@@ -599,7 +599,7 @@ int main(int argc, char **argv) {
              voids[iVoid].densCon, 
              voids[iVoid].voidProb);
 
-     fprintf(fpBarycenter, "%d %e %e %e\n", 
+     fprintf(fpBarycenter, "%d  %e %e %e\n", 
              voids[iVoid].voidID,
              voids[iVoid].barycenter[0],
              voids[iVoid].barycenter[1],

external/external_python_build.cmake

@@ -101,6 +101,7 @@ IF(INTERNAL_HEALPY)
            "-DNETCDF4_DIR=${NETCDF_BIN_DIR}"
            "-DPYTHON_LDFLAGS:STRING=${PYTHON_LDFLAGS}"
            "-DPYTHON_LOCAL_SITE_PACKAGE=${PYTHON_LOCAL_SITE_PACKAGE}"
+           "-DSUPPORT_ARCH_NATIVE=${SUPPORT_ARCH_NATIVE}"
            "-DTARGET_PATH=${CMAKE_BINARY_DIR}/ext_build/python" "-P")
 
   ExternalProject_Add(healpy

external/python_build.cmake

@@ -7,6 +7,9 @@ SET(ENV{PYTHONPATH} ${PYTHON_LOCAL_SITE_PACKAGE}:$ENV{PYTHONPATH})
 SET(ENV{CFITSIO_EXT_INC} ${CFITSIO_EXT_INC})
 SET(ENV{CFITSIO_EXT_LIB} ${CFITSIO_EXT_LIB})
 SET(ENV{CFITSIO_EXT_PREFIX} ${CFITSIO_EXT_PREFIX})
+IF (NOT SUPPORT_ARCH_NATIVE)
+  SET(ENV{HEALPY_WITHOUT_NATIVE} 1)
+ENDIF(NOT SUPPORT_ARCH_NATIVE)
 
 SET(PYTHON_BUILD_COMMAND ${PYTHON_EXECUTABLE} setup.py build)
 MESSAGE(STATUS "Running ${PYTHON_BUILD_COMMAND}")

pipeline/generateCatalog.py

@@ -20,6 +20,7 @@ if (len(sys.argv) > 1):
     print "  Cannot find parameter file %s!" % filename
     exit(-1)
   parms = imp.load_source("name", filename)
+  regenerateFlag = False
   globals().update(vars(parms))
 else:
   print " Using default parameters"
@@ -90,7 +91,7 @@ for sample in dataSampleList:
     launchGenerate(sample, GENERATE_PATH, workDir=workDir, 
                    inputDataDir=inputDataDir, zobovDir=zobovDir,
                    figDir=figDir, logFile=logFile, useLCDM=useLCDM,
-                   continueRun=continueRun)
+                   continueRun=continueRun, regenerate=regenerateFlag)
 
   # --------------------------------------------------------------------------
   if (startCatalogStage <= 2) and (endCatalogStage >= 2) and not sample.isCombo:

python_tools/pipeline_source/prepareCatalogs.in.py

@@ -81,6 +81,7 @@ endCatalogStage   = 4
 startAPStage = 1
 endAPStage = 1
 
+regenerateFlag = False
 ZOBOV_PATH = "@CMAKE_BINARY_DIR@/zobov/"
 CTOOLS_PATH = "@CMAKE_BINARY_DIR@/c_tools/"
 freshStack = True

python_tools/void_python_tools/backend/launchers.py

@@ -24,11 +24,18 @@ ncFloat = 'f8' # Double precision
 # -----------------------------------------------------------------------------
 def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, 
                    zobovDir=None, figDir=None, logFile=None, useLCDM=False,
-                   continueRun=None):
+                   continueRun=None,regenerate=False):
 
   if sample.dataType == "observation":
     sampleName = sample.fullName
 
+    if regenerate:
+      inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par"
+      outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName
+    else:
+      inputParameterFlag = ""
+      outputFile = zobovDir+"/zobov_slice_" + sampleName
+
     if sample.dataFile == "":
       datafile = inputDataDir+"/"+sampleName
     else:
@@ -50,16 +57,17 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
       zMax %g
       density_fake %g
       %s
-      """ % (datafile, maskFile, zobovDir+"/zobov_slice_"+sampleName,
+      %s
+      """ % (datafile, maskFile, outputFile,
              zobovDir+"/zobov_slice_"+sampleName+".par",
              sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity,
-             useLCDMFlag)
+             useLCDMFlag, inputParameterFlag)
 
     parmFile = os.getcwd()+"/generate_"+sample.fullName+".par"
 
     file(parmFile, mode="w").write(conf)
 
-    if not (continueRun and jobSuccessful(logFile, "Done!\n")):
+    if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")):
       cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile)
       os.system(cmd)
       if jobSuccessful(logFile, "Done!\n"):
@@ -97,6 +105,13 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
 
     datafile = inputDataDir+"/"+sample.dataFile
 
+    if regenerate:
+      inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par"
+      outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName 
+    else:
+      inputParameterFlag = ""
+      outputFile = zobovDir+"/zobov_slice_" + sampleName
+
     if sample.usePecVel:
       includePecVelString = "peculiarVelocities"
     else:
@@ -134,20 +149,21 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
       rangeZ_min %g
       rangeZ_max %g
       subsample %g
-      """ % (dataFileLine, zobovDir+"/zobov_slice_"+sampleName,
+      %s
+      """ % (dataFileLine, outputFile,
              zobovDir+"/zobov_slice_"+sampleName+".par",
              includePecVelString,
              useLightConeString,
              sample.dataUnit,
              xMin, xMax, yMin, yMax,
              sample.zBoundaryMpc[0], sample.zBoundaryMpc[1],
-             sample.subsample)
+             sample.subsample,inputParameterFlag)
 
     parmFile = os.getcwd()+"/generate_"+sample.fullName+".par"
 
     file(parmFile, mode="w").write(conf)
 
-    if not (continueRun and jobSuccessful(logFile, "Done!\n")):
+    if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")):
       cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile)
       os.system(cmd)
       if jobSuccessful(logFile, "Done!\n"):
@@ -258,6 +274,7 @@ def launchPrune(sample, binPath, thisDataPortion=None,
        sample.boxLen <= 1.e-1: 
       periodicLine += "z"
     periodicLine += "' "
+    periodicLine = ""
 
   if not (continueRun and jobSuccessful(logFile, "NetCDF: Not a valid ID\n")):
     cmd = binPath
@@ -300,6 +317,9 @@ def launchPrune(sample, binPath, thisDataPortion=None,
                           str(thisDataPortion)+"_"+\
                           str(sampleName)+".out"
     cmd += " &> " + logFile
+    f=file("run_prune.sh",mode="w")
+    f.write(cmd)
+    f.close()
     os.system(cmd)
 
     if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or \
@@ -383,14 +403,14 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
                 voidDir=None, freshStack=True, runSuffix=None,
                 zobovDir=None,
                 INCOHERENT=False, ranSeed=None, summaryFile=None, 
-                continueRun=None, dataType=None):
+                continueRun=None, dataType=None, prefixRun=""):
 
   sampleName = sample.fullName
 
   runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin,
                              stack.rMax, thisDataPortion)
 
-  logFile = logDir+"/stack_"+sampleName+"_"+runSuffix+".out"
+  logFile = logDir+("/%sstack_"%prefixRun)+sampleName+"_"+runSuffix+".out"
  
   treeFile = voidDir+"/tree.data"
 
@@ -417,7 +437,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
     maxDen = 0.2*float(maskIndex)/float(totalPart)
   else:
     maskIndex = 999999999
-    maxDen = 0.2
+    maxDen = -0.8
 
   if INCOHERENT:
     nullTestFlag = "INCOHERENT"
@@ -458,7 +478,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
    zobovDir+"/vol_"+sampleName+".dat",
    stack.rMin,
    stack.rMax,
-   zobovDir+"/zobov_slice_"+sampleName,
+   zobovDir+("/%szobov_slice_"%prefixRun)+sampleName,
    zobovDir+"/zobov_slice_"+sampleName+".par",
    maxDen,
    centralRadius,
@@ -474,7 +494,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
    zobovDir+"/boundaryDistances_"+thisDataPortion+"_"+sampleName+".out",
    rescaleFlag)
 
-  parmFile = os.getcwd()+"/stack_"+sample.fullName+".par"
+  parmFile = os.getcwd()+("/%sstack_"%prefixRun)+sample.fullName+".par"
 
   fp = file(parmFile, mode="w")
   fp.write(conf)

zobov/voztie.c

@@ -8,6 +8,7 @@ int main(int argc, char *argv[]) {
   FILE *part, *adj, *vol;
   char partfile[200], *suffix, adjfile[200], volfile[200], *outDir;
   float *vols, volstemp;
+  int *cnt_adj;
   
   PARTADJ *adjs;
 
@@ -78,13 +79,17 @@ int main(int argc, char *argv[]) {
   if (mockIndex == -1) mockIndex = np;
   // END PMS
 
+  cnt_adj = (int *)malloc(np*sizeof(int));
+  if (cnt_adj == NULL)
+    printf("Could not allocate cnt_adj.\n");
+
   adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ));
   if (adjs == NULL) printf("Couldn't allocate adjs.\n");
   vols = (float *)malloc(np*sizeof(float));
   if (vols == NULL) printf("Couldn't allocate vols.\n");
   orig = (int *)malloc(nvpmax*sizeof(int));
   if (orig == NULL) printf("Couldn't allocate orig.\n");
-  if ((vols == NULL) || (orig == NULL) || (adjs == NULL)) {
+  if ((cnt_adj==NULL) || (vols == NULL) || (orig == NULL) || (adjs == NULL)) {
     printf("Not enough memory to allocate. Exiting.\n");
     exit(0);
   }
@@ -238,11 +243,24 @@ printf("\n");
 // END OMS
   /* Adjacencies: first the numbers of adjacencies, 
      and the number we're actually going to write per particle */
+
+  // Recount the number of adjacencies after merge
+  for(i=0;i<mockIndex;i++)
+    cnt_adj[i] = 0;
+  for(i=0;i<mockIndex;i++)
+    {
+      for (j=0;j<adjs[i].nadj;j++)
+        {
+          cnt_adj[adjs[i].adj[j]]++;
+          cnt_adj[i]++;
+        }
+    }
+
 // PMS
   for (i=0;i<mockIndex;i++)
   //for (i=0;i<np;i++)
 // END PMS
-    fwrite(&adjs[i].nadj,1,sizeof(int),adj);
+    fwrite(&cnt_adj[i],1,sizeof(int),adj);
     
   /* Now the lists of adjacencies (without double counting) */
   // PMS
@@ -253,9 +271,11 @@ printf("\n");
       nout = 0;
       for (j=0;j<adjs[i].nadj; j++) if (adjs[i].adj[j] > i) nout++;
       fwrite(&nout,1,sizeof(int),adj);      
-      for (j=0;j<adjs[i].nadj; j++) 
+      for (j=0;j<adjs[i].nadj; j++) {
-	if (adjs[i].adj[j] > i) 
+        int id = adjs[i].adj[j];
+	if (id > i) 
 	  fwrite(&(adjs[i].adj[j]),1,sizeof(int),adj);
+      }
     }
 
   fclose(adj);
@@ -275,5 +295,10 @@ printf("\n");
 
   fclose(vol);
 
+  free(vols);
+  free(cnt_adj);
+  free(adjs);
+  free(orig);
+
   return(0);
 }