diff --git a/python_tools/void_python_tools/voidUtil/__init__.py b/python_tools/void_python_tools/voidUtil/__init__.py index 189dfa2..f90294d 100644 --- a/python_tools/void_python_tools/voidUtil/__init__.py +++ b/python_tools/void_python_tools/voidUtil/__init__.py @@ -21,3 +21,6 @@ from catalogUtil import * from plotDefs import * from plotUtil import * +from matchUtil import * +from xcorUtil import * +from profileUtil import * diff --git a/python_tools/void_python_tools/voidUtil/catalogUtil.py b/python_tools/void_python_tools/voidUtil/catalogUtil.py index 2c676cb..3374105 100644 --- a/python_tools/void_python_tools/voidUtil/catalogUtil.py +++ b/python_tools/void_python_tools/voidUtil/catalogUtil.py @@ -215,6 +215,7 @@ class Catalog: numVoids = 0 numPartTot = 0 numZonesTot = 0 + volNorm = 0 boxLen = np.zeros((3)) part = None zones2Parts = None @@ -327,6 +328,7 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadPart=True): partData, boxLen, volNorm, isObservationData, ranges = loadPart(sampleDir) numPartTot = len(partData) catalog.numPartTot = numPartTot + catalog.volNorm = volNorm catalog.partPos = partData catalog.part = [] for i in xrange(len(partData)): diff --git a/python_tools/void_python_tools/voidUtil/matchUtil.py b/python_tools/void_python_tools/voidUtil/matchUtil.py index 8222f21..ef39dae 100755 --- a/python_tools/void_python_tools/voidUtil/matchUtil.py +++ b/python_tools/void_python_tools/voidUtil/matchUtil.py @@ -18,7 +18,7 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. #+ - +__all__=['compareCatalogs',] from void_python_tools.backend import * import imp diff --git a/python_tools/void_python_tools/voidUtil/profileUtil.py b/python_tools/void_python_tools/voidUtil/profileUtil.py new file mode 100644 index 0000000..d0cb31f --- /dev/null +++ b/python_tools/void_python_tools/voidUtil/profileUtil.py @@ -0,0 +1,118 @@ +#+ +# VIDE -- Void IDentification and Examination -- ./python_tools/void_python_tools/plotting/plotTools.py +# Copyright (C) 2010-2013 Guilhem Lavaux +# Copyright (C) 2011-2013 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. +#+ +__all__=['buildProfile','fitHSWProfile','getHSWProfile',] + +from void_python_tools.backend.classes import * +from plotDefs import * +import numpy as np +import os +import void_python_tools.apTools as vp +from scipy.optimize import curve_fit +from scipy.interpolate import interp1d + +def HamausProfile(r, rs, dc): + alpha = -2.0*rs + 4.0 + if rs < 0.91: + beta = 17.5*rs - 6.5 + else: + beta = -9.8*rs + 18.4 + return dc * (1 - (r/rs)**alpha) / (1+ (r)**beta) + 1 + +# ----------------------------------------------------------------------------- +def buildProfile(catalog, rMin, rMax): + +# builds a stacked radial density profile from the given catalog +# catalog: void catalog +# rMin: minimum void radius, in Mpc/h +# rMax: maximum void radius, in Mpc/h +# +# returns: +# binCenters: array of radii in binned profile +# stackedProfile: the stacked density profile +# sigmas: 1-sigma uncertainty in each bin + + rMaxProfile = rMin*3 + 2 + periodicLine = getPeriodic(catalog.sampleInfo) + + print " Building particle tree..." + partTree = getPartTree(catalog) + + print " Selecting voids to stack..." + accepted = (catalog.voids[:].radius > rMin) & (catalog.voids[:].radius < rMax) + voidsToStack = catalog.voids[accepted] + + print " Stacking voids..." + allProfiles = [] + for void in voidsToStack: + center = void.barycenter + + localPart = catalog.partData[ getBall(partTree, center, rMaxProfile) ] + shiftedPart = shiftPart(localPart, center, periodicLine, catalog.ranges) + + dist = np.sqrt(np.sum(shiftedPart[:,:]**2, axis=1)) + thisProfile, radii = np.histogram(dist, bins=nBins, range=(0,rMaxProfile)) + deltaV = 4*np.pi/3*(radii[1:]**3-radii[0:(radii.size-1)]**3) + thisProfile = np.float32(thisProfile) + thisProfile /= deltaV + thisProfile /= catalog.volNorm + allProfiles.append(thisProfile) + binCenters = 0.5*(radii[1:] + radii[:-1]) + + stackedProfile = np.std(allProfiles, axis=0) / np.sqrt(nVoids) + sigmas = np.std(allProfiles, axis=0) / np.sqrt(nVoids) + + return binCenters, stackedProfile, sigmas + + +# ----------------------------------------------------------------------------- +def fitHSWProfile(radii, densities, sigmas): + +# fits the given density profile to the HSW function +# radii: array of radii in r/rV +# densities: array of densities +# sigmas: array of uncertainties +# +# returns: +# popt: best-fit values of dc and rs +# pcov: covariance matrix + + popt, pcov = curve_fit(HamausProfile, radii, densities, + sigma=sigmas) + maxfev=10000, xtol=5.e-3, + p0=[1.0,-1.0]) + + return popt, pcov + + +# ----------------------------------------------------------------------------- +def getHSWProfile(den, radius): + +# returns the HSW profile for the given sample density and void size +# (interpolated from best-fit values) +# density: density of sample +# radius: void size in Mpc/h + +# returns: +# binCenters: array of radii in binned profile +# stackedProfile: the density profile + + sample = catalog.sampleInfo + data = catalog.voids[:].radius + diff --git a/python_tools/void_python_tools/voidUtil/xcorUtil.py b/python_tools/void_python_tools/voidUtil/xcorUtil.py index c0e8609..bad91fe 100644 --- a/python_tools/void_python_tools/voidUtil/xcorUtil.py +++ b/python_tools/void_python_tools/voidUtil/xcorUtil.py @@ -39,8 +39,6 @@ #+ from void_python_tools.backend import * -from void_python_tools.plotting import * -import void_python_tools.xcor as xcorlib import imp import pickle import argparse @@ -55,6 +53,7 @@ from matplotlib import rc from matplotlib.ticker import NullFormatter import random import sys +__all__=['computeCrossCor',] # ------------------------------------------------------------------------------