diff --git a/c_source/pruning/pruneVoids.cpp b/c_source/pruning/pruneVoids.cpp index a667a95..2371043 100644 --- a/c_source/pruning/pruneVoids.cpp +++ b/c_source/pruning/pruneVoids.cpp @@ -113,23 +113,14 @@ double expanFun (double z, void * p) { } void openFiles(string outputDir, string sampleName, - string prefix, string dataPortion, int numPartTot, int numKept, - FILE** fpZobov, FILE** fpCenters, - FILE** fpCentersNoCut, - FILE** fpBarycenter, FILE** fpShapes, - FILE** fpSkyPositions); + FILE** fpOutput); -void closeFiles(FILE* fpZobov, FILE* fpCenters, - FILE* fpCentersNoCut, - FILE* fpBarycenter, FILE* fpShapes, - FILE* fpSkyPositions); +void closeFiles(FILE* fpOutput); -void outputVoids(string outputDir, string sampleName, string prefix, - string dataPortion, int numPartTot, +void outputVoids(string outputDir, string sampleName, int numPartTot, vector voids, - bool isObservation, double *boxLen, - bool doTrim, bool doCentralDenCut); + bool isObservation, double *boxLen); int main(int argc, char **argv) { @@ -145,8 +136,7 @@ int main(int argc, char **argv) { return 1; if (!args.configFile_given) { - if (pruneVoids_conf_required (&args, - PRUNEVOIDS_CONF_PACKAGE)) + if (pruneVoids_conf_required (&args, PRUNEVOIDS_CONF_PACKAGE)) return 1; } else { args_params.check_required = 1; @@ -933,25 +923,11 @@ int main(int argc, char **argv) { dataPortions[0] = "central"; dataPortions[1] = "all"; - printf(" Output fully trimmed catalog...\n"); - prefix = ""; - for (int i = 0; i < 2; i++) { - dataPortion = dataPortions[i]; - outputVoids(outputDir, sampleName, prefix, dataPortion, - numPartTot, - voids, - args.isObservation_flag, boxLen, true, true); - } - printf(" Output fully untrimmed catalog...\n"); - prefix = "untrimmed_"; - for (int i = 0; i < 2; i++) { - dataPortion = dataPortions[i]; - outputVoids(outputDir, sampleName, prefix, dataPortion, + outputVoids(outputDir, sampleName, numPartTot, voids, - args.isObservation_flag, boxLen, false, false); - } + args.isObservation_flag, boxLen); clock2 = clock(); printf(" Time: %f sec (for %d voids)\n", @@ -965,131 +941,47 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- void openFiles(string outputDir, string sampleName, - string prefix, string dataPortion, int numPartTot, int numKept, - FILE** fpZobov, FILE** fpCenters, - FILE** fpBarycenter, FILE** fpShapes, - FILE** fpExtra, - FILE** fpSkyPositions) { + FILE** fpOutput) { - *fpZobov = fopen((outputDir+"/"+prefix+"voidDesc_"+dataPortion+"_"+sampleName).c_str(), "w"); - fprintf(*fpZobov, "%d particles, %d voids.\n", numPartTot, numKept); - fprintf(*fpZobov, "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n"); + *fpOutput = fopen((outputDir+"/voidDatabase_"+sampleName).c_str(), "w"); + //fprintf(*fpZobov, "%d particles, %d voids.\n", numPartTot, numKept); + fprintf(*fpOutput, "Void ID, void type, center (x,y,z) (Mpc/h), volume (normalized), volume(Mpc/h^3), radius (Mpc/h), redshift, RA, Dec, density contrast, max extent, nearest edge, num part, parent ID, tree level, num children, central density, core particle, core density, zone vol, zone num part, num zones, void probability, ellipticity, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n"); - *fpBarycenter = fopen((outputDir+"/"+prefix+"macrocenters_"+dataPortion+"_"+sampleName).c_str(), "w"); - - *fpCenters = fopen((outputDir+"/"+prefix+"centers_"+dataPortion+"_"+sampleName).c_str(), "w"); - fprintf(*fpCenters, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part, parent ID, tree level, number of children, central density\n"); - - *fpSkyPositions = fopen((outputDir+"/"+prefix+"sky_positions_"+dataPortion+"_"+sampleName).c_str(), "w"); - fprintf(*fpSkyPositions, "# RA, dec, redshift, radius (Mpc/h), void ID\n"); - - *fpShapes = fopen((outputDir+"/"+prefix+"shapes_"+dataPortion+"_"+sampleName).c_str(), "w"); - fprintf(*fpShapes, "# void ID, ellip, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n"); - - *fpExtra = fopen((outputDir+"/"+prefix+"extraInfo_"+dataPortion+"_"+sampleName).c_str(), "w"); - fprintf(*fpExtra, "# void type, max radius, nearest edge\n"); - } // end openFiles // ---------------------------------------------------------------------------- -void closeFiles(FILE* fpZobov, FILE* fpCenters, - FILE* fpBarycenter, FILE* fpShapes, - FILE* fpExtra, - FILE* fpSkyPositions) { +void closeFiles(FILE* fpOutput) { - fclose(fpZobov); - fclose(fpCenters); - fclose(fpBarycenter); + fclose(fpOutput); + //fclose(fpCenters); + //fclose(fpBarycenter); //fclose(fpDistances); - fclose(fpShapes); - fclose(fpExtra); - fclose(fpSkyPositions); + //fclose(fpShapes); + //fclose(fpExtra); + //fclose(fpSkyPositions); } // end closeFile // ---------------------------------------------------------------------------- -void outputVoids(string outputDir, string sampleName, string prefix, - string dataPortion, int numPartTot, +void outputVoids(string outputDir, string sampleName, int numPartTot, vector voids, - bool isObservation, double *boxLen, bool doTrim, - bool doCentralDenCut) { + bool isObservation, double *boxLen) { - int iVoid; - VOID outVoid; - FILE *fp, *fpZobov, *fpCenters, *fpCentersNoCut, *fpBarycenter, - *fpShapes, *fpExtra, *fpSkyPositions; + int iVoid; + VOID outVoid; + FILE *fp, *fpOutput; - - openFiles(outputDir, sampleName, prefix, dataPortion, + openFiles(outputDir, sampleName, numPartTot, voids.size(), - &fpZobov, &fpCenters, &fpBarycenter, - &fpShapes, &fpExtra, &fpSkyPositions); + &fpOutput); - for (iVoid = 0; iVoid < voids.size(); iVoid++) { - outVoid = voids[iVoid]; + for (iVoid = 0; iVoid < voids.size(); iVoid++) { + outVoid = voids[iVoid]; - if (dataPortion == "central" && outVoid.voidType == EDGE_VOID) { - continue; - } - - if (doTrim && outVoid.isLeaf) { - continue; - } - - if (doCentralDenCut && outVoid.hasHighCentralDen) { - continue; - } - - double outCenter[3]; - outCenter[0] = outVoid.macrocenter[0]; - outCenter[1] = outVoid.macrocenter[1]; - outCenter[2] = outVoid.macrocenter[2]; - - //if (isObservation) { - // outCenter[0] = (outVoid.macrocenter[0]-boxLen[0]/2.)*100.; - // outCenter[1] = (outVoid.macrocenter[1]-boxLen[1]/2.)*100.; - // outCenter[2] = (outVoid.macrocenter[2]-boxLen[2]/2.)*100.; - //} - - fprintf(fpZobov, "%d %d %d %f %f %d %d %f %d %f %f\n", - iVoid, - outVoid.voidID, - outVoid.coreParticle, - outVoid.coreDens, - outVoid.zoneVol, - outVoid.zoneNumPart, - outVoid.numZones, - outVoid.vol, - outVoid.numPart, - outVoid.densCon, - outVoid.voidProb); - - fprintf(fpBarycenter, "%d %e %e %e\n", - outVoid.voidID, - outVoid.macrocenter[0], - outVoid.macrocenter[1], - outVoid.macrocenter[2]); - - fprintf(fpCenters, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d %d %d %d %.2f\n", - outCenter[0], - outCenter[1], - outCenter[2], - outVoid.vol, - outVoid.radius, - outVoid.redshift, - 4./3.*M_PI*pow(outVoid.radius, 3), - outVoid.voidID, - outVoid.densCon, - outVoid.numPart, - outVoid.parentID, - outVoid.level, - outVoid.numChildren, - outVoid.centralDen); - double phi = atan2(outVoid.macrocenter[1]-boxLen[1]/2., outVoid.macrocenter[0]-boxLen[0]/2.); if (phi < 0) phi += 2.*M_PI; @@ -1098,16 +990,33 @@ void outputVoids(string outputDir, string sampleName, string prefix, double theta = acos((outVoid.macrocenter[2]-boxLen[2]/2.) / outVoid.redshiftInMpc); double dec = (M_PI/2. - theta) * 180./M_PI; - - fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n", + + fprintf(fpOutput, "%d %d %e %e %e %e %e %e %e %e %e %e %e %e %d %d %d %d %e %d %e %e %d %d %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", + outVoid.voidID, + outVoid.voidType, + outVoid.macrocenter[0], + outVoid.macrocenter[1], + outVoid.macrocenter[2], + outVoid.vol, + 4./3.*M_PI*pow(outVoid.radius, 3), + outVoid.radius, + outVoid.redshift, RA, dec, - outVoid.redshift, - outVoid.radius, - outVoid.voidID); - - fprintf(fpShapes, "%d %.6f %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e\n", - outVoid.voidID, + outVoid.densCon, + outVoid.maxRadius, + outVoid.nearestFlag, + outVoid.numPart, + outVoid.parentID, + outVoid.level, + outVoid.numChildren, + outVoid.centralDen, + outVoid.coreParticle, + outVoid.coreDens, + outVoid.zoneVol, + outVoid.zoneNumPart, + outVoid.numZones, + outVoid.voidProb, outVoid.ellip, gsl_vector_get(outVoid.eval, 0), gsl_vector_get(outVoid.eval, 1), @@ -1123,14 +1032,8 @@ void outputVoids(string outputDir, string sampleName, string prefix, gsl_matrix_get(outVoid.evec, 2 ,2) ); - fprintf(fpExtra, "%d %.5f %.5f\n", - outVoid.voidType, - outVoid.maxRadius, - outVoid.nearestFlag - ); } // end iVoid - closeFiles(fpZobov, fpCenters, fpBarycenter, - fpShapes, fpExtra, fpSkyPositions); + closeFiles(fpOutput); } // end outputVoids diff --git a/python_source/backend/surveyTools.py b/python_source/backend/surveyTools.py index df63f79..9ecf962 100644 --- a/python_source/backend/surveyTools.py +++ b/python_source/backend/surveyTools.py @@ -299,35 +299,35 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile, log.write(" Saving galaxy flags to file...\n") np.savetxt(edgeGalFile, flagList, fmt="%d") - # paint galaxy flags onto healpix map for diagnostics - # TODO - drop this when done testing - log.write(" Saving diagnostic maps to file...\n") - flagMap = np.zeros(len(contourMap)) - justEdgeRA = RA[flagList == 1] - justEdgeDec = Dec[flagList == 1] + ## paint galaxy flags onto healpix map for diagnostics + ## TODO - drop this when done testing + #log.write(" Saving diagnostic maps to file...\n") + #flagMap = np.zeros(len(contourMap)) + #justEdgeRA = RA[flagList == 1] + #justEdgeDec = Dec[flagList == 1] - phi, theta = convertAngle(justEdgeRA, justEdgeDec) + #phi, theta = convertAngle(justEdgeRA, justEdgeDec) - ipix = healpy.ang2pix(nside, theta, phi) - for i in ipix: - flagMap[i] += 1 - #np.put(flagMap, ipix, 1) + #ipix = healpy.ang2pix(nside, theta, phi) + #for i in ipix: + # flagMap[i] += 1 + ##np.put(flagMap, ipix, 1) - healpy.write_map(outputDir+"/flagged_galaxies.fits", flagMap, - overwrite=True, - dtype=np.dtype('float64')) + #healpy.write_map(outputDir+"/flagged_galaxies.fits", flagMap, + # overwrite=True, + # dtype=np.dtype('float64')) - allGalMap = np.zeros(len(contourMap)) - phi, theta = convertAngle(RA, Dec) + #allGalMap = np.zeros(len(contourMap)) + #phi, theta = convertAngle(RA, Dec) - ipix = healpy.ang2pix(nside, theta, phi) - for i in ipix: - allGalMap[i] += 1 - #np.put(allGalMap, ipix, 1) + #ipix = healpy.ang2pix(nside, theta, phi) + #for i in ipix: + # allGalMap[i] += 1 + ##np.put(allGalMap, ipix, 1) - healpy.write_map(outputDir+"/all_galaxies.fits", allGalMap, - overwrite=True, - dtype=np.dtype('float64')) + #healpy.write_map(outputDir+"/all_galaxies.fits", allGalMap, + # overwrite=True, + # dtype=np.dtype('float64')) return diff --git a/python_source/voidUtil/catalogUtil.py b/python_source/voidUtil/catalogUtil.py index fe2deb4..e8cb0d9 100644 --- a/python_source/voidUtil/catalogUtil.py +++ b/python_source/voidUtil/catalogUtil.py @@ -31,7 +31,10 @@ import os NetCDFFile = Dataset ncFloat = 'f8' -# ----------------------------------------------------------------------------- +CATALOG_V1 = 1 +CATALOG_V2 = 2 + +# ----------------------------------------------------------------------- def loadPart(sampleDir): print(" Loading particle data...") sys.stdout.flush() @@ -270,11 +273,14 @@ def loadVoidCatalog(sampleDir, print("Loading catalog from ", sampleDir) - isOldCatalog = os.path.exists(sampleDir+"/mask_index.txt") + if os.path.exists(sampleDir+"/mask_index.txt"): + version = CATALOG_V1 + else: + version = CATALOG_V2 - if isOldCatalog and clearNearBoundaries: + if version == CATALOG_V1 and clearNearBoundaries: print("WARNING: Old catalog. Unable to clear near boundaries.") - if isOldCatalog and maxCentralDen != -1: + if version == CATALOG_V1 and maxCentralDen != -1: print("WARNING: Old catalog. Central density cuts already applied.") if replicateOldCentralVoids: @@ -310,19 +316,24 @@ def loadVoidCatalog(sampleDir, # for new catalogs, we will load by default the whole shebang, then # apply filters later. for old catalogs, we need to pick the right file - prefix = "untrimmed_" - if isOldCatalog and clearTree: prefix = "" + if version == CATALOG_V1: + if clearTree: + prefix = "" + else: + prefix = "untrimmed_" - dataPortion = "all" - if isOldCatalog and clearEdges: dataPortion = "central" + if clearEdges: + dataPortion = "central" + else: + dataPortion = "all" - print("Loading voids...") - fileName = sampleDir+"/"+prefix+"voidDesc_"+dataPortion+"_"+sample.fullName+".out" + print("Loading version-1 voids...") + fileName = sampleDir+"/"+prefix+"voidDesc_"+dataPortion+"_"+sample.fullName+".out" - catData = np.loadtxt(fileName, comments="#", skiprows=2) - catalog.voids = [] - for line in catData: - catalog.voids.append(Bunch(iVoid = int(line[0]), + catData = np.loadtxt(fileName, comments="#", skiprows=2) + catalog.voids = [] + for line in catData: + catalog.voids.append(Bunch(iVoid = int(line[0]), voidID = int(line[1]), coreParticle = line[2], coreDens = line[3], @@ -335,7 +346,6 @@ def loadVoidCatalog(sampleDir, voidProb = line[10], # below values to be read in or computed later radius = 0., - macrocenter = np.zeros((3)), redshift = 0, RA = 0, Dec = 0, @@ -351,78 +361,130 @@ def loadVoidCatalog(sampleDir, nearestEdge = 0. )) - catalog.numVoids = len(catalog.voids) - print(" Read %d voids" % catalog.numVoids) + catalog.numVoids = len(catalog.voids) + print(" Read %d voids" % catalog.numVoids) - print("Loading macrocenters...") - iLine = 0 - for line in open(sampleDir+"/"+prefix+"macrocenters_"+dataPortion+"_"+sample.fullName+".out"): - line = line.split() - catalog.voids[iLine].macrocenter[0] = float(line[1]) - catalog.voids[iLine].macrocenter[1] = float(line[2]) - catalog.voids[iLine].macrocenter[2] = float(line[3]) - iLine += 1 + print("Loading macrocenters...") + iLine = 0 + for line in open(sampleDir+"/"+prefix+"macrocenters_"+dataPortion+"_"+sample.fullName+".out"): + line = line.split() + catalog.voids[iLine].macrocenter[0] = float(line[1]) + catalog.voids[iLine].macrocenter[1] = float(line[2]) + catalog.voids[iLine].macrocenter[2] = float(line[3]) + iLine += 1 - iLine = 0 - fileName = sampleDir+"/"+prefix+"sky_positions_"+dataPortion+"_"+sample.fullName+".out" - catData = np.loadtxt(fileName, comments="#") - for line in catData: - catalog.voids[iLine].RA = float(line[0]) - catalog.voids[iLine].Dec = float(line[1]) - iLine += 1 + iLine = 0 + fileName = sampleDir+"/"+prefix+"sky_positions_"+dataPortion+"_"+sample.fullName+".out" + catData = np.loadtxt(fileName, comments="#") + for line in catData: + catalog.voids[iLine].RA = float(line[0]) + catalog.voids[iLine].Dec = float(line[1]) + iLine += 1 - print("Loading derived void information...") - fileName = sampleDir+"/"+prefix+"centers_"+dataPortion+"_"+sample.fullName+".out" - catData = np.loadtxt(fileName, comments="#") - for (iLine,line) in enumerate(catData): - catalog.voids[iLine].volume = float(line[6]) - catalog.voids[iLine].radius = float(line[4]) - catalog.voids[iLine].redshift = float(line[5]) - catalog.voids[iLine].parentID = float(line[10]) - catalog.voids[iLine].treeLevel = float(line[11]) - catalog.voids[iLine].numChildren = float(line[12]) - catalog.voids[iLine].centralDen = float(line[13]) - iLine += 1 - - fileName = sampleDir+"/"+prefix+"shapes_"+dataPortion+"_"+sample.fullName+".out" - catData = np.loadtxt(fileName, comments="#") - for (iLine,line) in enumerate(catData): - catalog.voids[iLine].ellipticity = float(line[1]) - - catalog.voids[iLine].eigenVals[0] = float(line[2]) - catalog.voids[iLine].eigenVals[1] = float(line[3]) - catalog.voids[iLine].eigenVals[2] = float(line[4]) - - catalog.voids[iLine].eigenVecs[0][0] = float(line[5]) - catalog.voids[iLine].eigenVecs[0][1] = float(line[6]) - catalog.voids[iLine].eigenVecs[0][2] = float(line[7]) - - catalog.voids[iLine].eigenVecs[1][0] = float(line[8]) - catalog.voids[iLine].eigenVecs[1][1] = float(line[9]) - catalog.voids[iLine].eigenVecs[1][2] = float(line[10]) - - catalog.voids[iLine].eigenVecs[2][0] = float(line[11]) - catalog.voids[iLine].eigenVecs[2][1] = float(line[12]) - catalog.voids[iLine].eigenVecs[2][2] = float(line[13]) - - iLine += 1 - - fileName = sampleDir+"/"+prefix+"extraInfo_"+dataPortion+"_"+sample.fullName+".out" - if os.path.exists(fileName): + print("Loading derived void information...") + fileName = sampleDir+"/"+prefix+"centers_"+dataPortion+"_"+sample.fullName+".out" catData = np.loadtxt(fileName, comments="#") for (iLine,line) in enumerate(catData): - catalog.voids[iLine].voidType = int(line[0]) - catalog.voids[iLine].maxRadius = float(line[1]) - catalog.voids[iLine].nearestEdge = float(line[2]) + catalog.voids[iLine].volume = float(line[6]) + catalog.voids[iLine].radius = float(line[4]) + catalog.voids[iLine].redshift = float(line[5]) + catalog.voids[iLine].parentID = float(line[10]) + catalog.voids[iLine].treeLevel = float(line[11]) + catalog.voids[iLine].numChildren = float(line[12]) + catalog.voids[iLine].centralDen = float(line[13]) + iLine += 1 + + fileName = sampleDir+"/"+prefix+"shapes_"+dataPortion+"_"+sample.fullName+".out" + catData = np.loadtxt(fileName, comments="#") + for (iLine,line) in enumerate(catData): + catalog.voids[iLine].ellipticity = float(line[1]) + + catalog.voids[iLine].eigenVals[0] = float(line[2]) + catalog.voids[iLine].eigenVals[1] = float(line[3]) + catalog.voids[iLine].eigenVals[2] = float(line[4]) + + catalog.voids[iLine].eigenVecs[0][0] = float(line[5]) + catalog.voids[iLine].eigenVecs[0][1] = float(line[6]) + catalog.voids[iLine].eigenVecs[0][2] = float(line[7]) + + catalog.voids[iLine].eigenVecs[1][0] = float(line[8]) + catalog.voids[iLine].eigenVecs[1][1] = float(line[9]) + catalog.voids[iLine].eigenVecs[1][2] = float(line[10]) + + catalog.voids[iLine].eigenVecs[2][0] = float(line[11]) + catalog.voids[iLine].eigenVecs[2][1] = float(line[12]) + catalog.voids[iLine].eigenVecs[2][2] = float(line[13]) + + iLine += 1 - iLine += 1 else: - print(" Old catalog: extra info file not found") + print("Loading version-2 voids...") + fileName = sampleDir+"/"+prefix+"voidDatabase_"+sample.fullName+".out" + + catalog.voids = [] + for line in catData: + + macrocenter = np.zeros((3)) + macrocenter[0] = float(line[2]) + macrocenter[1] = float(line[3]) + macrocenter[2] = float(line[4]) + + eigenVals = np.zeros((3)) + eigenVecs = np.zeros((3,3)) + + eigenVals[0] = float(line[26]) + eigenVals[1] = float(line[27]) + eigenVals[2] = float(line[28]) + + eigenVecs[0][0] = float(line[29]) + eigenVecs[0][1] = float(line[30]) + eigenVecs[0][2] = float(line[31]) + eigenVecs[1][0] = float(line[32]) + eigenVecs[1][1] = float(line[33]) + eigenVecs[1][2] = float(line[34]) + eigenVecs[2][0] = float(line[35]) + eigenVecs[2][1] = float(line[36]) + eigenVecs[2][2] = float(line[37]) + + catalog.voids.append(Bunch( + voidID = int(line[0]), + voidType = int(line[1]), + macrocenter = macrocenter, + voidVol = float(line[5]), + volume = float(line[6]), + radius = float(line[7]), + redshift = float(line[8]), + RA = float(line[9]), + Dec = float(line[10]), + densCon = float(line[11]), + maxRadius = float(line[12]), + nearestEdge = float(line[13]), + numPart = int(line[14]), + parentID = int(line[15]), + treeLevel = int(line[16]), + numChildren = int(line[17]), + centralDen = float(line[18]), + coreParticle = int(line[19]), + coreDens = float(line[20]), + zoneVol = float(line[21]), + zoneNumPart = int(line[22]), + numZones = int(line[23]), + voidProb = float(line[24]), + ellipticity = float(line[25]), + eigenVals = eigenVals, + eigenVecs = eigenVecs + )) + + catalog.numVoids = len(catalog.voids) + print(" Read %d voids" % catalog.numVoids) + + + # apply filters to new catalogs - if not isOldCatalog: + if version != CATALOG_V1: print("Filtering catalog...") if clearEdges: catalog = filterOnType(catalog, CENTRAL_VOID) if clearTree: catalog = filterOnTreeLevel(catalog, level=-1)