diff --git a/c_source/prep/prepObservation.cpp b/c_source/prep/prepObservation.cpp index 15f3e43..2110e50 100644 --- a/c_source/prep/prepObservation.cpp +++ b/c_source/prep/prepObservation.cpp @@ -229,6 +229,7 @@ void flagEdgeGalaxies(prepObservation_info& args , output_data.mask_index = output_data.id_gal.size(); // write a small text file with galaxy position (for diagnostic purposes) + // TODO - remove this FILE *fp; fp = fopen("galaxies.txt", "w"); for (int i = 0; i < data.size(); i++) { @@ -240,7 +241,7 @@ void flagEdgeGalaxies(prepObservation_info& args , } fclose(fp); -/* NOTE: temporarily moved to python for quick debugging. Will move back to +/* NOTE: temporarily moved to python for quick prototyping. Will move back to here once it's all sorted // convert redshift boundaries to covmoving if necessary @@ -463,8 +464,8 @@ int main(int argc, char **argv) fp = fopen("total_particles.txt", "w"); fprintf(fp, "%d", output_data.pos.size()); fclose(fp); - printf("Done!\n"); // END PMS - TODO REMOVE THIS + printf("Done!\n"); return 0; } diff --git a/c_source/prep/prepObservation.ggo b/c_source/prep/prepObservation.ggo index 756a9b4..548ccd9 100644 --- a/c_source/prep/prepObservation.ggo +++ b/c_source/prep/prepObservation.ggo @@ -19,5 +19,3 @@ option "useComoving" - "Convert to real space using LCDM cosmology" flag off option "omegaM" - "Omega Matter for fiducial cosmology" double optional default="0.27" option "nsideForContour" - "HEALPix NSIDE resolution for figuring out mask contours" int optional default="-1" - -option "meanPartSep" - "Estimated mean tracer seperation in h^3 / Mpc^3" double optional default="1" diff --git a/examples/example_observation/example_observation.py b/examples/example_observation/example_observation.py index 8064484..11e5221 100644 --- a/examples/example_observation/example_observation.py +++ b/examples/example_observation/example_observation.py @@ -30,7 +30,7 @@ continueRun = False # 1 : extract redshift slices from data # 2 : void extraction using zobov # 3 : removal of small voids and voids near the edge -startCatalogStage = 2 +startCatalogStage = 1 endCatalogStage = 3 basePath = os.path.dirname(os.path.abspath(__file__)) @@ -51,7 +51,7 @@ figDir = os.path.join(workDir,"figs","example_observation") numZobovThreads = 2 # optimization: number of subdivisions of the volume -numZobovDivisions = 2 +numZobovDivisions = 1 # Maximum density for merging voids # 0 (equivalent to infinitely large value) -> Merge everything (no threshold) diff --git a/python_source/backend/launchers.py b/python_source/backend/launchers.py index 4ab01e2..84c8606 100644 --- a/python_source/backend/launchers.py +++ b/python_source/backend/launchers.py @@ -76,24 +76,6 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, sample.maskFile = outputDir + "/constructed_mask.fits" figureOutMask(datafile, sample.nsideForContour, sample.maskFile) - # compute mean particle separation - (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0], - sample.zRange[1], sample.zRange[0], sample.zRange[1], "all", - sample.omegaM, useComoving=useComoving) - - numTracers = int(open(outputDir+"/mask_index.txt", "r").read()) - sample.meanPartSep = (1.*numTracers/boxVol/nbar)**(-1/3.) - - - # flag edge galaxies - galFile = outputDir + "galaxies.txt" - edgeGalFile = outputDir + "/galaxy_edge_flags.txt" - #edgeMaskFile = outputDir + "/mask_edge_map.fits" - contourFile = outputDir + "/contour_map.fits" - findEdgeGalaxies(galFile, sample.maskFile, edgeGalFile, contourFile, - sample.zBoundary[0], sample.zBoundary[1], sample.omegaM, - useComoving, sample.boundaryWidth, sample.meanPartSep) - if useComoving: useComovingFlag = "useComoving" else: @@ -111,12 +93,11 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, %s omegaM %g nsideForContour %g - meanPartSep %g """ % (datafile, sample.maskFile, outputFile, outputDir+"/zobov_slice_"+sampleName+".par", sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity, useComovingFlag, inputParameterFlag, sample.omegaM, - sample.nsideForContour, sample.meanPartSep) + sample.nsideForContour) parmFile = os.getcwd()+"/prep_"+sample.fullName+".par" @@ -126,6 +107,24 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, arg1 = "--configFile=%s" % parmFile with open(logFile, 'wt') as log: subprocess.call([binPath, arg1], stdout=log, stderr=log) + + # compute mean particle separation + (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0], + sample.zRange[1], sample.zRange[0], sample.zRange[1], "all", + sample.omegaM, useComoving=useComoving) + + numTracers = int(open(outputDir+"/mask_index.txt", "r").read()) + sample.meanPartSep = (1.*numTracers/boxVol/nbar)**(-1/3.) + + # flag edge galaxies with python routine for now + galFile = outputDir+"/zobov_slice_"+sampleName + #galFile = outputDir + "galaxies.txt" + edgeGalFile = outputDir + "/galaxy_edge_flags.txt" + contourFile = outputDir + "/contour_map.fits" + findEdgeGalaxies(galFile, sample.maskFile, edgeGalFile, contourFile, + sample.zBoundary[0], sample.zBoundary[1], sample.omegaM, + useComoving, sample.boundaryWidth, sample.meanPartSep, outputDir) + if jobSuccessful(logFile, "Done!\n"): endTime = time.time() walltime = endTime - startTime @@ -137,6 +136,16 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, else: print("already done!") + + # compute mean particle separation + # TODO - clean up this duplicate to cover all continueRun options + (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0], + sample.zRange[1], sample.zRange[0], sample.zRange[1], "all", + sample.omegaM, useComoving=useComoving) + + numTracers = int(open(outputDir+"/mask_index.txt", "r").read()) + sample.meanPartSep = (1.*numTracers/boxVol/nbar)**(-1/3.) + if os.access(parmFile, os.F_OK): os.unlink(parmFile) if os.access("contour_map.fits", os.F_OK): diff --git a/python_source/backend/surveyTools.py b/python_source/backend/surveyTools.py index 8c63070..fac8bfe 100644 --- a/python_source/backend/surveyTools.py +++ b/python_source/backend/surveyTools.py @@ -131,7 +131,7 @@ def figureOutMask(galFile, nside, outMaskFile): # figures out which galaxies live on a mask or redshift edge def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile, zmin, zmax, omegaM, useComoving, boundaryWidth, - meanPartSep): + meanPartSep, outputDir): if useComoving: zmin = comovingDistance(zmin, Om=omegaM)*LIGHT_SPEED @@ -145,15 +145,50 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile, npix = healpy.nside2npix(nside) # load in galaxies - galPos = np.genfromtxt(galFile) + #galPos = np.genfromtxt(galFile) + galPos = np.genfromtxt(outputDir+"/galaxies.txt") + with open(galFile, 'rb') as File: + chk = np.fromfile(File, dtype=np.int32,count=1) + Np = np.fromfile(File, dtype=np.int32,count=1)[0] + chk = np.fromfile(File, dtype=np.int32,count=1) + + chk = np.fromfile(File, dtype=np.int32,count=1) + x = np.fromfile(File, dtype=np.float32,count=Np) + chk = np.fromfile(File, dtype=np.int32,count=1) + + chk = np.fromfile(File, dtype=np.int32,count=1) + y = np.fromfile(File, dtype=np.float32,count=Np) + chk = np.fromfile(File, dtype=np.int32,count=1) + + chk = np.fromfile(File, dtype=np.int32,count=1) + z = np.fromfile(File, dtype=np.float32,count=Np) + chk = np.fromfile(File, dtype=np.int32,count=1) + + chk = np.fromfile(File, dtype=np.int32,count=1) + RA = np.fromfile(File, dtype=np.float32,count=Np) + chk = np.fromfile(File, dtype=np.int32,count=1) + + chk = np.fromfile(File, dtype=np.int32,count=1) + Dec = np.fromfile(File, dtype=np.float32,count=Np) + chk = np.fromfile(File, dtype=np.int32,count=1) + + chk = np.fromfile(File, dtype=np.int32,count=1) + redshift = np.fromfile(File, dtype=np.float32,count=Np) + chk = np.fromfile(File, dtype=np.int32,count=1) + + #print(galPos.shape) + #galPos = np.column_stack((x,y,z)) + #print(galPos.shape) + flagList = np.zeros(len(galPos[:,0])) galTree = scipy.spatial.cKDTree(galPos) # flag galaxies near mask edges # using the "ray marching" algorithm: follow rays along lines of sight # of all mask edges, flagging nearest neighbor galaxies as we go - + raySteps = np.arange(zmin, zmax, meanPartSep) + #print(meanPartSep, len(raySteps)) contourPixels = np.nonzero(contourMap)[0] #print(contourPixels) @@ -164,34 +199,42 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile, y = raySteps * vec[1] z = raySteps * vec[2] ray = np.array((x,y,z)).T - #print(ray) dist, nearest = galTree.query(ray) flagList[nearest] = 1 - #print(nearest) # flag galaxies near redsfhit boundaries - # TODO - save time by only covering portion of sphere with data - ds = np.sqrt(healpy.nside2pixarea(nside)) / 1000. - phi = np.arange(0, 2*np.pi, ds*2) - theta = np.arange(0, np.pi, ds) - vec = healpy.ang2vec(theta, phi) + # TODO - save time by only covering portion of sphere that has data + sphereIndices = np.arange(len(contourMap)) + vec = healpy.pix2vec(nside, sphereIndices) + vec = np.asarray(vec) maxEdge = zmax * vec + maxEdge = maxEdge.T dist, nearest = galTree.query(maxEdge) - #print(nearest) - #print(galPos[nearest]) flagList[nearest] = 2 minEdge = zmin * vec + minEdge = minEdge.T dist, nearest = galTree.query(minEdge) - #print(nearest) - #print(galPos[nearest]) flagList[nearest] = 3 # output flag information np.savetxt(edgeGalFile, flagList, fmt="%d") + # paint galaxy flags onto healpix map for diagnostics + # TODO - drop this when done testing + flagMap = np.zeros(len(contourMap)) + justEdgeRA = RA[flagList == 1] + justEdgeDec = Dec[flagList == 1] + + phi, theta = convertAngle(justEdgeRA, justEdgeDec) + + ipix = healpy.ang2pix(nside, theta, phi) + np.put(flagMap, ipix, 1) + + healpy.write_map("./flagged_galaxies.fits", flagMap, overwrite=True, + dtype=np.dtype('float64')) # # output galaxy edge flags diff --git a/python_source/vide_pipeline/__main__.py b/python_source/vide_pipeline/__main__.py index 840b4ca..c833ca2 100644 --- a/python_source/vide_pipeline/__main__.py +++ b/python_source/vide_pipeline/__main__.py @@ -112,7 +112,7 @@ for sample in dataSampleList: # ------------------------------------------------------------------------- if (startCatalogStage <= 3) and (endCatalogStage >= 3) and not sample.isCombo: - print(" Pruning void catalogs", "...", end='') + print(" Pruning void catalogs...", end='') sys.stdout.flush() logFile = logDir+"/pruneVoids_"+sampleName+".out"