From 02f756e8a5bc644fc5d059fdde3cc04a47c9e265 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 22 Mar 2013 08:51:45 -0500 Subject: [PATCH] integrated nico's cross-correlation backend and enabled plotting of the void distribution --- pipeline/generateCatalog.py | 2 + python_tools/setup.py | 2 +- python_tools/void_python_tools/__init__.py | 5 +- .../void_python_tools/plotting/plotTools.py | 50 ++++++++++++- .../void_python_tools/xcor/__init__.py | 20 +++++ .../void_python_tools/xcor/xcorlib.py | 73 +++++++++++++++++++ 6 files changed, 147 insertions(+), 5 deletions(-) create mode 100644 python_tools/void_python_tools/xcor/__init__.py create mode 100644 python_tools/void_python_tools/xcor/xcorlib.py diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index a786b99..5ef2771 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -150,5 +150,7 @@ if (startCatalogStage <= 4) and (endCatalogStage >= 4): dataPortion=thisDataPortion, setName=setName) plotNumberDistribution(workDir, dataSampleList, figDir, showPlot=False, dataPortion=thisDataPortion, setName=setName) + plotVoidDistribution(workDir, dataSampleList, figDir, showPlot=False, + dataPortion=thisDataPortion, setName=setName) print "\n Done!" diff --git a/python_tools/setup.py b/python_tools/setup.py index b0bc165..ad0518b 100644 --- a/python_tools/setup.py +++ b/python_tools/setup.py @@ -31,7 +31,7 @@ setup( cmdclass = {'build_ext': build_ext}, include_dirs = [np.get_include()], packages= - ['void_python_tools','void_python_tools.backend','void_python_tools.apTools', + ['void_python_tools','void_python_tools.backend','void_python_tools.apTools', 'void_python_tools.xcor', 'void_python_tools.apTools.profiles','void_python_tools.apTools.chi2', 'void_python_tools.plotting'], #ext_modules = [Extension("void_python_tools.chi2.velocityProfileFitNative", ["void_python_tools/chi2/velocityProfileFitNative.pyx"], libraries=["gsl", "gslcblas"]), Extension("void_python_tools.chi2.likelihoo", ["void_python_tools/chi2/likelihood.pyx"], libraries=["gsl", "gslcblas"])] ext_modules = [ diff --git a/python_tools/void_python_tools/__init__.py b/python_tools/void_python_tools/__init__.py index 5f07ae6..a84eedd 100644 --- a/python_tools/void_python_tools/__init__.py +++ b/python_tools/void_python_tools/__init__.py @@ -17,6 +17,7 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. #+ -from void_python_tools.backend import * -from void_python_tools.apTools import * +from void_python_tools.backend import * +from void_python_tools.apTools import * from void_python_tools.plotting import * +from void_python_tools.xcor import * diff --git a/python_tools/void_python_tools/plotting/plotTools.py b/python_tools/void_python_tools/plotting/plotTools.py index 5bc617c..317efef 100644 --- a/python_tools/void_python_tools/plotting/plotTools.py +++ b/python_tools/void_python_tools/plotting/plotTools.py @@ -18,7 +18,7 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. #+ __all__=['plotRedshiftDistribution', 'plotSizeDistribution', 'plot1dProfiles', - 'plotMarg1d', 'plotNumberDistribution'] + 'plotMarg1d', 'plotNumberDistribution', 'plotVoidDistribution'] from void_python_tools.backend.classes import * from plotDefs import * @@ -26,6 +26,7 @@ import numpy as np import os import pylab as plt import void_python_tools.apTools as vp +import void_python_tools.xcor as xcor # ----------------------------------------------------------------------------- def plotRedshiftDistribution(workDir=None, sampleList=None, figDir=None, @@ -266,7 +267,7 @@ def plotNumberDistribution(workDir=None, sampleList=None, figDir=None, boxVol *= 1.e-9 - filename = workDir+"/sample_"+sampleName+"/centers_"+dataPortion+"_"+\ + filename = workDir+"/sample_"+sampleName+"/centers_nocut_"+dataPortion+"_"+\ sampleName+".out" if not os.access(filename, os.F_OK): print "File not found: ", filename @@ -294,5 +295,50 @@ def plotNumberDistribution(workDir=None, sampleList=None, figDir=None, if showPlot: os.system("display %s" % figDir+"/fig_"+plotName+".png") +# ----------------------------------------------------------------------------- +def plotVoidDistribution(workDir=None, sampleList=None, figDir=None, + plotNameBase="dv", + showPlot=False, dataPortion=None, setName=None): + Nmesh = 256 + + for (iSample,sample) in enumerate(sampleList): + plt.clf() + plt.xlabel("x (Mpc/h)") + plt.xlabel("y (Mpc/h)") + + sampleName = sample.fullName + + plotName = plotNameBase+"_"+sampleName + + filename = workDir+"/sample_"+sampleName+"/centers_nocut_"+dataPortion+"_"+\ + sampleName+".out" + + if not os.access(filename, os.F_OK): + print "File not found: ", filename + continue + + void_file = open(filename,'r') + void_header1 = void_file.readline().split(",") + void_data1 = np.reshape(void_file.read().split(), + (-1,len(void_header1))).astype(np.float32) + void_file.close() + + xv = void_data1[:,0:3] + rv = void_data1[:,4] + vv = void_data1[:,6] + dv, wm, ws = xcor.cic(xv, sample.boxLen, Lboxcut = 0., Nmesh = Nmesh) + plt.imshow(np.sum(dv+1,2)/Nmesh,extent=[0,sample.boxLen,0,sample.boxLen], + aspect='equal', + cmap='YlGnBu_r',interpolation='gaussian') + plt.colorbar() + + plt.title("Voids: "+sampleName) + + plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight") + plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") + plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight") + + if showPlot: + os.system("display %s" % figDir+"/fig_"+plotName+".png") diff --git a/python_tools/void_python_tools/xcor/__init__.py b/python_tools/void_python_tools/xcor/__init__.py new file mode 100644 index 0000000..26f152d --- /dev/null +++ b/python_tools/void_python_tools/xcor/__init__.py @@ -0,0 +1,20 @@ +#+ +# VIDE -- Void IDEntification pipeline -- ./python_tools/void_python_tools/plotting/__init__.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. +#+ +from xcorlib import * diff --git a/python_tools/void_python_tools/xcor/xcorlib.py b/python_tools/void_python_tools/xcor/xcorlib.py new file mode 100644 index 0000000..d5b489a --- /dev/null +++ b/python_tools/void_python_tools/xcor/xcorlib.py @@ -0,0 +1,73 @@ +import numpy as np + +def cic( x, Lbox, Lboxcut = 0, Nmesh = 128, weights = None ): + + if weights == None: weights = 1 + wm = np.mean(weights) + ws = np.mean(weights**2) + + d = np.mod(x/(Lbox+2*Lboxcut)*Nmesh,1) + + box = ([Lboxcut,Lbox+Lboxcut],[Lboxcut,Lbox+Lboxcut],[Lboxcut,Lbox+Lboxcut]) + + rho = np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*(1-d[:,1])*(1-d[:,2]))[0] \ + + np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*(1-d[:,1])*(1-d[:,2]))[0],1,0) \ + + np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*d[:,1]*(1-d[:,2]))[0],1,1) \ + + np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*(1-d[:,1])*d[:,2])[0],1,2) \ + + np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*d[:,1]*(1-d[:,2]))[0],1,0),1,1) \ + + np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*(1-d[:,1])*d[:,2])[0],1,0),1,2) \ + + np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*d[:,1]*d[:,2])[0],1,1),1,2) \ + + np.roll(np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*d[:,1]*d[:,2])[0],1,0),1,1),1,2) + + rho /= wm + + rho = rho/rho.mean() - 1. + + return (rho, wm, ws) + + +def powcor( d1, d2, Lbox, Nmesh = 128, Nbin = 100, scale = 'lin', cor = False ): + + # Fourier transform + d1 = np.fft.fftn(d1) + d2 = np.fft.fftn(d2) + + # CIC correction + wid = np.indices(np.shape(d1)) - Nmesh/2 + #wid[np.where(wid >= Nmesh/2)] -= Nmesh + wid = wid*np.pi/Nmesh + 1e-100 + wcic = np.prod(np.sin(wid)/wid,0)**2 + + # Shell average power spectrum + dk = 2*np.pi/Lbox + Pk = np.conj(d1)*d2*(Lbox/Nmesh**2)**3 + Pk = np.fft.fftshift(Pk)/wcic**2 + + (Nm, km, Pkm, SPkm) = shellavg(np.real(Pk), dk, Nmesh, Nbin = Nbin, xmin = 0., xmax = Nmesh*dk/2, scale = scale) + + # Inverse Fourier transform and shell average correlation function + if cor == True: + dx = Lbox/Nmesh + Xr = np.fft.ifftshift(np.fft.ifftn(np.fft.ifftshift(Pk)))*(Nmesh/Lbox)**3 + + (Nmx, rm, Xrm, SXrm) = shellavg(np.real(Xr), dx, Nmesh, Nbin = Nbin, xmin = dx, xmax = 140., scale = scale) + + return ((Nm, km, Pkm, SPkm),(Nmx, rm, Xrm, SXrm)) + + else: return (Nm, km, Pkm, SPkm) + + +def shellavg( f, dx, Nmesh, Nbin = 100, xmin = 0., xmax = 1., scale = 'lin' ): + + x = np.indices(np.shape(f)) - Nmesh/2 + #x[np.where(x >= Nmesh/2)] -= Nmesh + x = dx*np.sqrt(np.sum(x**2,0)) + if scale == 'lin': bins = xmin+(xmax-xmin)* (np.arange(Nbin+1)/float(Nbin)) + if scale == 'log': bins = xmin*(xmax/xmin)**(np.arange(Nbin+1)/float(Nbin)) + + Nm = np.histogram(x, bins = bins)[0] + xm = np.histogram(x, bins = bins, weights = x)[0]/Nm + fm = np.histogram(x, bins = bins, weights = f)[0]/Nm + fs = np.sqrt((np.histogram(x, bins = bins, weights = f**2)[0]/Nm - fm**2)/(Nm-1)) + + return (Nm, xm, fm, fs) \ No newline at end of file