integrated nico's cross-correlation backend and enabled plotting of the void distribution

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!"

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 = [

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 *

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")
 

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 *

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)