From 1218ba3bdd896ad4b7a02c4089b168ed11234d5a Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Tue, 6 May 2014 17:22:26 -0500 Subject: [PATCH] removed extra cross-comparison files --- crossCompare/analysis/makeCocenterProfiles.py | 176 -------- crossCompare/analysis/overlapVoids.py | 103 ----- crossCompare/parm/sampleParm.py | 51 --- crossCompare/plotting/datasetsToPlot.py | 37 -- crossCompare/plotting/plotCocenterProfiles.py | 193 --------- crossCompare/plotting/plotDenMaps.py | 404 ------------------ crossCompare/plotting/plotMatchEllipRatio.py | 154 ------- crossCompare/plotting/plotMatchSizeRatio.py | 163 ------- crossCompare/plotting/plotNumberFunc.py | 305 ------------- 9 files changed, 1586 deletions(-) delete mode 100644 crossCompare/analysis/makeCocenterProfiles.py delete mode 100644 crossCompare/analysis/overlapVoids.py delete mode 100644 crossCompare/parm/sampleParm.py delete mode 100644 crossCompare/plotting/datasetsToPlot.py delete mode 100644 crossCompare/plotting/plotCocenterProfiles.py delete mode 100644 crossCompare/plotting/plotDenMaps.py delete mode 100644 crossCompare/plotting/plotMatchEllipRatio.py delete mode 100644 crossCompare/plotting/plotMatchSizeRatio.py delete mode 100644 crossCompare/plotting/plotNumberFunc.py diff --git a/crossCompare/analysis/makeCocenterProfiles.py b/crossCompare/analysis/makeCocenterProfiles.py deleted file mode 100644 index 4274936..0000000 --- a/crossCompare/analysis/makeCocenterProfiles.py +++ /dev/null @@ -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!" - diff --git a/crossCompare/analysis/overlapVoids.py b/crossCompare/analysis/overlapVoids.py deleted file mode 100644 index 269e9b1..0000000 --- a/crossCompare/analysis/overlapVoids.py +++ /dev/null @@ -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!" diff --git a/crossCompare/parm/sampleParm.py b/crossCompare/parm/sampleParm.py deleted file mode 100644 index 7aa74b2..0000000 --- a/crossCompare/parm/sampleParm.py +++ /dev/null @@ -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 diff --git a/crossCompare/plotting/datasetsToPlot.py b/crossCompare/plotting/datasetsToPlot.py deleted file mode 100644 index c1ebf1a..0000000 --- a/crossCompare/plotting/datasetsToPlot.py +++ /dev/null @@ -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" - diff --git a/crossCompare/plotting/plotCocenterProfiles.py b/crossCompare/plotting/plotCocenterProfiles.py deleted file mode 100644 index a14e7b0..0000000 --- a/crossCompare/plotting/plotCocenterProfiles.py +++ /dev/null @@ -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!" - diff --git a/crossCompare/plotting/plotDenMaps.py b/crossCompare/plotting/plotDenMaps.py deleted file mode 100644 index a12c1a2..0000000 --- a/crossCompare/plotting/plotDenMaps.py +++ /dev/null @@ -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!" - diff --git a/crossCompare/plotting/plotMatchEllipRatio.py b/crossCompare/plotting/plotMatchEllipRatio.py deleted file mode 100644 index 7847809..0000000 --- a/crossCompare/plotting/plotMatchEllipRatio.py +++ /dev/null @@ -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") - - diff --git a/crossCompare/plotting/plotMatchSizeRatio.py b/crossCompare/plotting/plotMatchSizeRatio.py deleted file mode 100644 index 315fed1..0000000 --- a/crossCompare/plotting/plotMatchSizeRatio.py +++ /dev/null @@ -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") - - diff --git a/crossCompare/plotting/plotNumberFunc.py b/crossCompare/plotting/plotNumberFunc.py deleted file mode 100644 index 12c022c..0000000 --- a/crossCompare/plotting/plotNumberFunc.py +++ /dev/null @@ -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") - -