From 022eec19bb82a582b15308c610aff83435fcaae0 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 29 Nov 2012 11:14:59 -0600 Subject: [PATCH] some bug fixes to the a-p analysis scripts; fixed major bug in prunevoids; updated catalog release scripts --- c_tools/mock/generateMock.cpp | 7 +- c_tools/stacking/pruneVoids.cpp | 5 +- pipeline/applyMaskToMock.py | 407 ++++++++++-------- pipeline/datasets/multidark.py | 7 +- pipeline/generateCatalog.py | 7 - pipeline/prepareCatalogs.py | 24 +- plotting/datasetsToPlot.py | 6 +- plotting/plotCompareDensCon.py | 17 +- plotting/plotCompareDist.py | 16 +- plotting/plotDensConVsR.py | 99 +++++ .../void_python_tools/backend/classes.py | 7 + .../void_python_tools/backend/launchers.py | 56 ++- .../void_python_tools/plotting/plotDefs.py | 18 +- 13 files changed, 437 insertions(+), 239 deletions(-) create mode 100755 plotting/plotDensConVsR.py diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index 93e7e72..8d4dff8 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -188,12 +188,13 @@ SimuData *doLoadMultidark(const char *multidarkname) &outd->Pos[0][i], &outd->Pos[1][i], &outd->Pos[2][i], &outd->Vel[2][i]); - outd->uniqueID[i] = 1.0; - //outd->uniqueID[i] = 1.0 * outd->Id[i]; + outd->uniqueID[i] = 1.0 * outd->Id[i]; + if (i < 10) printf("TEST %d %d\n", i, outd->Id[i]); if (outd->Id[i] == -99 && outd->Pos[0][i] == -99 && outd->Pos[1][i] == -99 && outd->Pos[2][i] == -99 && outd->Vel[2][i] == -99) { + printf("FOUND END\n"); break; } else { actualNumPart++; @@ -396,8 +397,6 @@ void generateOutput(SimuData *data, int axis, f.beginCheckpoint(); for (uint32_t i = 0; i < data->NumPart; i++) { - //printf("HELLO %d %d\n", i, data->Id[i]); - //f.writeReal32(data->Id[i]); f.writeReal32(data->uniqueID[i]); } f.endCheckpoint(); diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 86f850b..fbaf8d4 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -441,9 +441,6 @@ int main(int argc, char **argv) { } voids[iVoid].nearestEdge = nearestEdge; - //if (nearestEdge < voids[iVoid].nearestMock) { - // voids[iVoid].nearestMock = nearestEdge; - //} } // iVoid printf(" Picking winners and losers...\n"); @@ -453,7 +450,7 @@ int main(int argc, char **argv) { for (iVoid = 0; iVoid < numVoids; iVoid++) { if (voids[iVoid].densCon < 1.5) { - //voids[iVoid].accepted = 0; +// voids[iVoid].accepted = -4; } if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) { diff --git a/pipeline/applyMaskToMock.py b/pipeline/applyMaskToMock.py index c3fcd18..af7723c 100755 --- a/pipeline/applyMaskToMock.py +++ b/pipeline/applyMaskToMock.py @@ -1,8 +1,6 @@ #!/usr/bin/env python -# prepares input catalogs based on multidark simulations -# (borrows heavily from generateMock, but doesn't hold much in memory) -# also creates necessary analyzeVoids input files +# applies a mask to a given dataset import numpy as np import os @@ -32,6 +30,18 @@ parser.add_argument('--scripts', dest='script', action='store_const', parser.add_argument('--parmFile', dest='parmFile', default="", help='path to parameter file') +parser.add_argument('--subsamples', dest='subsample', action='store_const', + const=True, default=False, + help='write subsamples') +parser.add_argument('--halos', dest='halos', action='store_const', + const=True, default=False, + help='write halos') +parser.add_argument('--hod', dest='hod', action='store_const', + const=True, default=False, + help='write hod') +parser.add_argument('--all', dest='all', action='store_const', + const=True, default=False, + help='write everything') args = parser.parse_args() @@ -57,9 +67,9 @@ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): #------------------------------------------------------------------------------ # for given dataset parameters, outputs a script for use with analyzeVoids -def writeScript(setName, dataFileNameBase, - scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, - numSlices, useVel, lbox, minRadius, omegaM, subsample=1.0, +def writeObservationScript(setName, dataFileNameBase, maskFileName, + scriptDir, catalogDir, fileNums, redshifts, + useVel, minRadius, omegaM, suffix=".dat"): @@ -88,7 +98,7 @@ ranSeed = 101010 useLCDM = False bias = 1.16 -dataPortions = ["central"] +dataPortions = ["central", "all"] dataSampleList = [] """) @@ -112,25 +122,17 @@ numZobovThreads = {numZobovThreads} sampleInfo = """ newSample = Sample(dataFile = "{dataFile}", - dataFormat = "{dataFormat}", - dataUnit = {dataUnit}, fullName = "{sampleName}", nickName = "{sampleName}", - dataType = "simulation", + dataType = "observation", + maskFile = "{maskFileName}", zBoundary = ({zMin}, {zMax}), zRange = ({zMin}, {zMax}), - zBoundaryMpc = ({zMinMpc}, {zMaxMpc}), - omegaM = {omegaM}, minVoidRadius = {minRadius}, includeInHubble = True, partOfCombo = False, isCombo = False, - boxLen = {boxLen}, - usePecVel = {usePecVel}, - numSubvolumes = {numSubvolumes}, - mySubvolume = "{mySubvolume}", - useLightCone = {useLightCone}, - subsample = {subsample}) + usePecVel = {usePecVel}) dataSampleList.append(newSample) newSample.addStack({zMin}, {zMax}, {minRadius} , {minRadius}+2, True, False) newSample.addStack({zMin}, {zMax}, {minRadius} , {minRadius}+4, True, False) @@ -142,180 +144,249 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False) for (iFile, redshift) in enumerate(redshifts): fileNum = fileNums[iFile] zBox = float(redshift) - Om = float(omegaM) - zBoxMpc = LIGHT_SPEED/100.*vp.angularDiameter(zBox, Om=Om) - boxMaxMpc = zBoxMpc + lbox + + zMin = zBox + zMax = zBox + 0.1 + zMin, zMax = getMaskedZRange(lbox, zBox, zMin, zMax, omegaM) + dataFileName = dataFileNameBase + fileNum + suffix - # converter from redshift to comoving distance - zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000) - zVsDX = np.zeros(len(zVsDY)) - for i in xrange(len(zVsDY)): - zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om) + sampleName = getSampleName(setName, redshift, useVel, + iSlice=-1, iVol=-1) - if useLightCone: - boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \ - LIGHT_SPEED*lbox, zVsDX, zVsDY)-zBox - dzSafe = 0.03 - else: - boxWidthZ = lbox*100./LIGHT_SPEED - dzSafe = 0.0 - - for iSlice in xrange(numSlices): - sliceMin = zBox + dzSafe + iSlice*(boxWidthZ-dzSafe)/numSlices - sliceMax = zBox + dzSafe + (iSlice+1)*(boxWidthZ-dzSafe)/numSlices - - sliceMinMpc = sliceMin*LIGHT_SPEED/100. - sliceMaxMpc = sliceMax*LIGHT_SPEED/100. - - sliceMin = "%0.2f" % sliceMin - sliceMax = "%0.2f" % sliceMax - sliceMinMpc = "%0.1f" % sliceMinMpc - sliceMaxMpc = "%0.1f" % sliceMaxMpc - - dataFileName = dataFileNameBase + fileNum + suffix - - for iX in xrange(numSubvolumes): - for iY in xrange(numSubvolumes): - - mySubvolume = "%d%d" % (iX, iY) - - sampleName = getSampleName(setName, sliceMin, useVel, - iSlice=iSlice, iVol=mySubvolume) - - scriptFile.write(sampleInfo.format(dataFile=dataFileName, + scriptFile.write(sampleInfo.format(dataFile=dataFileName, dataFormat=dataFormat, - dataUnit=dataUnit, sampleName=sampleName, - zMin=sliceMin, - zMax=sliceMax, - zMinMpc=sliceMinMpc, - zMaxMpc=sliceMaxMpc, - omegaM=Om, - boxLen=lbox, + maskFileName=maskFileName, + zMin=zMin, + zMax=zMax, usePecVel=useVel, - minRadius=minRadius, - numSubvolumes=numSubvolumes, - mySubvolume=mySubvolume, - useLightCone=useLightCone, - subsample=subsample)) + minRadius=minRadius)) scriptFile.close() return #------------------------------------------------------------------------------ -#------------------------------------------------------------------------------ -if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir) -if not os.access(catalogDir, os.F_OK): os.mkdir(catalogDir) - -# ----------------------------------------------------------------------------- -# now the SDSS HOD -parFileText = """ -% cosmology -OMEGA_M {omegaM} -HUBBLE {hubble} -OMEGA_B 0.0469 -SIGMA_8 0.82 -SPECTRAL_INDX 0.95 -ITRANS 5 -REDSHIFT {redshift} - -% halo definition -%DELTA_HALO 200 -DELTA_HALO 740.74 % 200/Om_m -M_max 1.00E+16 - -% fit function types -pdfs 11 -pdfc 2 -EXCLUSION 4 - -% hod parameters -M_min {Mmin} -GALAXY_DENSITY 0.0111134 % computed automatically if M_min set, use for sanity -M1 {M1} -sigma_logM {sigma_logM} -alpha {alpha} -M_cut {Mcut} - -% simulation info -real_space_xi 1 -HOD 1 -populate_sim 1 -HaloFile {haloFile} -RESOLUTION {numPartPerSide} -BOX_SIZE {boxSize} - -% output -root_filename hod - """ - -print " Doing DR9 HOD" - - # these parameters come from Manera et al 2012, eq. 26 -parFileName = "./hod.par" -parFile = open(parFileName, 'w') -haloFile = catalogDir+haloFileBase+fileNums[iRedshift] -parFile.write(parFileText.format(omegaM=omegaM, - hubble=hubble, - redshift=redshift, - Mmin=1.23e13, - M1=1.e14, - sigma_logM=0.596, - alpha=1.0127, - Mcut=1.19399e13, - haloFile=haloFile, - numPartPerSide=numPart**(1/3.), - boxSize=lbox)) -parFile.close() - -os.system(hodPath+" "+parFileName+">& /dev/null") - # now place these particles on a lightcone, restrict redshift range, apply mask -mask = hp.read_map(maskFile) -nside = hp.get_nside(mask) +def applyMask(inFileName, outFileName, maskFileName, lbox, zBox, zMin, zMax, omegaM): + mask = hp.read_map(maskFileName) + nside = hp.get_nside(mask) -inFile = open('hod.mock', 'r') -outFile = open(catalogDir+"/mock.out")) + inFile = open(inFileName, 'r') + outFile = open(outFileName, 'w') -zBox = float(redshiftRange[0]) -Om = float(omegaM) + Om = float(omegaM) -# converter from redshift to comoving distance -zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000) -zVsDX = np.zeros(len(zVsDY)) -for i in xrange(len(zVsDY)): - zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om) - -for line in inFile: - line = line.split(',') - x = float(line[0]) - lbox/2. - y = float(line[1]) - lbox/2. - z = float(line[2]) - lbox/2. - vz = float(line[5]) + # converter from redshift to comoving distance + zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000) + zVsDX = np.zeros(len(zVsDY)) + for i in xrange(len(zVsDY)): + zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om) zBoxInMpc = vp.angularDiameter(zBox, Om=Om) - redshift = np.sqrt(x*x + y*y + z*z) - redshift = np.interp(zBoxInMpc+100./LIGHT_SPEED*redshift, zVsDX, zVsDY) + for (iLine,line) in enumerate(inFile): + if iLine < 5: + print >> outFile, line.rstrip() + continue - if redshift < redshiftRange[0] or redshift > redshiftRange[1]: continue + line = line.split(' ') + uniqueID = int(line[0]) + x = float(line[1]) - lbox/2. + y = float(line[2]) - lbox/2. + z = float(line[3]) - lbox/2. + vz = float(line[4]) - RA = np.atan((y-lbox/2.)/(x-lbox/2.)) * 100/np.pi + 180. - Dec = np.asin((z-lboc/2.)/(redshift*LIGHT_SPEED/100.)) * 180/np.pi + # TODO if usePecvel - phi = np.pi/180. * RA - theta = np.pi/2. - Dec*np.pi/180. - pos = np.zeros((3)) + redshift = np.sqrt(x*x + y*y + z*z) + redshift = np.interp(zBoxInMpc+100./LIGHT_SPEED*redshift, zVsDX, zVsDY) - pix = hp.ang2pix(nside, theta, phi) - if mask[pix] <= 0: continue + if redshift < zMin or redshift > zMax: continue + + vec = np.array((x,y,z)) + theta, phi = hp.vec2ang(vec) + theta = theta[0] + phi = phi[0] + RA = phi*180./np.pi + Dec = 180./np.pi*(np.pi/2.-theta) + + pix = hp.vec2pix(nside, x, y, z) + if mask[pix] <= 0.2: continue - print >> outFile, RA, Dec, redshift*LIGHT_SPEED, 0., x, y, z + print >> outFile, RA, Dec, redshift*LIGHT_SPEED, uniqueID, x, y, z -inFile.close() -outFile.close() + inFile.close() + outFile.close() -os.system("rm ./hod.*") +#------------------------------------------------------------------------------ +def getMaskedZRange(lbox, zBox, zMin, zMax, omegaM): + + if zMin < zBox: zMin = zBox + Om = float(omegaM) + + # converter from redshift to comoving distance + zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000) + zVsDX = np.zeros(len(zVsDY)) + for i in xrange(len(zVsDY)): + zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om) + + zBoxInMpc = vp.angularDiameter(zBox, Om=Om) + + boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \ + LIGHT_SPEED*lbox/2., zVsDX, zVsDY)-zBox + + print "RANGE", zBox, zBox+boxWidthZ + if zMax > zBox+boxWidthZ: zMax = zBox+boxWidthZ + + return zMin, zMax + +#------------------------------------------------------------------------------ +#------------------------------------------------------------------------------ + +inCatalogDir = catalogDir + +prefix = "masked_" + prefix +voidOutputDir += "masked/" +figDir += "masked/" +logDir += "masked/" +scriptDir += "masked/" +catalogDir += "masked/" +dataType = "observation" + +if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir) +if not os.access(catalogDir, os.F_OK): os.mkdir(catalogDir) + +#------------------------------------------------------------------------------ +# first the directly downsampled runs +# Note: ss0.002 ~ SDSS DR7 dim2 +# ss0.0004 ~ SDSS DR9 mid +baseResolution = float(numPart)/lbox/lbox/lbox # particles/Mpc^3 +for thisSubSample in subSamples: + + keepFraction = float(thisSubSample) / baseResolution + maxKeep = keepFraction * numPart + minRadius = int(np.ceil(lbox/maxKeep**(1./3))) + + if args.script or args.all: + print " Doing subsample", thisSubSample, " scripts" + setName = prefix+"ss"+str(thisSubSample) + writeObservationScript(setName, "masked_md.ss"+str(thisSubSample)+"_z", + maskFileName, + scriptDir, catalogDir, fileNums, + redshifts, False, minRadius, omegaM) + writeObservationScript(setName, "masked_md.ss"+str(thisSubSample)+"_z", + maskFileName, + scriptDir, catalogDir, fileNums, + redshifts, True, minRadius, omegaM) + + if args.subsample or args.all: + print " Doing subsample", thisSubSample + + for (iRedshift, redshift) in enumerate(redshifts): + print " redshift", redshift + + zMin = float(redshift) + zMax = zMin + 0.1 + zMin, zMax = getMaskedZRange(lbox, float(redshift), zMin, zMax, omegaM) + + if dataFormat == "multidark" or dataFormat == "random": + inFileName = inCatalogDir+"/md.ss"+str(thisSubSample)+"_z"+redshift+".dat" + outFileName = catalogDir+"/masked_md.ss"+str(thisSubSample)+"_z"+redshift+".dat" + + applyMask(inFileName, outFileName, maskFileName, lbox, float(redshift), + zMin, zMax, omegaM) +# ----------------------------------------------------------------------------- +# now halos +if (args.script or args.all) and dataFormat == "multidark": + print " Doing halo scripts" + + for minHaloMass in minHaloMasses: + + if dataFormat == "multidark": + setName = prefix+"halos_min"+str(minHaloMass) + writeScript(setName, "md.halos_min"+str(minHaloMass)+"_z", + scriptDir, catalogDir, fileNums, + redshifts, + numSubvolumes, numSlices, False, lbox, minRadius, omegaM) + writeScript(setName, "md.halos_min"+str(minHaloMass)+"_z", + scriptDir, catalogDir, fileNums, + redshifts, + numSubvolumes, numSlices, True, lbox, minRadius, omegaM) + +if args.halos or args.all: + print " Doing halos" + + for minHaloMass in minHaloMasses: + print " min halo mass = ", minHaloMass + + for (iRedshift, redshift) in enumerate(redshifts): + print " z = ", redshift + + dataFile = catalogDir+haloFileBase+fileNums[iRedshift] + inFile = open(dataFile, 'r') + numPart = 0 + for line in inFile: + line = line.split(',') + if minHaloMass == "none" or float(line[6]) > minHaloMass: + numPart += 1 + inFile.close() + + sampleName = "md.halos_min"+str(minHaloMass)+"_z"+redshift + outFile = open(catalogDir+"/"+sampleName+".dat", 'w') + +# ----------------------------------------------------------------------------- +# now the SDSS HOD + +if (args.script or args.all) and dataFormat == "multidark": + print " Doing DR7 HOD scripts" + if dataFormat == "multidark": + setName = prefix+"hod_dr72dim2" + writeScript(setName, "md.hod_dr72dim2_z", + scriptDir, catalogDir, fileNums, redshifts, + numSubvolumes, numSlices, False, lbox, 5, omegaM) + writeScript(setName, "md.hod_dr72dim2_z", + scriptDir, catalogDir, fileNums, redshifts, + numSubvolumes, numSlices, True, lbox, 5, omegaM) + +if args.hod or args.all: + print " Doing DR7 HOD" + for (iRedshift, redshift) in enumerate(redshifts): + print " z = ", redshift + + sampleName = getSampleName("md.hod_dr72dim2", redshift, False) + outFileName = catalogDir+"/"+sampleName+".dat" + +# ----------------------------------------------------------------------------- +# now the BOSS HOD +if (args.script or args.all) and dataFormat == "multidark": + print " Doing DR9 HOD scripts" + if dataFormat == "multidark": + setName = prefix+"hod_dr9mid" + writeScript(setName, "md.hod_dr9mid_z", + scriptDir, catalogDir, fileNums, redshifts, + numSubvolumes, numSlices, False, lbox, 15, omegaM) + writeScript(setName, "md.hod_dr9mid_z", + scriptDir, catalogDir, fileNums, redshifts, + numSubvolumes, numSlices, True, lbox, 15, omegaM) + +if args.hod or args.all: + print " Doing DR9 HOD" + for (iRedshift, redshift) in enumerate(redshifts): + print " z = ", redshift + + outFileName = catalogDir+"/"+sampleName+".dat" + + if dataFormat == "multidark" or dataFormat == "random": + sampleName = getSampleName("md.hod_dr9mid", redshift, False) + inFileName = inCatalogDir+"/"+sampleName+".dat" + outFileName = catalogDir+"/masked_"+sampleName+".dat" + + zMin = float(redshift) + zMax = zMin + 0.1 + zMin, zMax = getMaskedZRange(lbox, float(redshift), zMin, zMax, omegaM) + + applyMask(inFileName, outFileName, maskFileName, lbox, float(redshift), + zMin, zMax, omegaM) diff --git a/pipeline/datasets/multidark.py b/pipeline/datasets/multidark.py index a6a2563..38f081c 100644 --- a/pipeline/datasets/multidark.py +++ b/pipeline/datasets/multidark.py @@ -49,7 +49,7 @@ numSubvolumes = 1 prefix = "md_" # list of desired subsamples - these are in unts of h Mpc^-3! -#subSamples = [ 1.0 ] +#subSamples = [0.0004] subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002)) # common filename of halo files @@ -62,7 +62,7 @@ minHaloMasses = (("none", 2e12, 1.23e13)) # adjust these two parameters given the memory contraints on your system: # numZobovDivisions: how many sub-volumes per dimension will zobov process # numZobovThreads: how many sub-volumes to process at once? -numZobovDivisions = 4 +numZobovDivisions = 2 numZobovThreads = 2 # simulation information @@ -71,6 +71,9 @@ lbox = 1000 # Mpc/h omegaM = 0.27 hubble = 0.70 +# the mask file which is used by applyMaskToMock +maskFileName = os.getenv("HOME")+"/workspace/Voids/catalogs/boss/boss_mask_final.fits" + # END CONFIGURATION # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 5abacb5..b69ca74 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -9,13 +9,6 @@ import pickle # ------------------------------------------------------------------------------ -def my_import(name): - mod = __import__(name) - components = name.split('.') - for comp in components[1:]: - mod = getattr(mod, comp) - return mod - if (len(sys.argv) == 0): print "Usage: ./anylyzeVoids.py parameter_file.py" exit(-1) diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index 8d766b4..66f9c25 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -11,15 +11,6 @@ import void_python_tools as vp import argparse import imp -# ------------------------------------------------------------------------------ - -def my_import(name): - mod = __import__(name) - components = name.split('.') - for comp in components[1:]: - mod = getattr(mod, comp) - return mod - # ----------------------------------------------------------------------------- LIGHT_SPEED = 299792.458 @@ -83,7 +74,7 @@ def writeScript(setName, dataFileNameBase, import os from void_python_tools.backend.classes import * -continueRun = True # set to True to enable restarting aborted jobs +continueRun = False # set to True to enable restarting aborted jobs startCatalogStage = 1 endCatalogStage = 4 @@ -133,6 +124,7 @@ newSample = Sample(dataFile = "{dataFile}", zBoundaryMpc = ({zMinMpc}, {zMaxMpc}), omegaM = {omegaM}, minVoidRadius = {minRadius}, + profileBinSize = 1.0, includeInHubble = True, partOfCombo = False, isCombo = False, @@ -143,12 +135,12 @@ newSample = Sample(dataFile = "{dataFile}", useLightCone = {useLightCone}, subsample = {subsample}) dataSampleList.append(newSample) -newSample.addStack({zMin}, {zMax}, {minRadius} , {minRadius}+2, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius} , {minRadius}+4, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius}+2, {minRadius}+6, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius}+6, {minRadius}+10, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius}+10, {minRadius}+18, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False) +newSample.addStack({zMin}, {zMax}, 2*{minRadius} , 2*{minRadius}+2, True, False) +newSample.addStack({zMin}, {zMax}, 2*{minRadius} , 2*{minRadius}+4, True, False) +newSample.addStack({zMin}, {zMax}, 2*{minRadius}+2, 2*{minRadius}+6, True, False) +newSample.addStack({zMin}, {zMax}, 2*{minRadius}+6, 2*{minRadius}+10, True, False) +newSample.addStack({zMin}, {zMax}, 2*{minRadius}+10, 2*{minRadius}+18, True, False) +newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, False) """ for (iFile, redshift) in enumerate(redshifts): fileNum = fileNums[iFile] diff --git a/plotting/datasetsToPlot.py b/plotting/datasetsToPlot.py index e4d6ac3..476bd8d 100755 --- a/plotting/datasetsToPlot.py +++ b/plotting/datasetsToPlot.py @@ -4,9 +4,9 @@ workDir = "/home/psutter2/workspace/Voids/" figDir = "./figs" -sampleDirList = [ "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/", - #"multidark/md_ss01.0_pv/sample_md_ss1.0_pv_z0.56_d00/", - "multidark/md_halos_min1.23e13_pv/sample_md_halos_min1.23e13_pv_z0.56_d00/", +sampleDirList = [ + "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/", + "multidark/md_halos_min1.23e+13_pv/sample_md_halos_min1.23e+13_pv_z0.56_d00/", "random/ran_ss0.0004/sample_ran_ss0.0004_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/", diff --git a/plotting/plotCompareDensCon.py b/plotting/plotCompareDensCon.py index 8e21a8a..068be26 100755 --- a/plotting/plotCompareDensCon.py +++ b/plotting/plotCompareDensCon.py @@ -14,8 +14,6 @@ import argparse # ------------------------------------------------------------------------------ -from datasetsToPlot import * - plotNameBase = "compdenscon" obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing? @@ -24,10 +22,20 @@ 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) @@ -38,10 +46,10 @@ for sampleDir in sampleDirList: dataSampleList.append(pickle.load(input)) plt.clf() -plt.xlabel("Void Radius (Mpc/h)") +plt.xlabel("Void Density Contrast") plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]") plt.yscale('log') -plt.xlim(xmax=80.) +plt.xlim(xmax=5.) plotName = plotNameBase @@ -62,6 +70,7 @@ for (iSample,sample) in enumerate(dataSampleList): 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): diff --git a/plotting/plotCompareDist.py b/plotting/plotCompareDist.py index 77804fb..84273c0 100755 --- a/plotting/plotCompareDist.py +++ b/plotting/plotCompareDist.py @@ -14,8 +14,6 @@ import argparse # ------------------------------------------------------------------------------ -from datasetsToPlot import * - plotNameBase = "compdist" obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing? @@ -24,10 +22,20 @@ 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) @@ -38,10 +46,10 @@ for sampleDir in sampleDirList: dataSampleList.append(pickle.load(input)) plt.clf() -plt.xlabel("Void Radius (Mpc/h)") +plt.xlabel("Void Radius [Mpc/h]") plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]") plt.yscale('log') -plt.xlim(xmax=80.) +plt.xlim(xmax=120.) plotName = plotNameBase diff --git a/plotting/plotDensConVsR.py b/plotting/plotDensConVsR.py new file mode 100755 index 0000000..b9e15f7 --- /dev/null +++ b/plotting/plotDensConVsR.py @@ -0,0 +1,99 @@ +#!/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/python_tools/void_python_tools/backend/classes.py b/python_tools/void_python_tools/backend/classes.py index 3eb8b0a..18c2ffc 100755 --- a/python_tools/void_python_tools/backend/classes.py +++ b/python_tools/void_python_tools/backend/classes.py @@ -167,3 +167,10 @@ def jobSuccessful(logFile, doneString): def getStackSuffix(zMin, zMax, rMin, rMax, dataPortion): return "z"+str(zMin)+"-"+str(zMax)+"_"+str(rMin)+"-"+str(rMax)+\ "Mpc"+"_"+dataPortion + +def my_import(name): + mod = __import__(name) + components = name.split('.') + for comp in components[1:]: + mod = getattr(mod, comp) + return mod diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index f9d0ae2..53fb81c 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -764,12 +764,10 @@ def launchProfile(sample, stack, voidDir=None, logFile=None, continueRun=None): fp.close() return - # TEST if sample.profileBinSize == "auto": - density = 0.5 * 50 / Rcircular + density = 0.5 * 50 / Rcircular / 2 else: - #density = 0.5 * 50 / 60 - density = 0.5 * 50 / Rcircular / sample.profileBinSize + density = sample.profileBinSize sys.stdout = open(logFile, 'w') sys.stderr = open(logFile, 'a') @@ -836,10 +834,16 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None, while badChain: Rexpect = (stack.rMin+stack.rMax)/2 Rtruncate = stack.rMin*3. + 1 # TEST - ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, - Niter=300000, - Nburn=100000, - Rextracut=Rtruncate) + if sample.dataType == "observation": + ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, + Niter=300000, + Nburn=100000, + Rextracut=Rtruncate) + else: + ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, + Niter=300000, + Nburn=100000, + Rextracut=Rtruncate) badChain = (args[0][0] > 0.5 or args[0][1] > stack.rMax or \ args[0][2] > stack.rMax) and \ (ntries < maxtries) @@ -862,8 +866,8 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None, rescaleFactor = 1.0 # TEST - plotting raw galaxy counts - vp.draw_fit(*args, delta_tick=50, vmax=1.0, delta_v=0.01, - #vp.draw_fit(*args, delta_tick=0.2, vmax=1.0, delta_v=0.01, + #vp.draw_fit(*args, delta_tick=50, vmax=1.0, delta_v=0.01, + vp.draw_fit(*args, delta_tick=0.2, vmax=1.0, delta_v=0.01, nocontour=True, model_max=1.0, delta_model=0.1, plotTitle=plotTitle, show_ratio=showRatio, showContours=showContours, @@ -965,8 +969,12 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None, # (not not fully accurate) plots for (iZBin, stack) in enumerate(sample.getUniqueZBins()): - aveDist = vp.aveStretchCone(stack.zMin, stack.zMax, - skyFrac = sample.skyFraction) + if sample.dataType == "observation": + aveDist = vp.aveStretchCone(stack.zMin, stack.zMax, + skyFrac = sample.skyFraction) + else: + aveDist = vp.aveStretch(stack.zMin, stack.zMax) + aveDistList[iZBin, 0] = stack.zMin aveDistList[iZBin, 1] = aveDist aveDistList[iZBin, 2] = 0.00125 @@ -1010,16 +1018,24 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None, # skyFrac = sample.skyFraction, # voidRedshifts=voidRedshifts) - aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, - skyFrac = sample.skyFraction) + if sample.dataType == "observation": + aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, + skyFrac = sample.skyFraction) + else: + aveDist = vp.aveStretch(zBin.zMin, zBin.zMax) + expList[0, iR, iZBin, 2] = aveDist for (iZCheck,zBinCheck) in enumerate(allZBins): for (iRCheck,rBinCheck) in enumerate(allRBins): if zBin.zMin == zBinCheck.zMin and zBin.zMax == zBinCheck.zMax: if rBin.rMin == rBinCheck.rMin and rBin.rMax == rBinCheck.rMax: - aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, - skyFrac = sample.skyFraction) + if sample.dataType == "observation": + aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, + skyFrac = sample.skyFraction) + else: + aveDist = vp.aveStretch(zBin.zMin, zBin.zMax) + allExpList[iSample, iRCheck, iZCheck, 0] = exp allExpList[iSample, iRCheck, iZCheck, 1] = expError @@ -1112,8 +1128,12 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None, if zBin.zMaxPlot > plotZmax: plotZmax = zBin.zMaxPlot if zBin.zMinPlot < plotZmin: plotZmin = zBin.zMinPlot - aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, - skyFrac = sample.skyFraction) + if sample.dataType == "observation": + aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, + skyFrac = sample.skyFraction) + else: + aveDist = vp.aveStretch(zBin.zMin, zBin.zMax) + aveDistList[iZ, 0] = zBin.zMin aveDistList[iZ, 1] = aveDist aveDistList[iZ, 2] = 0.00125 diff --git a/python_tools/void_python_tools/plotting/plotDefs.py b/python_tools/void_python_tools/plotting/plotDefs.py index a999d41..d1da991 100644 --- a/python_tools/void_python_tools/plotting/plotDefs.py +++ b/python_tools/void_python_tools/plotting/plotDefs.py @@ -1,15 +1,15 @@ LIGHT_SPEED = 299792.458 -colorList = ['r', 'b', 'g', 'y', 'c', 'm', 'y', +colorList = ['r', 'b', 'g', 'y', 'c', 'm', 'brown', 'grey', - 'darkred', 'orange', 'pink', 'darkblue', - 'lightblue', 'chocolate', - 'indigo', 'lightseagreen', 'maroon', 'olive', - 'royalblue', 'palevioletred', 'seagreen', 'tomato', - 'aquamarine', 'darkslateblue', - 'khaki', 'lawngreen', 'mediumorchid', - 'orangered', 'thistle' - 'yellowgreen'] + 'darkred', 'orange', 'pink', 'darkblue', + 'lightblue', 'chocolate', + 'indigo', 'lightseagreen', 'maroon', 'olive', + 'royalblue', 'palevioletred', 'seagreen', 'tomato', + 'aquamarine', 'darkslateblue', + 'khaki', 'lawngreen', 'mediumorchid', + 'orangered', 'thistle' + 'yellowgreen'] linewidth = 4 fontsize = 12