removed extra cross-comparison files

This commit is contained in:
P.M. Sutter 2014-05-06 17:22:26 -05:00
parent b441758a64
commit 1218ba3bdd
9 changed files with 0 additions and 1586 deletions

View file

@ -1,176 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/analysis/makeCocenterProfiles.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.
#+
#!/usr/bin/env python
#+
# VIDE -- Void IDentification and Examination -- ./pipeline/apAnalysis.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.
#+
# takes voids of given radii, computes 1 profiles,
# then computes 1d profiles for higher-resolution catalogs using
# same positions
# computes radial density profiles centered on baseSample
import imp
import pickle
import os
import numpy as np
import argparse
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from void_python_tools.backend import *
from util import *
# ------------------------------------------------------------------------------
parser = argparse.ArgumentParser(description='Analyze.')
parser.add_argument('--parm', dest='parm', default='datasetsToAnalyze.py', help='path to parameter file')
parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False,
help='display the plot (default: just write eps)')
args = parser.parse_args()
# -----------------------------------------------------------------------------
# plot a slice of the density around the void in baseIDList,
# with any voids in the slice shown and any voids in baseIDList flagged
def saveProfiles(baseSample, stack, sampleList, profileList, radii,
figDir, showPlot, outputDir):
thisRadius = str(stack.rMin) + "-" + str(stack.rMax)
plotName = "1dprofile_cocenter_" + baseSample.fullName+"_"+thisRadius
np.savez(outputDir+"/1dprofile_cocentered_"+plotName+".dat",
profileList, radii)
return
# -----------------------------------------------------------------------------
filename = args.parm
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(outputDir, os.F_OK):
os.makedirs(outputDir)
if not os.access(logDir, os.F_OK):
os.makedirs(logDir)
if not os.access(figDir, os.F_OK):
os.makedirs(figDir)
# get list of base voids
with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
baseSample = pickle.load(input)
baseSampleName = baseSample.fullName
baseVoidList = np.loadtxt(workDir+baseSampleDir+"/centers_central_"+\
baseSampleName+".out")
sampleList = []
for sampleDir in sampleDirList:
if compareSampleTag in sampleDir: continue
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
sampleList.append(pickle.load(input))
sampleDirList.insert(0,baseSampleDir)
sampleList.insert(0,baseSample)
# pick our void sample
for stack in baseSample.stacks:
print " Stack:", stack.rMin, "-", stack.rMax
accepted = (baseVoidList[:,4] > stack.rMin) & (baseVoidList[:,4] < stack.rMax)
stackVoidList = baseVoidList[accepted]
print " We have", len(stackVoidList), "voids here"
profileList = []
radii = []
rMaxProfile = stack.rMin*3 + 2
if baseSample.profileBinSize == "auto":
density = 0.5 * 50 / rMaxProfile / 2
else:
density = baseSample.profileBinSize
nBins = rMaxProfile*density
for (iSample, sampleDir) in enumerate(sampleDirList):
if compareSampleTag in sampleDir: continue
sample = sampleList[iSample]
print " Working with", sample.fullName, "..."
sys.stdout.flush()
sampleName = sample.fullName
print " Loading particle data..."
partData, boxLen, volNorm = loadPart(workDir, sampleDir, sample)
stackedProfile = np.zeros((nBins))
print " Stacking voids..."
binCenters = []
for void in stackVoidList:
periodicLine = getPeriodic(sample)
center = void[0:3]
shiftedPart = shiftPart(partData, center, periodicLine, boxLen)
dist = np.sqrt(shiftedPart[:,0]**2 + shiftedPart[:,1]**2 + \
shiftedPart[:,2]**2)
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
stackedProfile += thisProfile
binCenters = 0.5*(radii[1:]+radii[:-1])
stackedProfile /= volNorm
stackedProfile /= len(stackVoidList)
profileList.append(stackedProfile)
# plot these profiles
print " Plotting..."
sys.stdout.flush()
#binCenters = 0.5*(radii[1:] + radii[:-1])
saveProfiles(baseSample, stack, sampleList, profileList, binCenters,
figDir, args.showPlot, outputDir)
print " Done!"

View file

@ -1,103 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/analysis/overlapVoids.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.
#+
#!/usr/bin/env python
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/analysis/mergerTree.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.
#+
# computes the overlap between two void catalogs
from void_python_tools.backend import *
from void_python_tools.plotting import *
import imp
import pickle
import os
import matplotlib.pyplot as plt
import numpy as np
import argparse
# ------------------------------------------------------------------------------
parser = argparse.ArgumentParser(description='Analyze.')
parser.add_argument('--parm', dest='parm', default='datasetsToAnalyze.py',
help='path to parameter file')
args = parser.parse_args()
# ------------------------------------------------------------------------------
filename = args.parm
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(outputDir, os.F_OK):
os.makedirs(outputDir)
if not os.access(logDir, os.F_OK):
os.makedirs(logDir)
outFileName = outputDir + "/" + "voidOverlap" #+ ".dat"
with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
baseSample = pickle.load(input)
for (iSample, sampleDir) in enumerate(sampleDirList):
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
sample = pickle.load(input)
print " Working with", sample.fullName, "...",
sys.stdout.flush()
sampleName = sample.fullName
binPath = CTOOLS_PATH+"/analysis/voidOverlap"
logFile = logDir+"/mergertree_"+baseSample.fullName+"_"+sampleName+".out"
stepOutputFileName = outFileName + "_" + baseSample.fullName + "_" + \
sampleName + "_"
launchVoidOverlap(baseSample, sample, workDir+baseSampleDir,
workDir+sampleDir, binPath,
thisDataPortion="central", logFile=logFile,
continueRun=False, workDir=workDir,
outputFile=stepOutputFileName,
#matchMethod="useID")
matchMethod="prox")
print " Done!"

View file

@ -1,51 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/parm/sampleParm.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.
#+
#!/usr/bin/env python
workDir = "" # base directory for all samples
outputDir = ""
logDir = "./logs/"
figDir = "./figs/"
# path to c_tools directory in VIDE
CTOOLS_PATH = "/home/psutter2/projects/Voids/vide/c_tools/"
# the path under workDir/ which which holds the sample you want to compare againt (e.g., the fiducial case)
baseSampleDir = "sample_/"
# comma-separated list of samples to compare against baseSampleDir
sampleDirList = [
"sample1_/",
"sample2_/",
"sample3_/",
]
dataPortions = [ "all" ]
# this name gets appended to all output filenames and plots
plotLabel = "test"
# title to place on plots
plotTitle = "Test"
# don't touch this for now; it will be fully implemented and explained later
compareSampleTag = ""
doTheory = False

View file

@ -1,37 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/plotting/datasetsToPlot.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.
#+
#!/usr/bin/env python
workDir = "/home/psutter2/workspace/Voids/"
figDir = "./figs"
sampleDirList = [
# "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/",
# "multidark/md_halos_min1.393e12_pv/sample_md_halos_min1.393e12_pv_z0.56_d00/",
# "random/ran_ss0.000175/sample_ran_ss0.000175_z0.56_d00/",
# "random/ran_ss0.1/sample_ran_ss0.1_z0.56_d00/",
# "multidark/md_hod_dr9mid_pv/sample_md_hod_dr9mid_pv_z0.56_d00/",
# "multidark/md_ss0.000175_pv/sample_md_ss0.000175_pv_z0.56_d00/",
"sdss_dr9/sample_lss.dr9cmassmid.dat/",
"lanl/masked/masked_lanl_hod_dr9mid_pv/sample_masked_lanl_hod_dr9mid_pv_z0.5/" ]
dataPortion = "central"

View file

@ -1,193 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/plotting/plotCocenterProfiles.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.
#+
#!/usr/bin/env python
#+
# VIDE -- Void IDentification and Examination -- ./pipeline/apAnalysis.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.
#+
# plots radial density profiles centered on baseVoid.
# requires makeCocenteredProfiles to be run first!
import imp
import pickle
import os
import numpy as np
import argparse
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from void_python_tools.backend import *
from util import *
from globalOptions import *
from scipy.optimize import curve_fit
matplotlib.rcParams.update({'font.size': 20})
# ------------------------------------------------------------------------------
parser = argparse.ArgumentParser(description='Analyze.')
parser.add_argument('--parm', dest='parm', default='datasetsToAnalyze.py', help='path to parameter file')
parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False,
help='display the plot (default: just write eps)')
args = parser.parse_args()
# -----------------------------------------------------------------------------
# Lavaux & Wandelt (2012) profile
def LWProfile(r, A0, A3, alpha):
return A0 + A3*r**3.
#return A0 + A3*r**alpha
# -----------------------------------------------------------------------------
# http://arxiv.org/pdf/astro-ph/0508297v1.pdf eq. 5
def PadillaProfile(r, A1, A2, alpha):
return 1.5-A1*np.exp(-(A2*r)**alpha)
# -----------------------------------------------------------------------------
# plot a slice of the density around the void in baseIDList,
# with any voids in the slice shown and any voids in baseIDList flagged
def plotProfiles(baseSample, stack, sampleList,
figDir, showPlot, outputDir, doTheory):
thisRadius = str(stack.rMin) + "-" + str(stack.rMax)
plotName = "1dprofile_cocenter_" + baseSample.fullName+"_"+thisRadius
filename = "1dprofile_cocenter_" + baseSample.fullName+"_"+thisRadius
npzfile = np.load(outputDir+"/1dprofile_cocentered_"+plotName+".dat.npz")
profileList = npzfile['arr_0']
radii = npzfile['arr_1']
plt.clf()
plt.xlabel(r"$R/R_{eff}$")
#plt.xlabel(r"$R/R_{v,\mathrm{max}}$")
plt.ylabel(r"$n / \bar n$")
plt.xlim(xmin=0.0, xmax=2.5)
plt.ylim(ymin=0.0, ymax=1.7)
#plt.xscale('log')
for (iSample, sample) in enumerate(sampleList):
lineTitle = sample.nickName[:-10]
if "DM LowDen" in lineTitle or "DM HighDen" in lineTitle: continue
thisPlotTitle = r"$R_{eff}$ = "+thisRadius+ r" $h^{-1}$Mpc"
legendTitle = "Fixed Center Samples"
thisProfile = profileList[iSample]
if np.all(thisProfile == 0.): continue
rV = (stack.rMin + stack.rMax)/2.
if len(radii) > 0:
scaledRadii = radii/rV
plt.plot(scaledRadii, thisProfile, label=lineTitle,
color=colorList[iSample],
linewidth=linewidth)
if doTheory:
nBins = len(scaledRadii)/2
try:
popt, pcov = curve_fit(PadillaProfile, scaledRadii[0:nBins],
thisProfile[0:nBins], maxfev=10000, xtol=5.e-3)
except RuntimeError:
print "Warning: no convergence reached"
label = r"A1=%.2f, A2=%.2f, $\alpha$=%.2f" % (popt[0], popt[1], popt[2])
rho = PadillaProfile(scaledRadii, popt[0], popt[1], popt[2])
plt.plot(scaledRadii, rho, '--', label=label,
color=colorList[iSample], linewidth=2)
plt.title(thisPlotTitle + " (center from " + plotTitle + ")", fontsize=18)
plt.legend(title = legendTitle, loc = "lower right", prop={'size':16})
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")
return
# -----------------------------------------------------------------------------
filename = args.parm
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(outputDir, os.F_OK):
os.makedirs(outputDir)
if not os.access(logDir, os.F_OK):
os.makedirs(logDir)
if not os.access(figDir, os.F_OK):
os.makedirs(figDir)
# get list of base voids
with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
baseSample = pickle.load(input)
baseSampleName = baseSample.fullName
baseVoidList = np.loadtxt(workDir+baseSampleDir+"/centers_central_"+\
baseSampleName+".out")
sampleList = []
for sampleDir in sampleDirList:
if compareSampleTag in sampleDir: continue
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
sampleList.append(pickle.load(input))
sampleDirList.insert(0,baseSampleDir)
sampleList.insert(0,baseSample)
# pick our void sample
for stack in baseSample.stacks:
print " Stack:", stack.rMin, "-", stack.rMax
# plot these profiles
print " Plotting..."
sys.stdout.flush()
plotProfiles(baseSample, stack, sampleList,
figDir, args.showPlot, outputDir, doTheory)
print " Done!"

View file

@ -1,404 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/plotting/plotDenMaps.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.
#+
#!/usr/bin/env python
#+
# VIDE -- Void IDentification and Examination -- ./pipeline/apAnalysis.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.
#+
# takes nVoids evenly distributed, plots a slice of the local density and
# overlays the voids
import imp
import pickle
import os
import numpy as np
import argparse
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from void_python_tools.backend import *
import void_python_tools.xcor as xcor
from netCDF4 import Dataset
#import pylab as plt
NetCDFFile = Dataset
ncFloat = 'f8'
matplotlib.rcParams.update({'font.size': 16})
# ------------------------------------------------------------------------------
mergerNameBase = "voidOverlap"
parser = argparse.ArgumentParser(description='Analyze.')
parser.add_argument('--parm', dest='parm', default='datasetsToAnalyze.py', help='path to parameter file')
parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False,
help='display the plot (default: just write eps)')
args = parser.parse_args()
nVoids = 10
# ------------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# plot a slice of the density around the void in baseIDList,
# with any voids in the slice shown and any voids in baseIDList flagged
def plotVoidAndDen(idList, voidList, partData, boxLen, figDir,
sliceCenter=None, sliceWidth=200,
baseIDList=None, baseRadius=0, nickName=None,
baseNickName=None,
baseIndex=0,
periodic=None, showPlot=False, plotName=None):
if len(voidList) <= 1: return
plt.clf()
#sliceWidth = 220
sliceWidth = max(220, sliceWidth)
# make an appropriate box
xwidth = sliceWidth
ywidth = sliceWidth
zwidth = sliceWidth/4.
#zwidth = max(sliceWidth/4., 50)
# get mean density
part = 1.*partData
totalNumPart = len(part)
totalVol = (part[:,0].max() - part[:,0].min()) * \
(part[:,1].max() - part[:,1].min()) * \
(part[:,2].max() - part[:,2].min())
meanDen = totalNumPart/totalVol
# single out the matched void
keepVoid = []
if len(np.atleast_1d(idList)) > 0:
keepVoid = voidList[voidList[:,7] == idList]
if len(np.shape(keepVoid)) > 1: keepVoid = keepVoid[0,:]
filter = voidList[:,7] != idList
voidList = voidList[filter,:]
# convert everything to relative coordinates
part[:,0] -= sliceCenter[0]
part[:,1] -= sliceCenter[1]
part[:,2] -= sliceCenter[2]
shiftUs = np.abs(part[:,0]) > boxLen[0]/2.
if ("x" in periodicLine): part[shiftUs,0] -= \
np.copysign(boxLen[0],part[shiftUs,0])
shiftUs = np.abs(part[:,1]) > boxLen[1]/2.
if ("y" in periodicLine): part[shiftUs,1] -= \
np.copysign(boxLen[1],part[shiftUs,1])
shiftUs = np.abs(part[:,2]) > boxLen[2]/2.
if ("z" in periodicLine): part[shiftUs,2] -= \
np.copysign(boxLen[2],part[shiftUs,2])
voidList = np.atleast_2d(voidList)
np.atleast_2d(voidList)[:,0] -= sliceCenter[0]
np.atleast_2d(voidList)[:,1] -= sliceCenter[1]
np.atleast_2d(voidList)[:,2] -= sliceCenter[2]
shiftUs = np.abs(voidList[:,0]) > boxLen[0]/2.
if ("x" in periodicLine):
voidList[shiftUs,0] -= \
np.copysign(boxLen[0],voidList[shiftUs,0])
shiftUs = np.abs(voidList[:,1]) > boxLen[1]/2.
if ("y" in periodicLine): voidList[shiftUs,1] -= \
np.copysign(boxLen[1],voidList[shiftUs,1])
shiftUs = np.abs(voidList[:,2]) > boxLen[2]/2.
if ("z" in periodicLine): voidList[shiftUs,2] -= \
np.copysign(boxLen[2],voidList[shiftUs,2])
if len(np.atleast_1d(keepVoid)) >= 1:
keepVoid[0] -= sliceCenter[0]
keepVoid[1] -= sliceCenter[1]
keepVoid[2] -= sliceCenter[2]
shiftUs = np.abs(keepVoid[0]) > boxLen[0]/2.
if ("x" in periodicLine) and shiftUs: keepVoid[0] -= \
np.copysign(boxLen[0],keepVoid[0])
shiftUs = np.abs(keepVoid[1]) > boxLen[1]/2.
if ("y" in periodicLine) and shiftUs: keepVoid[1] -= \
np.copysign(boxLen[1],keepVoid[1])
shiftUs = np.abs(keepVoid[2]) > boxLen[2]/2.
if ("z" in periodicLine) and shiftUs: keepVoid[2] -= \
np.copysign(boxLen[2],keepVoid[2])
xmin = -xwidth/2.
xmax = xwidth/2.
ymin = -ywidth/2.
ymax = ywidth/2.
zmin = -zwidth/2.
zmax = zwidth/2.
# pull out voids that were potential matches
filter = np.sqrt(voidList[:,0]**2 + voidList[:,1]**2 + voidList[:,2]**2) <=\
baseRadius*1.5
potentialMatches = voidList[filter]
# get centers and radii of any other voids in slice
zminVoid = -zwidth/16.
zmaxVoid = zwidth/16.
filter = (voidList[:,0] > xmin) & (voidList[:,0] < xmax) & \
(voidList[:,1] > ymin) & (voidList[:,1] < ymax) & \
(voidList[:,2] > zminVoid) & (voidList[:,2] < zmaxVoid)
voidList = voidList[filter,:]
# slice particles
filter = (part[:,0] > xmin) & (part[:,0] < xmax) & \
(part[:,1] > ymin) & (part[:,1] < ymax) & \
(part[:,2] > zmin) & (part[:,2] < zmax)
part = part[filter]
# plot density
extent = [xmin, xmax, ymin, ymax]
hist, xedges, yedges = np.histogram2d(part[:,0], part[:,1], normed=False,
bins=64)
#hist /= meanDen
hist = np.log10(hist+1)
plt.imshow(hist,
aspect='equal',
extent=extent,
interpolation='gaussian',
cmap='YlGnBu_r')
#plt.colorbar()
# overlay voids as circles
fig = plt.gcf()
ax = fig.add_subplot(1,1,1)
# the original void
circle = plt.Circle((0,0), baseRadius,
edgecolor='orange', facecolor=None, fill=False,
linewidth=5)
fig.gca().add_artist(circle)
# our matched void
if len(np.atleast_1d(keepVoid)) > 0:
if idList == baseIDList:
edgecolor = 'orange'
else:
edgecolor = 'red'
circle = plt.Circle((keepVoid[0], keepVoid[1]), keepVoid[4],
edgecolor=edgecolor, facecolor=None, fill=False,
linewidth=5)
fig.gca().add_artist(circle)
# other voids in the slice
for void in voidList:
if np.any(void[7] == idList):
continue
else:
color = 'white'
circle = plt.Circle((void[0], void[1]), void[4],
edgecolor=color, facecolor=None, fill=False,
linewidth=5)
fig.gca().add_artist(circle)
# potential match voids
#for void in potentialMatches:
# if np.any(void[7] == idList):
# continue
# else:
# color = 'green'
# circle = plt.Circle((void[0], void[1]), void[4],
# edgecolor=color, facecolor=None, fill=False,
# linewidth=5)
# fig.gca().add_artist(circle)
baseNickName = baseNickName[:-10].lstrip()
nickName = nickName[:-10].lstrip()
if idList == baseIDList:
title = "%d $h^{-1}$Mpc" % int(baseRadius)
title += " (" + baseNickName + ")"
else:
title = r"$\rightarrow$ "
if len(np.atleast_1d(keepVoid)) > 0:
title += "%d $h^{-1}$Mpc" % int(keepVoid[4]) + " (" + nickName + ")"
else:
title += "No match"
#title += "(" + str(int(baseIDList)) + ")"
plt.title(title, fontsize=20)
#plt.xlabel("x [$h^{-1}$Mpc]", fontsize=14)
#plt.ylabel("y [$h^{-1}$Mpc]", fontsize=14)
plotName += "_" + str(int(baseIndex))
#plotName += "_" + str(int(baseRadius))
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")
return
# -----------------------------------------------------------------------------
filename = args.parm
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(outputDir, os.F_OK):
os.makedirs(outputDir)
if not os.access(logDir, os.F_OK):
os.makedirs(logDir)
if not os.access(figDir, os.F_OK):
os.makedirs(figDir)
mergerFileBase = outputDir + "/" + mergerNameBase
# get list of base voids
with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
baseSample = pickle.load(input)
baseSampleName = baseSample.fullName
baseVoidList = np.loadtxt(workDir+baseSampleDir+"/centers_central_"+\
baseSampleName+".out")
# sort by size
radii = baseVoidList[:,4]
indices = np.argsort(radii)[::-1]
baseVoidList = baseVoidList[indices,:]
setName = baseSampleDir.split('/')[0]
# pick our void sample
bigVoidList = baseVoidList[0:10,:]
stride = len(baseVoidList)/nVoids
baseVoidList = baseVoidList[::stride]
baseVoidList = np.vstack((bigVoidList,baseVoidList))
sampleDirList.insert(0,baseSampleDir)
for (iSample, sampleDir) in enumerate(sampleDirList):
if compareSampleTag in sampleDir: continue
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
sample = pickle.load(input)
print " Working with", sample.fullName, "..."
sys.stdout.flush()
sampleName = sample.fullName
print " Loading particle data..."
sys.stdout.flush()
infoFile = workDir+"/"+sampleDir+"/zobov_slice_"+sample.fullName+".par"
File = NetCDFFile(infoFile, 'r')
ranges = np.zeros((3,2))
ranges[0][0] = getattr(File, 'range_x_min')
ranges[0][1] = getattr(File, 'range_x_max')
ranges[1][0] = getattr(File, 'range_y_min')
ranges[1][1] = getattr(File, 'range_y_max')
ranges[2][0] = getattr(File, 'range_z_min')
ranges[2][1] = getattr(File, 'range_z_max')
File.close()
mul = np.zeros((3))
mul[:] = ranges[:,1] - ranges[:,0]
boxLen = mul
partFile = workDir+"/"+sampleDir+"/zobov_slice_"+sample.fullName
#partFile = catalogDir+"/"+sample.dataFile
iLine = 0
partData = []
part = np.zeros((3))
File = file(partFile)
chk = np.fromfile(File, dtype=np.int32,count=1)
Np = np.fromfile(File, dtype=np.int32,count=1)
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)
x *= mul[0]
x += ranges[0][0]
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)
y *= mul[1]
y += ranges[1][0]
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)
z *= mul[2]
z += ranges[2][0]
chk = np.fromfile(File, dtype=np.int32,count=1)
File.close()
partData = np.column_stack((x,y,z))#.transpose()
for (iBaseVoid,baseVoid) in enumerate(baseVoidList):
print " Void:", int(baseVoid[7]), "(", int(baseVoid[4]), ")"
baseIDList = baseVoid[7]
sliceCenter = baseVoid[0:3]
sliceWidth = baseVoid[4]*4
# get matched void
idList = []
if sample.fullName == baseSample.fullName:
idList = baseIDList
else:
matchFile=mergerFileBase+"_"+baseSampleName+"_"+sampleName+"_summary.out"
if os.access(matchFile, os.F_OK):
matchList = np.loadtxt(matchFile)
for i,testID in enumerate(matchList[:,0]):
if testID == baseIDList:
if (matchList[i,8] > 0): idList.append(matchList[i,8])
idList = np.array(idList)
idList = idList.astype(int)
voidList = np.loadtxt(workDir+sampleDir+"/trimmed_nodencut_centers_central_"+\
sampleName+".out")
periodicLine = getPeriodic(sample)
plotVoidAndDen(idList, voidList, partData, boxLen, figDir,
sliceCenter=sliceCenter, sliceWidth=sliceWidth,
baseIDList=baseIDList, baseRadius=baseVoid[4],
baseIndex=iBaseVoid,
nickName=sample.nickName, periodic=periodicLine,
baseNickName=baseSample.nickName,
showPlot=args.showPlot,
plotName="denmap_"+setName+"_"+baseSampleName+"_"+sampleName)
print " Done!"

View file

@ -1,154 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/plotting/plotMatchEllipRatio.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.
#+
#!/usr/bin/env python
# plots cumulative distributions of number counts
import matplotlib
matplotlib.use('Agg')
from void_python_tools.backend import *
from void_python_tools.plotting import *
import void_python_tools.apTools as vp
import imp
import pickle
import os
import matplotlib.pyplot as plt
import numpy as np
import argparse
from globalOptions import *
# plots the ratio of ellipticities for matched voids
# ------------------------------------------------------------------------------
parser = argparse.ArgumentParser(description='Plot.')
parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False,
help='display the plot (default: just write eps)')
parser.add_argument('--parm', dest='parm', default='datasetsToPlot.py',
help='path to parameter file')
args = parser.parse_args()
# ------------------------------------------------------------------------------
print "Plotting ellipticity ratio"
filename = args.parm
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(figDir, os.F_OK):
os.makedirs(figDir)
dataSampleList = []
compareSampleList = []
with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
baseSample = pickle.load(input)
for sampleDir in sampleDirList:
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
thisSample = pickle.load(input)
if compareSampleTag in thisSample.fullName:
compareSampleList.append(thisSample)
else:
dataSampleList.append(thisSample)
plt.clf()
numSubPlots = len(dataSampleList)
fig, axesList = plt.subplots(numSubPlots, sharex=True, sharey=True)
axesList = np.atleast_1d(axesList)
for (iSample,sample) in enumerate(dataSampleList):
if sample.fullName == baseSample.fullName: continue
sampleName = sample.fullName
lineTitle = sample.nickName[:-10]
# plt.xlabel("Void Radius [Mpc/h]")
# plt.ylabel(r"1st Progenitor Relative Ellipticity")
#plt.yscale('log')
# plt.xlim(xmax=rMax)
plt.xlim(rMin, rMax)
plotNameBase = "matchrelellip"
plotName = plotNameBase + "_" + plotLabel# + "_" + sampleName
filename = outputDir+"/voidOverlap_"+baseSample.fullName+"_"+sampleName+"_summary.out"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
continue
data = np.loadtxt(filename, comments="#")
if data.ndim == 1:
print " Too few!"
continue
# find a sample to compare it to
for compareSample in compareSampleList:
if compareSample.nickName[:10] == sample.nickName[:10]:
filename = outputDir+"/voidOverlap_"+baseSample.fullName+"_"+compareSample.fullName+"_summary.out"
compareData = np.loadtxt(filename, comments="#")
axesList[iSample].scatter(compareData[:,1], compareData[:,10],
color='blue', alpha=alpha, s=pointsize)
#plt.scatter(data[:,1], data[:,10],
# label=lineTitle, color=colorList[iSample])
axesList[iSample].scatter(data[:,1], data[:,10],
label=lineTitle, color='red', alpha=alpha, s=pointsize)
axesList[iSample].legend(loc = "upper left", prop={'size':10})
plt.ylim(0., 4.0)
yticks = axesList[iSample].yaxis.get_major_ticks()
yticks[-1].label1.set_visible(False)
yticks[0].label1.set_visible(False)
#axesList[iSample].set_xlim([20,rMax])
fig.subplots_adjust(hspace=0)
plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
axesList[0].set_title(plotTitle, fontsize=14)
#axesList[0].set_title("Match ellipticity ratio - "+plotTitle)
fig.text(0.5, 0.04, r'$R_{eff}$ [$h^{-1}$Mpc]', ha='center', va='center', fontsize=14)
fig.text(0.06, 0.5, 'Match Ellipticity Ratio', ha='center', va='center', rotation='vertical', fontsize=14)
#plt.legend(title = "Samples", loc = "upper left", prop={'size':8})
#plt.title("Match ellipticity ratio - "+plotTitle+" - "+sample.nickName)
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 args.showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")

View file

@ -1,163 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/plotting/plotMatchSizeRatio.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.
#+
#!/usr/bin/env python
# plots cumulative distributions of number counts
import matplotlib
matplotlib.use('Agg')
from void_python_tools.backend import *
from void_python_tools.plotting import *
import void_python_tools.apTools as vp
import imp
import pickle
import os
import matplotlib.pyplot as plt
import numpy as np
import argparse
from globalOptions import *
from scipy.optimize import curve_fit
# plots the ratio of sizes for matched voids
# ------------------------------------------------------------------------------
parser = argparse.ArgumentParser(description='Plot.')
parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False,
help='display the plot (default: just write eps)')
parser.add_argument('--parm', dest='parm', default='datasetsToPlot.py',
help='path to parameter file')
args = parser.parse_args()
# -----------------------------------------------------------------------------
def Linear(x, a, b):
return a*x+b
# ------------------------------------------------------------------------------
print "Plotting match size ratio"
filename = args.parm
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(figDir, os.F_OK):
os.makedirs(figDir)
dataSampleList = []
compareSampleList = []
with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
baseSample = pickle.load(input)
for sampleDir in sampleDirList:
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
thisSample = pickle.load(input)
if compareSampleTag in thisSample.fullName:
compareSampleList.append(thisSample)
else:
dataSampleList.append(thisSample)
plt.clf()
#plt.yscale('log')
numSubPlots = len(dataSampleList)
fig, axesList = plt.subplots(numSubPlots, sharex=True, sharey=True)
axesList = np.atleast_1d(axesList)
for (iSample,sample) in enumerate(dataSampleList):
if sample.fullName == baseSample.fullName: continue
sampleName = sample.fullName
lineTitle = sample.nickName[:-10]
plotNameBase = "matchvolrelradius"
plotName = plotNameBase + "_" + plotLabel# + "_" + sampleName
filename = outputDir+"/voidOverlap_"+baseSample.fullName+"_"+sampleName+"_summary.out"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
continue
data = np.loadtxt(filename, comments="#")
if data.ndim == 1:
print " Too few!"
continue
# find a sample to compare it to
for compareSample in compareSampleList:
if compareSample.nickName[:10] == sample.nickName[:10]:
filename = outputDir+"/voidOverlap_"+baseSample.fullName+"_"+compareSample.fullName+"_summary.out"
compareData = np.loadtxt(filename, comments="#")
axesList[iSample].scatter(compareData[:,1], compareData[:,2],
color='blue', alpha=alpha, s=pointsize)
axesList[iSample].scatter(data[:,1], data[:,2],
label=lineTitle, color='red', alpha=alpha, s=pointsize)
try:
popt, pcov = curve_fit(Linear, data[:,1],
data[:,2], maxfev=10000, xtol=5.e-3)
except RuntimeError:
print "Warning: no convergence reached"
label = r"$a$=%.2f, $b$=%.2f" % (popt[0], popt[1])
radii = np.arange(rMin, rMax, 1)
relRad = Linear(radii, popt[0], popt[1])
#axesList[iSample].plot(radii, relRad, '-', label=label,
# color='black', linewidth=2)
axesList[iSample].legend(loc = "best", fancybox=True,prop={'size':10})
plt.ylim(0., 2.5)
yticks = axesList[iSample].yaxis.get_major_ticks()
yticks[-1].label1.set_visible(False)
yticks[0].label1.set_visible(False)
axesList[iSample].set_xlim([rMin,rMax])
fig.subplots_adjust(hspace=0)
plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
axesList[0].set_title(plotTitle, fontsize=14)
#axesList[0].set_title("Match size ratio - "+plotTitle)
fig.text(0.5, 0.04, r'$R_{eff}$ [$h^{-1}$Mpc]', ha='center', va='center', fontsize=14)
fig.text(0.06, 0.5, 'Match Relative Radius', ha='center', va='center', rotation='vertical', fontsize=14)
#plt.legend(title = "Samples", loc = "upper left", prop={'size':8})
#plt.title("Match size ratio - "+plotTitle+" - "+sample.nickName)
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 args.showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")

View file

@ -1,305 +0,0 @@
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/plotting/plotNumberFunc.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.
#+
#!/usr/bin/env python
# plots cumulative distributions of number counts
import matplotlib
matplotlib.use('Agg')
from void_python_tools.backend import *
from void_python_tools.plotting import *
import void_python_tools.apTools as vp
import imp
import pickle
import os
import matplotlib.pyplot as plt
import numpy as np
import argparse
import svdw
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d
from globalOptions import *
# ------------------------------------------------------------------------------
histBinWidth = 1 # Mpc
parser = argparse.ArgumentParser(description='Plot.')
parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False,
help='display the plot (default: just write eps)')
parser.add_argument('--binned', dest='binned', action='store_const',
const=True, default=False,
help='plot binned function (default: cumulative)')
parser.add_argument('--parm', dest='parm', default='datasetsToPlot.py',
help='path to parameter file')
parser.add_argument('--xmax', dest='xmax', default=120.,
help='x limit of plot')
parser.add_argument('--xmin', dest='xmin', default=20.,
help='x limit of plot')
args = parser.parse_args()
# ------------------------------------------------------------------------------
def svdwFunc(r, scaleFactor):
radius, cumu_ps = svdw.getSvdW(.01, 100, 100, scaleFactor=scaleFactor)
cumu_ps += 0.1
cumu_ps = np.log10(cumu_ps)
interped = interp1d(radius, cumu_ps)
interpVal = interped(r)
interpVal[interpVal<1.0] = 1.0
#print "HELLO", r, interpVal
return interpVal
# ------------------------------------------------------------------------------
def loadData(sampleDir, dataPortion, treePortion='all'):
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
sample = pickle.load(input)
filename = workDir+"/"+sampleDir+"/centers_"+dataPortion+"_"+sample.fullName+".out"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
return -1, -1, -1
data = np.loadtxt(filename, comments="#")
if data.ndim == 1:
print " Too few!"
return -1, -1, -1
if treePortion == "parents":
filter = data[:,10] == -1
data = data[filter]
elif treePortion == "children":
filter = data[:,10] != -1
data = data[filter]
data = data[:,4]
indices = np.arange(0, len(data), 1)
sorted = np.sort(data)
if sample.dataType == "observation":
boxVol = vp.getSurveyProps(sample.maskFile,
sample.zBoundary[0], sample.zBoundary[1],
sample.zRange[0], sample.zRange[1], "all",
selectionFuncFile=None)[0]
#selectionFuncFile=sample.selFunFile)[0]
boxVol *= obsFudgeFactor
else:
boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] -
sample.zBoundaryMpc[0])
boxVol *= 1.e-9 # Mpc->Gpc
indices /= boxVol
#xmin = sorted[0]
#xmax = sorted[-1]
#bins = int((xmax-xmin)/histBinWidth)
bins = args.xmax/histBinWidth
hist, binEdges = np.histogram(sorted, bins=bins, range=(0., args.xmax))
#hist, binEdges = np.histogram(sorted, bins=bins, range=(xmin,xmax))
binCenters = 0.5*(binEdges[1:] + binEdges[:-1])
if not args.binned:
foundStart = False
for iBin in xrange(len(hist)):
if not foundStart and hist[iBin] == 0:
continue
foundStart = True
hist[iBin] = np.sum(hist[iBin:])
hist /= boxVol
hist = np.log10(hist)
lineTitle = sample.nickName[:-10]
return hist, binCenters, lineTitle
def fill_between(x, y1, y2=0, ax=None, **kwargs):
"""Plot filled region between `y1` and `y2`.
This function works exactly the same as matplotlib's fill_between, except
that it also plots a proxy artist (specifically, a rectangle of 0 size)
so that it can be added it appears on a legend.
"""
ax = ax if ax is not None else plt.gca()
ax.fill_between(x, y1, y2, interpolate=True, **kwargs)
p = plt.Rectangle((0, 0), 0, 0, **kwargs)
ax.add_patch(p)
# ------------------------------------------------------------------------------
print "Plotting number function"
filename = args.parm
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(figDir, os.F_OK):
os.makedirs(figDir)
plt.clf()
plt.xlabel(r"$R_{eff}$ [$h^{-1}$Mpc]", fontsize=14)
plt.ylabel(r"log ($n$ > $R_{eff}$ [$h^3$ Gpc$^{-3}$])", fontsize=14)
#plt.yscale('log')
plt.xlim(xmin=5.)
plt.xlim(xmax=100.)
plt.ylim(ymin=1)
plt.ylim(ymax=5)
plotNameBase = "numberfunc"
plotName = plotNameBase + "_" + plotLabel
sampleDirList.append(baseSampleDir)
for (iSample,sampleDir) in enumerate(sampleDirList):
for dataPortion in dataPortions:
# get all the data
allHist = []
if "ZZZZ" in sampleDir:
for fileZ in fileList:
thisSampleDir = sampleDir.replace("ZZZZ", fileZ)
hist, binCenters, lineTitle = loadData(thisSampleDir, dataPortion)
if lineTitle == -1: continue
allHist.append(hist)
lineLabel = lineTitle.replace(fileZ, "all")
if dataPortion != 'all': lineLabel += ", " + dataPortion
maxHist = 1.*allHist[-1]
minHist = 1.*allHist[-1]
for iHist in xrange(len(allHist)-1):
maxHist = np.maximum(maxHist, allHist[iHist])
minHist = np.minimum(minHist, allHist[iHist])
trim = (maxHist > 1)
minHist = minHist[trim]
maxHist = maxHist[trim]
binCentersToUse = binCenters[trim]
alpha = 0.75
if dataPortion == "central":
hatch = '//'
else:
hatch = None
fill_between(binCentersToUse, minHist, maxHist,
label=lineLabel, color=colorList[iSample],
alpha=alpha,
hatch=hatch
)
else:
#treeList = ["children", "parents", "all"]
treeList = ["all"]
for (iTree,treeItem) in enumerate(treeList):
hist, binCenters, lineLabel = loadData(sampleDir, dataPortion, treePortion=treeItem)
trim = (hist > 1)
hist = hist[trim]
binCentersToUse = binCenters[trim]
if lineLabel == -1: continue
if dataPortion != 'all': lineLabel += ", " + dataPortion
if treeItem != "all": lineLabel += ", " + treeItem
iColor = iSample + iTree
if dataPortion == "central":
lineStyle = '--'
else:
lineStyle = '-'
if "DM" in lineLabel: lineColor = colorList[0]
if "Halos" in lineLabel: lineColor = colorList[1]
if "HOD" in lineLabel: lineColor = colorList[2]
if "FullDen" in lineLabel:
lineColor = colorList[4]
lineStyle = '-'
linewidth = 5
if "HighDen" in lineLabel or "HighRes" in lineLabel or \
"All" in lineLabel:
lineStyle = "-"
linewidth = 5
if "LowDen" in lineLabel or "LowRes" in lineLabel or \
"1.20" in lineLabel:
lineStyle = "--"
linewidth = 5
plt.plot(binCentersToUse, hist, lineStyle,
label=lineLabel, color=lineColor,
linewidth=linewidth)
#if doTheory and "LowRes" in sampleDir:
# popt, pcov = curve_fit(svdwFunc, binCentersToUse, hist, p0=25.)
# radius, cumu_ps = svdw.getSvdW(.01, 100, 100, scaleFactor=popt[0])
# cumu_ps = np.log10(cumu_ps)
# print "HELLO", binCentersToUse, hist, cumu_ps
# plt.plot(radius, cumu_ps, color='black', label="SVdW/%g" % popt[0])
# and now the theoretical curve
if doTheory:
#radius, cumu_ps = svdw.getSvdW(.01, 100, 100)
#cumu_ps = np.log10(cumu_ps)
#plt.plot(radius, cumu_ps, color='black', label="SVdW" )
iLine = 0
scaleFactorList = []
#scaleFactorList = [10, 5]
for scaleFactor in scaleFactorList:
#radius, cumu_ps = svdw.getSvdW(.01, 100, 100, scaleFactor=scaleFactor)
#cumu_ps = np.log10(cumu_ps)
#plt.plot(radius, cumu_ps, lineList[iLine], color='black', label="SVdW/%g" % scaleFactor)
iLine += 1
dvList = [-.07]
iLine = 0
for dv in dvList:
radius, cumu_ps = svdw.getSvdW(.01, 100, 1000, d_v=dv)
cumu_ps = np.log10(cumu_ps)
plt.plot(radius, cumu_ps, lineList[iLine], color='black', label=r"SVdW $\delta_v$=%g" % dv, linewidth=3)
iLine += 1
dvList = [-.015]
for dv in dvList:
radius, cumu_ps = svdw.getSvdW(.01, 100, 1000, d_v=dv)
#radius, cumu_ps = svdw.getSvdW(.01, 100, 100, d_v=dv)
cumu_ps = np.log10(cumu_ps)
plt.plot(radius, cumu_ps, lineList[iLine], color='black', label="SVdW $\delta_v$=%g" % dv, linewidth=3)
iLine += 1
plt.legend(loc = "best", fancybox=True, prop={'size':12})
#plt.legend(title = "Samples", loc = "upper right", prop={'size':8})
#plt.title("Number func - "+plotTitle)
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 args.showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png")