From aded7a7c2cd87fd1456208d99acfa4bb8bcd04e4 Mon Sep 17 00:00:00 2001 From: "Paul M. Sutter" Date: Wed, 5 Jun 2024 01:06:36 +0200 Subject: [PATCH] setting the stage for on-the-fly mask and edge calculation --- c_source/pruning/pruneVoids.cpp | 2 +- .../example_observation.py | 13 +- python_source/backend/__init__.py | 1 + python_source/backend/cosmologyTools.py | 62 +-------- python_source/backend/launchers.py | 3 +- python_source/backend/surveyTools.py | 129 ++++++++++++++++++ 6 files changed, 143 insertions(+), 67 deletions(-) create mode 100644 python_source/backend/surveyTools.py diff --git a/c_source/pruning/pruneVoids.cpp b/c_source/pruning/pruneVoids.cpp index a42f3bb..34be238 100644 --- a/c_source/pruning/pruneVoids.cpp +++ b/c_source/pruning/pruneVoids.cpp @@ -506,7 +506,7 @@ int main(int argc, char **argv) { for (iVoid = 0; iVoid < numVoids; iVoid++) { voidID = voids[iVoid].voidID; - printf(" DOING %d (of %d) %d %d %f\n", iVoid, numVoids, voidID, + printf(" DOING %d (of %d) %d %d %f\n", iVoid+1, numVoids, voidID, voids[iVoid].numPart, voids[iVoid].radius); diff --git a/examples/example_observation/example_observation.py b/examples/example_observation/example_observation.py index 72469d0..219cb1a 100644 --- a/examples/example_observation/example_observation.py +++ b/examples/example_observation/example_observation.py @@ -82,8 +82,13 @@ newSample = Sample( # assume sample is volume-limited? volumeLimited = True, - # HEALpix mask file - maskFile = inputDataDir+"/example_observation_mask.fits", + # HEALpix mask file - set to None to auto-compute + maskFile = None, + #maskFile = inputDataDir+"/example_observation_mask.fits", + + # if maskFile blank, desired resolution for HEALpix + # mask mapping, otherwise pulled from maskFile + nsideForMask = 128, # radial selection function (if not volume limited) selFunFile = None, @@ -100,12 +105,12 @@ newSample = Sample( # density of mock particles in cubic Mpc/h # (make this as high as you can afford) - fakeDensity = 0.0000001, + fakeDensity = 0.05, # if true, convert to comoving space using LCDM cosmology useComoving = True, - # cosmology + # cosmology assuming flat universe omegaM = 0.3, ) diff --git a/python_source/backend/__init__.py b/python_source/backend/__init__.py index 430e3cc..fc6f056 100644 --- a/python_source/backend/__init__.py +++ b/python_source/backend/__init__.py @@ -20,3 +20,4 @@ from .classes import * from .launchers import * from .cosmologyTools import * +from .surveyTools import * diff --git a/python_source/backend/cosmologyTools.py b/python_source/backend/cosmologyTools.py index 326dbac..d677e63 100644 --- a/python_source/backend/cosmologyTools.py +++ b/python_source/backend/cosmologyTools.py @@ -22,11 +22,10 @@ import numpy as np import scipy.integrate as integrate -import healpy as healpy import os from backend import * -__all__=['expansion', 'angularDiameter', 'aveExpansion', 'getSurveyProps'] +__all__=['expansion', 'angularDiameter', 'aveExpansion'] # returns 1/E(z) for the given cosmology def expansion(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): @@ -51,62 +50,3 @@ def aveExpansion(zStart, zEnd, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): ave = integrate.quad(expansion, zStart, zEnd, args=(Om, Ot, w0, wa))[0] ave = (zEnd-zStart)/ave return ave - -# ----------------------------------------------------------------------------- -# returns the volume and galaxy density for a given redshit slice -def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectionFuncFile=None, useComoving=False): - - LIGHT_SPEED = 299792.458 - - mask = healpy.read_map(maskFile) - area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi - - if useComoving: - zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=0.27) - zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=0.27) - selFunMin = LIGHT_SPEED/100.*angularDiameter(selFunMin, Om=0.27) - selFunMax = LIGHT_SPEED/100.*angularDiameter(selFunMax, Om=0.27) - else: - zmin = zmin * 3000 - zmax = zmax * 3000 - selFunMin *= 3000 - selFunMax *= 3000 - volume = area * (zmax**3 - zmin**3) / 3 - - if selectionFuncFile == None: - nbar = 1.0 - - elif not os.access(selectionFuncFile, os.F_OK): - print(" Warning, selection function file %s not found, using default of uniform selection." % selectionFuncFile) - nbar = 1.0 - - else: - selfunc = np.genfromtxt(selectionFuncFile) - selfunc = np.array(selfunc) - selfunc[:,0] = selfunc[:,0]/100. - selfuncUnity = selfunc - selfuncUnity[:,1] = 1.0 - selfuncMin = selfunc[0,0] - selfuncMax = selfunc[-1,0] - selfuncDx = selfunc[1,0] - selfunc[0,0] - selfuncN = np.size(selfunc[:,0]) - - selFunMin = max(selFunMin, selfuncMin) - selFunMax = min(selFunMax, selfuncMax) - - - def f(z): return selfunc[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 - def fTotal(z): return selfuncUnity[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 - - zrange = np.linspace(selFunMin, selFunMax) - - nbar = scipy.integrate.quad(f, selFunMin, selFunMax) - nbar = nbar[0] - - ntotal = scipy.integrate.quad(fTotal, 0.0, max(selfuncUnity[:,0])) - ntotal = ntotal[0] - - nbar = ntotal / area / nbar - - - return (volume, nbar) diff --git a/python_source/backend/launchers.py b/python_source/backend/launchers.py index da235cc..7ddcf29 100644 --- a/python_source/backend/launchers.py +++ b/python_source/backend/launchers.py @@ -36,6 +36,7 @@ from pylab import figure from netCDF4 import Dataset from backend.classes import * from backend.cosmologyTools import * +from backend.surveyTools import * import pickle import scipy.interpolate as interpolate @@ -531,7 +532,7 @@ def launchPrune(sample, binPath, cmd += " --outputDir=" + zobovDir cmd += " --sampleName=" + str(sampleName) log = open(logFile, 'w') - log.write(f"Command is {cmd}\n") + #log.write(f"Command is {cmd}\n") subprocess.call(cmd, stdout=log, stderr=log, shell=True) log.close() diff --git a/python_source/backend/surveyTools.py b/python_source/backend/surveyTools.py new file mode 100644 index 0000000..90a5ab0 --- /dev/null +++ b/python_source/backend/surveyTools.py @@ -0,0 +1,129 @@ +#+ +# VIDE -- Void IDentification and Examination -- ./python_tools/vide/apTools/chi2/cosmologyTools.py +# Copyright (C) 2010-2014 Guilhem Lavaux +# Copyright (C) 2011-2014 P. M. Sutter +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; version 2 of the License. +# +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +#+ +# a suite of functions to compute expansion rates, angular diameter +# distances, and expected void stretching + +import numpy as np +import healpy as healpy +import os +from backend import * + +__all__=['getSurveyProps', 'getNside', 'figureOutMask', 'findEdgeGalaxies'] + +# ----------------------------------------------------------------------------- +# returns the volume and galaxy density for a given redshit slice +def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectionFuncFile=None, useComoving=False): + + LIGHT_SPEED = 299792.458 + + mask = healpy.read_map(maskFile) + area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi + + if useComoving: + zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=0.27) + zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=0.27) + selFunMin = LIGHT_SPEED/100.*angularDiameter(selFunMin, Om=0.27) + selFunMax = LIGHT_SPEED/100.*angularDiameter(selFunMax, Om=0.27) + else: + zmin = zmin * 3000 + zmax = zmax * 3000 + selFunMin *= 3000 + selFunMax *= 3000 + volume = area * (zmax**3 - zmin**3) / 3 + + if selectionFuncFile == None: + nbar = 1.0 + + elif not os.access(selectionFuncFile, os.F_OK): + print(" Warning, selection function file %s not found, using default of uniform selection." % selectionFuncFile) + nbar = 1.0 + + else: + selfunc = np.genfromtxt(selectionFuncFile) + selfunc = np.array(selfunc) + selfunc[:,0] = selfunc[:,0]/100. + selfuncUnity = selfunc + selfuncUnity[:,1] = 1.0 + selfuncMin = selfunc[0,0] + selfuncMax = selfunc[-1,0] + selfuncDx = selfunc[1,0] - selfunc[0,0] + selfuncN = np.size(selfunc[:,0]) + + selFunMin = max(selFunMin, selfuncMin) + selFunMax = min(selFunMax, selfuncMax) + + + def f(z): return selfunc[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 + def fTotal(z): return selfuncUnity[np.ceil((z-selfuncMin)/selfuncDx), 1]*z**2 + + zrange = np.linspace(selFunMin, selFunMax) + + nbar = scipy.integrate.quad(f, selFunMin, selFunMax) + nbar = nbar[0] + + ntotal = scipy.integrate.quad(fTotal, 0.0, max(selfuncUnity[:,0])) + ntotal = ntotal[0] + + nbar = ntotal / area / nbar + + + return (volume, nbar) + + +# ----------------------------------------------------------------------------- +# returns the nside resolution from the given maskfile +def getNside(maskFile): + nside = 1.0 + + return nside + + +# ----------------------------------------------------------------------------- +# computes the mask from a given datafile and writes it to a file +def figureOutMask(galFile, nside, outMaskFile): + + npix = healpy.nside2npix(nside) + mask = np.zeros((npix)) + +for line in open(galFile): + line = line.split() + RA = np.float(line[3]) + Dec = np.float(line[4]) + z = np.float(line[5]) + + phi = np.pi/180.*RA + theta = Dec*np.pi/180. + theta = np.pi/2. - Dec*np.pi/180. + pos = np.zeros((3)) + + pix = healpy.ang2pix(nside, theta, phi) + mask[pix] = 1. + + healpy.write_map(outMaskFile, mask) + + return mask + +# ----------------------------------------------------------------------------- +# figures out which galaxies live on a mask edge, and also writes the edge +# map to an auxillary file +def findEdgeGalaxies(galFile, maskFile): + + + return