From 21431f9a315261ce0e63d10d661a5ac296618914 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 2 Nov 2012 08:19:23 -0500 Subject: [PATCH 01/32] added a-p analysis script to repo; some minor bug fixes to adjust tools to work with simulation data --- GetQhull.cmake | 18 - c_tools/mock/generateMock.cpp | 1 + pipeline/generateCatalog.py | 3 +- pipeline/prepareCatalogs.py | 6 +- .../void_python_tools/backend/backend.py | 1453 ----------------- .../void_python_tools/backend/launchers.py | 60 +- 6 files changed, 41 insertions(+), 1500 deletions(-) delete mode 100644 GetQhull.cmake delete mode 100755 python_tools/void_python_tools/backend/backend.py diff --git a/GetQhull.cmake b/GetQhull.cmake deleted file mode 100644 index ac54728..0000000 --- a/GetQhull.cmake +++ /dev/null @@ -1,18 +0,0 @@ -SET(QHULL_BASE_PATH CACHE PATH "Qhull base path") - -find_path(QHULL_INCLUDE_PATH qhull_a.h HINTS ${QHULL_BASE_PATH}/src/libqhull) -find_path(QHULL_CPP_INCLUDE_PATH Qhull.h HINTS ${QHULL_BASE_PATH}/src/libqhullcpp) -find_library(QHULL_LIBRARY qhull_p HINTS ${QHULL_BASE_PATH}/lib) -find_library(QHULL_CPP_LIBRARY qhullcpp HINTS ${QHULL_BASE_PATH}/lib) -find_library(QHULL_P_LIBRARY qhullstatic_p HINTS ${QHULL_BASE_PATH}/lib) - -if ((NOT QHULL_INCLUDE_PATH) OR (NOT QHULL_CPP_LIBRARY)) - message(SEND_ERROR "Qhull library not found") -endif((NOT QHULL_INCLUDE_PATH) OR (NOT QHULL_CPP_LIBRARY)) - -set(QHULL_INCLUDES ${QHULL_INCLUDE_PATH} ${QHULL_INCLUDE_PATH}/.. ${QHULL_CPP_INCLUDE_PATH} ${QHULL_BASE_PATH}/src) -set(QHULL_LIBRARIES ${QHULL_CPP_LIBRARY} ${QHULL_P_LIBRARY}) - -add_definitions(-Dqh_QHpointer) - -mark_as_advanced(QHULL_INCLUDE_PATH QHULL_CPP_INCLUDE_PATH QHULL_LIBRARY QHULL_CPP_LIBRARY QHULL_P_LIBRARY) diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index 3316436..797ca65 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -391,6 +391,7 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& f.add_att("range_y_max", ranges[1][1]); f.add_att("range_z_min", ranges[2][0]); f.add_att("range_z_max", ranges[2][1]); + f.add_att("mask_index", -1); NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart); NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim); diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index eec1402..025bbc2 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -23,9 +23,8 @@ if (len(sys.argv) > 1): filename = sys.argv[1] print " Loading parameters from", filename if not os.access(filename, os.F_OK): - print " Cannot find parameter file!" + print " Cannot find parameter file %s!" % filename exit(-1) - #parms = __import__(filename[:-3], globals(), locals(), ['*']) parms = imp.load_source("name", filename) globals().update(vars(parms)) else: diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index c2bb424..f45e31b 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -196,8 +196,10 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False) mySubvolume = "%d%d" % (iX, iY) - sampleName = getSampleName(prefix, base, redshift, useVel, - iSlice=iSlice, iVol=mySubvolume) + sampleName = getSampleName(prefix, base, sliceMin, useVel, + iSlice=-1, iVol=mySubvolume) + #sampleName = getSampleName(prefix, base, redshift, useVel, + # iSlice=iSlice, iVol=mySubvolume) scriptFile.write(sampleInfo.format(dataFile=dataFileName, dataFormat=dataFormat, diff --git a/python_tools/void_python_tools/backend/backend.py b/python_tools/void_python_tools/backend/backend.py deleted file mode 100755 index 5209d25..0000000 --- a/python_tools/void_python_tools/backend/backend.py +++ /dev/null @@ -1,1453 +0,0 @@ -#!/usr/bin/env python - -# classes and routines used to support the analyzeVoids script - -import os -import glob -import voidProject as vp -import numpy as np -import numpy.ma as ma -import os -import shutil -import glob -import subprocess -import sys -from pylab import figure -from netCDF4 import Dataset - -NetCDFFile = Dataset -ncFloat='f8' # Double precision - -class Stack: - zMin = 0.0 - zMax = 0.1 - rMin = 5 - rMax = 15 - zMinPlot = 0.0 - zMaxPlot = 0.1 - includeInHubble = True - partOfCombo = False - needProfile = True - rescaleMode = "rmax" # options: "rmax" to scale to largest void in stack - # "rv" normalize each void - - def __init__(self, zMin, zMax, rMin, rMax, includeInHubble, partOfCombo, - zMinPlot=None, needProfile=True, rescaleMode="rmax"): - self.zMin = zMin - self.zMax = zMax - self.rMin = rMin - self.rMax = rMax - self.zMaxPlot = zMax - self.includeInHubble = includeInHubble - self.partOfCombo = partOfCombo - self.needProfile = needProfile - self.rescaleMode = rescaleMode - - if zMinPlot == None: - self.zMinPlot = self.zMin - else: - self.zMinPlot = zMinPlot - -class Sample: - dataType = "observation" - dataFile = "lss.dr72dim.dat" - fullName = "lss.dr72dim.dat" - nickName = "dim" - zobovDir = "" - maskFile = "rast_window_512.fits" - selFunFile = "czselfunc.all.dr72dim.dat" - skyFraction = 0.19 - zBoundary = (0.0, 0.1) - zRange = (0.0, 0.1) - minVoidRadius = 5 - fakeDensity = 0.01 - profileBinSize = 2 # Mpc - volumeLimited = True - includeInHubble = True - partOfCombo = False - isCombo = False - comboList = [] - numSubVolumes = 2 - - # applies to simulations only - boxLen = 1024 # Mpc/h - usePecVel = False - - stacks = [] - - def __init__(self, dataFile="", fullName="", - nickName="", maskFile="", selFunFile="", - zBoundary=(), zRange=(), - minVoidRadius=0, fakeDensity=0.01, volumeLimited=True, - includeInHubble=True, partOfCombo=False, isCombo=False, - comboList=(), profileBinSize=2.0, skyFraction=0.19, - dataType="observation", numSubVolumes=2, - boxLen=1024, usePecVel=False): - self.dataFile = dataFile - self.fullName = fullName - self.nickName = nickName - self.maskFile = maskFile - self.selFunFile = selFunFile - self.zBoundary = zBoundary - self.zRange = zRange - self.minVoidRadius = minVoidRadius - self.fakeDensity = fakeDensity - self.volumeLimited = volumeLimited - self.includeInHubble = includeInHubble - self.partOfCombo = partOfCombo - self.isCombo = isCombo - self.comboList = comboList - self.zobovDir = None - self.profileBinSize = profileBinSize - self.skyFraction = skyFraction - self.dataType = dataType - self.numSubVolumes = numSubVolumes - self.boxLen = boxLen - self.usePecVel = usePecVel - - self.stacks = [] - - def addStack(self, zMin, zMax, rMin, rMax, - includeInHubble, partOfCombo,zMinPlot=None, - needProfile=True, rescaleMode="rmax"): - self.stacks.append(Stack(zMin, zMax, rMin, rMax, - includeInHubble, partOfCombo, - zMinPlot=zMinPlot, needProfile=needProfile, - rescaleMode=rescaleMode)) - - def getHubbleStacks(self): - stacksForHubble = [] - for stack in self.stacks: - if stack.includeInHubble: - stacksForHubble.append(stack) - return stacksForHubble - - def getUniqueRBins(self): - uniqueRStacks = [] - for stack in self.stacks: - if not stack.includeInHubble: - continue - alreadyHere = False - for stackCheck in uniqueRStacks: - if stack.rMin == stackCheck.rMin and stack.rMax == stackCheck.rMax: - alreadyHere = True - break - if not alreadyHere: - uniqueRStacks.append(stack) - return uniqueRStacks - - def getUniqueZBins(self): - uniqueZStacks = [] - for stack in self.stacks: - if not stack.includeInHubble: - continue - alreadyHere = False - for stackCheck in uniqueZStacks: - if stack.zMin == stackCheck.zMin and stack.zMax == stackCheck.zMax: - alreadyHere = True - break - if not alreadyHere: - uniqueZStacks.append(stack) - return uniqueZStacks - -# ----------------------------------------------------------------------------- -# ----------------------------------------------------------------------------- -def jobSuccessful(logFile, doneString): - jobDone = False - checkLine = "" - if os.access(logFile, os.F_OK): - filelines = file(logFile, "r").readlines() - if len(filelines) >= 1: - checkLine = filelines[-1] - jobDone = (checkLine == doneString) - return jobDone - -def getStackSuffix(zMin, zMax, rMin, rMax, dataPortion): - return "z"+str(zMin)+"-"+str(zMax)+"_"+str(rMin)+"-"+str(rMax)+\ - "Mpc"+"_"+dataPortion - - -# ----------------------------------------------------------------------------- -# ----------------------------------------------------------------------------- -# routines which communicate with individual data analysis portions - -# they make the analyzeVoids.py script easier to read - -# ----------------------------------------------------------------------------- -def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, - zobovDir=None, figDir=None, logFile=None, useLCDM=False, - continueRun=None, dataType=None): - - if dataType == "observation": - sampleName = sample.fullName - - if sample.dataFile == "": - datafile = inputDataDir+"/"+sampleName - else: - dataFile = sample.dataFile - - maskFile = sample.maskFile - - if useLCDM: - useLCDMFlag = "useLCDM" - else: - useLCDMFlag = "" - - conf=""" - catalog %s - mask %s - output %s - params %s - zMin %g - zMax %g - density_fake %g - %s - """ % (datafile, maskFile, zobovDir+"/zobov_slice_"+sampleName, - zobovDir+"/zobov_slice_"+sampleName+".par", - sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity, - useLCDMFlag) - - parmFile = os.getcwd()+"/generate.par" - - file(parmFile, mode="w").write(conf) - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - cmd = "%s --configFile=%s >& %s" % (binPath,parmFile,logFile) - os.system(cmd) - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - if os.access(parmFile, os.F_OK): - os.unlink(parmFile) - - if os.access("contour_map.fits", os.F_OK): - os.system("mv %s %s" % ("contour_map.fits", zobovDir)) - os.system("mv %s %s" % ("mask_map.fits", zobovDir)) - - if os.access("comoving_distance.txt", os.F_OK): - os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) - - if os.access("mask_index.txt", os.F_OK): - os.system("mv %s %s" % ("mask_index.txt", zobovDir)) - os.system("mv %s %s" % ("total_particles.txt", zobovDir)) - os.system("mv %s %s" % ("sample_info.txt", zobovDir)) - - if os.access("galaxies.txt", os.F_OK): - os.system("mv %s %s" % ("galaxies.txt", zobovDir)) - os.system("mv %s %s" % ("mock_galaxies.txt", zobovDir)) - os.system("mv %s %s" % ("mock_boundary.txt", zobovDir)) - os.system("mv %s %s" % ("mock_sphere.txt", zobovDir)) - - else: # simulation - sampleName = sample.fullName - - datafile = inputDataDir+"/"+sample.dataFile - - if sample.includePecVel: - includePecVelString = "peculiarVelocities" - else: - includePecVelString = "" - - conf=""" - multidark %s - output %s - outputParameter %s - %s - rangeX_min %g - rangeX_max %g - rangey_min %g - rangeY_max %g - rangeZ_min %g - rangeZ_max %g - """ % (datafile, zobovDir+"/zobov_slice_"+sampleName, - zobovDir+"/zobov_slice_"+sampleName+".par", - includePecVelString, - 0, sample.boxLen, 0, sample.boxLen, - sample.zBoundary[0], sample.zBoundary[1]) - - parmFile = os.getcwd()+"/generate.par" - - file(parmFile, mode="w").write(conf) - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - cmd = "%s --configFile=%s >& %s" % (binPath,parmFile,logFile) - os.system(cmd) - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - if os.access(parmFile, os.F_OK): - os.unlink(parmFile) - -# ----------------------------------------------------------------------------- -def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, - dataType=None, numSubVolumes=None): - - sampleName = sample.fullName - - datafile = zobovDir+"zobov_slice_"+sampleName - - logFile = logDir+"/zobov_"+sampleName+".out" - - vozScript = "./scr_"+sampleName - - if os.access(vozScript, os.F_OK): - os.unlink(vozScript) - - if dataType == "observation": - maskIndex = open(zobovDir+"/mask_index.txt", "r").read() - totalPart = open(zobovDir+"/total_particles.txt", "r").read() - maxDen = 0.2*float(maskIndex)/float(totalPart) - else: - maskIndex = 999999999 - maxDen = 0.2 - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - for fileName in glob.glob(zobovDir+"/part._"+sampleName+".*"): - os.unlink(fileName) - - if os.access(zobovDir+"/adj_"+sampleName+".dat", os.F_OK): - os.unlink(zobovDir+"adj_"+sampleName+".dat") - - if os.access(zobovDir+"/voidDesc_"+sampleName+".out", os.F_OK): - os.unlink(zobovDir+"/voidDesc_"+sampleName+".out") - - cmd = "%s/vozinit %s 0.1 1.0 2 %s %g %s %s %s >& %s" % \ - (binPath, datafile, \ - "_"+sampleName, numSubVolumes, \ - binPath, zobovDir, maskIndex, logFile) - os.system(cmd) - - cmd = "./%s >> %s 2>&1" % (vozScript, logFile) - os.system(cmd) - - cmd = "%s/jozov %s %s %s %s %s %g %s >> %s 2>&1" % \ - (binPath, \ - zobovDir+"/adj_"+sampleName+".dat", \ - zobovDir+"/vol_"+sampleName+".dat", \ - zobovDir+"/voidPart_"+sampleName+".dat", \ - zobovDir+"/voidZone_"+sampleName+".dat", \ - zobovDir+"/voidDesc_"+sampleName+".out", \ - maxDen, maskIndex, logFile) - os.system(cmd) - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - if os.access(vozScript, os.F_OK): - os.unlink(vozScript) - -# ----------------------------------------------------------------------------- -def launchPrune(sample, binPath, thisDataPortion=None, - summaryFile=None, logFile=None, zobovDir=None, - continueRun=None, dataType=None): - - sampleName = sample.fullName - - numVoids = sum(1 for line in \ - open(zobovDir+"/voidDesc_"+sampleName+".out")) - numVoids -= 2 - - if dataType == "observation": - mockIndex = open(zobovDir+"/mask_index.txt", "r").read() - totalPart = open(zobovDir+"/total_particles.txt", "r").read() - maxDen = 0.2*float(mockIndex)/float(totalPart) - else: - mockIndex = "9999999999" - maxDen = 0.2 - - if not (continueRun and jobSuccessful(logFile, "NetCDF: Not a valid ID\n")): - cmd = binPath - cmd += " --partFile=" + zobovDir+"/zobov_slice_"+str(sampleName) - cmd += " --voidDesc=" + zobovDir+"/voidDesc_"+str(sampleName)+".out" - cmd += " --void2Zone="+zobovDir+"/voidZone_"+str(sampleName)+".dat" - cmd += " --zone2Part=" + zobovDir+"/voidPart_"+str(sampleName)+".dat" - cmd += " --partVol=" + zobovDir+"/vol_"+str(sampleName)+".dat" - cmd += " --extraInfo=" + zobovDir+"/zobov_slice_"+str(sampleName)+\ - ".par" - cmd += " --tolerance=1.0" - cmd += " --dataPortion=" + thisDataPortion - cmd += " --mockIndex=" + str(mockIndex) - cmd += " --maxCentralDen=" + str(maxDen) - cmd += " --zMin=" + str(sample.zBoundary[0]) - cmd += " --zMax=" + str(sample.zBoundary[1]) - cmd += " --rMin=" + str(sample.minVoidRadius) - cmd += " --numVoids=" + str(numVoids) - cmd += " --output=" + zobovDir+"/voidDesc_"+\ - str(thisDataPortion)+"_"+\ - str(sampleName)+".out" - cmd += " --outCenters=" + zobovDir+"/barycenters_"+\ - str(thisDataPortion)+"_"+\ - str(sampleName)+".out" - cmd += " --outInfo=" + zobovDir+"/centers_"+\ - str(thisDataPortion)+"_"+\ - str(sampleName)+".out" - cmd += " --outNoCutInfo=" + zobovDir+"/centers_nocut_"+\ - str(thisDataPortion)+"_"+\ - str(sampleName)+".out" - cmd += " --outSkyPositions=" + zobovDir+"/sky_positions_"+\ - str(thisDataPortion)+"_"+\ - str(sampleName)+".out" - cmd += " --outDistances=" + zobovDir+"/boundaryDistances_"+\ - str(thisDataPortion)+"_"+\ - str(sampleName)+".out" - cmd += " >& " + logFile - os.system(cmd) - - if jobSuccessful(logFile, "NetCDF: Not a valid ID\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - -# ----------------------------------------------------------------------------- -def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, - voidDir=None, freshStack=True, runSuffix=None, - zobovDir=None, - INCOHERENT=False, ranSeed=None, summaryFile=None, - continueRun=None, dataType=None): - - sampleName = sample.fullName - - runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - stack.rMax, thisDataPortion) - - logFile = logDir+"/stack_"+sampleName+"_"+runSuffix+".out" - - treeFile = voidDir+"/tree.data" - - if (freshStack) and os.access(treeFile, os.F_OK): - os.unlink(treeFile) - - centralRadius = stack.rMin * 0.25 - - # restrict to relavent ranges of sample - zMin = max(sample.zRange[0],stack.zMin) * 3000 - zMax = min(sample.zRange[1],stack.zMax) * 3000 - - if dataType == "observation": - obsFlag = "observation" - else: - obsFlag = "" - - Rextracut = stack.rMin*3 + 1 - Rcircular = stack.rMin*3 + 2 - - if dataType == "observation": - maskIndex = open(zobovDir+"/mask_index.txt", "r").read() - totalPart = open(zobovDir+"/total_particles.txt", "r").read() - maxDen = 0.2*float(maskIndex)/float(totalPart) - else: - maskIndex = 999999999 - maxDen = 0.2 - - if INCOHERENT: - nullTestFlag = "INCOHERENT" - else: - nullTestFlag = "" - - if stack.rescaleMode == "rmax": - rescaleFlag = "rescale" - else: - rescaleFlag = "" - - conf=""" - desc %s - partzone %s - zonevoid %s - volumefile %s - Rmin %g - Rmax %g - particles %s - extraInfo %s - densityThreshold %g - centralRadius %g - edgeAvoidance %g - Circular %g - Zmin %g - Zmax %g - %s - %s - ranSeed %d - dataPortion %s - barycenters %s - boundaryDistances %s - %s - """ % \ - (zobovDir+"/voidDesc_"+thisDataPortion+"_"+sampleName+".out", - zobovDir+"/voidPart_"+sampleName+".dat", - zobovDir+"/voidZone_"+sampleName+".dat", - zobovDir+"/vol_"+sampleName+".dat", - stack.rMin, - stack.rMax, - zobovDir+"/zobov_slice_"+sampleName, - zobovDir+"/zobov_slice_"+sampleName+".par", - maxDen, - centralRadius, - Rextracut, - Rcircular, - zMin, - zMax, - obsFlag, - nullTestFlag, - ranSeed, - thisDataPortion, - zobovDir+"/barycenters_"+thisDataPortion+"_"+sampleName+".out", - zobovDir+"/boundaryDistances_"+thisDataPortion+"_"+sampleName+".out", - rescaleFlag) - - parmFile = os.getcwd()+"/stack.par" - - fp = file(parmFile, mode="w") - fp.write(conf) - - # continue stacking if requested; pull files to local directory - if os.access(treeFile, os.F_OK): - os.system("mv %s %s" % (treeFile, "./"+tree.data)) - fp.write("loadSearchTree\n") - fp.write("getTree\n") - fp.write("doExtraction\n") - else: - fp.write("dumpSearchTree\n") - fp.write("dumpTree\n") - fp.write("doExtraction\n") - fp.close() - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - cmd = "%s --configFile=%s >& %s" % \ - (binPath, parmFile, logFile) - os.system(cmd) - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - return - - minDist = None - maxDist = None - numVoids = None - numParticles = None - for line in open(logFile): - if "Minimum void distance is" in line: - line = line.split() - minDist = float(line[4])/3000 - - if "Maximum void distance is" in line: - line = line.split() - maxDist = float(line[4])/3000 - - if "Selected" in line: - line = line.split() - numVoids = line[1] - - if "Number of particles =" in line: - line = line.split() - numParticles = line[4] - - open(voidDir+"/num_voids.txt", "w").write(numVoids) - open(voidDir+"/num_particles.txt", "w").write(numParticles) - - if os.access(voidDir+"/NOVOID", os.F_OK): - os.unlink(voidDir+"/NOVOID") - - if (numVoids == "0"): - print " No voids found; skipping!" - fp = open(voidDir+"/NOVOID", "w") - fp.write("no voids found\n") - fp.close() - return - - emptyStack = False - for line in open(logFile): - if "EMPTY STACK" in line: - emptyStack = True - if emptyStack: - print " Stack is empty; skipping!" - fp = open(voidDir+"/NOVOID", "w") - fp.write("empty stack\n") - fp.close() - return - - # figure out box volume and average density - maskFile = sample.maskFile - sulFunFile = sample.selFunFile - - if not os.access(sample.selFunFile, os.F_OK) and not volumeLimited: - print " Cannot find", selFunFile, "!" - exit(-1) - - sys.stdout = open(logFile, 'a') - sys.stderr = open(logFile, 'a') - zMin = sample.zRange[0] - zMax = sample.zRange[1] - if not sample.volumeLimited: - props = vp.getSurveyProps(maskFile, stack.zMin, - stack.zMax, zMin, zMax, "all", - selectionFuncFile=sample.selFunFile) - else: - zMinForVol = sample.zBoundary[0] - zMaxForVol = sample.zBoundary[1] - props = vp.getSurveyProps(maskFile, zMinForVol, - zMaxForVol, zMin, zMax, "all") - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - - boxVol = props[0] - nbar = props[1] - - if sample.volumeLimited: - nbar = 1.0 - - summaryLine = runSuffix + " " + \ - thisDataPortion + " " + \ - str(stack.zMin) + " " + \ - str(stack.zMax) + \ - " " + str(numVoids) + " " + str(minDist) + \ - " " + str(maxDist) - if summaryFile != None: - open(summaryFile, "a").write(summaryLine+"\n") - - # move local outputs to storage directory - if os.access("tree.data", os.F_OK): - normalization = float(numParticles)/float(boxVol)/nbar - dataTemp = file("centers.txt", "r").readlines() - fp = file("normalizations.txt", "w") - for iVoid in xrange(len(dataTemp)): - fp.write(str(normalization)+"\n") - fp.close() - - os.system("mv %s %s" % ("tree.data", treeFile)) - os.system("mv %s %s" % ("void_indexes.txt", voidDir+"/")) - os.system("mv %s %s" % ("posx.nc", voidDir+"/posx.nc")) - os.system("mv %s %s" % ("posy.nc", voidDir+"/posy.nc")) - os.system("mv %s %s" % ("posz.nc", voidDir+"/posz.nc")) - os.system("mv %s %s" % ("redshifts.nc", voidDir+"/redshifts.nc")) - os.system("mv %s %s" % ("indexes.nc", voidDir+"/")) - os.system("mv %s %s" % ("kdtree_stackvoids.dat", voidDir+"/")) - os.system("mv %s %s" % ("centers.txt", voidDir+"/")) - os.system("mv %s %s" % ("sky_positions.txt", voidDir+"/")) - os.system("mv %s %s" % ("check.txt", voidDir+"/")) - os.system("mv %s %s" % ("tracer.txt", voidDir+"/")) - os.system("mv %s %s" % ("normalizations.txt", voidDir+"/")) - os.system("mv %s %s" % ("boundaryDistances.txt", voidDir+"/")) - - if os.access(parmFile, os.F_OK): - os.unlink(parmFile) - - return - -# ----------------------------------------------------------------------------- -def launchCombine(sample, stack, voidDir=None, logFile=None, - zobovDir=None, workDir=None, thisDataPortion=None): - - sampleName = sample.fullName - - runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - stack.rMax, thisDataPortion) - - sys.stdout = open(logFile, 'w') - sys.stderr = open(logFile, 'a') - - if os.access(voidDir+"/num_voids.txt", os.F_OK): - os.unlink(voidDir+"/num_voids.txt") - - doneGalUpdate = os.access(zobovDir+"donegalupdate", os.F_OK) - - for comboName in sample.comboList: - if not os.access(zobovDir+"/galaxies.txt", os.F_OK): - shutil.copy(workDir+"/sample_"+comboName+"/galaxies.txt", zobovDir) - elif not doneGalUpdate: - dataTemp = file(workDir+"/sample_"+comboName+"/galaxies.txt", -"r").read() - file(zobovDir+"/galaxies.txt", "a").write(dataTemp) - - sourceStackDir = workDir+"/sample_"+comboName+"/stacks_"+\ - runSuffix - - if os.access(sourceStackDir+"/NOVOID", os.F_OK): - continue - - if not os.access(voidDir+"/num_voids.txt", os.F_OK): - # just copy - shutil.copy(sourceStackDir+"/num_voids.txt", voidDir) - shutil.copy(sourceStackDir+"/num_particles.txt", voidDir) - shutil.copy(sourceStackDir+"/posx.nc", voidDir) - shutil.copy(sourceStackDir+"/posy.nc", voidDir) - shutil.copy(sourceStackDir+"/posz.nc", voidDir) - shutil.copy(sourceStackDir+"/indexes.nc", voidDir) - shutil.copy(sourceStackDir+"/redshifts.nc", voidDir) - shutil.copy(sourceStackDir+"/centers.txt", voidDir) - shutil.copy(sourceStackDir+"/void_indexes.txt", voidDir) - shutil.copy(sourceStackDir+"/sky_positions.txt", voidDir) - shutil.copy(sourceStackDir+"/normalizations.txt", voidDir) - shutil.copy(sourceStackDir+"/boundaryDistances.txt", voidDir) - - else: - # append data - dataTemp = file(voidDir+"/num_voids.txt", "r").readlines() - dataTemp2 = file(sourceStackDir+"/num_voids.txt", "r").\ - readlines() - dataTemp = np.array(dataTemp, dtype='i') - dataTemp2 = np.array(dataTemp2, dtype='i') - dataTemp += dataTemp2 - dataTemp = str(dataTemp[0]) - file(voidDir+"/num_voids.txt", "w").write(str(dataTemp)) - - dataTemp = file(voidDir+"/num_particles.txt", "r").readlines() - dataTemp2 = file(sourceStackDir+"/num_particles.txt", "r").\ - readlines() - dataTemp = np.array(dataTemp, dtype='i') - dataTemp2 = np.array(dataTemp2, dtype='i') - dataTemp += dataTemp2 - dataTemp = str(dataTemp[0]) - file(voidDir+"/num_particles.txt", "w").write(str(dataTemp)) - - dataTemp = file(sourceStackDir+"/centers.txt", "r").read() - file(voidDir+"/centers.txt", "a").write(dataTemp) - dataTemp = file(sourceStackDir+"/normalizations.txt", "r").\ - read() - file(voidDir+"/normalizations.txt", "a").write(dataTemp) - - dataTemp = file(sourceStackDir+"/boundaryDistances.txt","r").\ - read() - file(voidDir+"/boundaryDistances.txt", "a").write(dataTemp) - - dataTemp = file(sourceStackDir+"/sky_positions.txt", "r").\ - read() - file(voidDir+"/sky_positions.txt", "a").write(dataTemp) - - idxTemp = file(sourceStackDir+"/void_indexes.txt", "r").\ - readlines() - idxTemp = np.array(idxTemp, dtype='i') - dataTemp = (NetCDFFile(voidDir+"/posx.nc").\ - variables['array'])[0:] - idxTemp[:] += len(dataTemp) - fp = open(voidDir+"/void_indexes.txt", "a") - for idx in idxTemp: - fp.write(str(idx)+"\n") - fp.close() - dataTemp = () - - fp = NetCDFFile(voidDir+"/posx.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/posx.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/posx.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/posy.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/posy.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/posy.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/posz.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/posz.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/posz.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/redshifts.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/redshifts.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/redshifts.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/indexes.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/indexes.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/indexes.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - # add NOVOID if necessary - if not os.access(voidDir+"/centers.txt", os.F_OK): - fp = open(voidDir+"/NOVOID", "w") - fp.write("no voids in combo\n") - fp.close() - - fp = open(zobovDir+"/donegalupdate", "w") - fp.close() - - print "Done!" - - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - #else: - # print "already done!" - -# ----------------------------------------------------------------------------- -def launchProfile(sample, stack, voidDir=None, logFile=None, continueRun=None): - - sampleName = sample.fullName - - if os.access(voidDir+"/NOVOID", os.F_OK): - print "no stack here; skipping!" - return - - numVoids = open(voidDir+"/num_voids.txt", "r").readline() - numParticles = open(voidDir+"/num_particles.txt", "r").readline() - - Rextracut = stack.rMin*3 + 1 - Rcircular = stack.rMin*3 + 2 - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - profileParmFile = voidDir+"/params.txt" - fp = open(profileParmFile, "w") - fp.write(str(stack.rMin)+"\n") - fp.write(str(stack.rMax)+"\n") - fp.write(str(Rcircular)+"\n") - fp.write(numVoids+"\n") - fp.close() - - if sample.zRange[0] > stack.zMax or sample.zRange[1] < stack.zMin: - print "outside sample redshift range; skipping!" - fp = open(voidDir+"/NOVOID", "w") - fp.write("outside z range: profile\n") - fp.close() - return - - # TEST - if sample.profileBinSize == "auto": - density = 0.5 * 50 / Rcircular - else: - #density = 0.5 * 50 / 60 - density = 0.5 * 50 / Rcircular / sample.profileBinSize - - sys.stdout = open(logFile, 'w') - sys.stderr = open(logFile, 'a') - vp.build_1d_profile(base_dir=voidDir, density=density, - rescaleMode=stack.rescaleMode) - vp.build_2d_profile(base_dir=voidDir, density=density) - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - - if jobSuccessful(logFile, "Done!\n"): - print "done", numVoids - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - -# ----------------------------------------------------------------------------- -def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None, - continueRun=None, thisDataPortion=None): - - sampleName = sample.fullName - - runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - stack.rMax, thisDataPortion) - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - if os.access(voidDir+"/NOVOID", os.F_OK): - print "no voids here; skipping!" - return - - if stack.zMin < sample.zRange[0] or stack.zMax > sample.zRange[1]: - print "outside sample redshift range; skipping!" - return - - if sample.partOfCombo or not sample.includeInHubble: - print "sample not needed for further analysis; skipping!" - fp = open(voidDir+"/NOFIT", "w") - fp.write("sample not needed for hubble\n") - fp.close() - return - - if not stack.includeInHubble: - print "radius not needed for further analysis; skipping!" - return - - if not stack.includeInHubble: - print "redshift not needed for further analysis; skipping!" - return - - if os.access(figDir+"/stackedVoid_"+sampleName+"_"+\ - runSuffix+".eps", os.F_OK): - os.unlink(figDir+"/stackedVoid_"+sampleName+"_"+\ - runSuffix+".eps") - - sys.stdout = open(logFile, 'w') - sys.stderr = open(logFile, 'a') - - badChain = True - ntries = 0 - maxtries = 5 - while badChain: - Rexpect = (stack.rMin+stack.rMax)/2 - Rtruncate = stack.rMax*3. + 1 # TEST - 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) - ntries += 1 - - np.save(voidDir+"/chain.npy", ret) - np.savetxt(voidDir+"/fits.out", fits) - - plotTitle = "Sample: "+sample.nickName+\ - ", z = "+str(stack.zMin)+"-"+\ - str(stack.zMax)+", R = "+str(stack.rMin)+"-"+\ - str(stack.rMax)+r" $h^{-1}$ Mpc" - - showRatio = (ntries <= maxtries) - showContours = (ntries <= maxtries) - - if stack.rescaleMode == "rmax": - rescaleFactor = stack.rMax - else: - 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, - nocontour=True, model_max=1.0, delta_model=0.1, - plotTitle=plotTitle, show_ratio=showRatio, - showContours=showContours, - rescaleFactor=rescaleFactor) - figure(1).savefig(figDir+"/stackedVoid_"+sampleName+\ - "_"+runSuffix+".eps") - figure(1).savefig(figDir+"/stackedVoid_"+sampleName+\ - "_"+runSuffix+".pdf") - figure(1).savefig(figDir+"/stackedVoid_"+sampleName+\ - "_"+runSuffix+".png") - - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - if jobSuccessful(logFile, "Done!\n"): - print "done (", ntries, " tries)" - else: - print "FAILED!" - exit(-1) - - # record the measured stretch - exp = args[1] - expError = args[2] - outline = str(exp) + " " + str(expError) - open(voidDir+"/expansion.out", "w").write(outline) - - if os.access(voidDir+"/NOFIT", os.F_OK): - os.unlink(voidDir+"/NOFIT") - if ntries > maxtries: - print " No reliable fit found; skipping!" - fp = open(voidDir+"/NOFIT", "w") - fp.write("bad ellipticity fit\n") - fp.close() - return - - else: - print "already done!" - - -# ----------------------------------------------------------------------------- -def launchHubble(dataPortions=None, dataSampleList=None, logDir=None, - INCOHERENT=None, workDir=None, figDir=None, errorBars=None, - ZOBOV_PATH=None, continueRun=None, voidDir=None, - doPlot = True): - - for thisDataPortion in dataPortions: - print " For data portion", thisDataPortion - sys.stdout.flush() - - # merge the bins from all samples - numSamplesForHubble = 0 - maxRBins = 0 - maxZBins = 0 - allRBins = [] - allZBins = [] - - for sample in dataSampleList: - if sample.includeInHubble: - numSamplesForHubble += 1 - - for stack in sample.getHubbleStacks(): - alreadyHere = False - for stackCheck in allRBins: - if stack.rMin == stackCheck.rMin and stack.rMax==stackCheck.rMax: - alreadyHere = True - break - if not alreadyHere: - allRBins.append(stack) - - alreadyHere = False - for stackCheck in allZBins: - if stack.zMin == stackCheck.zMin and stack.zMax==stackCheck.zMax: - alreadyHere = True - break - if not alreadyHere: - allZBins.append(stack) - - allExpList = np.zeros((numSamplesForHubble,len(allRBins),len(allZBins),3)) - allExpList[:] = -1 - aveDistList = np.zeros((len(allZBins),3)) - - iSample = -1 - for sample in dataSampleList: - if not sample.includeInHubble: - continue - - iSample += 1 - - zobovDir = sample.zobovDir - sampleName = sample.fullName - - print " For data set", sampleName, "...", - sys.stdout.flush() - - logFile = logDir+"/hubble_"+sampleName+"_"+thisDataPortion+".out" - - #AVEDIST_PATH = ZOBOV_PATH+"/mytools/computeAverageDistortionPMS" - - # compute the average stretch in each redshift bin for cleaner-looking - # (not not fully accurate) plots - for (iZBin, stack) in enumerate(sample.getUniqueZBins()): - - aveDist = vp.aveStretchCone(stack.zMin, stack.zMax, - skyFrac = sample.skyFraction) - aveDistList[iZBin, 0] = stack.zMin - aveDistList[iZBin, 1] = aveDist - aveDistList[iZBin, 2] = 0.00125 - - numFits = 0 - fp = file(zobovDir+'/calculatedExpansions_'+thisDataPortion+'.txt', - mode="w") - - expList = np.zeros((1, len(sample.getUniqueRBins()), - len(sample.getUniqueZBins()), 3)) - - for (iZBin, zBin) in enumerate(sample.getUniqueZBins()): - for (iR, rBin) in enumerate(sample.getUniqueRBins()): - - runSuffix = getStackSuffix(zBin.zMin, zBin.zMax, rBin.rMin, - rBin.rMax, thisDataPortion) - - expFile = zobovDir+"/stacks_"+runSuffix+"/expansion.out" - if not (os.access(zobovDir+"/stacks_"+runSuffix+"/NOVOID", \ - os.F_OK) or - os.access(zobovDir+"/stacks_"+runSuffix+"/NOFIT", \ - os.F_OK) or - not os.access(expFile, os.F_OK)): - exp = float(open(expFile, "r").readline().split()[0]) - expError = float(open(expFile, "r").readline().split()[1]) - else: - exp = -1 - expError = 0 - - expList[0, iR, iZBin, 0] = exp - expList[0, iR, iZBin, 1] = expError - - # TODO - # compute the per-stack stretch for liklihood calculation - #runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - # stack.rMax, thisDataPortion) - #voidDir = zobovDir+"/stacks_" + runSuffix - #centersFile = voidDir+"/centers.txt" - #voidRedshifts = np.loadtxt(centersFile)[:,5] - #aveDist = vp.aveWeightedStretchCone(stack.zMin, stack.zMax, - # skyFrac = sample.skyFraction, - # voidRedshifts=voidRedshifts) - - aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, - skyFrac = sample.skyFraction) - 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) - allExpList[iSample, iRCheck, iZCheck, 0] = exp - allExpList[iSample, iRCheck, iZCheck, 1] = expError - - fp.write(str(exp) + " " + str(expError) + " "+ " "+ \ - str(aveDist) + " -1 " + runSuffix+"\n") - numFits += 1 - - fp.close() - - rlist = np.zeros((len(sample.getUniqueRBins()),2)) - for (iR,rBin) in enumerate(sample.getUniqueRBins()): - rlist[iR, 0] = rBin.rMin - rlist[iR, 1] = rBin.rMax - - zbase = np.zeros((len(sample.getUniqueZBins()),2)) - for (iZ,zBin) in enumerate(sample.getUniqueZBins()): - zbase[iZBin,0] = zBin.zMinPlot - zbase[iZBin,1] = zBin.zMaxPlot - #np.savetxt(zobovDir+'/zbase_'+thisDataPortion+'.txt', zbase) - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - sys.stdout = open(logFile, 'w') - sys.stderr = open(logFile, 'a') - plotTitle = "Sample: "+sample.nickName+", "+thisDataPortion+" voids" - if doPlot: - #vp.do_all_obs(zbase, expList, workDir+"/avedistortion_", - vp.do_all_obs(zbase, expList, aveDistList, - rlist, plotTitle=plotTitle, plotAve=True) - figure(1).savefig(figDir+"/hubble_"+sampleName+"_"+thisDataPortion+\ - ".eps",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_"+sampleName+"_"+thisDataPortion+\ - ".pdf",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_"+sampleName+"_"+thisDataPortion+\ - ".png",bbox_inches='tight') - else: - print "Skipping plot" - print "Done!" - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - # now do all samples in one plot - print " For data set combined...", - sys.stdout.flush() - - logFile = logDir+"/hubble_combined_"+thisDataPortion+".out" - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - - if errorBars == "ESTIMATED": - # replace calculated errors with estimated ones - errorBarFile = workDir+"/calculatedErrorBars_"+thisDataPortion+".txt" - if os.access(errorBarFile, os.F_OK): - print "REPLACING ERROR BARS!", - sys.stdout.flush() - fp = file(errorBarFile, mode="r") - ignoreLine = fp.readline() - for iZ in xrange(len(allExpList[0,0,:,0])): - ignoreLine = fp.readline() - - for iZ in xrange(len(allExpList[0,0,:,0])): - for iSample in xrange(len(allExpList[:,0,0,0])): - for iR in xrange(len(allExpList[0,:,0,0])): - line = fp.readline().split() - allExpList[iSample,iR,iZ,2] = allExpList[iSample,iR,iZ,1] - allExpList[iSample,iR,iZ,1] = float(line[1]) - - fp.close() - - rlist = np.zeros((len(allRBins),2)) - for (iR,rBin) in enumerate(allRBins): - rlist[iR, 0] = rBin.rMin - rlist[iR, 1] = rBin.rMax - - zbase = np.zeros((len(allZBins),2)) - plotZmax = 0.0 - plotZmin = 1.e99 - for (iZ,zBin) in enumerate(allZBins): - zbase[iZ,0] = zBin.zMinPlot - zbase[iZ,1] = zBin.zMaxPlot - - if zBin.zMaxPlot > plotZmax: plotZmax = zBin.zMaxPlot - if zBin.zMinPlot < plotZmin: plotZmin = zBin.zMinPlot - - aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, - skyFrac = sample.skyFraction) - aveDistList[iZ, 0] = zBin.zMin - aveDistList[iZ, 1] = aveDist - aveDistList[iZ, 2] = 0.00125 - - - shortSampleNames = list() - for sample in dataSampleList: - if sample.includeInHubble: - shortSampleNames.append(sample.nickName) - sys.stdout = open(logFile, 'w') - sys.stderr = open(logFile, 'a') - if doPlot: - if INCOHERENT: - #plotTitle = "all samples, incoherent "+\ - # thisDataPortion+" voids" - plotTitle = '' - else: - #plotTitle = "all samples, "+thisDataPortion+" voids" - plotTitle = '' - vp.do_all_obs(zbase, allExpList, aveDistList, - rlist, plotTitle=plotTitle, sampleNames=shortSampleNames, - plotAve=True, mulfac = 1.0, biasLine = 1.16, - plotZmin=plotZmin, plotZmax=plotZmax+0.2) - figure(1).savefig(figDir+"/hubble_combined_"+thisDataPortion+\ - ".eps",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_combined_"+thisDataPortion+\ - ".pdf",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_combined_"+thisDataPortion+\ - ".png",bbox_inches='tight') - - if INCOHERENT: - plotTitle = "all samples, incoherent "+\ - thisDataPortion+" voids (systematics corrected)" - else: - plotTitle = "all samples, "+thisDataPortion+\ - " voids (systematics corrected)" - vp.do_all_obs(zbase, allExpList, aveDistList, - rlist, plotTitle=plotTitle, sampleNames=shortSampleNames, - plotAve=True, mulfac = 1.16, - plotZmin=plotZmin, plotZmax=plotZmax+0.2) - figure(1).savefig(figDir+"/hubble_combined_"+thisDataPortion+\ - "_debiased.eps",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_combined_"+thisDataPortion+\ - "_debiased.pdf",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_combined_"+thisDataPortion+\ - "_debiased.png",bbox_inches='tight') - else: - print "Skipping plot" - print "Done!" - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - - # save all expansion data to a single file - fp = file(workDir+'/calculatedExpansions_'+thisDataPortion+'.txt', - mode="w") - fp.write(str(numSamplesForHubble) + " " + - str(len(allRBins)) + " " + str(len(allZBins)) + "\n") - for zBin in allZBins: - fp.write(str(zBin.zMin) + " " + str(zBin.zMax) + "\n") - - for iZ in xrange(len(allExpList[0,0,:,0])): - for iSample in xrange(len(allExpList[:,0,0,0])): - for iR in xrange(len(allExpList[0,:,0,0])): - fp.write(str(allExpList[iSample,iR,iZ,0]) + " " +\ - str(allExpList[iSample,iR,iZ,1]) + " " +\ - str(allExpList[iSample,iR,iZ,2]) + "\n") - fp.close() - - # save all void distribution data to a single file - fp = file(workDir+'/voidDistributions_'+thisDataPortion+'.txt', - mode="w") - fp.write(str(numSamplesForHubble) + " " + - str(len(allRBins)) + " " + str(len(allZBins)) + "\n") - for zBin in allZBins: - fp.write(str(zBin.zMin) + " " + str(zBin.zMax) + "\n") - - for zBin in allZBins: - for sample in dataSampleList: - if not sample.includeInHubble: continue - for rBin in allRBins: - runSuffix = getStackSuffix(zBin.zMin, zBin.zMax, rBin.rMin, - rBin.rMax, thisDataPortion) - voidDir = sample.zobovDir+"/stacks_" + runSuffix - centersFile = voidDir+"/centers.txt" - if os.access(centersFile, os.F_OK): - voidRedshifts = np.loadtxt(centersFile)[:,5] - #fp.write(str(len(voidRedshifts))+" ") - np.savetxt(fp, voidRedshifts[None]) - else: - fp.write("-1\n") - fp.close() - - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - -# ----------------------------------------------------------------------------- -def launchLikelihood(dataPortions=None, logDir=None, workDir=None, - continueRun=False): - - for thisDataPortion in dataPortions: - print " For data portion", thisDataPortion, "...", - sys.stdout.flush() - - logFile = logDir+"/likelihood_"+thisDataPortion+".out" - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - - sys.stdout = open(logFile, 'w') - sys.stderr = open(logFile, 'a') - - vp.build1dLikelihood(workDir+"/calculatedExpansions_"+\ - thisDataPortion+".txt", - workDir+"/voidDistributions_" + \ - thisDataPortion+".txt", - nbins = 20, - OmStart = 0.0, - OmEnd = 1.0, - biasStart = 1.0, - biasEnd = 1.2, - outputBase = workDir+"/1dlikelihoods_"+thisDataPortion+"_", - useBinAve = False) - - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - -# ----------------------------------------------------------------------------- -def launchEstimateErrorBars(workDir=None, nullDir=None, numNulls=None, - dataPortion=None): - - IVALUE = 0 - IERROR = 1 - - # open first null file to get sizes - nullFile = nullDir + "/" + str(1) + "/calculatedExpansions_" + \ - str(dataPortion) + ".txt" - fp = open(nullFile, "r") - binLine = fp.readline() - numSamples = int(binLine.split(" ")[0]) - numRBins = int(binLine.split(" ")[1]) - numZBins = int(binLine.split(" ")[2]) - numBins = numRBins * numSamples * numZBins - nRTotal = numRBins * numSamples - fp.close() - - # load null data - nullData = np.zeros([numNulls, numBins, 2]) - - for iNull in xrange(numNulls): - nullFile = nullDir + "/" + str(iNull) + "/calculatedExpansions_" + \ - str(dataPortion) + ".txt" - fp = open(nullFile, "r") - - binLine = fp.readline() - numSamples = int(binLine.split(" ")[0]) - numRBins = int(binLine.split(" ")[1]) - numZBins = int(binLine.split(" ")[2]) - numBins = numRBins * numSamples * numZBins - - zList = np.zeros( [numZBins,2] ) - for iZ in xrange(int(numZBins)): - line = fp.readline().split() - zList[iZ,0] = float(line[0]) - zList[iZ,1] = float(line[1]) - - for iZ in xrange(numZBins): - iRTotal = 0 - for iSample in xrange(numSamples): - for iR in xrange(numRBins): - line = fp.readline().split() - nullData[iNull,iRTotal*numZBins+iZ,IVALUE] = float(line[0]) - nullData[iNull,iRTotal*numZBins+iZ,IERROR] = float(line[1]) - iRTotal += 1 - fp.close() - - # use 68% limits to get error bars - errorBars = np.zeros([numBins, 2]) - for iBin in xrange(numBins): - binData = nullData[:,iBin,IVALUE] - validData = binData[:] > -1 - binData = binData[validData] - binData = np.sort(binData) - - if np.size(binData) > 0: - errorBars[iBin, 0] = binData[int(0.16*np.size(binData))] - errorBars[iBin, 1] = binData[int(0.84*np.size(binData))] - - mean = (errorBars[iBin, 1] + errorBars[iBin, 0])/2 - sig = (errorBars[iBin, 1] - errorBars[iBin, 0])/2 - errorBars[iBin, 0] = mean - errorBars[iBin, 1] = sig - else: - errorBars[iBin, :] = -1 - - # writes out error bars - fp = file(workDir+"/calculatedErrorBars_" + \ - str(dataPortion) + ".txt", mode="w") - fp.write(str(numSamples) + " " + str(numRBins) + " " + str(numZBins) + "\n") - for iZ in xrange(numZBins): - fp.write(str(zList[iZ,0]) + " " + str(zList[iZ,1]) + "\n") - - for iZ in xrange(numZBins): - iRTotal = 0 - for iSample in xrange(numSamples): - for iR in xrange(numRBins): - fp.write(str(errorBars[iRTotal*numZBins+iZ,0]) + " " + \ - str(errorBars[iRTotal*numZBins+iZ,1]) + "\n") - iRTotal += 1 - fp.close() - - - - diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 9368c33..92536cd 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -150,6 +150,9 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, else: print "already done!" + if os.access("comoving_distance.txt", os.F_OK): + os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) + if os.access(parmFile, os.F_OK): os.unlink(parmFile) @@ -460,34 +463,37 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, return # figure out box volume and average density - maskFile = sample.maskFile - sulFunFile = sample.selFunFile + if sample.dataType == "observation": + maskFile = sample.maskFile + sulFunFile = sample.selFunFile - if not os.access(sample.selFunFile, os.F_OK) and not volumeLimited: - print " Cannot find", selFunFile, "!" - exit(-1) + if not os.access(sample.selFunFile, os.F_OK) and not sample.volumeLimited: + print " Cannot find", selFunFile, "!" + exit(-1) - sys.stdout = open(logFile, 'a') - sys.stderr = open(logFile, 'a') - zMin = sample.zRange[0] - zMax = sample.zRange[1] - if not sample.volumeLimited: - props = vp.getSurveyProps(maskFile, stack.zMin, - stack.zMax, zMin, zMax, "all", - selectionFuncFile=sample.selFunFile) + sys.stdout = open(logFile, 'a') + sys.stderr = open(logFile, 'a') + zMin = sample.zRange[0] + zMax = sample.zRange[1] + if not sample.volumeLimited: + props = vp.getSurveyProps(maskFile, stack.zMin, + stack.zMax, zMin, zMax, "all", + selectionFuncFile=sample.selFunFile) + else: + zMinForVol = sample.zBoundary[0] + zMaxForVol = sample.zBoundary[1] + props = vp.getSurveyProps(maskFile, zMinForVol, + zMaxForVol, zMin, zMax, "all") + sys.stdout = sys.__stdout__ + sys.stderr = sys.__stderr__ + + boxVol = props[0] + nbar = props[1] + if sample.volumeLimited: + nbar = 1.0 else: - zMinForVol = sample.zBoundary[0] - zMaxForVol = sample.zBoundary[1] - props = vp.getSurveyProps(maskFile, zMinForVol, - zMaxForVol, zMin, zMax, "all") - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - - boxVol = props[0] - nbar = props[1] - - if sample.volumeLimited: nbar = 1.0 + boxVol = sample.boxLen**3 summaryLine = runSuffix + " " + \ thisDataPortion + " " + \ @@ -1173,7 +1179,11 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None, voidDir = sample.zobovDir+"/stacks_" + runSuffix centersFile = voidDir+"/centers.txt" if os.access(centersFile, os.F_OK): - voidRedshifts = np.loadtxt(centersFile)[:,5] + voidRedshifts = np.loadtxt(centersFile) + if voidRedshifts.ndim > 1: + voidRedshifts = voidRedshifts[:,5] + else: + voidRedshifts = voidRedshifts[5] #fp.write(str(len(voidRedshifts))+" ") np.savetxt(fp, voidRedshifts[None]) else: From e2ee7b22dd66525c4e532742e07da091fa029ec5 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Tue, 6 Nov 2012 08:01:08 -0600 Subject: [PATCH 02/32] removed periodic boundary handing in stacking, since we may be dealing with subsections From 0e2f551f7a6f7e05ca87234cc6a64eb62634d8cb Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 6 Nov 2012 10:58:16 -0500 Subject: [PATCH 03/32] Added Healpy to the dependencies to build --- external/external_build.cmake | 1 + external/external_python_build.cmake | 27 +++++++++++++++++++++++++++ external/python_build.cmake | 3 +++ external/python_install.cmake | 4 ++++ 4 files changed, 35 insertions(+) diff --git a/external/external_build.cmake b/external/external_build.cmake index 555d3dd..c28c6db 100644 --- a/external/external_build.cmake +++ b/external/external_build.cmake @@ -243,6 +243,7 @@ ExternalProject_Add(cfitsio INSTALL_COMMAND make install ) SET(CFITSIO_LIBRARY ${CMAKE_BINARY_DIR}/ext_build/cfitsio/lib/libcfitsio.a) +SET(CFITSIO_INCLUDE_PATH ${CMAKE_BINARY_DIR}/ext_build/cfitsio/include) ################# # Build Healpix diff --git a/external/external_python_build.cmake b/external/external_python_build.cmake index 92ca65f..e4580ae 100644 --- a/external/external_python_build.cmake +++ b/external/external_python_build.cmake @@ -2,6 +2,7 @@ INCLUDE(FindPythonInterp) SET(INTERNAL_NETCDF4_PYTHON ON) SET(INTERNAL_CYTHON ON) +SET(INTERNAL_HEALPY ON) IF(INTERNAL_CYTHON) @@ -12,6 +13,10 @@ IF(INTERNAL_NETCDF4_PYTHON) SET(NETCDF4_PYTHON_URL "http://netcdf4-python.googlecode.com/files/netCDF4-1.0.1.tar.gz" CACHE STRING "URL to download NetCDF4-python from") ENDIF(INTERNAL_NETCDF4_PYTHON) +IF (INTERNAL_HEALPY) + SET(HEALPY_URL "http://github.com/healpy/healpy/archive/1.4.1.tar.gz" CACHE STRING "URL to download Healpy from") +ENDIF(INTERNAL_HEALPY) + execute_process( COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_SOURCE_DIR}/external/detect_site.py ${CMAKE_BINARY_DIR}/ext_build/python RESULT_VARIABLE RET_VALUE @@ -67,6 +72,28 @@ IF(INTERNAL_NETCDF4_PYTHON) BUILD_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_build.cmake INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake ) + SET(PREV_PYTHON_BUILD ${PREV_PYTHON_BUILD} netcdf4-python) ENDIF(INTERNAL_NETCDF4_PYTHON) +IF(INTERNAL_HEALPY) + SET(BUILD_ENVIRONMENT + ${CMAKE_COMMAND} + "-DPYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}" + "-DPYTHON_CPPFLAGS:STRING=${PYTHON_CPPFLAGS}" + "-DCFITSIO_EXT_LIB=${CFITSIO_LIBRARY}" + "-DCFITSIO_EXT_INC=${CFITSIO_INCLUDE_PATH}" + "-DNETCDF4_DIR=${NETCDF_BIN_DIR}" + "-DPYTHON_LDFLAGS:STRING=${PYTHON_LDFLAGS}" + "-DPYTHON_LOCAL_SITE_PACKAGE=${PYTHON_LOCAL_SITE_PACKAGE}" + "-DTARGET_PATH=${CMAKE_BINARY_DIR}/ext_build/python" "-P") + ExternalProject_Add(healpy + DEPENDS ${PREV_PYTHON_BUILD} + URL ${HEALPY_URL} + PREFIX ${BUILD_PREFIX}/healpy-prefix + CONFIGURE_COMMAND echo "No configure" + BUILD_IN_SOURCE 1 + BUILD_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_build.cmake + INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake + ) +ENDIF(INTERNAL_HEALPY) diff --git a/external/python_build.cmake b/external/python_build.cmake index f3de8de..ca9970e 100644 --- a/external/python_build.cmake +++ b/external/python_build.cmake @@ -4,6 +4,9 @@ SET(ENV{CPPFLAGS} ${PYTHON_CPPFLAGS}) SET(ENV{LDFLAGS} ${PYTHON_LDFLAGS}) SET(ENV{VOID_GSL} ${VOID_GSL}) SET(ENV{PYTHONPATH} ${PYTHON_LOCAL_SITE_PACKAGE}:$ENV{PYTHONPATH}) +SET(ENV{CFITSIO_EXT_INC} ${CFITSIO_EXT_INC}) +SET(ENV{CFITSIO_EXT_LIB} ${CFITSIO_EXT_LIB}) + SET(PYTHON_BUILD_COMMAND ${PYTHON_EXECUTABLE} setup.py build) MESSAGE(STATUS "Running ${PYTHON_BUILD_COMMAND}") execute_process( diff --git a/external/python_install.cmake b/external/python_install.cmake index 5afa3cb..e33e377 100644 --- a/external/python_install.cmake +++ b/external/python_install.cmake @@ -3,8 +3,12 @@ SET(ENV{NETCDF4_DIR} ${NETCDF4_DIR}) SET(ENV{CPPFLAGS} ${PYTHON_CPPFLAGS}) SET(ENV{LDFLAGS} ${PYTHON_LDFLAGS}) SET(ENV{VOID_GSL} ${VOID_GSL}) +SET(ENV{CFITSIO_EXT_INC} ${CFITSIO_EXT_INC}) +SET(ENV{CFITSIO_EXT_LIB} ${CFITSIO_EXT_LIB}) SET(ENV{PYTHONPATH} ${PYTHON_LOCAL_SITE_PACKAGE}:$ENV{PYTHONPATH}) + SET(PYTHON_INSTALL_COMMAND ${PYTHON_EXECUTABLE} setup.py install --prefix=${TARGET_PATH} --install-lib=${PYTHON_LOCAL_SITE_PACKAGE}) + message(STATUS "Running ${PYTHON_INSTALL_COMMAND}") execute_process( COMMAND ${PYTHON_INSTALL_COMMAND} From 4764b42a63c95b94bf86d89cd4dda4d82424b825 Mon Sep 17 00:00:00 2001 From: nhamaus Date: Tue, 6 Nov 2012 17:25:29 +0100 Subject: [PATCH 04/32] changed baseResolution to float --- pipeline/prepareCatalogs.py | 2 +- pipeline/prepareGadgetCatalog.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index c2bb424..359826f 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -227,7 +227,7 @@ if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir) # first the directly downsampled runs # Note: ss0.002 ~ SDSS DR7 dim2 # ss0.0004 ~ SDSS DR9 mid -baseResolution = numPart/lbox/lbox/lbox # particles/Mpc^3 +baseResolution = float(numPart)/lbox/lbox/lbox # particles/Mpc^3 for thisSubSample in subSamples: keepFraction = float(thisSubSample) / baseResolution diff --git a/pipeline/prepareGadgetCatalog.py b/pipeline/prepareGadgetCatalog.py index 2fc83f3..1f8342f 100755 --- a/pipeline/prepareGadgetCatalog.py +++ b/pipeline/prepareGadgetCatalog.py @@ -220,7 +220,7 @@ if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir) # first the directly downsampled runs # Note: ss0.002 ~ SDSS DR7 dim2 # ss0.0004 ~ SDSS DR9 mid -baseResolution = numPart/lbox/lbox/lbox # particles/Mpc^3 +baseResolution = float(numPart)/lbox/lbox/lbox # particles/Mpc^3 for thisSubSample in subSamples: keepFraction = float(thisSubSample) / baseResolution From c9d3a5540eee72536a6ab23e0ca06245ca8ffd70 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 6 Nov 2012 13:25:25 -0500 Subject: [PATCH 05/32] Added argparse/setuptools. It is temporarily set to always on while I figure a way of detection python version --- external/external_python_build.cmake | 39 +++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/external/external_python_build.cmake b/external/external_python_build.cmake index e4580ae..d934e17 100644 --- a/external/external_python_build.cmake +++ b/external/external_python_build.cmake @@ -3,20 +3,29 @@ INCLUDE(FindPythonInterp) SET(INTERNAL_NETCDF4_PYTHON ON) SET(INTERNAL_CYTHON ON) SET(INTERNAL_HEALPY ON) - +SET(INTERNAL_ARGPARSE ON) IF(INTERNAL_CYTHON) SET(CYTHON_URL "http://cython.org/release/Cython-0.17.1.tar.gz" CACHE STRING "URL to download Cython from") + mark_as_advanced(CYTHON_URL) ENDIF(INTERNAL_CYTHON) IF(INTERNAL_NETCDF4_PYTHON) SET(NETCDF4_PYTHON_URL "http://netcdf4-python.googlecode.com/files/netCDF4-1.0.1.tar.gz" CACHE STRING "URL to download NetCDF4-python from") + mark_as_advanced(NETCDF4_PYTHON_URL) ENDIF(INTERNAL_NETCDF4_PYTHON) IF (INTERNAL_HEALPY) SET(HEALPY_URL "http://github.com/healpy/healpy/archive/1.4.1.tar.gz" CACHE STRING "URL to download Healpy from") + mark_as_advanced(HEALPY_URL) ENDIF(INTERNAL_HEALPY) +IF(INTERNAL_ARGPARSE) + SET(SETUPTOOLS_URL "http://pypi.python.org/packages/source/s/setuptools/setuptools-0.6c11.tar.gz" CACHE STRING "URL to download setuptools from") + SET(ARGPARSE_URL "http://argparse.googlecode.com/files/argparse-1.2.1.tar.gz" CACHE STRING "URL to download argparse from") + mark_as_advanced(ARGPARSE_URL SETUPTOOLS_URL) +ENDIF(INTERNAL_ARGPARSE) + execute_process( COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_SOURCE_DIR}/external/detect_site.py ${CMAKE_BINARY_DIR}/ext_build/python RESULT_VARIABLE RET_VALUE @@ -97,3 +106,31 @@ IF(INTERNAL_HEALPY) INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake ) ENDIF(INTERNAL_HEALPY) + +IF(INTERNAL_ARGPARSE) + SET(BUILD_ENVIRONMENT + ${CMAKE_COMMAND} + "-DPYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}" + "-DPYTHON_LOCAL_SITE_PACKAGE=${PYTHON_LOCAL_SITE_PACKAGE}" + "-DTARGET_PATH=${CMAKE_BINARY_DIR}/ext_build/python" "-P") + + ExternalProject_Add(setuptools + URL ${SETUPTOOLS_URL} + PREFIX ${BUILD_PREFIX}/setuptools-prefix + CONFIGURE_COMMAND echo "No configure" + BUILD_IN_SOURCE 1 + BUILD_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_build.cmake + INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake + ) + + ExternalProject_Add(argparse + DEPENDS setuptools + URL ${ARGPARSE_URL} + PREFIX ${BUILD_PREFIX}/argparse-prefix + CONFIGURE_COMMAND echo "No configure" + BUILD_IN_SOURCE 1 + BUILD_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_build.cmake + INSTALL_COMMAND ${BUILD_ENVIRONMENT} ${CMAKE_SOURCE_DIR}/external/python_install.cmake + ) + SET(AUXILIARY_PYTHON_DEPEND ${AUXILIARY_PYTHON_DEPEND} argparse) +ENDIF(INTERNAL_ARGPARSE) From 0f53c385b8efa97cadce68b100eff56f1e188ad7 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 6 Nov 2012 15:36:15 -0500 Subject: [PATCH 06/32] Cleaned up shown variables by marking most URLs as 'advanced variables' --- external/external_build.cmake | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/external/external_build.cmake b/external/external_build.cmake index c28c6db..fd407e4 100644 --- a/external/external_build.cmake +++ b/external/external_build.cmake @@ -22,28 +22,34 @@ SET(INTERNAL_QHULL ON) IF(INTERNAL_GENGETOPT) SET(GENGETOPT_URL "ftp://ftp.gnu.org/gnu/gengetopt/gengetopt-2.22.5.tar.gz" CACHE STRING "URL to download gengetopt from") + mark_as_advanced(GENGETOPT_URL) ENDIF(INTERNAL_GENGETOPT) IF(INTERNAL_HDF5) SET(HDF5_URL "http://www.hdfgroup.org/ftp/HDF5/current/src/hdf5-1.8.9.tar.gz" CACHE STRING "URL to download HDF5 from") + mark_as_advanced(HDF5_URL) ENDIF(INTERNAL_HDF5) IF(INTERNAL_NETCDF) SET(NETCDF_URL "http://www.unidata.ucar.edu/downloads/netcdf/ftp/netcdf-4.1.3.tar.gz" CACHE STRING "URL to download NetCDF from") + mark_as_advanced(NETCDF_URL) ENDIF(INTERNAL_NETCDF) IF(INTERNAL_BOOST) SET(BOOST_URL "http://sourceforge.net/projects/boost/files/boost/1.49.0/boost_1_49_0.tar.gz/download" CACHE STRING "URL to download Boost from") + mark_as_advanced(BOOST_URL) ELSE(INTERNAL_BOOST) find_package(Boost 1.49.0 COMPONENTS format spirit phoenix python FATAL_ERROR) ENDIF(INTERNAL_BOOST) IF(INTERNAL_GSL) SET(GSL_URL "ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz" CACHE STRING "URL to download GSL from ") + mark_as_advanced(GSL_URL) ENDIF(INTERNAL_GSL) IF(INTERNAL_QHULL) SET(QHULL_URL "http://www.qhull.org/download/qhull-2012.1-src.tgz" CACHE STRING "URL to download QHull from") + mark_as_advanced(QHULL_URL) ENDIF(INTERNAL_QHULL) @@ -72,10 +78,11 @@ if (INTERNAL_GENGETOPT) BUILD_IN_SOURCE 1 INSTALL_COMMAND make install ) - SET(GENGETOPT ${GENGETOPT_BIN_DIR}/bin/gengetopt) + SET(GENGETOPT ${GENGETOPT_BIN_DIR}/bin/gengetopt CACHE FILEPATH "Path GenGetOpt binary") else(INTERNAL_GENGETOPT) find_program(GENGETOPT gengetopt) endif(INTERNAL_GENGETOPT) +mark_as_advanced(GENGETOPT) ############### # Build HDF5 @@ -113,6 +120,7 @@ else(INTERNAL_HDF5) find_library(HDF5HL_LIBRARY hdf5_hl) endif (INTERNAL_HDF5) SET(CONFIGURE_CPP_FLAGS "${CONFIGURE_CPP_FLAGS} -I${HDF5_INCLUDE_PATH}") +mark_as_advanced(HDF5_INCLUDE_PATH HDF5_LIBRARY HDF5_CPP_LIBRARY HDF5HL_LIBRARY HDF5HL_CPP_LIBRARY) ############### # Build NetCDF @@ -154,6 +162,7 @@ ELSE(INTERNAL_NETCDF) SET(CONFIGURE_CPP_FLAGS ${CONFIGURE_CPP_FLAGS} -I${NETCDF_INCLUDE_PATH} -I${NETCDFCPP_INCLUDE_PATH}) endif (INTERNAL_NETCDF) +mark_as_advanced(NETCDF_LIBRARY NETCDFCPP_LIBRARY NETCDF_INCLUDE_PATH NETCDFCPP_INCLUDE_PATH) ################## # Build BOOST @@ -171,8 +180,9 @@ if (INTERNAL_BOOST) INSTALL_COMMAND echo "No install" ) set(Boost_INCLUDE_DIRS ${BOOST_SOURCE_DIR} CACHE STRING "Boost path" FORCE) - set(Boost_LIBRARIES ${BOOST_SOURCE_DIR}/stage/lib/libboost_python.a) + set(Boost_LIBRARIES ${BOOST_SOURCE_DIR}/stage/lib/libboost_python.a CACHE STRING "Boost libraries" FORCE) endif (INTERNAL_BOOST) +mark_as_advanced(Boost_INCLUDE_DIRS Boost_LIBRARIES) ################## # Build GSl @@ -201,6 +211,7 @@ ELSE(INTERNAL_GSL) find_library(GSLCBLAS_LIBRARY gslcblas) find_path(GSL_INCLUDE_PATH NAMES gsl/gsl_blas.h) ENDIF(INTERNAL_GSL) +mark_as_advanced(GSL_LIBRARY GSLCBLAS_LIBRARY GSL_INCLUDE_PATH) ################## # Build CosmoTool @@ -298,6 +309,7 @@ if (INTERNAL_QHULL) add_definitions(-Dqh_QHpointer) else(INTERNAL_QHULL) + message(FATAL_ERROR "Only packaged QHull is supported") endif(INTERNAL_QHULL) SET(QHULL_LIBRARIES ${QHULL_CPP_LIBRARY} ${QHULL_LIBRARY} ) From 61acd5938330ff3f9ec6e6114f0b00022427b5ff Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Tue, 6 Nov 2012 14:41:09 -0600 Subject: [PATCH 07/32] futher fixes for handling sub-volumes in analysis --- c_tools/libzobov/voidTree.hpp | 4 +- c_tools/mock/generateFromCatalog.cpp | 1 + c_tools/stacking/pruneVoids.cpp | 29 ++-- pipeline/generateCatalog.py | 5 +- pipeline/prepareCatalogs.py | 7 +- pipeline/prepareGadgetCatalog.py | 5 +- .../void_python_tools/backend/launchers.py | 3 +- .../void_python_tools/plotting/plotTools.py | 154 +++++++++++++++++- 8 files changed, 182 insertions(+), 26 deletions(-) diff --git a/c_tools/libzobov/voidTree.hpp b/c_tools/libzobov/voidTree.hpp index a316027..bb825dc 100644 --- a/c_tools/libzobov/voidTree.hpp +++ b/c_tools/libzobov/voidTree.hpp @@ -126,8 +126,8 @@ public: // Go bigger, though I would say we should not to. } while (iter_candidate != candidateList.begin()) ; - if (vid_good_candidate < 0) - std::cout << "Failure to lookup parent (2)" << std::endl; + //if (vid_good_candidate < 0) + // std::cout << "Failure to lookup parent (2)" << std::endl; return vid_good_candidate; } diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index f7efcc7..febc6d3 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/c_tools/mock/generateFromCatalog.cpp @@ -414,6 +414,7 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn for (uint32_t i = 0; i < pdata.pos.size(); i++) { f.writeReal32((pdata.pos[i].xyz[j]+Lmax)/(2*Lmax)); +if (i < 10) printf("TEST WRITE %d %e\n", (pdata.pos[i].xyz[j]+Lmax)/(2*Lmax)); } f.endCheckpoint(); } diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 4580e60..d6d5890 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -254,8 +254,9 @@ int main(int argc, char **argv) { for (iVoid = 0; iVoid < numVoids; iVoid++) { voidID = voids[iVoid].voidID; - //printf(" DOING %d (of %d) %d %d\n", iVoid, numVoids, voidID, - // voids[iVoid].numPart); + printf(" DOING %d (of %d) %d %d %f\n", iVoid, numVoids, voidID, + voids[iVoid].numPart, + voids[iVoid].radius); voids[iVoid].center[0] = part[voids[iVoid].coreParticle].x; voids[iVoid].center[1] = part[voids[iVoid].coreParticle].y; @@ -411,30 +412,36 @@ int main(int argc, char **argv) { } // iVoid printf(" Picking winners and losers...\n"); - numKept = numVoids; + for (iVoid = 0; iVoid < numVoids; iVoid++) { + voids[iVoid].accepted = 1; + } + for (iVoid = 0; iVoid < numVoids; iVoid++) { if (strcmp(args_info.dataPortion_arg, "edge") == 0 && tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) { - numKept--; voids[iVoid].accepted = 0; } if (strcmp(args_info.dataPortion_arg, "central") == 0 && tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestMock) { - numKept--; + voids[iVoid].accepted = 0; + } + + if (voids[iVoid].radius < args_info.rMin_arg) { voids[iVoid].accepted = 0; } if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) { - numKept--; voids[iVoid].accepted = -1; } - if (voids[iVoid].radius < args_info.rMin_arg) { - numKept--; - voids[iVoid].accepted = 0; - } } + + numKept = 0; + for (iVoid = 0; iVoid < numVoids; iVoid++) { + if (voids[iVoid].accepted == 1) numKept++; + } + printf(" Number kept: %d (out of %d)\n", numKept, numVoids); printf(" Output...\n"); @@ -452,7 +459,7 @@ int main(int argc, char **argv) { if (voids[iVoid].accepted != 1) continue; fprintf(fp, "%d %d %d %f %f %d %d %f %d %f %f\n", - i, + iVoid, voids[iVoid].voidID, voids[iVoid].coreParticle, voids[iVoid].coreDens, diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 025bbc2..7f95f66 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -107,6 +107,9 @@ if (startCatalogStage <= 4) and (endCatalogStage >= 4): sys.stdout.flush() for thisDataPortion in dataPortions: - plotNumberCounts(workDir, dataSampleList, figDir, showPlot=True, dataPortion=thisDataPortion, setName=setName) + plotRedshiftDistribution(workDir, dataSampleList, figDir, showPlot=False, + dataPortion=thisDataPortion, setName=catalogName) + plotSizeDistribution(workDir, dataSampleList, figDir, showPlot=False, + dataPortion=thisDataPortion, setName=catalogName) print "\n Done!" diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index f45e31b..af9d997 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -7,7 +7,6 @@ import numpy as np import os import sys -from Scientific.IO.NetCDF import NetCDFFile import void_python_tools as vp import argparse @@ -95,10 +94,10 @@ from void_python_tools.backend.classes import * continueRun = False startCatalogStage = 1 -endCatalogStage = 3 +endCatalogStage = 4 startAPStage = 1 -endAPStage = 6 +endAPStage = 7 ZOBOV_PATH = os.getenv("PWD")+"/../zobov/" CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/" @@ -297,7 +296,7 @@ if args.script or args.all: for line in inFile: numPart += 1 inFile.close() - minRadius = int(np.ceil(lbox/numPart**(1./3.))) + minRadius = 2*int(np.ceil(lbox/numPart**(1./3.))) if dataFormat == "multidark": writeScript(prefix, "halos", scriptDir, catalogDir, redshifts, diff --git a/pipeline/prepareGadgetCatalog.py b/pipeline/prepareGadgetCatalog.py index 2fc83f3..f0548fa 100755 --- a/pipeline/prepareGadgetCatalog.py +++ b/pipeline/prepareGadgetCatalog.py @@ -5,7 +5,6 @@ import numpy as np import os import sys -from Scientific.IO.NetCDF import NetCDFFile import void_python_tools as vp import argparse @@ -96,10 +95,10 @@ from void_python_tools.backend.classes import * continueRun = False startCatalogStage = 1 -endCatalogStage = 3 +endCatalogStage = 4 startAPStage = 1 -endAPStage = 6 +endAPStage = 7 ZOBOV_PATH = os.getenv("PWD")+"/../zobov/" CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/" diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 92536cd..e75343e 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -76,7 +76,6 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, if os.access("contour_map.fits", os.F_OK): os.system("mv %s %s" % ("contour_map.fits", zobovDir)) - os.system("mv %s %s" % ("mask_map.fits", zobovDir)) if os.access("comoving_distance.txt", os.F_OK): os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) @@ -815,7 +814,7 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None, maxtries = 5 while badChain: Rexpect = (stack.rMin+stack.rMax)/2 - Rtruncate = stack.rMax*3. + 1 # TEST + Rtruncate = stack.rMin*3. + 1 # TEST ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, Niter=300000, Nburn=100000, diff --git a/python_tools/void_python_tools/plotting/plotTools.py b/python_tools/void_python_tools/plotting/plotTools.py index 4fad79a..a7c5058 100644 --- a/python_tools/void_python_tools/plotting/plotTools.py +++ b/python_tools/void_python_tools/plotting/plotTools.py @@ -1,4 +1,5 @@ -__all__=['plotNumberCounts'] +__all__=['plotRedshiftDistribution', 'plotSizeDistribution', 'plot1dProfiles', + 'plotMarg1d'] from void_python_tools.backend.classes import * from plotDefs import * @@ -7,11 +8,12 @@ import os import pylab as plt # ----------------------------------------------------------------------------- -def plotNumberCounts(workDir=None, sampleList=None, figDir=None, plotNameBase="numbercount", +def plotRedshiftDistribution(workDir=None, sampleList=None, figDir=None, + plotNameBase="zdist", showPlot=False, dataPortion=None, setName=None): plt.clf() - plt.xlabel("Comoving Distance (Mpc/h)") + plt.xlabel("Redshift") plt.ylabel("Number of Voids") plotTitle = setName @@ -67,3 +69,149 @@ def plotNumberCounts(workDir=None, sampleList=None, figDir=None, plotNameBase="n os.system("display %s" % figDir+"/fig_"+plotName+".png") +# ----------------------------------------------------------------------------- +def plotSizeDistribution(workDir=None, sampleList=None, figDir=None, + plotNameBase="sizedist", + showPlot=False, dataPortion=None, setName=None): + + plt.clf() + plt.xlabel("Void Radius (Mpc/h)") + plt.ylabel("Number of Voids") + + plotTitle = setName + + plotName = plotNameBase + + xMin = 1.e00 + xMax = 0 + + for (iSample,sample) in enumerate(sampleList): + + sampleName = sample.fullName + lineTitle = sampleName + + filename = workDir+"/sample_"+sampleName+"/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 + + xMin = 5 + xMax = 140 + + range = (xMin, xMax) + nbins = np.ceil((xMax-xMin)/5) + + #thisMax = np.max(data[:,5]) + #thisMin = np.min(data[:,5]) + #if thisMax > xMax: xMax = thisMax + #if thisMin < xMin: xMin = thisMin + + plt.hist(data[:,4], bins=nbins, + label=lineTitle, color=colorList[iSample], + histtype = "step", range=range, + linewidth=linewidth) + + plt.legend(title = "Samples", loc = "upper right") + plt.title(plotTitle) + + plt.xlim(xMin, xMax) + #plt.xlim(xMin, xMax*1.4) # make room for legend + + 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") + + +# ----------------------------------------------------------------------------- +def plot1dProfiles(workDir=None, sampleList=None, figDir=None, + plotNameBase="1dprofile", + showPlot=False, dataPortion=None, setName=None): + + plt.clf() + plt.xlabel(r"$R/R_{v,\mathrm{max}}$") + plt.ylabel(r"$n / \bar n$") + + for (iSample,sample) in enumerate(sampleList): + + sampleName = sample.fullName + + for (iStack,stack) in enumerate(sample.stacks): + + plotTitle = setName + plotName = plotNameBase + + runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, + stack.rMax, dataPortion) + + plotTitle = sampleName + ", z = "+str(stack.zMin)+"-"+str(stack.zMax)+", R = "+str(stack.rMin)+"-"+str(stack.rMax)+ r" $h^{-1}$ Mpc" + + filename = workDir+"/sample_"+sampleName+"/stacks_"+runSuffix+\ + "profile_1d.txt" + + 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[:,1] /= stack.rMax + plt.ylim(ymin=0.0, ymax=np.amax(data[:,2])+0.1) + plt.xlim(xmin=0.0, xmax=2.1) + plt.plot(data[:,1], data[:,2], label=lineTitle, color=colorList[0], + linewidth=linewidth) + + 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 showPlot: + os.system("display %s" % figDir+"/fig_"+plotName+".png") + + +# ----------------------------------------------------------------------------- +def plotMarg1d(workDir=None, sampleList=None, figDir=None, + plotNameBase="marg1d", + showPlot=False, dataPortion=None, setName=None): + + plotNames = ("Om", "w0", "wa") + plotTitles = ("$\Omega_M$", "$w_0$", "$w_a$") + files = ("Om", "w0", "wa") + + for iPlot in range(len(plotNames)): + + plt.clf() + + plotName = plotNameBase+"_"+plotNames[iPlot]+"_"+dataPortion + plotTitle = plotTitles[iPlot] + dataFile = workDir + "/likelihoods_"+dataPortion+"_"+files[iPlot]+".dat" + + plt.xlabel(plotTitle, fontsize="20") + plt.ylabel("Likelihood", fontsize="20") + plt.ylim(0.0, 1.0) + + data = np.loadtxt(dataFile, comments="#") + plt.plot(data[:,0], data[:,1], color='k', linewidth=linewidth) + + 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") + + + From 15f03654d72e367d6f84339d4550032eaae31686 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 6 Nov 2012 15:42:11 -0500 Subject: [PATCH 08/32] Check if python < 2.7, and build argparse if necessary --- external/external_python_build.cmake | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/external/external_python_build.cmake b/external/external_python_build.cmake index d934e17..d6b0a2f 100644 --- a/external/external_python_build.cmake +++ b/external/external_python_build.cmake @@ -3,7 +3,13 @@ INCLUDE(FindPythonInterp) SET(INTERNAL_NETCDF4_PYTHON ON) SET(INTERNAL_CYTHON ON) SET(INTERNAL_HEALPY ON) -SET(INTERNAL_ARGPARSE ON) + +IF (PYTHON_VERSION_STRING VERSION_LESS 2.7) + MESSAGE(STATUS "Python version is less than 2.7, argparse is needed.") + SET(INTERNAL_ARGPARSE ON) +ELSE (PYTHON_VERSION_STRING VERSION_LESS 2.7) + MESSAGE(STATUS "Python version is greater than 2.7, argparse is already bundled.") +ENDIF (PYTHON_VERSION_STRING VERSION_LESS 2.7) IF(INTERNAL_CYTHON) SET(CYTHON_URL "http://cython.org/release/Cython-0.17.1.tar.gz" CACHE STRING "URL to download Cython from") From 058b12bb3f52cd182c3d395e48f705cf29cc95e8 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 6 Nov 2012 15:56:03 -0500 Subject: [PATCH 09/32] healpy setup.py also needs to know where cfitsio is located --- external/external_build.cmake | 5 +++-- external/external_python_build.cmake | 1 + external/python_build.cmake | 1 + external/python_install.cmake | 1 + 4 files changed, 6 insertions(+), 2 deletions(-) diff --git a/external/external_build.cmake b/external/external_build.cmake index fd407e4..8835789 100644 --- a/external/external_build.cmake +++ b/external/external_build.cmake @@ -253,8 +253,9 @@ ExternalProject_Add(cfitsio BUILD_IN_SOURCE 1 INSTALL_COMMAND make install ) -SET(CFITSIO_LIBRARY ${CMAKE_BINARY_DIR}/ext_build/cfitsio/lib/libcfitsio.a) -SET(CFITSIO_INCLUDE_PATH ${CMAKE_BINARY_DIR}/ext_build/cfitsio/include) +SET(CFITSIO_PREFIX ${CMAKE_BINARY_DIR}/ext_build/cfitsio) +SET(CFITSIO_LIBRARY ${CFITSIO_PREFIX}/lib/libcfitsio.a) +SET(CFITSIO_INCLUDE_PATH ${CFITSIO_PREFIX}/include) ################# # Build Healpix diff --git a/external/external_python_build.cmake b/external/external_python_build.cmake index d6b0a2f..6e78b42 100644 --- a/external/external_python_build.cmake +++ b/external/external_python_build.cmake @@ -97,6 +97,7 @@ IF(INTERNAL_HEALPY) "-DPYTHON_CPPFLAGS:STRING=${PYTHON_CPPFLAGS}" "-DCFITSIO_EXT_LIB=${CFITSIO_LIBRARY}" "-DCFITSIO_EXT_INC=${CFITSIO_INCLUDE_PATH}" + "-DCFITSIO_EXT_PREFIX=${CFITSIO_PREFIX}" "-DNETCDF4_DIR=${NETCDF_BIN_DIR}" "-DPYTHON_LDFLAGS:STRING=${PYTHON_LDFLAGS}" "-DPYTHON_LOCAL_SITE_PACKAGE=${PYTHON_LOCAL_SITE_PACKAGE}" diff --git a/external/python_build.cmake b/external/python_build.cmake index ca9970e..b834e9b 100644 --- a/external/python_build.cmake +++ b/external/python_build.cmake @@ -6,6 +6,7 @@ SET(ENV{VOID_GSL} ${VOID_GSL}) SET(ENV{PYTHONPATH} ${PYTHON_LOCAL_SITE_PACKAGE}:$ENV{PYTHONPATH}) SET(ENV{CFITSIO_EXT_INC} ${CFITSIO_EXT_INC}) SET(ENV{CFITSIO_EXT_LIB} ${CFITSIO_EXT_LIB}) +SET(ENV{CFITSIO_EXT_PREFIX} ${CFITSIO_EXT_PREFIX}) SET(PYTHON_BUILD_COMMAND ${PYTHON_EXECUTABLE} setup.py build) MESSAGE(STATUS "Running ${PYTHON_BUILD_COMMAND}") diff --git a/external/python_install.cmake b/external/python_install.cmake index e33e377..ee04a3c 100644 --- a/external/python_install.cmake +++ b/external/python_install.cmake @@ -4,6 +4,7 @@ SET(ENV{CPPFLAGS} ${PYTHON_CPPFLAGS}) SET(ENV{LDFLAGS} ${PYTHON_LDFLAGS}) SET(ENV{VOID_GSL} ${VOID_GSL}) SET(ENV{CFITSIO_EXT_INC} ${CFITSIO_EXT_INC}) +SET(ENV{CFITSIO_EXT_PREFIX} ${CFITSIO_EXT_PREFIX}) SET(ENV{CFITSIO_EXT_LIB} ${CFITSIO_EXT_LIB}) SET(ENV{PYTHONPATH} ${PYTHON_LOCAL_SITE_PACKAGE}:$ENV{PYTHONPATH}) From de5a7cb354a1da4d4c7e5e102e5e568be5c81407 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Wed, 7 Nov 2012 09:27:11 -0600 Subject: [PATCH 10/32] fix to have separate file name and redshift lists in script generation --- pipeline/prepareGadgetCatalog.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/pipeline/prepareGadgetCatalog.py b/pipeline/prepareGadgetCatalog.py index e3e1ffa..0d4ac57 100755 --- a/pipeline/prepareGadgetCatalog.py +++ b/pipeline/prepareGadgetCatalog.py @@ -33,10 +33,13 @@ dataFormat = "gadget" useLightCone = False # common filename of particle files -redshiftFileBase = "mdr1_particles_z" +particleFileBase = "mdr1_particles_z" -# list of redshifts for the particle files -# to get particle file name, we take redshiftFileBase+redshift +# list of file numbers for the particle files +# to get particle file name, we take particleFileBase+fileNum +fileNums = (("0.0", "0.53", "1.0")) + +# redshift of each file in the above list redshifts = (("0.0", "0.53", "1.0")) # how many independent slices along the z-axis? @@ -80,7 +83,7 @@ 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, redshifts, numSubvolumes, + scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, useVel, lbox, minRadius, omegaM, subsample=1.0, suffix=".dat"): @@ -149,7 +152,8 @@ newSample = Sample(dataFile = "{dataFile}", dataSampleList.append(newSample) """ - for redshift in redshifts: + for (iFile, redshift) in enumerate(redshifts): + fileNum = fileNums[iFile] zBox = float(redshift) Om = float(omegaM) zBoxMpc = LIGHT_SPEED/100.*vp.angularDiameter(zBox, Om=Om) @@ -181,7 +185,7 @@ dataSampleList.append(newSample) sliceMinMpc = "%0.1f" % sliceMinMpc sliceMaxMpc = "%0.1f" % sliceMaxMpc - dataFileName = dataFileNameBase + redshift + suffix + dataFileName = dataFileNameBase + fileNum + suffix for iX in xrange(numSubvolumes): for iY in xrange(numSubvolumes): @@ -228,10 +232,12 @@ for thisSubSample in subSamples: print " Doing subsample", thisSubSample, " scripts" setName = prefix+"ss"+str(thisSubSample) - writeScript(setName, redshiftFileBase, scriptDir, catalogDir, redshifts, + writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums, + redshifts, numSubvolumes, numSlices, True, lbox, minRadius, omegaM, subsample=thisSubSample, suffix="") - writeScript(setName, redshiftFileBase, scriptDir, catalogDir, redshifts, + writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums, + redshifts, numSubvolumes, numSlices, False, lbox, minRadius, omegaM, subsample=thisSubSample, suffix="") From 053c46e57b094c94fe0c2fa0430e934be34c2aed Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Wed, 7 Nov 2012 11:51:51 -0600 Subject: [PATCH 11/32] fix to allow switching on/off cosmological distortions --- c_tools/mock/generateMock.ggo | 2 +- python_tools/void_python_tools/backend/launchers.py | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/c_tools/mock/generateMock.ggo b/c_tools/mock/generateMock.ggo index 7148f48..7ca7158 100644 --- a/c_tools/mock/generateMock.ggo +++ b/c_tools/mock/generateMock.ggo @@ -26,6 +26,6 @@ option "rangeZ_max" - "Maximum range in Z for making the box (after distorti option "preReShift" - "Reshift the zero of the Z axis" flag off option "peculiarVelocities" - "Added peculiar velocities distortion" flag off -option "cosmo" - "Apply cosmological redshift" flag on +option "cosmo" - "Apply cosmological redshift" flag off option "subsample" - "Subsample the input simulation by the specified amount" double optional diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index e75343e..8f9e397 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -101,6 +101,11 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, else: includePecVelString = "" + if sample.useLightCone: + useLightConeString = "cosmo" + else: + useLightConeString = "" + if sample.dataFormat == "multidark": dataFileLine = "multidark " + datafile elif sample.dataFormat == "gadget": @@ -119,6 +124,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, output %s outputParameter %s %s + %s rangeX_min %g rangeX_max %g rangeY_min %g @@ -129,6 +135,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, """ % (dataFileLine, zobovDir+"/zobov_slice_"+sampleName, zobovDir+"/zobov_slice_"+sampleName+".par", includePecVelString, + useLightConeString, xMin, xMax, yMin, yMax, sample.zBoundaryMpc[0], sample.zBoundaryMpc[1], sample.subsample) From 9b8ce57a97aa731192f25dae58a716b83ca6571c Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 8 Nov 2012 10:00:08 -0500 Subject: [PATCH 12/32] Imported code for flexible gadget units. Imported load ramses code. Always invoke metricTransform. Add the option to import box parameters from an already existing mock catalog --- c_tools/mock/generateMock.cpp | 166 +++++++++++++++++++++++++++++----- c_tools/mock/generateMock.ggo | 4 + 2 files changed, 149 insertions(+), 21 deletions(-) diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index 797ca65..bfdf82f 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -18,9 +18,67 @@ using namespace CosmoTool; #define LIGHT_SPEED 299792.458 +static double gadgetUnit=1e-3; + +SimuData *doLoadRamses(const char *basename, int baseid, int velAxis, bool goRedshift) +{ + SimuData *d, *outd; + + d = loadRamsesSimu(basename, baseid, -1, true, 0); + outd = new SimuData; + + outd->NumPart = d->TotalNumPart; + outd->BoxSize = d->BoxSize; + outd->TotalNumPart = outd->NumPart; + outd->Hubble = d->Hubble; + outd->Omega_Lambda = d->Omega_Lambda; + outd->Omega_M = d->Omega_M; + outd->time = d->time; + + for (int k = 0; k < 3; k++) + outd->Pos[k] = new float[outd->NumPart]; + outd->Vel[2] = new float[outd->NumPart]; + delete d; + + int curCpu = 0; + cout << "loading cpu 0 " << endl; + while (d = loadRamsesSimu(basename, baseid, curCpu, true, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID)) + { + for (int k = 0; k < 3; k++) + for (int i = 0; i < d->NumPart; i++) + { + assert(d->Id[i] >= 1); + assert(d->Id[i] <= outd->TotalNumPart); + outd->Pos[k][d->Id[i]-1] = d->Pos[k][i]; + outd->Vel[2][d->Id[i]-1] = d->Vel[velAxis][i]; + } + + if (goRedshift) + for (int i = 0; i < d->NumPart; i++) + outd->Pos[velAxis][d->Id[i]-1] += d->Vel[velAxis][i]/100.; + + delete d; + curCpu++; + cout << "loading cpu " << curCpu << endl; + } + + return outd; +} + + + SimuData *myLoadGadget(const char *fname, int id, int flags) { - return loadGadgetMulti(fname, id, flags); + SimuData *sim = loadGadgetMulti(fname, id, flags); + sim->BoxSize *= gadgetUnit*1000; + for (int j = 0; j < 3; j++) + { + if (sim->Pos[j] != 0) { + for (long i = 0; i < sim->NumPart; i++) + sim->Pos[j][i] *= gadgetUnit*1000; + } + } + return sim; } SimuData *doLoadSimulation(const char *gadgetname, int velAxis, bool goRedshift, SimuData *(*loadFunction)(const char *fname, int id, int flags)) @@ -111,7 +169,7 @@ SimuData *doLoadMultidark(const char *multidarkname) fscanf(fp, "%f\n", &outd->Omega_M); fscanf(fp, "%f\n", &outd->Hubble); fscanf(fp, "%f\n", &outd->time); - fscanf(fp, "%d\n", &outd->NumPart); + fscanf(fp, "%ld\n", &outd->NumPart); outd->time = 1./(1.+outd->time); // convert to scale factor outd->TotalNumPart = outd->NumPart; @@ -181,7 +239,7 @@ Interpolate make_cosmological_redshift(double OM, double OL, double z0, double z return buildFromVector(pairs); } -void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact) +void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact, bool cosmo_flag) { int x0, x1, x2; @@ -200,7 +258,11 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double } double z0 = 1/data->time - 1; - Interpolate z_vs_D = make_cosmological_redshift(data->Omega_M, data->Omega_Lambda, 0., z0+8*data->BoxSize*100/LIGHT_SPEED); // Redshift 2*z0 should be sufficient ? + Interpolate z_vs_D = + make_cosmological_redshift(data->Omega_M, data->Omega_Lambda, + 0., z0+8*data->BoxSize*100/LIGHT_SPEED); + // Redshift 2*z0 should be sufficient ? This is fragile. + //A proper solver is needed here. double z_base = reshift ? z0 : 0; TotalExpansion e_computer; double baseComovingDistance; @@ -229,13 +291,15 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double // Distorted redshift if (reduced_red == 0) z = 0; - else { - z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.; - } + else if (cosmo_flag) + z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.; + else + z = reduced_red*LIGHT_SPEED/100.0; + expfact[i] = z / z_old; - // Add peculiar velocity - if (pecvel) - z += v/100; + // Add peculiar velocity + if (pecvel) + z += v/100; } catch(const InvalidRangeException& e) { cout << "Trying to interpolate out of the tabulated range." << endl; @@ -404,6 +468,63 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& delete[] expansion_fac; } +void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, generateMock_info& args_info) +{ + NcFile f(args_info.inputParameter_arg); + NcVar *v; + int *particle_id; + double *expansion_fac; + + boxed = new SimuData; + boxed->Hubble = simu->Hubble; + boxed->Omega_M = simu->Omega_M; + boxed->Omega_Lambda = simu->Omega_Lambda; + boxed->time = simu->time; + boxed->BoxSize = simu->BoxSize; + + NcVar *v_id = f.get_var("particle_ids"); + long *edges1; + double ranges[3][2]; + double mul[3]; + + edges1 = v_id->edges(); + assert(v_id->num_dims()==1); + + boxed->NumPart = edges1[0]; + delete[] edges1; + + particle_id = new int[boxed->NumPart]; + + v_id->get(particle_id, boxed->NumPart); + + ranges[0][0] = f.get_att("range_x_min")->as_double(0); + ranges[0][1] = f.get_att("range_x_max")->as_double(0); + ranges[1][0] = f.get_att("range_y_min")->as_double(0); + ranges[1][1] = f.get_att("range_y_max")->as_double(0); + ranges[2][0] = f.get_att("range_z_min")->as_double(0); + ranges[2][1] = f.get_att("range_z_max")->as_double(0); + + for (int j = 0; j < 3; j++) + { + boxed->Pos[j] = new float[boxed->NumPart]; + boxed->Vel[j] = 0; + mul[j] = 1.0/(ranges[j][1] - ranges[j][0]); + } + + uint32_t k = 0; + for (uint32_t i = 0; i < boxed->NumPart; i++) + { + int id = particle_id[i]; + + for (int j = 0; j < 3; j++) + { + boxed->Pos[j][i] = (simu->Pos[j][id]-ranges[j][0])*mul[j]; + } + } + + delete[] particle_id; +} + int main(int argc, char **argv) { generateMock_info args_info; @@ -434,13 +555,15 @@ int main(int argc, char **argv) generateMock_conf_print_version(); + gadgetUnit=args_info.gadgetUnit_arg; + if (args_info.ramsesBase_given || args_info.ramsesId_given) { if (args_info.ramsesBase_given && args_info.ramsesId_given) { - //simu = doLoadRamses(args_info.ramsesBase_arg, - // args_info.ramsesId_arg, - // args_info.axis_arg, false); - } + simu = doLoadRamses(args_info.ramsesBase_arg, + args_info.ramsesId_arg, + args_info.axis_arg, false); + } else { cerr << "Both ramsesBase and ramsesId are required to be able to load snapshots" << endl; @@ -490,14 +613,15 @@ int main(int argc, char **argv) double *expfact; - if (args_info.cosmo_flag){ - metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag, args_info.peculiarVelocities_flag, expfact); - } else { - expfact = new double[simu->NumPart]; - for (int j = 0; j < simu->NumPart; j++) expfact[j] = 1.0; - } + metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag, + args_info.peculiarVelocities_flag, expfact, + args_info.cosmo_flag); + + if (args_info.inputParameter_given) + makeBoxFromParameter(simu, expfact, simuOut, args_info); + else + makeBox(simu, expfact, simuOut, args_info); - makeBox(simu, expfact, simuOut, args_info); delete simu; generateOutput(simuOut, args_info.axis_arg, args_info.output_arg); diff --git a/c_tools/mock/generateMock.ggo b/c_tools/mock/generateMock.ggo index 7148f48..9ef988f 100644 --- a/c_tools/mock/generateMock.ggo +++ b/c_tools/mock/generateMock.ggo @@ -29,3 +29,7 @@ option "peculiarVelocities" - "Added peculiar velocities distortion" flag off option "cosmo" - "Apply cosmological redshift" flag on option "subsample" - "Subsample the input simulation by the specified amount" double optional + +option "inputParameter" - "Input geometry (optional, warning!)" string optional +option "gadgetUnit" - "Unit of length in gadget file in Mpc/h" double optional default="0.001" + From 95db063b7b97259ecce24bae154f0da9f5a048ae Mon Sep 17 00:00:00 2001 From: nhamaus Date: Thu, 8 Nov 2012 18:03:16 +0100 Subject: [PATCH 13/32] fixed gadget input unit to Mpc/h --- python_tools/void_python_tools/backend/launchers.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 8f9e397..6a094c8 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -110,7 +110,8 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, dataFileLine = "multidark " + datafile elif sample.dataFormat == "gadget": dataFileLine = "gadget " + datafile - + useGadgetUnit = "gadgetUnit=1" + iX = float(sample.mySubvolume[0]) iY = float(sample.mySubvolume[1]) @@ -125,6 +126,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, outputParameter %s %s %s + %s rangeX_min %g rangeX_max %g rangeY_min %g @@ -136,6 +138,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, zobovDir+"/zobov_slice_"+sampleName+".par", includePecVelString, useLightConeString, + useGadgetUnit, xMin, xMax, yMin, yMax, sample.zBoundaryMpc[0], sample.zBoundaryMpc[1], sample.subsample) From ed3f18f6d0ecbfde66fac15dede1b102f5f52253 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 8 Nov 2012 16:06:06 -0600 Subject: [PATCH 14/32] can now select gadget units in preparation script --- pipeline/prepareGadgetCatalog.py | 3 +++ python_tools/void_python_tools/backend/launchers.py | 5 ++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/pipeline/prepareGadgetCatalog.py b/pipeline/prepareGadgetCatalog.py index 0d4ac57..51f67cc 100755 --- a/pipeline/prepareGadgetCatalog.py +++ b/pipeline/prepareGadgetCatalog.py @@ -28,6 +28,7 @@ dataType = "simulation" # available formats for simulation: gadget, multidark dataFormat = "gadget" +dataUnit = 1 # as multiple of Mpc/h # place particles on the lightcone? useLightCone = False @@ -131,6 +132,7 @@ logDir = "{logDir}/{setName}/" sampleInfo = """ newSample = Sample(dataFile = "{dataFile}", dataFormat = "{dataFormat}", + dataUnit = {dataUnit}, fullName = "{sampleName}", nickName = "{sampleName}", dataType = "simulation", @@ -197,6 +199,7 @@ dataSampleList.append(newSample) scriptFile.write(sampleInfo.format(dataFile=dataFileName, dataFormat=dataFormat, + dataUnit=dataUnit, sampleName=sampleName, zMin=sliceMin, zMax=sliceMax, diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 6a094c8..52beefb 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -110,7 +110,6 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, dataFileLine = "multidark " + datafile elif sample.dataFormat == "gadget": dataFileLine = "gadget " + datafile - useGadgetUnit = "gadgetUnit=1" iX = float(sample.mySubvolume[0]) iY = float(sample.mySubvolume[1]) @@ -126,7 +125,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, outputParameter %s %s %s - %s + gadgetUnit %g rangeX_min %g rangeX_max %g rangeY_min %g @@ -138,7 +137,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, zobovDir+"/zobov_slice_"+sampleName+".par", includePecVelString, useLightConeString, - useGadgetUnit, + sample.gadgetUnit, xMin, xMax, yMin, yMax, sample.zBoundaryMpc[0], sample.zBoundaryMpc[1], sample.subsample) From 9c10c3714f8118fb0f0b4dbde9e96a2c9fb724e0 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 9 Nov 2012 07:10:08 -0600 Subject: [PATCH 15/32] fixed bug for data unit, and fixed bug for detecting completion of pruning --- python_tools/void_python_tools/backend/classes.py | 4 +++- python_tools/void_python_tools/backend/launchers.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/python_tools/void_python_tools/backend/classes.py b/python_tools/void_python_tools/backend/classes.py index 7466acb..60de62f 100755 --- a/python_tools/void_python_tools/backend/classes.py +++ b/python_tools/void_python_tools/backend/classes.py @@ -38,6 +38,7 @@ class Sample: dataType = "observation" dataFormat = "sdss" dataFile = "lss.dr72dim.dat" + dataUNit = 1 fullName = "lss.dr72dim.dat" nickName = "dim" zobovDir = "" @@ -68,7 +69,7 @@ class Sample: stacks = [] - def __init__(self, dataFile="", fullName="", + def __init__(self, dataFile="", fullName="", dataUnit=1, nickName="", maskFile="", selFunFile="", zBoundary=(), zRange=(), zBoundaryMpc=(), minVoidRadius=0, fakeDensity=0.01, volumeLimited=True, @@ -106,6 +107,7 @@ class Sample: self.dataFormat = dataFormat self.subsample = subsample self.useLightCone = useLightCone + self.dataUnit = dataUnit self.stacks = [] diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 52beefb..71d6ea3 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -287,7 +287,7 @@ def launchPrune(sample, binPath, thisDataPortion=None, cmd += " >& " + logFile os.system(cmd) - if jobSuccessful(logFile, "NetCDF: Not a valid ID\n"): + if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or jobSuccessful(logFile, "Done!"): print "done" else: print "FAILED!" From 9cdc037946624f0fa90d78087643f77bad6b2c7b Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 9 Nov 2012 08:06:27 -0600 Subject: [PATCH 16/32] some bugs fixed related to gadget unit handling --- pipeline/prepareCatalogs.py | 5 ++++- python_tools/void_python_tools/backend/classes.py | 2 +- python_tools/void_python_tools/backend/launchers.py | 4 ++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index bbddc2d..d331d85 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -19,11 +19,12 @@ figDir = os.getenv("PWD")+"/../figs/multidark/" logDir = os.getenv("PWD")+"/../logs/multidark/" dataFormat = "multidark" +dataUnit = 1 dataType = "simulation" useLightCone = True -redshifts = (("0.0", "0.53", "1.0")) redshiftFileBase = "mdr1_particles_z" +redshifts = (("0.0", "0.53", "1.0")) numSlices = 4 numSubvolumes = 1 @@ -128,6 +129,7 @@ logDir = "{logDir}/{sampleName}/" sampleInfo = """ newSample = Sample(dataFile = "{dataFile}", dataFormat = "{dataFormat}", + dataUnit = {dataUnit}, fullName = "{sampleName}", nickName = "{sampleName}", dataType = "simulation", @@ -202,6 +204,7 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False) scriptFile.write(sampleInfo.format(dataFile=dataFileName, dataFormat=dataFormat, + dataUnit=dataUnit, sampleName=sampleName, zMin=sliceMin, zMax=sliceMax, diff --git a/python_tools/void_python_tools/backend/classes.py b/python_tools/void_python_tools/backend/classes.py index 60de62f..572e543 100755 --- a/python_tools/void_python_tools/backend/classes.py +++ b/python_tools/void_python_tools/backend/classes.py @@ -78,7 +78,7 @@ class Sample: dataType="observation", numSubDivisions=2, boxLen=1024, usePecVel=False, omegaM=0.27, numSubvolumes=1, mySubvolume=1, dataFormat="sdss", - subsample="1.0", useLightCone=True): + subsample=1.0, useLightCone=True): self.dataFile = dataFile self.fullName = fullName self.nickName = nickName diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 71d6ea3..b8cb57f 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -132,12 +132,12 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, rangeY_max %g rangeZ_min %g rangeZ_max %g - subsample %g + subsample %s """ % (dataFileLine, zobovDir+"/zobov_slice_"+sampleName, zobovDir+"/zobov_slice_"+sampleName+".par", includePecVelString, useLightConeString, - sample.gadgetUnit, + sample.dataUnit, xMin, xMax, yMin, yMax, sample.zBoundaryMpc[0], sample.zBoundaryMpc[1], sample.subsample) From 93ba823764428a5f4be2f8aa960b94acecc41da6 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Tue, 13 Nov 2012 13:55:40 -0600 Subject: [PATCH 17/32] plotting now supports cumulative number distributions --- pipeline/generateCatalog.py | 2 + .../void_python_tools/plotting/plotDefs.py | 2 + .../void_python_tools/plotting/plotTools.py | 65 ++++++++++++++++++- 3 files changed, 67 insertions(+), 2 deletions(-) diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 7f95f66..cf4dc23 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -111,5 +111,7 @@ if (startCatalogStage <= 4) and (endCatalogStage >= 4): dataPortion=thisDataPortion, setName=catalogName) plotSizeDistribution(workDir, dataSampleList, figDir, showPlot=False, dataPortion=thisDataPortion, setName=catalogName) + plotNumberDistribution(workDir, dataSampleList, figDir, showPlot=False, + dataPortion=thisDataPortion, setName=catalogName) print "\n Done!" diff --git a/python_tools/void_python_tools/plotting/plotDefs.py b/python_tools/void_python_tools/plotting/plotDefs.py index 980e442..a999d41 100644 --- a/python_tools/void_python_tools/plotting/plotDefs.py +++ b/python_tools/void_python_tools/plotting/plotDefs.py @@ -1,3 +1,5 @@ +LIGHT_SPEED = 299792.458 + colorList = ['r', 'b', 'g', 'y', 'c', 'm', 'y', 'brown', 'grey', 'darkred', 'orange', 'pink', 'darkblue', diff --git a/python_tools/void_python_tools/plotting/plotTools.py b/python_tools/void_python_tools/plotting/plotTools.py index a7c5058..1d42578 100644 --- a/python_tools/void_python_tools/plotting/plotTools.py +++ b/python_tools/void_python_tools/plotting/plotTools.py @@ -1,11 +1,12 @@ __all__=['plotRedshiftDistribution', 'plotSizeDistribution', 'plot1dProfiles', - 'plotMarg1d'] + 'plotMarg1d', 'plotNumberDistribution'] from void_python_tools.backend.classes import * from plotDefs import * import numpy as np import os import pylab as plt +import void_python_tools.apTools as vp # ----------------------------------------------------------------------------- def plotRedshiftDistribution(workDir=None, sampleList=None, figDir=None, @@ -43,7 +44,7 @@ def plotRedshiftDistribution(workDir=None, sampleList=None, figDir=None, zMax = sample.zRange[1] range = (zMin, zMax) - nbins = np.ceil((zMax-zMin)/0.1) + nbins = np.ceil((zMax-zMin)/0.02) thisMax = np.max(data[:,5]) thisMin = np.min(data[:,5]) @@ -215,3 +216,63 @@ def plotMarg1d(workDir=None, sampleList=None, figDir=None, +# ----------------------------------------------------------------------------- +def plotNumberDistribution(workDir=None, sampleList=None, figDir=None, + plotNameBase="numberdist", + showPlot=False, dataPortion=None, setName=None): + + plt.clf() + plt.xlabel("Void Radius (Mpc/h)") + plt.ylabel(r"N > R [h^3 Mpc^{-3}]") + + plotTitle = setName + + plotName = plotNameBase + + plt.yscale('log') + + for (iSample,sample) in enumerate(sampleList): + + sampleName = sample.fullName + lineTitle = sampleName + + if sample.dataType == "observation": + boxVol = vp.getSurveyProps(sample.maskFile, stack.zMin, stack.zMax, + sample.zRange[0], sample.zRange[1], "all", + selectionFuncFile=sample.selFunFile)[0] + else: + boxVol = sample.boxLen*sample.boxLen*(zBoundaryMpc[1]-zBoundaryMpc[0]) + + filename = workDir+"/sample_"+sampleName+"/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[:,4] + 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") + plt.title(plotTitle) + + plt.xlim(xMin, xMax) + #plt.xlim(xMin, xMax*1.4) # make room for legend + + 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") + + + From 08ea4ef0a2327e1cfe1e015501ab69cba55244c1 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Tue, 13 Nov 2012 13:56:44 -0600 Subject: [PATCH 18/32] more informative sample_info file --- c_tools/mock/generateFromCatalog.cpp | 19 +++++++++++++------ c_tools/mock/generateMock.cpp | 13 +++++++++++++ .../void_python_tools/backend/launchers.py | 1 + 3 files changed, 27 insertions(+), 6 deletions(-) diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index febc6d3..f0816b0 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/c_tools/mock/generateFromCatalog.cpp @@ -451,6 +451,19 @@ if (i < 10) printf("TEST WRITE %d %e\n", (pdata.pos[i].xyz[j]+Lmax)/(2*Lmax)); //v2->put(expansion_fac, pdata.pos.size()); //delete[] expansion_fac; + + FILE *infoFile = fopen("sample_info.txt", "w"); + fprintf(infoFile, "x_min = %f\n", -Lmax/100.); + fprintf(infoFile, "x_max = %f\n", Lmax/100.); + fprintf(infoFile, "y_min = %f\n", -Lmax/100.); + fprintf(infoFile, "y_max = %f\n", Lmax/100.); + fprintf(infoFile, "z_min = %f\n", -Lmax/100.); + fprintf(infoFile, "z_max = %f\n", Lmax/100.); + fprintf(infoFile, "mask_index = %d\n", pdata.mask_index); + fprintf(infoFile, "total_particles = %d\n", pdata.pos.size()); + fclose(infoFile); + + } int main(int argc, char **argv) @@ -517,12 +530,6 @@ int main(int argc, char **argv) fprintf(fp, "%d", output_data.mask_index); fclose(fp); - fp = fopen("sample_info.txt", "w"); - fprintf(fp, "Lmax = %f\n", output_data.Lmax); - fprintf(fp, "mask_index = %d\n", output_data.mask_index); - fprintf(fp, "total_particles = %d\n", output_data.pos.size()); - fclose(fp); - fp = fopen("total_particles.txt", "w"); fprintf(fp, "%d", output_data.pos.size()); fclose(fp); diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index bfdf82f..dcf1acb 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -360,6 +360,7 @@ void generateOutput(SimuData *data, int axis, f.writeReal32(data->Pos[x2][i]); } f.endCheckpoint(); + } void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& args_info) @@ -466,6 +467,18 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& delete[] particle_id; delete[] expansion_fac; + + + FILE *fp = fopen("sample_info.txt", "w"); + fprintf(fp, "x_min = %f\n", ranges[0][0]); + fprintf(fp, "x_max = %f\n", ranges[0][1]); + fprintf(fp, "y_min = %f\n", ranges[1][0]); + fprintf(fp, "y_max = %f\n", ranges[1][1]); + fprintf(fp, "z_min = %f\n", ranges[2][0]); + fprintf(fp, "z_max = %f\n", ranges[2][1]); + fprintf(fp, "mask_index = -1\n"); + fprintf(fp, "total_particles = %d\n", boxed->NumPart); + fclose(fp); } void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, generateMock_info& args_info) diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index b8cb57f..2686c5b 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -160,6 +160,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, if os.access("comoving_distance.txt", os.F_OK): os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) + os.system("mv %s %s" % ("sample_info.txt", zobovDir)) if os.access(parmFile, os.F_OK): os.unlink(parmFile) From afe22f71fa0088863707fc96e4a3aba90ef29efe Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Tue, 13 Nov 2012 15:23:37 -0600 Subject: [PATCH 19/32] added start of extra plotting scripts for comparison across data sets --- plotting/plotCompareDist.py | 90 +++++++++++++++++++ .../void_python_tools/plotting/plotTools.py | 3 +- 2 files changed, 92 insertions(+), 1 deletion(-) create mode 100755 plotting/plotCompareDist.py diff --git a/plotting/plotCompareDist.py b/plotting/plotCompareDist.py new file mode 100755 index 0000000..30038d2 --- /dev/null +++ b/plotting/plotCompareDist.py @@ -0,0 +1,90 @@ +#!/usr/bin/env python + +# plots cumulative distributions of number counts + +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 pylab as plt +import argparse + +# ------------------------------------------------------------------------------ + +workDir = "/home/psutter2/workspace/Voids/" +figDir = "./figs" + +sampleDirList = [ "multidark/md_halos/sample_md_halos_z0.03_d00/", + "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.03_d00/" ] + +plotNameBase = "compdist" + +dataPortion = "central" + +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)') +args = parser.parse_args() + +# ------------------------------------------------------------------------------ + +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 Radius (Mpc/h)") +plt.ylabel(r"N > R [h^3 Mpc^{-3}]") +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, stack.zMin, stack.zMax, + sample.zRange[0], sample.zRange[1], "all", + selectionFuncFile=sample.selFunFile)[0] + else: + boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1]-sample.zBoundaryMpc[0]) + + 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[:,4] + 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") +#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/plotting/plotTools.py b/python_tools/void_python_tools/plotting/plotTools.py index 1d42578..89a545b 100644 --- a/python_tools/void_python_tools/plotting/plotTools.py +++ b/python_tools/void_python_tools/plotting/plotTools.py @@ -241,7 +241,8 @@ def plotNumberDistribution(workDir=None, sampleList=None, figDir=None, sample.zRange[0], sample.zRange[1], "all", selectionFuncFile=sample.selFunFile)[0] else: - boxVol = sample.boxLen*sample.boxLen*(zBoundaryMpc[1]-zBoundaryMpc[0]) + boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] - + sample.zBoundaryMpc[0]) filename = workDir+"/sample_"+sampleName+"/centers_"+dataPortion+"_"+\ sampleName+".out" From 390fb2f4e2d9b2bbea8834b91f945e0015c84f72 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Wed, 14 Nov 2012 04:37:30 -0600 Subject: [PATCH 20/32] bug fixes to plotting. prepareCatalogs now handles gadget files; prepareGadgetCatalogs now deprecated --- pipeline/generateCatalog.py | 5 + pipeline/prepareCatalogs.py | 182 +++++++++++------- plotting/plotCompareDist.py | 17 +- .../void_python_tools/plotting/plotTools.py | 10 +- 4 files changed, 133 insertions(+), 81 deletions(-) diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index cf4dc23..55a05ce 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -5,6 +5,7 @@ from void_python_tools.backend import * from void_python_tools.plotting import * import imp +import pickle # ------------------------------------------------------------------------------ @@ -60,6 +61,10 @@ for sample in dataSampleList: if not os.access(zobovDir, os.F_OK): os.makedirs(zobovDir) + # save this sample's information + with open(zobovDir+"/sample_info.dat", 'wb') as output: + pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL) + # --------------------------------------------------------------------------- if (startCatalogStage <= 1) and (endCatalogStage >= 1) and not sample.isCombo: print " Extracting tracers from catalog...", diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index d331d85..451d2d4 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -10,33 +10,67 @@ import sys import void_python_tools as vp import argparse +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- +# CONFIGURATION + +# directory for the input simulation/observational particle files catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/multidark/" -scriptDir = os.getenv("PWD")+"/multidark/" + +# path to HOD code hodPath = os.getenv("HOME")+"/projects/Voids/hod/HOD.x" -outputDir = os.getenv("HOME")+"/workspace/Voids/multidark/" +# where to put the final void catalog, figures, and output logs +voidOutputDir = os.getenv("HOME")+"/workspace/Voids/multidark/" figDir = os.getenv("PWD")+"/../figs/multidark/" logDir = os.getenv("PWD")+"/../logs/multidark/" -dataFormat = "multidark" -dataUnit = 1 +# where to place the pipeline scripts +scriptDir = os.getenv("PWD")+"/multidark/" + +# simulation or observation? dataType = "simulation" + +# available formats for simulation: gadget, multidark +dataFormat = "multidark" +dataUnit = 1 # as multiple of Mpc/h + +# place particles on the lightcone? useLightCone = True -redshiftFileBase = "mdr1_particles_z" +# common filename of particle files +particleFileBase = "mdr1_particles_z" + +# list of file numbers for the particle files +# to get particle file name, we take particleFileBase+fileNum +fileNums = (("0.0", "0.53", "1.0")) + +# redshift of each file in the above list redshifts = (("0.0", "0.53", "1.0")) + +# how many independent slices along the z-axis? numSlices = 4 + +# how many subdivisions along the x- and y- axis? +# ( = 2 will make 4 subvolumes for each slice, = 3 will make 9, etc.) numSubvolumes = 1 +# prefix to give all outputs prefix = "md_" + +# list of desired subsamples subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002)) -numPart = 1024*1024*1024 # number of particles in base catalog -lbox = 1000 +# simulation information +numPart = 1024*1024*1024 +lbox = 1000 # Mpc/h omegaM = 0.27 hubble = 0.70 +# END CONFIGURATION # ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- + LIGHT_SPEED = 299792.458 parser = argparse.ArgumentParser(description='options') @@ -58,35 +92,31 @@ parser.add_argument('--all', dest='all', action='store_const', args = parser.parse_args() #------------------------------------------------------------------------------ -def getSampleName(prefix, base, redshift, useVel, iSlice=-1, iVol=-1): - useVelString = "" - if useVel: useVelString = "_pv" +def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): - sampleName = prefix + base + useVelString + "_z" + redshift + sampleName = setName - if iSlice != -1: - sampleName += "_s" + str(iSlice) + #if useVel: sampleName += "_pv" - if iVol != -1: - sampleName += "_d" + iVol + sampleName += "_z" + redshift + + #if iSlice != -1: sampleName += "_s" + str(iSlice) + + if iVol != -1: sampleName += "_d" + iVol return sampleName #------------------------------------------------------------------------------ # for given dataset parameters, outputs a script for use with analyzeVoids -def writeScript(prefix, base, scriptDir, catalogDir, redshifts, numSubvolumes, +def writeScript(setName, dataFileNameBase, + scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, useVel, lbox, minRadius, omegaM, subsample=1.0, suffix=".dat"): - setName = prefix + base - - dataFileNameBase = setName #+ ".dat" if useVel: setName += "_pv" - scriptFileName = scriptDir + "/" + "md_" + base - if useVel: scriptFileName += "_pv" - scriptFileName += ".py" + scriptFileName = scriptDir + "/" + setName + ".py" scriptFile = open(scriptFileName, 'w') scriptFile.write("""#!/usr/bin/env/python @@ -113,17 +143,17 @@ dataPortions = ["central"] dataSampleList = [] """) - dataInfo = """ -catalogName = "{sampleName}" -workDir = "{outputDir}/{sampleName}/" -inputDataDir = "{catalogDir}" -figDir = "{figDir}/{sampleName}/" -logDir = "{logDir}/{sampleName}/" + dataInfo = """ +setName = "{setName}" + +workDir = "{voidOutputDir}/{setName}/" +inputDataDir = "{inputDataDir}" +figDir = "{figDir}/{setName}/" +logDir = "{logDir}/{setName}/" """ - scriptFile.write(dataInfo.format(sampleName=setName, figDir=figDir, - logDir=logDir, outputDir=outputDir, - catalogDir=catalogDir, + scriptFile.write(dataInfo.format(setName=setName, figDir=figDir, + logDir=logDir, voidOutputDir=voidOutputDir, inputDataDir=catalogDir)) sampleInfo = """ @@ -156,15 +186,15 @@ 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) """ - - for redshift in redshifts: + 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 # converter from redshift to comoving distance - zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000) + 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) @@ -178,31 +208,28 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False) dzSafe = 0.0 for iSlice in xrange(numSlices): - # trim off the first and last 10,000 km/s 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 + "_z" + redshift + suffix + dataFileName = dataFileNameBase + fileNum + suffix for iX in xrange(numSubvolumes): for iY in xrange(numSubvolumes): - + mySubvolume = "%d%d" % (iX, iY) - sampleName = getSampleName(prefix, base, sliceMin, useVel, - iSlice=-1, iVol=mySubvolume) - #sampleName = getSampleName(prefix, base, redshift, useVel, - # iSlice=iSlice, iVol=mySubvolume) + 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, @@ -240,28 +267,37 @@ for thisSubSample in subSamples: if args.script or args.all: print " Doing subsample", thisSubSample, " scripts" + setName = prefix+"ss"+str(thisSubSample) if dataFormat == "multidark": - writeScript(prefix, - "ss"+str(thisSubSample), scriptDir, catalogDir, redshifts, + writeScript(setName, "md.ss"+str(thisSubSample)+"_z", + scriptDir, catalogDir, fileNums, + redshifts, numSubvolumes, numSlices, False, lbox, minRadius, omegaM, subsample=1.0) - writeScript(prefix, "ss"+str(thisSubSample), scriptDir, catalogDir, + writeScript(setName, "md.ss"+str(thisSubSample)+"_z", + scriptDir, catalogDir, + fileNums, redshifts, numSubvolumes, numSlices, True, lbox, minRadius, omegaM, subsample=1.0) elif dataFormat == "gadget": - writeScript("", redshiftFileBase, scriptDir, catalogDir, redshifts, + writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums, + redshifts, + numSubvolumes, numSlices, False, lbox, minRadius, omegaM, + subsample=thisSubSample, suffix="") + writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums, + redshifts, numSubvolumes, numSlices, True, lbox, minRadius, omegaM, subsample=thisSubSample, suffix="") if args.subsample or args.all: print " Doing subsample", thisSubSample - for redshift in redshifts: + for (iRedshift, redshift) in enumerate(redshifts): print " redshift", redshift - dataFile = catalogDir+"/"+redshiftFileBase+redshift + dataFile = catalogDir+"/"+particleFileBase+fileNums[iRedshift] inFile = open(dataFile, 'r') - sampleName = getSampleName("ss"+str(thisSubSample), redshift, False) + sampleName = "md.ss"+str(thisSubSample)+"_z"+redshift outFile = open(catalogDir+"/"+sampleName+".dat", 'w') outFile.write("%f\n" %(lbox)) @@ -293,7 +329,7 @@ for thisSubSample in subSamples: if args.script or args.all: print " Doing halo scripts" - dataFile = catalogDir+"/mdr1_halos_z"+redshifts[0] + dataFile = catalogDir+"/mdr1_halos_z"+fileNums[0] inFile = open(dataFile, 'r') numPart = 0 for line in inFile: numPart += 1 @@ -302,25 +338,27 @@ if args.script or args.all: minRadius = 2*int(np.ceil(lbox/numPart**(1./3.))) if dataFormat == "multidark": - writeScript(prefix, "halos", scriptDir, catalogDir, redshifts, + setName = prefix+"halos" + writeScript(setName, "md.halos_z", scriptDir, catalogDir, fileNums, + redshifts, numSubvolumes, numSlices, False, lbox, minRadius, omegaM) - writeScript(prefix, "halos", scriptDir, catalogDir, redshifts, numSubvolumes, + writeScript(setName, "md.halos_z", scriptDir, catalogDir, fileNums, + redshifts, numSubvolumes, numSlices, True, lbox, minRadius, omegaM) if args.halos or args.all: print " Doing halos" - for redshift in redshifts: + for (iRedshift, redshift) in enumerate(redshifts): print " z = ", redshift - dataFile = catalogDir+"/mdr1_halos_z"+redshift + dataFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift] inFile = open(dataFile, 'r') numPart = 0 for line in inFile: numPart += 1 inFile.close() - sampleName = getSampleName("halos", redshift, False) - + sampleName = "md.halos_z"+redshift outFile = open(catalogDir+"/"+sampleName+".dat", 'w') outFile.write("%f\n" %(lbox)) @@ -390,21 +428,22 @@ root_filename hod if args.script or args.all: print " Doing DR7 HOD scripts" if dataFormat == "multidark": - writeScript(prefix, - "hod_dr72dim2", scriptDir, catalogDir, redshifts, + setName = prefix+"hod_dr72dim2" + writeScript(setName, "md.hod_dr72dim2_z", + scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, False, lbox, 5, omegaM) - writeScript(prefix, - "hod_dr72dim2", scriptDir, catalogDir, redshifts, + 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 redshift in redshifts: + for (iRedshift, redshift) in enumerate(redshifts): print " z = ", redshift parFileName = "./hod.par" parFile = open(parFileName, 'w') - haloFile = catalogDir+"/mdr1_halos_z"+redshift + haloFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift] parFile.write(parFileText.format(omegaM=omegaM, hubble=hubble, redshift=redshift, @@ -420,7 +459,7 @@ if args.hod or args.all: os.system(hodPath+" "+parFileName+">& /dev/null") - sampleName = getSampleName("hod_dr72dim2", redshift, False) + sampleName = getSampleName("md.hod_dr72dim2", redshift, False) outFileName = catalogDir+"/"+sampleName+".dat" os.system("mv hod.mock" + " " + outFileName) @@ -431,21 +470,22 @@ if args.hod or args.all: if args.script or args.all: print " Doing DR9 HOD scripts" if dataFormat == "multidark": - writeScript(prefix, - "hod_dr9mid", scriptDir, catalogDir, redshifts, + setName = prefix+"hod_dr9mid" + writeScript(setName, "md.hod_dr9mid_z", + scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, False, lbox, 5, omegaM) - writeScript(prefix, - "hod_dr9mid", scriptDir, catalogDir, redshifts, + writeScript(prefix, "md.hod_dr9mid_z", + scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, True, lbox, 5, omegaM) if args.hod or args.all: print " Doing DR9 HOD" - for redshift in redshifts: + for (iRedshift, redshift) in enumerate(redshifts): print " z = ", redshift parFileName = "./hod.par" parFile = open(parFileName, 'w') - haloFile = catalogDir+"/mdr1_halos_z"+redshift + haloFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift] parFile.write(parFileText.format(omegaM=omegaM, hubble=hubble, redshift=redshift, @@ -461,7 +501,7 @@ if args.hod or args.all: os.system(hodPath+" "+parFileName+">& /dev/null") - sampleName = getSampleName("hod_dr9mid", redshift, False) + sampleName = getSampleName("md.hod_dr9mid", redshift, False) outFileName = catalogDir+"/"+sampleName+".dat" os.system("mv hod.mock" + " " + outFileName) diff --git a/plotting/plotCompareDist.py b/plotting/plotCompareDist.py index 30038d2..b2c47bd 100755 --- a/plotting/plotCompareDist.py +++ b/plotting/plotCompareDist.py @@ -8,7 +8,8 @@ import void_python_tools.apTools as vp import imp import pickle import os -import pylab as plt +import matplotlib.pyplot as plt +import numpy as np import argparse # ------------------------------------------------------------------------------ @@ -17,7 +18,9 @@ workDir = "/home/psutter2/workspace/Voids/" figDir = "./figs" sampleDirList = [ "multidark/md_halos/sample_md_halos_z0.03_d00/", - "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.03_d00/" ] + "multidark/md_ss0.1/sample_md_ss0.1_z0.03_d00/", + "multidark/md_hod_dr9mid/sample_md_hod_dr9mid_z0.03_d00/", + "sdss_dr9/sample_lss.dr9cmassmid.dat/" ] plotNameBase = "compdist" @@ -42,7 +45,7 @@ for sampleDir in sampleDirList: plt.clf() plt.xlabel("Void Radius (Mpc/h)") -plt.ylabel(r"N > R [h^3 Mpc^{-3}]") +plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]") plt.yscale('log') plotName = plotNameBase @@ -53,11 +56,15 @@ for (iSample,sample) in enumerate(dataSampleList): lineTitle = sampleName if sample.dataType == "observation": - boxVol = vp.getSurveyProps(sample.maskFile, stack.zMin, stack.zMax, + boxVol = vp.getSurveyProps(sample.maskFile, + sample.zBoundary[0], sample.zBoundary[1], sample.zRange[0], sample.zRange[1], "all", selectionFuncFile=sample.selFunFile)[0] else: - boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1]-sample.zBoundaryMpc[0]) + boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] - + sample.zBoundaryMpc[0]) + + boxVol *= 1.e-9 # Mpc->Gpc filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\ sampleName+".out" diff --git a/python_tools/void_python_tools/plotting/plotTools.py b/python_tools/void_python_tools/plotting/plotTools.py index 89a545b..c36ecad 100644 --- a/python_tools/void_python_tools/plotting/plotTools.py +++ b/python_tools/void_python_tools/plotting/plotTools.py @@ -223,7 +223,7 @@ def plotNumberDistribution(workDir=None, sampleList=None, figDir=None, plt.clf() plt.xlabel("Void Radius (Mpc/h)") - plt.ylabel(r"N > R [h^3 Mpc^{-3}]") + plt.ylabel(r"N > R [$h^3$ Mpc$^{-3}$]") plotTitle = setName @@ -237,13 +237,16 @@ def plotNumberDistribution(workDir=None, sampleList=None, figDir=None, lineTitle = sampleName if sample.dataType == "observation": - boxVol = vp.getSurveyProps(sample.maskFile, stack.zMin, stack.zMax, + boxVol = vp.getSurveyProps(sample.maskFile, + sample.zBoundary[0], sample.zBoundary[1], sample.zRange[0], sample.zRange[1], "all", selectionFuncFile=sample.selFunFile)[0] else: boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] - sample.zBoundaryMpc[0]) + boxVol *= 1.e-9 + filename = workDir+"/sample_"+sampleName+"/centers_"+dataPortion+"_"+\ sampleName+".out" if not os.access(filename, os.F_OK): @@ -265,9 +268,6 @@ def plotNumberDistribution(workDir=None, sampleList=None, figDir=None, plt.legend(title = "Samples", loc = "upper right") plt.title(plotTitle) - plt.xlim(xMin, xMax) - #plt.xlim(xMin, xMax*1.4) # make room for legend - 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") From 292afacaf53e23977e63c167dee0f92847e528ee Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Wed, 14 Nov 2012 06:31:53 -0600 Subject: [PATCH 21/32] now includes automatic detection and handling of periodic boxes; all logic for dealing with this moved to pruneVoids --- c_tools/stacking/pruneVoids.cpp | 124 +++++++++++------- c_tools/stacking/pruneVoids.ggo | 2 + .../void_python_tools/backend/launchers.py | 8 ++ 3 files changed, 87 insertions(+), 47 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index d6d5890..1390f3c 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -39,7 +39,7 @@ typedef struct voidZoneStruct { typedef struct voidStruct { float vol, coreDens, zoneVol, densCon, voidProb, radius; int voidID, numPart, numZones, coreParticle, zoneNumPart; - float maxRadius, nearestMock, centralDen, redshift; + float maxRadius, nearestMock, centralDen, redshift, redshiftInMpc; float center[3], barycenter[3]; int accepted; } VOID; @@ -89,16 +89,34 @@ int main(int argc, char **argv) { double ranges[2][3], boxLen[3], mul; double volNorm, radius; int clock1, clock2; + int periodicX=0, periodicY=0, periodicZ=0; numVoids = args_info.numVoids_arg; mockIndex = args_info.mockIndex_arg; tolerance = args_info.tolerance_arg; clock1 = clock(); - printf("Pruning parameters: %f %f %f\n", args_info.zMin_arg, - args_info.zMax_arg, - args_info.rMin_arg); - + printf("Pruning parameters: %f %f %f %s\n", args_info.zMin_arg, + args_info.zMax_arg, + args_info.rMin_arg, + args_info.periodic_arg); + + // check for periodic box + if (!args_info.isObservation_flag) { + if ( strchr(args_info.periodic_arg, 'x') != NULL) { + periodicX = 1; + printf("Will assume x-direction is periodic.\n"); + } + if ( strchr(args_info.periodic_arg, 'y') != NULL) { + periodicY = 1; + printf("Will assume y-direction is periodic.\n"); + } + if ( strchr(args_info.periodic_arg, 'z') != NULL) { + periodicZ = 1; + printf("Will assume z-direction is periodic.\n"); + } + } + // load box size printf("\n Getting info...\n"); NcFile f_info(args_info.extraInfo_arg); @@ -291,17 +309,14 @@ int main(int argc, char **argv) { voids[iVoid].barycenter[1] = 0.; voids[iVoid].barycenter[2] = 0.; -// TODO handle periodic boundaries? for (p = 0; p < voids[iVoid].numPart; p++) { dist[0] = voidPart[p].x - voids[iVoid].center[0]; dist[1] = voidPart[p].y - voids[iVoid].center[1]; dist[2] = voidPart[p].z - voids[iVoid].center[2]; - //if (!args_info.isObservation_flag) { - // dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0])); - // dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1])); - // dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2])); - //} + if (periodicX) dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0])); + if (periodicY) dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1])); + if (periodicZ) dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2])); voids[iVoid].barycenter[0] += voidPart[p].vol*(dist[0]); voids[iVoid].barycenter[1] += voidPart[p].vol*(dist[1]); @@ -324,11 +339,9 @@ int main(int argc, char **argv) { dist[1] = voidPart[p].y - voids[iVoid].barycenter[1]; dist[2] = voidPart[p].z - voids[iVoid].barycenter[2]; - //if (!args_info.isObservation_flag) { - // dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0])); - // dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1])); - // dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2])); - //} + if (periodicX) dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0])); + if (periodicY) dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1])); + if (periodicZ) dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2])); dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); if (dist2 < centralRad) centralDen += 1; @@ -336,21 +349,21 @@ int main(int argc, char **argv) { voids[iVoid].centralDen = centralDen / (4./3. * M_PI * pow(centralRad, 3./2.)); // compute maximum extent - if (args_info.isObservation_flag) { - maxDist = 0.; - for (p = 0; p < voids[iVoid].numPart; p++) { - for (p2 = p; p2 < voids[iVoid].numPart; p2++) { + //if (args_info.isObservation_flag) { + // maxDist = 0.; + // for (p = 0; p < voids[iVoid].numPart; p++) { + // for (p2 = p; p2 < voids[iVoid].numPart; p2++) { - dist[0] = voidPart[p].x - voidPart[p2].x; - dist[1] = voidPart[p].y - voidPart[p2].y; - dist[2] = voidPart[p].z - voidPart[p2].z; +// dist[0] = voidPart[p].x - voidPart[p2].x; +// dist[1] = voidPart[p].y - voidPart[p2].y; +// dist[2] = voidPart[p].z - voidPart[p2].z; - dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); - if (dist2 > maxDist) maxDist = dist2; - } - } - voids[iVoid].maxRadius = sqrt(maxDist)/2.; - } else { +// dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); +// if (dist2 > maxDist) maxDist = dist2; +// } +// } +// voids[iVoid].maxRadius = sqrt(maxDist)/2.; +// } else { maxDist = 0.; for (p = 0; p < voids[iVoid].numPart; p++) { @@ -358,11 +371,15 @@ int main(int argc, char **argv) { dist[0] = voidPart[p].y - voids[iVoid].barycenter[1]; dist[0] = voidPart[p].z - voids[iVoid].barycenter[2]; + if (periodicX) dist[0] = fmin(dist[0], abs(boxLen[0]-dist[0])); + if (periodicY) dist[1] = fmin(dist[1], abs(boxLen[1]-dist[1])); + if (periodicZ) dist[2] = fmin(dist[2], abs(boxLen[2]-dist[2])); + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); if (dist2 > maxDist) maxDist = dist2; } voids[iVoid].maxRadius = sqrt(maxDist); - } +// } if (args_info.isObservation_flag) { // compute distance from center to nearest mock @@ -382,28 +399,41 @@ int main(int argc, char **argv) { } if (args_info.isObservation_flag) { - voids[iVoid].redshift = + voids[iVoid].redshiftInMpc = sqrt(pow(voids[iVoid].barycenter[0] - boxLen[0]/2.,2) + pow(voids[iVoid].barycenter[1] - boxLen[1]/2.,2) + pow(voids[iVoid].barycenter[2] - boxLen[2]/2.,2)); - voids[iVoid].redshift = voids[iVoid].redshift; - redshift = voids[iVoid].redshift; - nearestEdge = fmin(fabs(redshift-args_info.zMin_arg*LIGHT_SPEED), - fabs(redshift-args_info.zMax_arg*LIGHT_SPEED)); + voids[iVoid].redshiftInMpc = voids[iVoid].redshiftInMpc; + redshift = voids[iVoid].redshiftInMpc; + nearestEdge = fmin(fabs(redshift-args_info.zMin_arg*LIGHT_SPEED/100.), + fabs(redshift-args_info.zMax_arg*LIGHT_SPEED/100.)); + voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.; + } else { + + voids[iVoid].redshiftInMpc = voids[iVoid].barycenter[2]; voids[iVoid].redshift = voids[iVoid].barycenter[2]/LIGHT_SPEED*100.; + + nearestEdge = 1.e99; - nearestEdge = fmin( - fabs(voids[iVoid].barycenter[0] - ranges[0][0]), - fabs(voids[iVoid].barycenter[0] - ranges[0][1])); - nearestEdge = fmin(nearestEdge, - fabs(voids[iVoid].barycenter[1] - ranges[1][0])); - nearestEdge = fmin(nearestEdge, - fabs(voids[iVoid].barycenter[1] - ranges[1][1])); - nearestEdge = fmin(nearestEdge, - fabs(voids[iVoid].barycenter[2] - ranges[2][0])); - nearestEdge = fmin(nearestEdge, - fabs(voids[iVoid].barycenter[2] - ranges[2][1])); + if (!periodicX) { + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].barycenter[0] - ranges[0][0])); + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].barycenter[0] - ranges[0][1])); + } + if (!periodicY) { + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].barycenter[1] - ranges[1][0])); + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].barycenter[1] - ranges[1][1])); + } + if (!periodicZ) { + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].barycenter[2] - ranges[2][0])); + nearestEdge = fmin(nearestEdge, + fabs(voids[iVoid].barycenter[2] - ranges[2][1])); + } } if (nearestEdge < voids[iVoid].nearestMock) { @@ -505,7 +535,7 @@ int main(int argc, char **argv) { fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n", atan((voids[iVoid].barycenter[1]-boxLen[1]/2.)/(voids[iVoid].barycenter[0]-boxLen[0]/2.)) * 180/M_PI + 180, - asin((voids[iVoid].barycenter[2]-boxLen[2]/2.)/voids[iVoid].redshift) * 180/M_PI, + asin((voids[iVoid].barycenter[2]-boxLen[2]/2.)/voids[iVoid].redshiftInMpc) * 180/M_PI, voids[iVoid].redshift, voids[iVoid].radius, voids[iVoid].voidID); diff --git a/c_tools/stacking/pruneVoids.ggo b/c_tools/stacking/pruneVoids.ggo index bd97e50..36e5143 100644 --- a/c_tools/stacking/pruneVoids.ggo +++ b/c_tools/stacking/pruneVoids.ggo @@ -40,6 +40,8 @@ option "outSkyPositions" - "output sky positions of voids" string required option "dataPortion" - "all, central, or edge" string required +option "periodic" - "Set of edges which are periodic" string optional default="xy" + option "tolerance" - "Fraction of void width to consider edge" double optional default="1.0" option "centralRadFrac" - "Fraction of void radii to consider central region" double optional default="4" diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 2686c5b..7237408 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -244,11 +244,18 @@ def launchPrune(sample, binPath, thisDataPortion=None, totalPart = open(zobovDir+"/total_particles.txt", "r").read() maxDen = 0.2*float(mockIndex)/float(totalPart) observationLine = " --isObservation" + periodicLine = "--periodic=''" else: mockIndex = -1 maxDen = 0.2 observationLine = "" + periodicLine = " --periodic='" + if sample.numSubvolumes == 1: periodicLine += "xy" + if sample.zBoundaryMpc[0] == 0 and \ + sample.zBoundaryMpc[1] == sample.boxLen : periodicLine += "z" + periodicLine += "' " + if not (continueRun and jobSuccessful(logFile, "NetCDF: Not a valid ID\n")): cmd = binPath cmd += " --partFile=" + zobovDir+"/zobov_slice_"+str(sampleName) @@ -267,6 +274,7 @@ def launchPrune(sample, binPath, thisDataPortion=None, cmd += " --rMin=" + str(sample.minVoidRadius) cmd += " --numVoids=" + str(numVoids) cmd += observationLine + cmd += periodicLine cmd += " --output=" + zobovDir+"/voidDesc_"+\ str(thisDataPortion)+"_"+\ str(sampleName)+".out" From 443b4bdeecda199f3a5b87929871ee1b41155ddb Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 15 Nov 2012 06:37:59 -0600 Subject: [PATCH 22/32] fixed bug in gadget generation script; prepareGadget now uses separate parameter files; support for random catalogs now enabled --- pipeline/generateCatalog.py | 6 +- pipeline/prepareCatalogs.py | 173 ++++++------ pipeline/prepareGadgetCatalog.py | 246 ------------------ .../void_python_tools/backend/launchers.py | 8 +- 4 files changed, 91 insertions(+), 342 deletions(-) delete mode 100755 pipeline/prepareGadgetCatalog.py diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 55a05ce..cd8408e 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -113,10 +113,10 @@ if (startCatalogStage <= 4) and (endCatalogStage >= 4): for thisDataPortion in dataPortions: plotRedshiftDistribution(workDir, dataSampleList, figDir, showPlot=False, - dataPortion=thisDataPortion, setName=catalogName) + dataPortion=thisDataPortion, setName=setName) plotSizeDistribution(workDir, dataSampleList, figDir, showPlot=False, - dataPortion=thisDataPortion, setName=catalogName) + dataPortion=thisDataPortion, setName=setName) plotNumberDistribution(workDir, dataSampleList, figDir, showPlot=False, - dataPortion=thisDataPortion, setName=catalogName) + dataPortion=thisDataPortion, setName=setName) print "\n Done!" diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index 451d2d4..d23d36e 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -9,66 +9,17 @@ import os import sys import void_python_tools as vp import argparse +import imp -# ----------------------------------------------------------------------------- -# ----------------------------------------------------------------------------- -# CONFIGURATION +# ------------------------------------------------------------------------------ -# directory for the input simulation/observational particle files -catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/multidark/" +def my_import(name): + mod = __import__(name) + components = name.split('.') + for comp in components[1:]: + mod = getattr(mod, comp) + return mod -# path to HOD code -hodPath = os.getenv("HOME")+"/projects/Voids/hod/HOD.x" - -# where to put the final void catalog, figures, and output logs -voidOutputDir = os.getenv("HOME")+"/workspace/Voids/multidark/" -figDir = os.getenv("PWD")+"/../figs/multidark/" -logDir = os.getenv("PWD")+"/../logs/multidark/" - -# where to place the pipeline scripts -scriptDir = os.getenv("PWD")+"/multidark/" - -# simulation or observation? -dataType = "simulation" - -# available formats for simulation: gadget, multidark -dataFormat = "multidark" -dataUnit = 1 # as multiple of Mpc/h - -# place particles on the lightcone? -useLightCone = True - -# common filename of particle files -particleFileBase = "mdr1_particles_z" - -# list of file numbers for the particle files -# to get particle file name, we take particleFileBase+fileNum -fileNums = (("0.0", "0.53", "1.0")) - -# redshift of each file in the above list -redshifts = (("0.0", "0.53", "1.0")) - -# how many independent slices along the z-axis? -numSlices = 4 - -# how many subdivisions along the x- and y- axis? -# ( = 2 will make 4 subvolumes for each slice, = 3 will make 9, etc.) -numSubvolumes = 1 - -# prefix to give all outputs -prefix = "md_" - -# list of desired subsamples -subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002)) - -# simulation information -numPart = 1024*1024*1024 -lbox = 1000 # Mpc/h -omegaM = 0.27 -hubble = 0.70 - -# END CONFIGURATION -# ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- LIGHT_SPEED = 299792.458 @@ -85,12 +36,25 @@ parser.add_argument('--halos', dest='halos', action='store_const', help='write halos') parser.add_argument('--hod', dest='hod', action='store_const', const=True, default=False, - help='write hos') + help='write hod') parser.add_argument('--all', dest='all', action='store_const', const=True, default=False, help='write everything') +parser.add_argument('--parmFile', dest='parmFile', + default="", + 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)) + + #------------------------------------------------------------------------------ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): @@ -253,6 +217,7 @@ newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False) #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ 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 @@ -288,45 +253,73 @@ for thisSubSample in subSamples: redshifts, numSubvolumes, numSlices, True, lbox, minRadius, omegaM, subsample=thisSubSample, suffix="") + elif dataFormat == "random": + writeScript(setName, "ran.ss"+str(thisSubSample)+"_z", + scriptDir, catalogDir, fileNums, + redshifts, + numSubvolumes, numSlices, False, lbox, minRadius, omegaM, + subsample=1.0) if args.subsample or args.all: print " Doing subsample", thisSubSample for (iRedshift, redshift) in enumerate(redshifts): print " redshift", redshift - dataFile = catalogDir+"/"+particleFileBase+fileNums[iRedshift] - inFile = open(dataFile, 'r') + + if dataFormat == "multidark": + dataFile = catalogDir+"/"+particleFileBase+fileNums[iRedshift] + inFile = open(dataFile, 'r') - sampleName = "md.ss"+str(thisSubSample)+"_z"+redshift - outFile = open(catalogDir+"/"+sampleName+".dat", 'w') + sampleName = "md.ss"+str(thisSubSample)+"_z"+redshift + outFile = open(catalogDir+"/"+sampleName+".dat", 'w') - outFile.write("%f\n" %(lbox)) - outFile.write("%s" %(omegaM)) - outFile.write("%s" %(hubble)) - outFile.write("%s\n" %(redshift)) - outFile.write("%d\n" %(maxKeep)) + outFile.write("%f\n" %(lbox)) + outFile.write("%s\n" %(omegaM)) + outFile.write("%s\n" %(hubble)) + outFile.write("%s\n" %(redshift)) + outFile.write("%d\n" %(maxKeep)) - numKept = 0 - for line in inFile: - if np.random.uniform() > keepFraction: continue - numKept += 1 - if numKept > maxKeep: break - line = line.split(',') - x = float(line[0]) - y = float(line[1]) - z = float(line[2]) - vz = float(line[3]) + numKept = 0 + for line in inFile: + if np.random.uniform() > keepFraction: continue + numKept += 1 + if numKept > maxKeep: break + line = line.split(',') + x = float(line[0]) + y = float(line[1]) + z = float(line[2]) + vz = float(line[3]) - outFile.write("%e %e %e %e\n" %(x,y,z,vz)) + outFile.write("%e %e %e %e\n" %(x,y,z,vz)) + + outFile.write("-99 -99 -99 -99\n") + inFile.close() + outFile.close() + + elif dataFormat == "random": + sampleName = "ran.ss"+str(thisSubSample)+"_z"+redshift + outFile = open(catalogDir+"/"+sampleName+".dat", 'w') + + outFile.write("%f\n" %(lbox)) + outFile.write("%s\n" %(omegaM)) + outFile.write("%s\n" %(hubble)) + outFile.write("%s\n" %(redshift)) + outFile.write("%d\n" %(maxKeep)) + + for i in xrange(int(maxKeep)): + x = np.random.uniform()*lbox + y = np.random.uniform()*lbox + z = np.random.uniform()*lbox + + outFile.write("%e %e %e %e\n" %(x,y,z, 0.)) + + outFile.write("-99 -99 -99 -99\n") + outFile.close() - outFile.write("-99 -99 -99 -99\n") - print "KEEPING:", numKept, "...predicted:", maxKeep - inFile.close() - outFile.close() # ----------------------------------------------------------------------------- # now halos -if args.script or args.all: +if (args.script or args.all) and dataFormat == "multidark": print " Doing halo scripts" dataFile = catalogDir+"/mdr1_halos_z"+fileNums[0] @@ -362,8 +355,8 @@ if args.halos or args.all: outFile = open(catalogDir+"/"+sampleName+".dat", 'w') outFile.write("%f\n" %(lbox)) - outFile.write("%s" %(omegaM)) - outFile.write("%s" %(hubble)) + outFile.write("%s\n" %(omegaM)) + outFile.write("%s\n" %(hubble)) outFile.write("%s\n" %(redshift)) outFile.write("%d\n" %(numPart)) @@ -425,7 +418,7 @@ BOX_SIZE {boxSize} root_filename hod """ -if args.script or args.all: +if (args.script or args.all) and dataFormat == "multidark": print " Doing DR7 HOD scripts" if dataFormat == "multidark": setName = prefix+"hod_dr72dim2" @@ -467,16 +460,16 @@ if args.hod or args.all: # ----------------------------------------------------------------------------- # now the BOSS HOD -if args.script or args.all: +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, 5, omegaM) - writeScript(prefix, "md.hod_dr9mid_z", + numSubvolumes, numSlices, False, lbox, 15, omegaM) + writeScript(setName, "md.hod_dr9mid_z", scriptDir, catalogDir, fileNums, redshifts, - numSubvolumes, numSlices, True, lbox, 5, omegaM) + numSubvolumes, numSlices, True, lbox, 15, omegaM) if args.hod or args.all: print " Doing DR9 HOD" diff --git a/pipeline/prepareGadgetCatalog.py b/pipeline/prepareGadgetCatalog.py deleted file mode 100755 index 51f67cc..0000000 --- a/pipeline/prepareGadgetCatalog.py +++ /dev/null @@ -1,246 +0,0 @@ -#!/usr/bin/env python - -# script which prepares inputs for gadget-based void run - -import numpy as np -import os -import sys -import void_python_tools as vp -import argparse - -# ----------------------------------------------------------------------------- -# ----------------------------------------------------------------------------- -# CONFIGURATION - -# directory for the input simulation/observational particle files -catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/gadget/" - -# where to put the final void catalog, figures, and output logs -voidOutputDir = os.getenv("HOME")+"/workspace/Voids/gadget/" -figDir = os.getenv("PWD")+"/../figs/gadget/" -logDir = os.getenv("PWD")+"/../logs/gadget/" - -# where to place the pipeline scripts -scriptDir = os.getenv("PWD")+"/gadget/" - -# simulation or observation? -dataType = "simulation" - -# available formats for simulation: gadget, multidark -dataFormat = "gadget" -dataUnit = 1 # as multiple of Mpc/h - -# place particles on the lightcone? -useLightCone = False - -# common filename of particle files -particleFileBase = "mdr1_particles_z" - -# list of file numbers for the particle files -# to get particle file name, we take particleFileBase+fileNum -fileNums = (("0.0", "0.53", "1.0")) - -# redshift of each file in the above list -redshifts = (("0.0", "0.53", "1.0")) - -# how many independent slices along the z-axis? -numSlices = 4 - -# how many subdivisions along the x- and y- axis? -# ( = 2 will make 4 subvolumes for each slice, = 3 will make 9, etc.) -numSubvolumes = 1 - -# prefix to give all outputs -prefix = "gadget_" - -# list of desired subsamples -subSamples = [ 1.0, 0.1 ] - -# simulation information -numPart = 1024*1024*1024 -lbox = 1000 # Mpc/h -omegaM = 0.27 -hubble = 0.70 - -# END CONFIGURATION -# ----------------------------------------------------------------------------- -# ----------------------------------------------------------------------------- - -#------------------------------------------------------------------------------ -LIGHT_SPEED = 299792.458 - -def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): - - sampleName = setName - - if useVel: setName += "_pv" - - if iSlice != -1: sampleName += "_s" + str(iSlice) - - if iVol != -1: sampleName += "_d" + iVol - - return sampleName - -#------------------------------------------------------------------------------ -# 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, - suffix=".dat"): - - if useVel: setName += "_pv" - - scriptFileName = scriptDir + "/" + setName + ".py" - scriptFile = open(scriptFileName, 'w') - - scriptFile.write("""#!/usr/bin/env/python -import os -from void_python_tools.backend.classes import * - -continueRun = False -startCatalogStage = 1 -endCatalogStage = 4 - -startAPStage = 1 -endAPStage = 7 - -ZOBOV_PATH = os.getenv("PWD")+"/../zobov/" -CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/" -freshStack = True -errorBars = "CALCULATED" -numIncoherentRuns = 100 -ranSeed = 101010 -useLCDM = False -bias = 1.16 - -dataPortions = ["central"] -dataSampleList = [] -""") - - dataInfo = """ -setName = "{setName}" - -workDir = "{voidOutputDir}/{setName}/" -inputDataDir = "{inputDataDir}" -figDir = "{figDir}/{setName}/" -logDir = "{logDir}/{setName}/" - """ - scriptFile.write(dataInfo.format(setName=setName, figDir=figDir, - logDir=logDir, voidOutputDir=voidOutputDir, - inputDataDir=catalogDir)) - - sampleInfo = """ -newSample = Sample(dataFile = "{dataFile}", - dataFormat = "{dataFormat}", - dataUnit = {dataUnit}, - fullName = "{sampleName}", - nickName = "{sampleName}", - dataType = "simulation", - 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}", - numSubDivisions = 4, - useLightCone = {useLightCone}, - subsample = {subsample}) -dataSampleList.append(newSample) - """ - - 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 - - # 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) - - 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, redshift, useVel, - iSlice=iSlice, iVol=mySubvolume) - - scriptFile.write(sampleInfo.format(dataFile=dataFileName, - dataFormat=dataFormat, - dataUnit=dataUnit, - sampleName=sampleName, - zMin=sliceMin, - zMax=sliceMax, - zMinMpc=sliceMinMpc, - zMaxMpc=sliceMaxMpc, - omegaM=Om, - boxLen=lbox, - usePecVel=useVel, - minRadius=minRadius, - numSubvolumes=numSubvolumes, - mySubvolume=mySubvolume, - useLightCone=useLightCone, - subsample=subsample)) - - scriptFile.close() - return - - -#------------------------------------------------------------------------------ -#------------------------------------------------------------------------------ -if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir) - -#------------------------------------------------------------------------------ -# 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))) - - print " Doing subsample", thisSubSample, " scripts" - setName = prefix+"ss"+str(thisSubSample) - writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums, - redshifts, - numSubvolumes, numSlices, True, lbox, minRadius, omegaM, - subsample=thisSubSample, suffix="") - writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums, - redshifts, - numSubvolumes, numSlices, False, lbox, minRadius, omegaM, - subsample=thisSubSample, suffix="") - diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 7237408..33ef904 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -106,7 +106,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, else: useLightConeString = "" - if sample.dataFormat == "multidark": + if sample.dataFormat == "multidark" or sample.dataFormat == "random": dataFileLine = "multidark " + datafile elif sample.dataFormat == "gadget": dataFileLine = "gadget " + datafile @@ -197,8 +197,10 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None): if os.access(zobovDir+"/voidDesc_"+sampleName+".out", os.F_OK): os.unlink(zobovDir+"/voidDesc_"+sampleName+".out") - cmd = "%s/vozinit %s 0.1 1.0 2 %s %g %s %s %s >& %s" % \ - (binPath, datafile, \ + numThreads = 2 + + cmd = "%s/vozinit %s 0.1 1.0 %g %s %g %s %s %s >& %s" % \ + (binPath, datafile, numThreads, \ "_"+sampleName, sample.numSubDivisions, \ binPath, zobovDir, maskIndex, logFile) os.system(cmd) From 6868a026f5ba237df76904b674b7da5dd65b6bc5 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 15 Nov 2012 06:47:06 -0600 Subject: [PATCH 23/32] fixed bug where zobov wasn't getting correct nmber of threads; fixed and now hard-coded to 2 parallel threads --- pipeline/datasets/multidark.py | 62 +++++++++++++++++++ .../void_python_tools/backend/launchers.py | 4 +- 2 files changed, 64 insertions(+), 2 deletions(-) create mode 100644 pipeline/datasets/multidark.py diff --git a/pipeline/datasets/multidark.py b/pipeline/datasets/multidark.py new file mode 100644 index 0000000..ae7d8c4 --- /dev/null +++ b/pipeline/datasets/multidark.py @@ -0,0 +1,62 @@ +import os + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- +# CONFIGURATION + +# directory for the input simulation/observational particle files +catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/multidark/" + +# path to HOD code +hodPath = os.getenv("HOME")+"/projects/Voids/hod/HOD.x" + +# where to put the final void catalog, figures, and output logs +voidOutputDir = os.getenv("HOME")+"/workspace/Voids/multidark/" +figDir = os.getenv("PWD")+"/../figs/multidark/" +logDir = os.getenv("PWD")+"/../logs/multidark/" + +# where to place the pipeline scripts +scriptDir = os.getenv("PWD")+"/multidark/" + +# simulation or observation? +dataType = "simulation" + +# available formats for simulation: gadget, multidark +dataFormat = "multidark" +dataUnit = 1 # as multiple of Mpc/h + +# place particles on the lightcone? +useLightCone = True + +# common filename of particle files +particleFileBase = "mdr1_particles_z" + +# list of file numbers for the particle files +# to get particle file name, we take particleFileBase+fileNum +fileNums = (("0.0", "0.53", "1.0")) + +# redshift of each file in the above list +redshifts = (("0.0", "0.53", "1.0")) + +# how many independent slices along the z-axis? +numSlices = 4 + +# how many subdivisions along the x- and y- axis? +# ( = 2 will make 4 subvolumes for each slice, = 3 will make 9, etc.) +numSubvolumes = 1 + +# prefix to give all outputs +prefix = "md_" + +# list of desired subsamples +subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002)) + +# simulation information +numPart = 1024*1024*1024 +lbox = 1000 # Mpc/h +omegaM = 0.27 +hubble = 0.70 + +# END CONFIGURATION +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 33ef904..e7a1fc6 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -200,8 +200,8 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None): numThreads = 2 cmd = "%s/vozinit %s 0.1 1.0 %g %s %g %s %s %s >& %s" % \ - (binPath, datafile, numThreads, \ - "_"+sampleName, sample.numSubDivisions, \ + (binPath, datafile, sample.numSubDivisions, \ + "_"+sampleName, numThreads, \ binPath, zobovDir, maskIndex, logFile) os.system(cmd) From b6359bb0a97855e67771572a2091dab7329076d1 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 15 Nov 2012 07:31:41 -0600 Subject: [PATCH 24/32] number of zobov threads and subdivisions now part of dataset definition --- pipeline/datasets/multidark.py | 9 ++++++++- pipeline/generateCatalog.py | 3 ++- pipeline/prepareCatalogs.py | 8 ++++++-- python_tools/void_python_tools/backend/classes.py | 4 +--- python_tools/void_python_tools/backend/launchers.py | 9 ++++----- 5 files changed, 21 insertions(+), 12 deletions(-) diff --git a/pipeline/datasets/multidark.py b/pipeline/datasets/multidark.py index ae7d8c4..24a5cf9 100644 --- a/pipeline/datasets/multidark.py +++ b/pipeline/datasets/multidark.py @@ -48,9 +48,16 @@ numSubvolumes = 1 # prefix to give all outputs prefix = "md_" -# list of desired subsamples +# list of desired subsamples - these are in unts of h Mpc^-3! +#subSamples = [ 1.0 ] subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002)) +# 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 +numZobovThreads = 2 + # simulation information numPart = 1024*1024*1024 lbox = 1000 # Mpc/h diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index cd8408e..2f200e7 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -88,7 +88,8 @@ for sample in dataSampleList: sys.stdout.flush() launchZobov(sample, ZOBOV_PATH, zobovDir=zobovDir, logDir=logDir, - continueRun=continueRun) + continueRun=continueRun, numZobovDivisions=numZobovDivisions, + numZobovThreads=numZobovThreads) # ------------------------------------------------------------------------- if (startCatalogStage <= 3) and (endCatalogStage >= 3) and not sample.isCombo: diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index d23d36e..9ef7697 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -115,10 +115,15 @@ workDir = "{voidOutputDir}/{setName}/" inputDataDir = "{inputDataDir}" figDir = "{figDir}/{setName}/" logDir = "{logDir}/{setName}/" + +numZobovDivisions = {numZobovDivisions} +numZobovThreads = {numZobovThreads} """ scriptFile.write(dataInfo.format(setName=setName, figDir=figDir, logDir=logDir, voidOutputDir=voidOutputDir, - inputDataDir=catalogDir)) + inputDataDir=catalogDir, + numZobovDivisions=numZobovDivisions, + numZobovThreads=numZobovThreads)) sampleInfo = """ newSample = Sample(dataFile = "{dataFile}", @@ -139,7 +144,6 @@ newSample = Sample(dataFile = "{dataFile}", usePecVel = {usePecVel}, numSubvolumes = {numSubvolumes}, mySubvolume = "{mySubvolume}", - numSubDivisions = 4, useLightCone = {useLightCone}, subsample = {subsample}) dataSampleList.append(newSample) diff --git a/python_tools/void_python_tools/backend/classes.py b/python_tools/void_python_tools/backend/classes.py index 572e543..3eb8b0a 100755 --- a/python_tools/void_python_tools/backend/classes.py +++ b/python_tools/void_python_tools/backend/classes.py @@ -63,7 +63,6 @@ class Sample: usePecVel = False subsample = 1.0 useLightCone = True - numSubDivisions = 1 numSubvolumes = 1 mySubvolume = 1 @@ -75,9 +74,9 @@ class Sample: minVoidRadius=0, fakeDensity=0.01, volumeLimited=True, includeInHubble=True, partOfCombo=False, isCombo=False, comboList=(), profileBinSize=2.0, skyFraction=0.19, - dataType="observation", numSubDivisions=2, boxLen=1024, usePecVel=False, omegaM=0.27, numSubvolumes=1, mySubvolume=1, dataFormat="sdss", + dataType="observation", subsample=1.0, useLightCone=True): self.dataFile = dataFile self.fullName = fullName @@ -98,7 +97,6 @@ class Sample: self.profileBinSize = profileBinSize self.skyFraction = skyFraction self.dataType = dataType - self.numSubDivisions = numSubDivisions self.boxLen = boxLen self.usePecVel = usePecVel self.omegaM = omegaM diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index e7a1fc6..92692f3 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -166,7 +166,8 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, os.unlink(parmFile) # ----------------------------------------------------------------------------- -def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None): +def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, + numZobovDivisions=None, numZobovThreads=None): sampleName = sample.fullName @@ -197,11 +198,9 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None): if os.access(zobovDir+"/voidDesc_"+sampleName+".out", os.F_OK): os.unlink(zobovDir+"/voidDesc_"+sampleName+".out") - numThreads = 2 - cmd = "%s/vozinit %s 0.1 1.0 %g %s %g %s %s %s >& %s" % \ - (binPath, datafile, sample.numSubDivisions, \ - "_"+sampleName, numThreads, \ + (binPath, datafile, numZobovDivisions, \ + "_"+sampleName, numZobovThreads, \ binPath, zobovDir, maskIndex, logFile) os.system(cmd) From 971cf21ea2b56f248bbc8e02a0b69149117d9f1d Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 15 Nov 2012 10:11:51 -0600 Subject: [PATCH 25/32] potential fix so that end of prunevoids is correctly detected --- python_tools/void_python_tools/backend/launchers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 92692f3..e15e8fb 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -297,7 +297,7 @@ def launchPrune(sample, binPath, thisDataPortion=None, cmd += " >& " + logFile os.system(cmd) - if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or jobSuccessful(logFile, "Done!"): + if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or jobSuccessful(logFile, "Done!") or jobSuccessful(logFile, ""): print "done" else: print "FAILED!" From 3ca6699603c2c8fc17b959cf77afd1c43326fd61 Mon Sep 17 00:00:00 2001 From: nhamaus Date: Fri, 16 Nov 2012 09:58:33 +0100 Subject: [PATCH 26/32] fix so that end of prunevoids is correctly detected --- python_tools/void_python_tools/backend/launchers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index e15e8fb..70699b3 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -297,7 +297,7 @@ def launchPrune(sample, binPath, thisDataPortion=None, cmd += " >& " + logFile os.system(cmd) - if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or jobSuccessful(logFile, "Done!") or jobSuccessful(logFile, ""): + if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or jobSuccessful(logFile, "Done!\n"): print "done" else: print "FAILED!" From 0ef211a1ccb450ddd422a6e62dc72b49bd770689 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 16 Nov 2012 09:48:02 -0600 Subject: [PATCH 27/32] some small bug fixes, also prepared prunveVoids to potentially remove voids based on density contrast criteria --- c_tools/stacking/pruneVoids.cpp | 9 ++++++++- pipeline/prepareCatalogs.py | 1 + plotting/plotCompareDist.py | 15 +++++++++++---- .../void_python_tools/backend/launchers.py | 5 +++-- 4 files changed, 23 insertions(+), 7 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 1390f3c..466151a 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -447,6 +447,13 @@ int main(int argc, char **argv) { } for (iVoid = 0; iVoid < numVoids; iVoid++) { +// TEST + //if (voids[iVoid].densCon > 1.5) { + // voids[iVoid].accepted = 0; + //} +// END TEST + + if (strcmp(args_info.dataPortion_arg, "edge") == 0 && tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) { voids[iVoid].accepted = 0; @@ -582,6 +589,6 @@ int main(int argc, char **argv) { clock1 = clock(); - printf("Done!\n"); + printf("Done!"); return 0; } // end main diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index 9ef7697..17c9dc5 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -480,6 +480,7 @@ if args.hod or args.all: for (iRedshift, redshift) in enumerate(redshifts): print " z = ", redshift + # these parameters come from Manera et al 2012, eq. 26 parFileName = "./hod.par" parFile = open(parFileName, 'w') haloFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift] diff --git a/plotting/plotCompareDist.py b/plotting/plotCompareDist.py index b2c47bd..b837173 100755 --- a/plotting/plotCompareDist.py +++ b/plotting/plotCompareDist.py @@ -17,15 +17,20 @@ import argparse workDir = "/home/psutter2/workspace/Voids/" figDir = "./figs" -sampleDirList = [ "multidark/md_halos/sample_md_halos_z0.03_d00/", - "multidark/md_ss0.1/sample_md_ss0.1_z0.03_d00/", - "multidark/md_hod_dr9mid/sample_md_hod_dr9mid_z0.03_d00/", +sampleDirList = [ "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/", + "multidark/md_halos_pv/sample_md_halos_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/", + "multidark/md_ss0.0004_pv/sample_md_ss0.0004_pv_z0.56_d00/", "sdss_dr9/sample_lss.dr9cmassmid.dat/" ] plotNameBase = "compdist" dataPortion = "central" +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, @@ -47,6 +52,7 @@ plt.clf() plt.xlabel("Void Radius (Mpc/h)") plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]") plt.yscale('log') +plt.xlim(xmax=80.) plotName = plotNameBase @@ -60,6 +66,7 @@ for (iSample,sample) in enumerate(dataSampleList): 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]) @@ -84,7 +91,7 @@ for (iSample,sample) in enumerate(dataSampleList): label=lineTitle, color=colorList[iSample], linewidth=linewidth) -plt.legend(title = "Samples", loc = "upper right") +plt.legend(title = "Samples", loc = "upper right", prop={'size':8}) #plt.title(plotTitle) plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight") diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index e15e8fb..be2de94 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -245,7 +245,7 @@ def launchPrune(sample, binPath, thisDataPortion=None, totalPart = open(zobovDir+"/total_particles.txt", "r").read() maxDen = 0.2*float(mockIndex)/float(totalPart) observationLine = " --isObservation" - periodicLine = "--periodic=''" + periodicLine = " --periodic=''" else: mockIndex = -1 maxDen = 0.2 @@ -297,7 +297,8 @@ def launchPrune(sample, binPath, thisDataPortion=None, cmd += " >& " + logFile os.system(cmd) - if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or jobSuccessful(logFile, "Done!") or jobSuccessful(logFile, ""): + if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or \ + jobSuccessful(logFile, "Done!"): print "done" else: print "FAILED!" From a53e3bf2904bec7c46b794ed8843d41794427f12 Mon Sep 17 00:00:00 2001 From: Paul Matt Sutter Date: Sat, 17 Nov 2012 02:49:02 -0500 Subject: [PATCH 28/32] reconfigured to put RA, Dec, redshift, and unique (catalog) ID in the binary ZOBOV file --- c_tools/mock/generateFromCatalog.cpp | 46 +++++++++++++++++++++++----- c_tools/mock/generateMock.cpp | 40 ++++++++++++++++++++++-- c_tools/stacking/pruneVoids.cpp | 10 +++--- external/cosmotool/src/loadSimu.hpp | 3 +- pipeline/prepareCatalogs.py | 12 ++++---- 5 files changed, 90 insertions(+), 21 deletions(-) diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index f0816b0..496f484 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/c_tools/mock/generateFromCatalog.cpp @@ -40,6 +40,7 @@ struct ParticleData vector ra; vector dec; vector redshift; + vector catalogID; int id_mask; // PMS int mask_index; @@ -85,6 +86,7 @@ void loadData(const string& fname, NYU_VData & data) { NYU_Data d; f >> d.index >> d.sector >> d.region >> d.ra >> d.dec >> d.cz >> d.fgotten >> d.phi_z; + d.uniqueID = d.index; data.push_back(d); } } @@ -125,6 +127,7 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data, output_data.ra.resize(data.size()); output_data.dec.resize(data.size()); output_data.redshift.resize(data.size()); + output_data.uniqueID.resize(data.size()); for (int j = 0; j < 3; j++) { @@ -143,7 +146,6 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data, 1.e-6, 1.e-6, &result, &error, &nEval); double Dc = result*LIGHT_SPEED; - cout << "HELLO " << data[i].cz << " " << Dc << endl; p.xyz[0] = Dc*cos(ra)*cos(dec); p.xyz[1] = Dc*sin(ra)*cos(dec); p.xyz[2] = Dc*sin(dec); @@ -157,6 +159,7 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data, output_data.ra[i] = ra; output_data.dec[i] = dec; output_data.redshift[i] = data[i].cz; + output_data.uniqueID[i] = data[i].uniqueID; for (int j = 0; j < 3; j++) { @@ -273,6 +276,7 @@ void generateSurfaceMask(generateFromCatalog_info& args , output_data.ra.push_back(-1); output_data.dec.push_back(-1); output_data.redshift.push_back(-1); + output_data.uniqueID.push_back(-1); //printf("INSERT MOCK %d %e %e %e\n", idx, p.xyz[0], p.xyz[1], p.xyz[2]); insertion++; } @@ -310,6 +314,7 @@ void generateSurfaceMask(generateFromCatalog_info& args , output_data.ra.push_back(-1); output_data.dec.push_back(-1); output_data.redshift.push_back(-1); + output_data.uniqueID.push_back(-1); insertion++; fprintf(fp, "%e %e %e\n", @@ -341,6 +346,7 @@ void generateSurfaceMask(generateFromCatalog_info& args , output_data.ra.push_back(-1); output_data.dec.push_back(-1); output_data.redshift.push_back(-1); + output_data.uniqueID.push_back(-1); insertion++; fprintf(fp, "%e %e %e\n", (p.xyz[0]), @@ -357,6 +363,7 @@ void generateSurfaceMask(generateFromCatalog_info& args , output_data.ra.push_back(-1); output_data.dec.push_back(-1); output_data.redshift.push_back(-1); + output_data.uniqueID.push_back(-1); insertion++; fprintf(fp, "%e %e %e\n", (p.xyz[0]), @@ -414,11 +421,40 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn for (uint32_t i = 0; i < pdata.pos.size(); i++) { f.writeReal32((pdata.pos[i].xyz[j]+Lmax)/(2*Lmax)); -if (i < 10) printf("TEST WRITE %d %e\n", (pdata.pos[i].xyz[j]+Lmax)/(2*Lmax)); } f.endCheckpoint(); } + cout << format("Writing RA...") << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < pdata.pos.size(); i++) { + f.writeReal32(pdata.pos[i].ra); + } + f.endCheckpoint(); + + cout << format("Writing Dec...") << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < pdata.pos.size(); i++) { + f.writeReal32(pdata.pos[i].Dec); + } + f.endCheckpoint(); + + cout << format("Writing Redshift...") << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < pdata.pos.size(); i++) { + f.writeReal32(pdata.pos[i].cz); + } + f.endCheckpoint(); + + cout << format("Writing Unique ID...") << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < pdata.pos.size(); i++) { + f.writeReal32(pdata.pos[i].uniqueID); + } + f.endCheckpoint(); + + NcFile fp(paramname.c_str(), NcFile::Replace); + NcFile fp(paramname.c_str(), NcFile::Replace); NcFile fp(paramname.c_str(), NcFile::Replace); fp.add_att("range_x_min", -Lmax/100.); @@ -434,9 +470,6 @@ if (i < 10) printf("TEST WRITE %d %e\n", (pdata.pos[i].xyz[j]+Lmax)/(2*Lmax)); NcDim *NumPart_dim = fp.add_dim("numpart_dim", nOutputPart); NcVar *v = fp.add_var("particle_ids", ncInt, NumPart_dim); - NcVar *vra = fp.add_var("RA", ncInt, NumPart_dim); - NcVar *vdec = fp.add_var("DEC", ncInt, NumPart_dim); - NcVar *vredshift = fp.add_var("z", ncInt, NumPart_dim); //NcVar *v2 = fp.add_var("expansion", ncDouble, NumPart_dim); //double *expansion_fac = new double[pdata.pos.size()]; @@ -445,9 +478,6 @@ if (i < 10) printf("TEST WRITE %d %e\n", (pdata.pos[i].xyz[j]+Lmax)/(2*Lmax)); // expansion_fac[i] = 1.0; v->put(&pdata.id_gal[0], nOutputPart); - vra->put(&pdata.ra[0], nOutputPart); - vdec->put(&pdata.dec[0], nOutputPart); - vredshift->put(&pdata.redshift[0], nOutputPart); //v2->put(expansion_fac, pdata.pos.size()); //delete[] expansion_fac; diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index dcf1acb..06954c0 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -178,14 +178,17 @@ SimuData *doLoadMultidark(const char *multidarkname) for (int k = 0; k < 3; k++) outd->Pos[k] = new float[outd->NumPart]; outd->Vel[2] = new float[outd->NumPart]; + outd->uniqueID = new float[outd->NumPart]; cout << "loading multidark particles" << endl; actualNumPart = 0; for (int i = 0; i < outd->NumPart; i++) { - fscanf(fp, "%f %f %f %f\n", &outd->Pos[0][i], &outd->Pos[1][i], + fscanf(fp, "%d %f %f %f %f\n", &outd->uniqueID[i], + &outd->Pos[0][i], &outd->Pos[1][i], &outd->Pos[2][i], &outd->Vel[2][i]); - if (outd->Pos[0][i] == -99 && outd->Pos[1][i] == -99 && + if (outd->uniqueID[i] == -99 && + outd->Pos[0][i] == -99 && outd->Pos[1][i] == -99 && outd->Pos[2][i] == -99 && outd->Vel[2][i] == -99) { break; } else { @@ -361,6 +364,39 @@ void generateOutput(SimuData *data, int axis, } f.endCheckpoint(); + cout << "Writing RA..." << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < data->NumPart; i++) + { + f.writeReal32(data->uniqueID[i]); + } + f.endCheckpoint(); + + cout << "Writing Dec..." << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < data->NumPart; i++) + { + f.writeReal32(data->uniqueID[i]); + } + f.endCheckpoint(); + + cout << "Writing redshift..." << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < data->NumPart; i++) + { + f.writeReal32(data->uniqueID[i]); + } + f.endCheckpoint(); + + cout << "Writing unique ID..." << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < data->NumPart; i++) + { + f.writeReal32(data->uniqueID[i]); + } + f.endCheckpoint(); + + } void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& args_info) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 466151a..4cb065e 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -529,7 +529,7 @@ int main(int argc, char **argv) { outCenter[2] = (voids[iVoid].barycenter[2]-boxLen[2]/2.)*100.; } - fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d\n", + fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f\n", outCenter[0], outCenter[1], outCenter[2], @@ -537,7 +537,8 @@ int main(int argc, char **argv) { voids[iVoid].radius, voids[iVoid].redshift, 4./3.*M_PI*pow(voids[iVoid].radius, 3), - voids[iVoid].voidID + voids[iVoid].voidID, + voids[iVoid].densCon, ); fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n", @@ -572,7 +573,7 @@ int main(int argc, char **argv) { } - fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d\n", + fprintf(fpInfo, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f\n", outCenter[0], outCenter[1], outCenter[2], @@ -580,7 +581,8 @@ int main(int argc, char **argv) { voids[iVoid].radius, voids[iVoid].redshift, 4./3.*M_PI*pow(voids[iVoid].radius, 3), - voids[iVoid].voidID); + voids[iVoid].voidID, + voids[iVoid].densCon); } fclose(fpInfo); diff --git a/external/cosmotool/src/loadSimu.hpp b/external/cosmotool/src/loadSimu.hpp index a7db941..1aa5c2a 100644 --- a/external/cosmotool/src/loadSimu.hpp +++ b/external/cosmotool/src/loadSimu.hpp @@ -24,9 +24,10 @@ namespace CosmoTool int *Id; float *Pos[3]; float *Vel[3]; + float *uniqueID; int *type; public: - SimuData() : Id(0),NumPart(0),type(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; } + SimuData() : Id(0),NumPart(0),type(0), uniqueID(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; uniqueID=0} ~SimuData() { for (int j = 0; j < 3; j++) diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index 17c9dc5..49e67b5 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -284,7 +284,7 @@ for thisSubSample in subSamples: outFile.write("%d\n" %(maxKeep)) numKept = 0 - for line in inFile: + for (i,line) in enumerate(inFile): if np.random.uniform() > keepFraction: continue numKept += 1 if numKept > maxKeep: break @@ -294,9 +294,9 @@ for thisSubSample in subSamples: z = float(line[2]) vz = float(line[3]) - outFile.write("%e %e %e %e\n" %(x,y,z,vz)) + outFile.write("%d %e %e %e %e\n" %(i,x,y,z,vz)) - outFile.write("-99 -99 -99 -99\n") + outFile.write("-99 -99 -99 -99 -99\n") inFile.close() outFile.close() @@ -315,9 +315,9 @@ for thisSubSample in subSamples: y = np.random.uniform()*lbox z = np.random.uniform()*lbox - outFile.write("%e %e %e %e\n" %(x,y,z, 0.)) + outFile.write("%d %e %e %e %e\n" % (i, x,y,z, 0.)) - outFile.write("-99 -99 -99 -99\n") + outFile.write("-99 -99 -99 -99 -99\n") outFile.close() @@ -375,7 +375,7 @@ if args.halos or args.all: vz = float(line[5]) # write to output file - outFile.write("%e %e %e %e\n" %(x,y,z,vz)) + outFile.write("%d %e %e %e %e\n" %(numKept,x,y,z,vz)) inFile.close() outFile.close() From 10dfe29a26119767c7cbf42f4124df089de46d85 Mon Sep 17 00:00:00 2001 From: Paul Matt Sutter Date: Sat, 17 Nov 2012 13:00:54 -0500 Subject: [PATCH 29/32] Supports minimum halo mass cuts. Start of scripts to generate masked mock sets, some files added to later support more general preparation scripts --- pipeline/applyMaskToMock.py | 321 ++++++++++++++++++ pipeline/datasets/mock_dr9mid.py | 62 ++++ pipeline/datasets/multidark.py | 7 + pipeline/generateCatalog.py | 3 + pipeline/prepareCatalogs.py | 98 +++--- plotting/datasetsToPlot.py | 17 + plotting/plotCompareDensCon.py | 93 +++++ plotting/plotCompareDist.py | 13 +- .../void_python_tools/backend/__init__.py | 1 + 9 files changed, 557 insertions(+), 58 deletions(-) create mode 100755 pipeline/applyMaskToMock.py create mode 100644 pipeline/datasets/mock_dr9mid.py create mode 100755 plotting/datasetsToPlot.py create mode 100755 plotting/plotCompareDensCon.py diff --git a/pipeline/applyMaskToMock.py b/pipeline/applyMaskToMock.py new file mode 100755 index 0000000..c3fcd18 --- /dev/null +++ b/pipeline/applyMaskToMock.py @@ -0,0 +1,321 @@ +#!/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 + +import numpy as np +import os +import sys +import void_python_tools as vp +import argparse +import imp +import healpy as hp + +# ------------------------------------------------------------------------------ + +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 + +parser = argparse.ArgumentParser(description='options') +parser.add_argument('--scripts', dest='script', action='store_const', + const=True, default=False, + help='write scripts') +parser.add_argument('--parmFile', dest='parmFile', + default="", + 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)) + + +#------------------------------------------------------------------------------ +def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): + + sampleName = setName + + sampleName += "_z" + redshift + + if iVol != -1: sampleName += "_d" + iVol + + return sampleName + +#------------------------------------------------------------------------------ +# 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, + suffix=".dat"): + + + if useVel: setName += "_pv" + + scriptFileName = scriptDir + "/" + setName + ".py" + scriptFile = open(scriptFileName, 'w') + + scriptFile.write("""#!/usr/bin/env/python +import os +from void_python_tools.backend.classes import * + +continueRun = True # set to True to enable restarting aborted jobs +startCatalogStage = 1 +endCatalogStage = 4 + +startAPStage = 1 +endAPStage = 7 + +ZOBOV_PATH = os.getenv("PWD")+"/../zobov/" +CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/" +freshStack = True +errorBars = "CALCULATED" +numIncoherentRuns = 100 +ranSeed = 101010 +useLCDM = False +bias = 1.16 + +dataPortions = ["central"] +dataSampleList = [] +""") + + + dataInfo = """ +setName = "{setName}" + +workDir = "{voidOutputDir}/{setName}/" +inputDataDir = "{inputDataDir}" +figDir = "{figDir}/{setName}/" +logDir = "{logDir}/{setName}/" + +numZobovDivisions = {numZobovDivisions} +numZobovThreads = {numZobovThreads} + """ + scriptFile.write(dataInfo.format(setName=setName, figDir=figDir, + logDir=logDir, voidOutputDir=voidOutputDir, + inputDataDir=catalogDir, + numZobovDivisions=numZobovDivisions, + numZobovThreads=numZobovThreads)) + + sampleInfo = """ +newSample = Sample(dataFile = "{dataFile}", + dataFormat = "{dataFormat}", + dataUnit = {dataUnit}, + fullName = "{sampleName}", + nickName = "{sampleName}", + dataType = "simulation", + 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}) +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) + """ + 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 + + # 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) + + 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, + dataFormat=dataFormat, + dataUnit=dataUnit, + sampleName=sampleName, + zMin=sliceMin, + zMax=sliceMax, + zMinMpc=sliceMinMpc, + zMaxMpc=sliceMaxMpc, + omegaM=Om, + boxLen=lbox, + usePecVel=useVel, + minRadius=minRadius, + numSubvolumes=numSubvolumes, + mySubvolume=mySubvolume, + useLightCone=useLightCone, + subsample=subsample)) + + 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) + +inFile = open('hod.mock', 'r') +outFile = open(catalogDir+"/mock.out")) + +zBox = float(redshiftRange[0]) +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]) + + 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) + + if redshift < redshiftRange[0] or redshift > redshiftRange[1]: continue + + 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 + + phi = np.pi/180. * RA + theta = np.pi/2. - Dec*np.pi/180. + pos = np.zeros((3)) + + pix = hp.ang2pix(nside, theta, phi) + if mask[pix] <= 0: continue + + print >> outFile, RA, Dec, redshift*LIGHT_SPEED, 0., x, y, z + +inFile.close() +outFile.close() + +os.system("rm ./hod.*") + + diff --git a/pipeline/datasets/mock_dr9mid.py b/pipeline/datasets/mock_dr9mid.py new file mode 100644 index 0000000..cd189f3 --- /dev/null +++ b/pipeline/datasets/mock_dr9mid.py @@ -0,0 +1,62 @@ +import os + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- +# CONFIGURATION + +# directory for the input simulation/observational particle files +catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/mock_dr9mid/" + +# path to HOD code +hodPath = os.getenv("HOME")+"/projects/Voids/hod/HOD.x" + +# path to mask +maskFile = os.getenv("HOME")+"/workspace/Voids/catalogs/boss/final_boss_mask.fits") + +# where to put the final void catalog, figures, and output logs +voidOutputDir = os.getenv("HOME")+"/workspace/Voids/mock_dr9mid/" +figDir = os.getenv("PWD")+"/../figs/mock_dr9mid/" +logDir = os.getenv("PWD")+"/../logs/mock_dr9mid/" + +# where to place the pipeline scripts +scriptDir = os.getenv("PWD")+"/mock_dr9mid/" + +# simulation or observation? +dataType = "observation" + +# available formats for simulation: gadget, multidark +dataFormat = "multidark" +dataUnit = 1 # as multiple of Mpc/h + +# place particles on the lightcone? +useLightCone = True + +# list of file numbers for the particle files +# to get particle file name, we take particleFileBase+fileNum +fileNums = (("0.53")) + +# redshift range of the mock +redshiftRange = (0.53, 0.6) + +# prefix to give all outputs +prefix = "mock_" + +# common filename of halo files +haloFileBase = "mdr1_halos_z" + +# 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 = 2 +numZobovThreads = 2 + +# simulation information +numPart = 1024*1024*1024 +lbox = 1000 # Mpc/h +omegaM = 0.27 +hubble = 0.70 + + +# END CONFIGURATION +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- diff --git a/pipeline/datasets/multidark.py b/pipeline/datasets/multidark.py index 24a5cf9..a6a2563 100644 --- a/pipeline/datasets/multidark.py +++ b/pipeline/datasets/multidark.py @@ -52,6 +52,13 @@ prefix = "md_" #subSamples = [ 1.0 ] subSamples = ((0.1, 0.05, 0.01, 0.002, 0.001, 0.0004, 0.0002)) +# common filename of halo files +haloFileBase = "mdr1_halos_z" + +# minimum halo mass cuts to apply for the halo catalog +# use "none" to get all halos +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? diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 2f200e7..de37da7 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -64,6 +64,9 @@ for sample in dataSampleList: # save this sample's information with open(zobovDir+"/sample_info.dat", 'wb') as output: pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL) + fp = open(zobovDir+"/sample_info.txt", 'w') + fp.write("Redshift range: %f - %f" %(sample.zBoundary[0], sample.zBoundary[1]) + fp.close() # --------------------------------------------------------------------------- if (startCatalogStage <= 1) and (endCatalogStage >= 1) and not sample.isCombo: diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index 49e67b5..b149401 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -60,12 +60,8 @@ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): sampleName = setName - #if useVel: sampleName += "_pv" - sampleName += "_z" + redshift - #if iSlice != -1: sampleName += "_s" + str(iSlice) - if iVol != -1: sampleName += "_d" + iVol return sampleName @@ -87,7 +83,7 @@ def writeScript(setName, dataFileNameBase, import os from void_python_tools.backend.classes import * -continueRun = False +continueRun = True # set to True to enable restarting aborted jobs startCatalogStage = 1 endCatalogStage = 4 @@ -326,59 +322,69 @@ for thisSubSample in subSamples: if (args.script or args.all) and dataFormat == "multidark": print " Doing halo scripts" - dataFile = catalogDir+"/mdr1_halos_z"+fileNums[0] - inFile = open(dataFile, 'r') - numPart = 0 - for line in inFile: numPart += 1 - inFile.close() + for minHaloMass in minHaloMasses: + dataFile = catalogDir+haloFileBase+fileNums[0] + inFile = open(dataFile, 'r') + numPart = 0 + for line in inFile: + line = line.split(',') + if numHaloMass == "none" or float(line[6]) > minHaloMass: + numPart += 1 + inFile.close() - minRadius = 2*int(np.ceil(lbox/numPart**(1./3.))) - - if dataFormat == "multidark": - setName = prefix+"halos" - writeScript(setName, "md.halos_z", scriptDir, catalogDir, fileNums, - redshifts, - numSubvolumes, numSlices, False, lbox, minRadius, omegaM) - writeScript(setName, "md.halos_z", scriptDir, catalogDir, fileNums, - redshifts, numSubvolumes, + minRadius = 2*int(np.ceil(lbox/numPart**(1./3.))) + + 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, numSlices, True, lbox, minRadius, omegaM) if args.halos or args.all: print " Doing halos" - for (iRedshift, redshift) in enumerate(redshifts): - print " z = ", redshift + for minHaloMass in minHaloMasses: + print " min halo mass = ", minHaloMass - dataFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift] - inFile = open(dataFile, 'r') - numPart = 0 - for line in inFile: numPart += 1 - inFile.close() + for (iRedshift, redshift) in enumerate(redshifts): + print " z = ", redshift - sampleName = "md.halos_z"+redshift - outFile = open(catalogDir+"/"+sampleName+".dat", 'w') + dataFile = catalogDir+haloFileBase+fileNums[iRedshift] + inFile = open(dataFile, 'r') + for line in inFile: + line = line.split(',') + if numHaloMass == "none" or float(line[6]) > minHaloMass: + inFile.close() - outFile.write("%f\n" %(lbox)) - outFile.write("%s\n" %(omegaM)) - outFile.write("%s\n" %(hubble)) - outFile.write("%s\n" %(redshift)) - outFile.write("%d\n" %(numPart)) + sampleName = "md.halos_z"+redshift + outFile = open(catalogDir+"/"+sampleName+".dat", 'w') - inFile = open(dataFile, 'r') - numKept = 0 - for line in inFile: - numKept += 1 - line = line.split(',') - x = float(line[0]) - y = float(line[1]) - z = float(line[2]) - vz = float(line[5]) + outFile.write("%f\n" %(lbox)) + outFile.write("%s\n" %(omegaM)) + outFile.write("%s\n" %(hubble)) + outFile.write("%s\n" %(redshift)) + outFile.write("%d\n" %(numPart)) - # write to output file - outFile.write("%d %e %e %e %e\n" %(numKept,x,y,z,vz)) + inFile = open(dataFile, 'r') + numKept = 0 + for line in inFile: + if numHaloMass == "none" or float(line[6]) > minHaloMass: + numKept += 1 + line = line.split(',') + x = float(line[0]) + y = float(line[1]) + z = float(line[2]) + vz = float(line[5]) - inFile.close() - outFile.close() + # write to output file + outFile.write("%d %e %e %e %e\n" %(numKept,x,y,z,vz)) + + inFile.close() + outFile.close() # ----------------------------------------------------------------------------- # now the SDSS HOD diff --git a/plotting/datasetsToPlot.py b/plotting/datasetsToPlot.py new file mode 100755 index 0000000..1217ed8 --- /dev/null +++ b/plotting/datasetsToPlot.py @@ -0,0 +1,17 @@ +#!/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_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/", + "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/", + "multidark/md_ss0.0004_pv/sample_md_ss0.0004_pv_z0.56_d00/", + "sdss_dr9/sample_lss.dr9cmassmid.dat/" ] + +dataPortion = "central" + diff --git a/plotting/plotCompareDensCon.py b/plotting/plotCompareDensCon.py new file mode 100755 index 0000000..8e21a8a --- /dev/null +++ b/plotting/plotCompareDensCon.py @@ -0,0 +1,93 @@ +#!/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 + +# ------------------------------------------------------------------------------ + +from datasetsToPlot import * + +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)') +args = parser.parse_args() + +# ------------------------------------------------------------------------------ + +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 Radius (Mpc/h)") +plt.ylabel(r"N > R [$h^3$ Gpc$^{-3}$]") +plt.yscale('log') +plt.xlim(xmax=80.) + +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 + + 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/plotting/plotCompareDist.py b/plotting/plotCompareDist.py index b837173..77804fb 100755 --- a/plotting/plotCompareDist.py +++ b/plotting/plotCompareDist.py @@ -14,21 +14,10 @@ import argparse # ------------------------------------------------------------------------------ -workDir = "/home/psutter2/workspace/Voids/" -figDir = "./figs" - -sampleDirList = [ "multidark/md_ss0.1_pv/sample_md_ss0.1_pv_z0.56_d00/", - "multidark/md_halos_pv/sample_md_halos_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/", - "multidark/md_ss0.0004_pv/sample_md_ss0.0004_pv_z0.56_d00/", - "sdss_dr9/sample_lss.dr9cmassmid.dat/" ] +from datasetsToPlot import * plotNameBase = "compdist" -dataPortion = "central" - obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing? parser = argparse.ArgumentParser(description='Plot.') diff --git a/python_tools/void_python_tools/backend/__init__.py b/python_tools/void_python_tools/backend/__init__.py index 045ef20..95cddfa 100644 --- a/python_tools/void_python_tools/backend/__init__.py +++ b/python_tools/void_python_tools/backend/__init__.py @@ -1,2 +1,3 @@ from classes import * from launchers import * +from catalogPrep import * From 90964d2be3273310214b7e0210d69035d6b600d7 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Sat, 17 Nov 2012 13:01:39 -0600 Subject: [PATCH 30/32] fixed bugs associated with saving RA and Dec --- c_tools/mock/generateFromCatalog.cpp | 13 +++++++------ c_tools/mock/generateMock.cpp | 14 +++++++------- c_tools/stacking/pruneVoids.cpp | 5 ++--- external/cosmotool/src/loadSimu.hpp | 2 ++ pipeline/generateCatalog.py | 8 ++++---- .../void_python_tools/backend/launchers.py | 6 +----- 6 files changed, 23 insertions(+), 25 deletions(-) diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index 496f484..1024f48 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/c_tools/mock/generateFromCatalog.cpp @@ -27,6 +27,7 @@ struct NYU_Data double cz; double fgotten; double phi_z; + double uniqueID; }; struct Position @@ -41,6 +42,7 @@ struct ParticleData vector dec; vector redshift; vector catalogID; + vector uniqueID; int id_mask; // PMS int mask_index; @@ -408,6 +410,7 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn UnformattedWrite f(fname); static const char axis[] = { 'X', 'Y', 'Z' }; double Lmax = pdata.Lmax; + double r2d = 180./M_PI; f.beginCheckpoint(); f.writeInt32(pdata.pos.size()); @@ -428,33 +431,31 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn cout << format("Writing RA...") << endl; f.beginCheckpoint(); for (uint32_t i = 0; i < pdata.pos.size(); i++) { - f.writeReal32(pdata.pos[i].ra); + f.writeReal32(pdata.ra[i]*r2d); } f.endCheckpoint(); cout << format("Writing Dec...") << endl; f.beginCheckpoint(); for (uint32_t i = 0; i < pdata.pos.size(); i++) { - f.writeReal32(pdata.pos[i].Dec); + f.writeReal32(pdata.dec[i]*r2d); } f.endCheckpoint(); cout << format("Writing Redshift...") << endl; f.beginCheckpoint(); for (uint32_t i = 0; i < pdata.pos.size(); i++) { - f.writeReal32(pdata.pos[i].cz); + f.writeReal32(pdata.redshift[i]); } f.endCheckpoint(); cout << format("Writing Unique ID...") << endl; f.beginCheckpoint(); for (uint32_t i = 0; i < pdata.pos.size(); i++) { - f.writeReal32(pdata.pos[i].uniqueID); + f.writeReal32(pdata.uniqueID[i]); } f.endCheckpoint(); - NcFile fp(paramname.c_str(), NcFile::Replace); - NcFile fp(paramname.c_str(), NcFile::Replace); NcFile fp(paramname.c_str(), NcFile::Replace); fp.add_att("range_x_min", -Lmax/100.); diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index 06954c0..06d06d3 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -178,16 +178,16 @@ SimuData *doLoadMultidark(const char *multidarkname) for (int k = 0; k < 3; k++) outd->Pos[k] = new float[outd->NumPart]; outd->Vel[2] = new float[outd->NumPart]; - outd->uniqueID = new float[outd->NumPart]; + outd->Id = new int[outd->NumPart]; cout << "loading multidark particles" << endl; actualNumPart = 0; for (int i = 0; i < outd->NumPart; i++) { - fscanf(fp, "%d %f %f %f %f\n", &outd->uniqueID[i], + fscanf(fp, "%d %d %f %f %f\n", &outd->Id[i], &outd->Pos[0][i], &outd->Pos[1][i], &outd->Pos[2][i], &outd->Vel[2][i]); - if (outd->uniqueID[i] == -99 && + 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) { break; @@ -368,7 +368,7 @@ void generateOutput(SimuData *data, int axis, f.beginCheckpoint(); for (uint32_t i = 0; i < data->NumPart; i++) { - f.writeReal32(data->uniqueID[i]); + f.writeReal32(data->Id[i]); } f.endCheckpoint(); @@ -376,7 +376,7 @@ void generateOutput(SimuData *data, int axis, f.beginCheckpoint(); for (uint32_t i = 0; i < data->NumPart; i++) { - f.writeReal32(data->uniqueID[i]); + f.writeReal32(data->Id[i]); } f.endCheckpoint(); @@ -384,7 +384,7 @@ void generateOutput(SimuData *data, int axis, f.beginCheckpoint(); for (uint32_t i = 0; i < data->NumPart; i++) { - f.writeReal32(data->uniqueID[i]); + f.writeReal32(data->Id[i]); } f.endCheckpoint(); @@ -392,7 +392,7 @@ void generateOutput(SimuData *data, int axis, f.beginCheckpoint(); for (uint32_t i = 0; i < data->NumPart; i++) { - f.writeReal32(data->uniqueID[i]); + f.writeReal32(data->Id[i]); } f.endCheckpoint(); diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 4cb065e..b50de5b 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -538,8 +538,7 @@ int main(int argc, char **argv) { voids[iVoid].redshift, 4./3.*M_PI*pow(voids[iVoid].radius, 3), voids[iVoid].voidID, - voids[iVoid].densCon, - ); + voids[iVoid].densCon); fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n", atan((voids[iVoid].barycenter[1]-boxLen[1]/2.)/(voids[iVoid].barycenter[0]-boxLen[0]/2.)) * 180/M_PI + 180, @@ -591,6 +590,6 @@ int main(int argc, char **argv) { clock1 = clock(); - printf("Done!"); + printf("Done!\n"); return 0; } // end main diff --git a/external/cosmotool/src/loadSimu.hpp b/external/cosmotool/src/loadSimu.hpp index 1aa5c2a..35567f2 100644 --- a/external/cosmotool/src/loadSimu.hpp +++ b/external/cosmotool/src/loadSimu.hpp @@ -41,6 +41,8 @@ namespace CosmoTool delete[] type; if (Id) delete[] Id; + if (uniqueID) + delete[] uniqueID; } }; diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index de37da7..9ebd693 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -62,11 +62,11 @@ for sample in dataSampleList: os.makedirs(zobovDir) # save this sample's information - with open(zobovDir+"/sample_info.dat", 'wb') as output: + with open(zobovDir+"/sample_info.dat", 'w') as output: pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL) - fp = open(zobovDir+"/sample_info.txt", 'w') - fp.write("Redshift range: %f - %f" %(sample.zBoundary[0], sample.zBoundary[1]) - fp.close() + #fp = open(zobovDir+"/sample_info.txt", 'w') + #fp.write("Redshift range: %f - %f" %(sample.zBoundary[0], sample.zBoundary[1]) + #fp.close() # --------------------------------------------------------------------------- if (startCatalogStage <= 1) and (endCatalogStage >= 1) and not sample.isCombo: diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index ae425c2..280f568 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -297,12 +297,8 @@ def launchPrune(sample, binPath, thisDataPortion=None, cmd += " >& " + logFile os.system(cmd) -<<<<<<< HEAD if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or \ - jobSuccessful(logFile, "Done!"): -======= - if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or jobSuccessful(logFile, "Done!\n"): ->>>>>>> 3d6fae9882fde44515a59e1144646d77bbd1a124 + jobSuccessful(logFile, "Done!\n"): print "done" else: print "FAILED!" From 0b54a31d0ca356917305baeba293aa17249cb769 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Sun, 18 Nov 2012 14:23:50 -0600 Subject: [PATCH 31/32] many bug fixes from the previous update; new changes appear to work now --- c_tools/libzobov/contour_pixels.cpp | 30 +++++++++++++++++++ c_tools/libzobov/contour_pixels.hpp | 1 + c_tools/mock/generateFromCatalog.cpp | 12 ++++---- c_tools/mock/generateMock.cpp | 2 ++ pipeline/generateCatalog.py | 20 +++++++++++-- pipeline/prepareCatalogs.py | 24 ++++++++------- .../void_python_tools/backend/launchers.py | 4 +-- 7 files changed, 72 insertions(+), 21 deletions(-) diff --git a/c_tools/libzobov/contour_pixels.cpp b/c_tools/libzobov/contour_pixels.cpp index 6486ae7..6aa7935 100644 --- a/c_tools/libzobov/contour_pixels.cpp +++ b/c_tools/libzobov/contour_pixels.cpp @@ -54,3 +54,33 @@ void computeContourPixels(Healpix_Map& m, vector& contour) write_Healpix_map_to_fits(h, contour_map, planckType()); } } + +void computeMaskPixels(Healpix_Map& m, vector& contour) +{ + for (int p = 0; p < m.Npix(); p++) + { + + if (m[p]>0) + { + contour.push_back(p); + // This is boundary go to next pixel + } + } + + if (DEBUG) + { + Healpix_Map contour_map; + + contour_map.SetNside(m.Nside(), RING); + contour_map.fill(0); + for (int p = 0; p < contour.size(); p++) + { + contour_map[contour[p]]=1; + } + + fitshandle h; + h.create("!mask_map.fits"); + write_Healpix_map_to_fits(h, contour_map, planckType()); + } +} + diff --git a/c_tools/libzobov/contour_pixels.hpp b/c_tools/libzobov/contour_pixels.hpp index d3a893b..c576477 100644 --- a/c_tools/libzobov/contour_pixels.hpp +++ b/c_tools/libzobov/contour_pixels.hpp @@ -7,5 +7,6 @@ void computeContourPixels(Healpix_Map& map, std::vector& contour); void computeFilledPixels(Healpix_Map& map, std::vector& contour); +void computeMaskPixels(Healpix_Map& map, std::vector& contour); #endif diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index 1024f48..1d49da3 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/c_tools/mock/generateFromCatalog.cpp @@ -295,10 +295,10 @@ void generateSurfaceMask(generateFromCatalog_info& args , int nPart = 100; // TEST - for (int iDir = 0; iDir < 0; iDir++) { - for (int iFace = 0; iFace < 0; iFace++) { - //for (int iDir = 0; iDir < 3; iDir++) { - //for (int iFace = 0; iFace < 2; iFace++) { + //for (int iDir = 0; iDir < 0; iDir++) { + //for (int iFace = 0; iFace < 0; iFace++) { + for (int iDir = 0; iDir < 3; iDir++) { + for (int iFace = 0; iFace < 2; iFace++) { int iy = (iDir + 1) % 3; int iz = (iDir + 2) % 3; @@ -483,6 +483,7 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn //delete[] expansion_fac; +/* FILE *infoFile = fopen("sample_info.txt", "w"); fprintf(infoFile, "x_min = %f\n", -Lmax/100.); fprintf(infoFile, "x_max = %f\n", Lmax/100.); @@ -493,7 +494,7 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn fprintf(infoFile, "mask_index = %d\n", pdata.mask_index); fprintf(infoFile, "total_particles = %d\n", pdata.pos.size()); fclose(infoFile); - +*/ } @@ -545,6 +546,7 @@ int main(int argc, char **argv) mask.Import(o_mask); computeContourPixels(mask,pixel_list); + computeMaskPixels(mask,full_mask_list); // We compute a cube holding all the galaxies + the survey surface mask diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index 06d06d3..7dacef6 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -505,6 +505,7 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& delete[] expansion_fac; +/* FILE *fp = fopen("sample_info.txt", "w"); fprintf(fp, "x_min = %f\n", ranges[0][0]); fprintf(fp, "x_max = %f\n", ranges[0][1]); @@ -515,6 +516,7 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& fprintf(fp, "mask_index = -1\n"); fprintf(fp, "total_particles = %d\n", boxed->NumPart); fclose(fp); +*/ } void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, generateMock_info& args_info) diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 9ebd693..b60a6db 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -64,9 +64,23 @@ for sample in dataSampleList: # save this sample's information with open(zobovDir+"/sample_info.dat", 'w') as output: pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL) - #fp = open(zobovDir+"/sample_info.txt", 'w') - #fp.write("Redshift range: %f - %f" %(sample.zBoundary[0], sample.zBoundary[1]) - #fp.close() + + fp = open(zobovDir+"/sample_info.txt", 'w') + fp.write("Sample name: %s\n" % sample.fullName) + fp.write("Sample nickname: %s\n" % sample.nickName) + fp.write("Data type: %s\n" % sample.dataType) + fp.write("Redshift range: %f - %f\n" %(sample.zBoundary[0],sample.zBoundary[1])) + fp.write("Estimated mean particle separation: %g\n" % sample.minVoidRadius) + + if (sample.dataType == "simulation"): + fp.write("Particles placed on lightcone: %g\n" % sample.useLightCone) + fp.write("Peculiar velocities included: %g\n" % sample.usePecVel) + fp.write("Additional subsampling fraction: %g\n" % sample.subsample) + fp.write("Simulation box length (Mpc/h): %g\n" % sample.boxLen) + fp.write("Simulation Omega_M: %g\n" % sample.omegaM) + fp.write("Number of simulation subvolumes: %g\n", sample.numSubvolumes) + fp.write("My subvolume index: %g\n", sample.mySubvolume) + fp.close() # --------------------------------------------------------------------------- if (startCatalogStage <= 1) and (endCatalogStage >= 1) and not sample.isCombo: diff --git a/pipeline/prepareCatalogs.py b/pipeline/prepareCatalogs.py index b149401..13a96e9 100755 --- a/pipeline/prepareCatalogs.py +++ b/pipeline/prepareCatalogs.py @@ -323,26 +323,28 @@ if (args.script or args.all) and dataFormat == "multidark": print " Doing halo scripts" for minHaloMass in minHaloMasses: + # estimate number of halos to get density dataFile = catalogDir+haloFileBase+fileNums[0] inFile = open(dataFile, 'r') numPart = 0 for line in inFile: line = line.split(',') - if numHaloMass == "none" or float(line[6]) > minHaloMass: - numPart += 1 + if minHaloMass == "none" or float(line[6]) > minHaloMass: + numPart += 1 inFile.close() minRadius = 2*int(np.ceil(lbox/numPart**(1./3.))) if dataFormat == "multidark": - setName = prefix+"halos_min"str(minHaloMass) + 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, - numSlices, True, lbox, minRadius, omegaM) + redshifts, + numSubvolumes, numSlices, True, lbox, minRadius, omegaM) if args.halos or args.all: print " Doing halos" @@ -355,9 +357,11 @@ if args.halos or args.all: dataFile = catalogDir+haloFileBase+fileNums[iRedshift] inFile = open(dataFile, 'r') + numPart = 0 for line in inFile: line = line.split(',') - if numHaloMass == "none" or float(line[6]) > minHaloMass: + if minHaloMass == "none" or float(line[6]) > minHaloMass: + numPart += 1 inFile.close() sampleName = "md.halos_z"+redshift @@ -370,18 +374,16 @@ if args.halos or args.all: outFile.write("%d\n" %(numPart)) inFile = open(dataFile, 'r') - numKept = 0 - for line in inFile: - if numHaloMass == "none" or float(line[6]) > minHaloMass: - numKept += 1 - line = line.split(',') + for (iHalo,line) in enumerate(inFile): + line = line.split(',') + if minHaloMass == "none" or float(line[6]) > minHaloMass: x = float(line[0]) y = float(line[1]) z = float(line[2]) vz = float(line[5]) # write to output file - outFile.write("%d %e %e %e %e\n" %(numKept,x,y,z,vz)) + outFile.write("%d %e %e %e %e\n" %(iHalo,x,y,z,vz)) inFile.close() outFile.close() diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 280f568..0d24ca7 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -83,7 +83,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, if os.access("mask_index.txt", os.F_OK): os.system("mv %s %s" % ("mask_index.txt", zobovDir)) os.system("mv %s %s" % ("total_particles.txt", zobovDir)) - os.system("mv %s %s" % ("sample_info.txt", zobovDir)) + #os.system("mv %s %s" % ("sample_info.txt", zobovDir)) if os.access("galaxies.txt", os.F_OK): os.system("mv %s %s" % ("galaxies.txt", zobovDir)) @@ -160,7 +160,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, if os.access("comoving_distance.txt", os.F_OK): os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) - os.system("mv %s %s" % ("sample_info.txt", zobovDir)) + #os.system("mv %s %s" % ("sample_info.txt", zobovDir)) if os.access(parmFile, os.F_OK): os.unlink(parmFile) From 9d975224728fddbfbefd895948eea9adc56a68d5 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Sun, 18 Nov 2012 18:09:18 -0600 Subject: [PATCH 32/32] more bug fixes --- c_tools/stacking/pruneVoids.cpp | 49 ++++++++++--------- .../void_python_tools/backend/launchers.py | 5 +- 2 files changed, 28 insertions(+), 26 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index b50de5b..88d077f 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -272,7 +272,7 @@ int main(int argc, char **argv) { for (iVoid = 0; iVoid < numVoids; iVoid++) { voidID = voids[iVoid].voidID; - printf(" DOING %d (of %d) %d %d %f\n", iVoid, numVoids, voidID, + printf(" DOING %d (of %d) %d %d %f\n", iVoid, numVoids, voidID, voids[iVoid].numPart, voids[iVoid].radius); @@ -349,22 +349,22 @@ int main(int argc, char **argv) { voids[iVoid].centralDen = centralDen / (4./3. * M_PI * pow(centralRad, 3./2.)); // compute maximum extent - //if (args_info.isObservation_flag) { - // maxDist = 0.; - // for (p = 0; p < voids[iVoid].numPart; p++) { - // for (p2 = p; p2 < voids[iVoid].numPart; p2++) { + if (args_info.isObservation_flag) { + maxDist = 0.; + for (p = 0; p < voids[iVoid].numPart; p++) { + for (p2 = p; p2 < voids[iVoid].numPart; p2++) { -// dist[0] = voidPart[p].x - voidPart[p2].x; -// dist[1] = voidPart[p].y - voidPart[p2].y; -// dist[2] = voidPart[p].z - voidPart[p2].z; + dist[0] = voidPart[p].x - voidPart[p2].x; + dist[1] = voidPart[p].y - voidPart[p2].y; + dist[2] = voidPart[p].z - voidPart[p2].z; -// dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); -// if (dist2 > maxDist) maxDist = dist2; -// } -// } -// voids[iVoid].maxRadius = sqrt(maxDist)/2.; -// } else { - maxDist = 0.; + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + if (dist2 > maxDist) maxDist = dist2; + } + } + voids[iVoid].maxRadius = sqrt(maxDist)/2.; + } else { + maxDist = 0.; for (p = 0; p < voids[iVoid].numPart; p++) { dist[0] = voidPart[p].x - voids[iVoid].barycenter[0]; @@ -379,7 +379,7 @@ int main(int argc, char **argv) { if (dist2 > maxDist) maxDist = dist2; } voids[iVoid].maxRadius = sqrt(maxDist); -// } + } if (args_info.isObservation_flag) { // compute distance from center to nearest mock @@ -447,12 +447,14 @@ int main(int argc, char **argv) { } for (iVoid = 0; iVoid < numVoids; iVoid++) { -// TEST - //if (voids[iVoid].densCon > 1.5) { - // voids[iVoid].accepted = 0; - //} -// END TEST - + if (voids[iVoid].densCon < 1.5) { + //voids[iVoid].accepted = 0; + } + + // toss out voids that are obviously wrong + if (voids[iVoid].densCon > 1.e3) { + voids[iVoid].accepted = 0; + } if (strcmp(args_info.dataPortion_arg, "edge") == 0 && tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) { @@ -471,7 +473,6 @@ int main(int argc, char **argv) { if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) { voids[iVoid].accepted = -1; } - } numKept = 0; @@ -489,7 +490,7 @@ int main(int argc, char **argv) { fpSkyPositions = fopen(args_info.outSkyPositions_arg, "w"); fprintf(fp, "%d particles, %d voids.\n", mockIndex, numKept); fprintf(fp, "see column in master void file\n"); - fprintf(fpInfo, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID\n"); + fprintf(fpInfo, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast\n"); fprintf(fpSkyPositions, "# RA, dec, redshift, radius (Mpc/h), void ID\n"); for (iVoid = 0; iVoid < numVoids; iVoid++) { diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 0d24ca7..be2c12b 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -76,6 +76,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, if os.access("contour_map.fits", os.F_OK): os.system("mv %s %s" % ("contour_map.fits", zobovDir)) + os.system("mv %s %s" % ("mask_map.fits", zobovDir)) if os.access("comoving_distance.txt", os.F_OK): os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) @@ -378,8 +379,8 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, %s ranSeed %d dataPortion %s - barycenters %s - boundaryDistances %s + #barycenters %s + #boundaryDistances %s %s """ % \ (zobovDir+"/voidDesc_"+thisDataPortion+"_"+sampleName+".out",