From b87a42650826286aba3f3f5bc8476612b41ff719 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 6 Sep 2013 20:02:46 -0500 Subject: [PATCH] added more tools and plotting for cross-comparison catalog analysis --- crossCompare/analysis/makeCocenterProfiles.py | 157 +++++++ crossCompare/analysis/overlapVoids.py | 81 ++++ crossCompare/parm/sampleParm.py | 32 ++ crossCompare/plotting/plotCocenterProfiles.py | 174 ++++++++ crossCompare/plotting/plotCompareDensCon.py | 121 ------ crossCompare/plotting/plotDenMaps.py | 385 ++++++++++++++++++ crossCompare/plotting/plotDensConVsR.py | 118 ------ crossCompare/plotting/plotMatchEllipRatio.py | 135 ++++++ crossCompare/plotting/plotMatchSizeRatio.py | 144 +++++++ crossCompare/plotting/plotNumberFunc.py | 144 +++++-- 10 files changed, 1224 insertions(+), 267 deletions(-) create mode 100755 crossCompare/analysis/makeCocenterProfiles.py create mode 100755 crossCompare/analysis/overlapVoids.py create mode 100755 crossCompare/parm/sampleParm.py create mode 100755 crossCompare/plotting/plotCocenterProfiles.py delete mode 100644 crossCompare/plotting/plotCompareDensCon.py create mode 100755 crossCompare/plotting/plotDenMaps.py delete mode 100644 crossCompare/plotting/plotDensConVsR.py create mode 100755 crossCompare/plotting/plotMatchEllipRatio.py create mode 100755 crossCompare/plotting/plotMatchSizeRatio.py diff --git a/crossCompare/analysis/makeCocenterProfiles.py b/crossCompare/analysis/makeCocenterProfiles.py new file mode 100755 index 0000000..53075c1 --- /dev/null +++ b/crossCompare/analysis/makeCocenterProfiles.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python +#+ +# VIDE -- Void IDEntification pipeline -- ./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!" + diff --git a/crossCompare/analysis/overlapVoids.py b/crossCompare/analysis/overlapVoids.py new file mode 100755 index 0000000..6fcd67d --- /dev/null +++ b/crossCompare/analysis/overlapVoids.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python +#+ +# VIDE -- Void IDEntification pipeline -- ./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) + +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!" diff --git a/crossCompare/parm/sampleParm.py b/crossCompare/parm/sampleParm.py new file mode 100755 index 0000000..e2381d0 --- /dev/null +++ b/crossCompare/parm/sampleParm.py @@ -0,0 +1,32 @@ +#!/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 diff --git a/crossCompare/plotting/plotCocenterProfiles.py b/crossCompare/plotting/plotCocenterProfiles.py new file mode 100755 index 0000000..8a2c0a8 --- /dev/null +++ b/crossCompare/plotting/plotCocenterProfiles.py @@ -0,0 +1,174 @@ +#!/usr/bin/env python +#+ +# VIDE -- Void IDEntification pipeline -- ./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!" + diff --git a/crossCompare/plotting/plotCompareDensCon.py b/crossCompare/plotting/plotCompareDensCon.py deleted file mode 100644 index 5e96dfb..0000000 --- a/crossCompare/plotting/plotCompareDensCon.py +++ /dev/null @@ -1,121 +0,0 @@ -#+ -# VIDE -- Void IDEntification pipeline -- ./crossCompare/plotting/plotCompareDensCon.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. -#+ -#!/usr/bin/env python - -# plots cumulative distributions of number counts versus density contrast - -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 - -# ------------------------------------------------------------------------------ - -plotNameBase = "compdenscon" - -obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing? - -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('--parmFile', dest='parmFile', default='datasetsToPlot.py', - help='path to parameter file') -args = parser.parse_args() - -# ------------------------------------------------------------------------------ - -filename = args.parmFile -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 = [] - -for sampleDir in sampleDirList: - with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input: - dataSampleList.append(pickle.load(input)) - -plt.clf() -plt.xlabel("Void Density Contrast") -plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]") -plt.yscale('log') -plt.xlim(xmax=5.) - -plotName = plotNameBase - -for (iSample,sample) in enumerate(dataSampleList): - - sampleName = sample.fullName - lineTitle = sampleName - - if sample.dataType == "observation": - boxVol = vp.getSurveyProps(sample.maskFile, - sample.zBoundary[0], sample.zBoundary[1], - sample.zRange[0], sample.zRange[1], "all", - selectionFuncFile=sample.selFunFile)[0] - boxVol *= obsFudgeFactor - else: - boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] - - sample.zBoundaryMpc[0]) - - boxVol *= 1.e-9 # Mpc->Gpc - - print " Loading", sampleName - filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\ - sampleName+".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 - data = data[:,8] - indices = np.arange(0, len(data), 1) - sorted = np.sort(data) - - plt.plot(sorted, indices[::-1]/boxVol, '-', - label=lineTitle, color=colorList[iSample], - linewidth=linewidth) - -plt.legend(title = "Samples", loc = "upper right", prop={'size':8}) -#plt.title(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") - - diff --git a/crossCompare/plotting/plotDenMaps.py b/crossCompare/plotting/plotDenMaps.py new file mode 100755 index 0000000..6eb2e11 --- /dev/null +++ b/crossCompare/plotting/plotDenMaps.py @@ -0,0 +1,385 @@ +#!/usr/bin/env python +#+ +# VIDE -- Void IDEntification pipeline -- ./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!" + diff --git a/crossCompare/plotting/plotDensConVsR.py b/crossCompare/plotting/plotDensConVsR.py deleted file mode 100644 index 87a44e8..0000000 --- a/crossCompare/plotting/plotDensConVsR.py +++ /dev/null @@ -1,118 +0,0 @@ -#+ -# VIDE -- Void IDEntification pipeline -- ./crossCompare/plotting/plotDensConVsR.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. -#+ -#!/usr/bin/env python - -# plots cumulative distributions of number counts versus density contrast - -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 - -# ------------------------------------------------------------------------------ - -plotNameBase = "densconvsr" - -obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing? - -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('--parmFile', dest='parmFile', default='datasetsToPlot.py', - help='path to parameter file') -args = parser.parse_args() - -# ------------------------------------------------------------------------------ - -filename = args.parmFile -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 = [] - -for sampleDir in sampleDirList: - with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input: - dataSampleList.append(pickle.load(input)) - -plt.clf() -plt.ylabel("Void Density Contrast") -plt.xlabel("Void Radius [Mpc/h]") -#plt.xlim(xmax=5.) -plt.yscale('log') - -plotName = plotNameBase - -for (iSample,sample) in enumerate(dataSampleList): - - sampleName = sample.fullName - lineTitle = sampleName - - if sample.dataType == "observation": - boxVol = vp.getSurveyProps(sample.maskFile, - sample.zBoundary[0], sample.zBoundary[1], - sample.zRange[0], sample.zRange[1], "all", - selectionFuncFile=sample.selFunFile)[0] - boxVol *= obsFudgeFactor - else: - boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] - - sample.zBoundaryMpc[0]) - - boxVol *= 1.e-9 # Mpc->Gpc - - print " Loading", sampleName - filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\ - sampleName+".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 - - plt.plot(data[:,4], data[:,8], '-', - label=lineTitle, color=colorList[iSample], - linewidth=linewidth) - -plt.legend(title = "Samples", loc = "upper right", prop={'size':8}) -#plt.title(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") - - diff --git a/crossCompare/plotting/plotMatchEllipRatio.py b/crossCompare/plotting/plotMatchEllipRatio.py new file mode 100755 index 0000000..6793127 --- /dev/null +++ b/crossCompare/plotting/plotMatchEllipRatio.py @@ -0,0 +1,135 @@ +#!/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") + + diff --git a/crossCompare/plotting/plotMatchSizeRatio.py b/crossCompare/plotting/plotMatchSizeRatio.py new file mode 100755 index 0000000..215b6d2 --- /dev/null +++ b/crossCompare/plotting/plotMatchSizeRatio.py @@ -0,0 +1,144 @@ +#!/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") + + diff --git a/crossCompare/plotting/plotNumberFunc.py b/crossCompare/plotting/plotNumberFunc.py index 1de232d..00ee6ec 100755 --- a/crossCompare/plotting/plotNumberFunc.py +++ b/crossCompare/plotting/plotNumberFunc.py @@ -2,6 +2,8 @@ # 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 @@ -11,12 +13,13 @@ 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 * # ------------------------------------------------------------------------------ -obsFudgeFactor = 1.0 # what fraction of the volume are we *reall* capturing? -#obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing? - histBinWidth = 1 # Mpc parser = argparse.ArgumentParser(description='Plot.') @@ -35,9 +38,20 @@ parser.add_argument('--xmin', dest='xmin', default=20., 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): +def loadData(sampleDir, dataPortion, treePortion='all'): with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input: sample = pickle.load(input) @@ -51,10 +65,18 @@ def loadData(sampleDir, dataPortion): 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], @@ -90,7 +112,7 @@ def loadData(sampleDir, dataPortion): hist = np.log10(hist) - lineTitle = sample.nickName + lineTitle = sample.nickName[:-10] return hist, binCenters, lineTitle @@ -122,17 +144,19 @@ if not os.access(figDir, os.F_OK): os.makedirs(figDir) plt.clf() -plt.xlabel("Void Radius [Mpc/h]") -plt.ylabel(r"log (N > R [$h^3$ Gpc$^{-3}$])") +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(xmax=args.xmax) -#plt.xlim(xmin=args.xmin) -#plt.ylim(ymin=-1) -#plt.ylim(ymax=6) +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 @@ -145,7 +169,7 @@ for (iSample,sampleDir) in enumerate(sampleDirList): allHist.append(hist) lineLabel = lineTitle.replace(fileZ, "all") - lineLabel += ", " + dataPortion + if dataPortion != 'all': lineLabel += ", " + dataPortion maxHist = 1.*allHist[-1] minHist = 1.*allHist[-1] @@ -170,23 +194,87 @@ for (iSample,sampleDir) in enumerate(sampleDirList): else: - hist, binCenters, lineLabel = loadData(sampleDir, dataPortion) - trim = (hist > 1) - hist = hist[trim] - binCentersToUse = binCenters[trim] - if lineLabel == -1: continue - lineLabel += ", " + dataPortion - if dataPortion == "central": - lineStyle = '--' - else: - lineStyle = '-' - plt.plot(binCentersToUse, hist, lineStyle, - label=lineLabel, color=colorList[iSample], - linewidth=linewidth) - + #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 -plt.legend(title = "Samples", loc = "upper right", prop={'size':8}) -plt.title("Number func - "+plotTitle) + 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")