diff --git a/c_tools/libzobov/loadZobov.cpp b/c_tools/libzobov/loadZobov.cpp index 7399874..2d5a2ff 100644 --- a/c_tools/libzobov/loadZobov.cpp +++ b/c_tools/libzobov/loadZobov.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -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(); } diff --git a/c_tools/libzobov/particleInfo.cpp b/c_tools/libzobov/particleInfo.cpp index fa6f467..d3eb2e8 100644 --- a/c_tools/libzobov/particleInfo.cpp +++ b/c_tools/libzobov/particleInfo.cpp @@ -1,3 +1,4 @@ +#include #include #include #include "particleInfo.hpp" @@ -5,6 +6,34 @@ using namespace std; using namespace CosmoTool; +template +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 +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][1] = f_info.get_att("range_x_max")->as_double(0); - info.ranges[1][0] = f_info.get_att("range_y_min")->as_double(0); - info.ranges[1][1] = f_info.get_att("range_y_max")->as_double(0); - info.ranges[2][0] = f_info.get_att("range_z_min")->as_double(0); - info.ranges[2][1] = f_info.get_att("range_z_max")->as_double(0); - info.mask_index = f_info.get_att("mask_index")->as_int(0); //PMS - isObservation = f_info.get_att("is_observation")->as_int(0); //PMS + info.ranges[0][0] = retrieve_attr_safe_double(f_info, "range_x_min", 0); + info.ranges[0][1] = retrieve_attr_safe_double(f_info, "range_x_max", 0); + info.ranges[1][0] = retrieve_attr_safe_double(f_info, "range_y_min", 0); + info.ranges[1][1] = retrieve_attr_safe_double(f_info, "range_y_max", 0); + info.ranges[2][0] = retrieve_attr_safe_double(f_info, "range_z_min", 0); + info.ranges[2][1] = retrieve_attr_safe_double(f_info, "range_z_max", 0); + info.mask_index = retrieve_attr_safe_int(f_info, "mask_index", 0); + isObservation = retrieve_attr_safe_int(f_info, "is_observation", 0); for (int i = 0; i < 3; i++) info.length[i] = info.ranges[i][1] - info.ranges[i][0]; diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index ad334e9..b3864c6 100644 --- a/c_tools/mock/generateMock.cpp +++ b/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]; - assert(boxed->Pos[j][loaded] > 0); - assert(boxed->Pos[j][loaded] < 1); + 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] <= 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); @@ -408,6 +409,7 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a boxed->new_attribute("particle_id", particle_id, delete_adaptor); boxed->new_attribute("num_snapshots", num_snapshots, delete_adaptor); boxed->new_attribute("snapshot_split", snapshot_split, delete_adaptor); + boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor); v_id->get(particle_id, boxed->NumPart); v_snap->get(snapshot_split, *num_snapshots); diff --git a/c_tools/mock/loaders/gadget_loader.cpp b/c_tools/mock/loaders/gadget_loader.cpp index 1c108ee..906f6c1 100644 --- a/c_tools/mock/loaders/gadget_loader.cpp +++ b/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; diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 52dfd4b..92d5090 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/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) - voids[iVoid].barycenter[1] = boxLen[1] - voids[iVoid].barycenter[1]; + if (voids[iVoid].barycenter[1] < 0) + 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) - voids[iVoid].barycenter[2] = boxLen[2] - voids[iVoid].barycenter[2]; + if (voids[iVoid].barycenter[2] < 0) + 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[0] = fabs(voidPart[p].z - voids[iVoid].barycenter[2]); + dist[1] = fabs(voidPart[p].y - voids[iVoid].barycenter[1]); + 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], diff --git a/external/external_python_build.cmake b/external/external_python_build.cmake index 6e78b42..ec5c4ee 100644 --- a/external/external_python_build.cmake +++ b/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 diff --git a/external/python_build.cmake b/external/python_build.cmake index b834e9b..f733405 100644 --- a/external/python_build.cmake +++ b/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}") diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 5941933..fdf3729 100755 --- a/pipeline/generateCatalog.py +++ b/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: diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index 2aa3b1e..08ab0ea 100755 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/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 diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 52bbaae..2671e5c 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/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) diff --git a/zobov/voztie.c b/zobov/voztie.c index 77530cb..cc8c006 100644 --- a/zobov/voztie.c +++ b/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 i) nout++; fwrite(&nout,1,sizeof(int),adj); - for (j=0;j i) + for (j=0;j 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); }