guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 23:31:12 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'master' of ssh://bitbucket.org/cosmicvoids/void_identification

Browse source
This commit is contained in:
Guilhem Lavaux 2013-03-31 17:01:47 -04:00
commit 5a0f2383b3
26 changed files with 1232 additions and 306 deletions
Download patch file Download diff file Expand all files Collapse all files

51
analysis/datasetsToAnalyze.py Normal file
View file

@ -0,0 +1,51 @@
#+
# VIDE -- Void IDEntification pipeline -- ./crossCompare/analysis/datasetsToAnalyze.py
# Copyright (C) 2010-2013 Guilhem Lavaux
# Copyright (C) 2011-2013 P. M. Sutter
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; version 2 of the License.
#
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+
#!/usr/bin/env python
outputDir = "/home/psutter2/workspace/Voids/analysis/xcor/"
# Sim parameters
Ni = 1000237 # Number of dark matter particles per dimension in simulation
ss = 0.1 # Subsampling fraction of dark matter particles to read
Mpart = 8.721e9 # Particle mass [M_sol]
Npart = int(ss*Ni) # Particle number of subsample
Lboxcut = 0. # Size of optional margin to be cut from the box [h^(-1)Mpc]
Nmesh = 256 # Interpolation meshlength
Nbin = 70 # Number of bins for power spectrum and correlation function
r_H = 3000. # Hubble scale [h^(-1)Mpc]
ns = 0.95 # Spectral index
sigma_8 = 0.82 # Sigma_8
h = 0.7 # Dimensionless Hubble parameter
# Input files
matterDir = '/home/psutter2/workspace/Voids/catalogs/mergertree512/'
haloDir = '/home/psutter2/workspace/Voids/catalogs/mergertree512/'
matterFilename = 'mf_4s_1G_512_1.000'
haloFilename = 'mf_4s_1G_512_bgc2_1.000.sdf'
voidBaseDir = "/home/psutter2/workspace/Voids/"
sampleDirList = [
"mergertree512/mt_ss0.01/sample_mt_ss0.01_z0.00_d00/",
]
dataPortion = "central"

352
analysis/xcor.py Executable file
View file

@ -0,0 +1,352 @@
#!/usr/bin/env python
#+
# VIDE -- Void IDEntification pipeline -- ./crossCompare/analysis/mergerTree.py
# Copyright (C) 2010-2013 Guilhem Lavaux
# Copyright (C) 2011-2013 P. M. Sutter
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; version 2 of the License.
#
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+
from void_python_tools.backend import *
from void_python_tools.plotting import *
import void_python_tools.xcor as xcorlib
import imp
import pickle
import argparse
import os
import string
import numpy as np
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib import rc
from matplotlib.ticker import NullFormatter
import random
import sys
# ------------------------------------------------------------------------------
dataNameBase = "xcor"
parser = argparse.ArgumentParser(description='Analyze.')
parser.add_argument('--parmFile', dest='parmFile', default='datasetsToAnalyze.py',
help='path to parameter file')
args = parser.parse_args()
# ------------------------------------------------------------------------------
filename = args.parmFile
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(outputDir, os.F_OK):
os.makedirs(outputDir)
for (iSample, sampleDir) in enumerate(sampleDirList):
with open(voidBaseDir+sampleDir+"/sample_info.dat", 'rb') as input:
sample = pickle.load(input)
print " Working with", sample.fullName, "...",
sys.stdout.flush()
sampleName = sample.fullName
# Sim parameters
Lbox = sample.boxLen # Boxlength [h^(-1)Mpc]
Lbox -= 2*Lboxcut # Reduced boxlength [h^(-1)Mpc]
Om = sample.omegaM # Omega_m
Ol = 1.-Om # Omega_l
z = sample.zRange[0] # Redshift
a = 1./(1.+z) # Scale factor
rho_m = Mpart*(Ni/Lbox)**3 # Background density [(h/Mpc)^3]
# Input files
voidDir = voidBaseDir+'/'+sampleDir+'/'
voidFilename1 = 'centers_central_'+sample.fullName+'.out'
voidFilename2 = 'voidDesc_central_'+sample.fullName+'.out'
# Read files
matter_file = open(matterDir+matterFilename,'r')
matter_data = []
count = 0
sep = 1e8
for (i,line) in enumerate(matter_file):
if i/int(sep) == i/sep: print str(round(i*100./Ni**3,1))+' percent of all particles read'
if random.random() > ss: continue
count += 1
matter_data.append(line.split(',')[0:3])
print str(count)+' particles read in total'
matter_data = np.asarray(matter_data,dtype=np.float32)
matter_file.close()
halo_file = open(haloDir+haloFilename,'r')
halo_data = np.reshape(halo_file.read().replace('\n',',').split(',')[0:-1],(-1,12)).astype(np.float32)
halo_file.close()
void_file = open(voidDir+voidFilename1,'r')
void_header1 = void_file.readline().split(",")
void_data1 = np.reshape(void_file.read().split(),(-1,len(void_header))).astype(np.float32)
void_file.close()
void_file = open(voidDir+voidFilename2,'r')
void_header2 = void_file.readline().split(",")+void_file.readline().split(",")
void_data2 = np.reshape(void_file.read().split(),(-1,11)).astype(np.float32)
void_file.close()
# Define arrays
xm = matter_data[:,0:3]
xh = halo_data[:,0:3]
mh = halo_data[:,6]
xv = void_data1[:,0:3]
rv = void_data1[:,4]
vv = void_data1[:,6]
mv = void_data2[:,8]*Mpart
# Interpolate to mesh
dm, wm, ws = xcorlib.cic(xm, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh)
dh, wm, ws = xcorlib.cic(xh, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh)
dv, wm, ws = xcorlib.cic(xv, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh)
# Load dark matter grid
#output = open('dm_'+str(Nmesh)+'_ss'+str(ss)+'_z'+str(z)+'.dat', 'rb')
#dm = pickle.load(output)
#output.close()
# Save dark matter grid
#output = open('dm_'+str(Nmesh)+'_ss'+str(ss)+'_z'+str(z)+'.dat', 'wb')
#pickle.dump(dm,output)
#output.close()
# Power spectra & correlation functions
((Nm, km, Pmm, SPmm),(Nmx, rm, Xmm, SXmm)) = xcorlib.powcor(dm, dm, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Pvm, SPvm),(Nmx, rm, Xvm, SXvm)) = xcorlib.powcor(dv, dm, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Phm, SPhm),(Nmx, rm, Xhm, SXhm)) = xcorlib.powcor(dh, dm, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Pvv, SPvv),(Nmx, rm, Xvv, SXvv)) = xcorlib.powcor(dv, dv, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Pvh, SPvh),(Nmx, rm, Xvh, SXvh)) = xcorlib.powcor(dv, dh, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Phh, SPhh),(Nmx, rm, Xhh, SXhh)) = xcorlib.powcor(dh, dh, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
# Number densities
nm = np.empty(len(km))
nh = np.empty(len(km))
nv = np.empty(len(km))
nm[:] = Npart/Lbox**3
nh[:] = len(xh)/Lbox**3
nv[:] = len(xv)/Lbox**3
# Number functions
Mbin = 40
Vbin = 40
Nh, Mh = np.histogram(mh, bins = mh.min()*(mh.max()/mh.min())**(np.arange(Mbin+1)/float(Mbin)))
Nvm, Mv = np.histogram(mv, bins = mv.min()*(mv.max()/mv.min())**(np.arange(Vbin+1)/float(Vbin)))
Nvv, Vv = np.histogram(vv, bins = vv.min()*(vv.max()/vv.min())**(np.arange(Vbin+1)/float(Vbin)))
# Bias
b_hh = np.sqrt(Phh/Pmm)
b_vv = np.sqrt(Pvv/Pmm)
b_hm = Phm/Pmm
b_vm = Pvm/Pmm
b_vh = Pvh/Phh
knl = 0.04 # Wavenumber above which nonlinearities kick in [h/Mpc]
idls = np.where(km <= knl)[0]
bm_hm = np.average(b_hm[idls],weights=Nm[idls])
bm_vm = np.average(b_vm[idls],weights=Nm[idls])
bm_vh = np.average(b_vh[idls],weights=Nm[idls])
# Shot Noise
sn_hh = Phh - Phm**2/Pmm
sn_vh = Pvh - Pvm*Phm/Pmm
sn_vv = Pvv - Pvm**2/Pmm
# Plots
ms = 4
mew = 0.2
margin = 1.2
kmin = km.min()/margin
kmax = km.max()*margin
rmin = rm.min()/margin
rmax = rm.max()*margin
plt.imshow(np.sum(dm+1,2)/Nmesh,extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian')
plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$')
plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$')
plt.title(r'Dark matter')
plt.colorbar()
plt.savefig(outputDir+'/dm_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
plt.imshow(np.sum(dv+1,2)/Nmesh,extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian')
plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$')
plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$')
plt.title(r'Voids')
plt.colorbar()
plt.savefig(outputDir+'/dv_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
plt.imshow(np.sum(dh+1,2)/Nmesh,extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian')
plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$')
plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$')
plt.title(r'Halos')
plt.colorbar()
plt.savefig(outputDir+'/dh_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
pa, = plt.plot(Mh[:-1], Nh, 'ro-', ms=ms, mew=mew)
pb, = plt.plot(Vv[:-1]*1e9, Nvv, 'bo-', ms=ms, mew=mew)
plt.xlabel(r'$M \;[h^{-1}M_{\odot}]$ , $V \;[h^{-3}\mathrm{kpc}^3]$')
plt.ylabel(r'$N(M,V)$')
plt.title(r'Number of halos and voids')
plt.xscale('log')
plt.yscale('log')
plt.xlim(min(10**np.floor(np.log10(Vv.min()))*1e9,10**np.floor(np.log10(Mh.min()))), max(10**np.ceil(np.log10(Mh.max())),10**np.ceil(np.log10(Vv.max()))*1e8))
plt.ylim(10**np.floor(np.log10(Nh.min())), 10**np.ceil(np.log10(Nh.max())))
plt.annotate(r'$\frac{4\pi}{3}\langle r_\mathrm{v}\rangle^3$', xy=(4./3.*np.pi*rv.mean()**3*1e9,1.1), xytext=(-50,235),textcoords='offset points',arrowprops=dict(fc='k',arrowstyle="->",connectionstyle="angle,angleA=0,angleB=90,rad=10"))
plt.legend([pa,pb],['halos','voids'],'best' )
plt.savefig(outputDir+'/number_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
plt.subplot(211)
plt.subplots_adjust(wspace=0,hspace=0)
plt.plot(np.sort(vv), mv[np.argsort(vv)], 'r-', ms=ms, mew=mew, lw=0.01)
plt.plot(np.sort(vv), np.sort(mv), 'k-', ms=ms, mew=mew)
plt.plot(np.sort(vv), np.sort(vv)*1e9, 'k--', ms=ms, mew=mew)
plt.title(r'Mass-volume relation of voids')
plt.ylabel(r'$M \;[h^{-1}M_\odot]$')
plt.xscale('log')
plt.yscale('log')
plt.subplot(211).xaxis.set_major_formatter(NullFormatter())
plt.subplot(212)
plt.plot(np.sort(vv), mv[np.argsort(vv)]/np.sort(vv)/rho_m-1., 'b-', ms=ms, mew=mew, lw=0.01)
plt.plot(np.sort(vv), np.sort(mv)/np.sort(vv)/rho_m-1., 'k-', ms=ms, mew=mew)
plt.xlabel(r'$V \;[h^{-3}\mathrm{Mpc}^3]$')
plt.ylabel(r'$\delta$')
plt.xscale('log')
plt.yscale('linear')
plt.ylim(-1.01,-0.861)
plt.savefig(outputDir+'/massvol_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
pa ,= plt.plot(km, Phh, 'r-', ms=ms, mew=mew)
plt.plot(km, Phh-sn_hh, 'r:', ms=ms, mew=mew)
plt.fill_between(km, Phh+SPhh, abs(Phh-SPhh), color='r', alpha=0.2)
pb ,= plt.plot(km, Phm, 'y-', ms=ms, mew=mew)
plt.fill_between(km, Phm+SPhm, abs(Phm-SPhm), color='y', alpha=0.2)
pc ,= plt.plot(km, Pmm, 'k-', ms=ms, mew=mew)
plt.plot(km, Pmm-1./nm, 'k:', ms=ms, mew=mew)
plt.fill_between(km, Pmm+SPmm, abs(Pmm-SPmm), color='k', alpha=0.2)
pd ,= plt.plot(km, Pvh, 'g-', ms=ms, mew=mew)
plt.plot(km, -Pvh, 'g--', ms=ms, mew=mew)
plt.plot(km, abs(Pvh-sn_vh), 'g:', ms=ms, mew=mew)
plt.fill_between(km, abs(Pvh+SPvh), abs(Pvh-SPvh), color='g', alpha=0.2)
pe ,= plt.plot(km, Pvm, 'm-', ms=ms, mew=mew)
plt.plot(km, -Pvm, 'm--', ms=ms, mew=mew)
plt.fill_between(km, abs(Pvm+SPvm), abs(Pvm-SPvm), color='m', alpha=0.2)
pf ,= plt.plot(km, Pvv, 'b-', ms=ms, mew=mew)
plt.plot(km, Pvv-sn_vv, 'b:', ms=ms, mew=mew)
plt.fill_between(km, Pvv+SPvv, abs(Pvv-SPvv), color='b', alpha=0.2)
plt.annotate(r'$\frac{\pi}{\langle r_\mathrm{v}\rangle}$', xy=(np.pi/(rv.mean()),1.01*10**np.floor(np.log10(abs(Pvh).min()))/margin), xytext=(10,280),textcoords='offset points',arrowprops=dict(fc='k',arrowstyle="->",connectionstyle="angle,angleA=0,angleB=90,rad=10"))
plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$')
plt.ylabel(r'$P(k) \;[h^{-3}\mathrm{Mpc}^3]$')
plt.title(r'Power spectra')
plt.xscale('log')
plt.yscale('log')
plt.xlim(kmin,kmax)
plt.ylim(10**np.floor(np.log10(abs(Pvh).min()))/margin, max(10**np.ceil(np.log10(Phh.max())),10**np.ceil(np.log10(Pvv.max())))*margin)
plt.legend([pa, pb, pc, pd, pe, pf],['hh', 'hm', 'mm', 'vh', 'vm', 'vv'],'lower left' )
plt.savefig(outputDir+'/power_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
pa ,= plt.plot(rm, Xhh, 'r-', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xhh+SXhh), abs(Xhh-SXhh), color='r', alpha=0.2)
pb ,= plt.plot(rm, Xhm, 'y-', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xhm+SXhm), abs(Xhm-SXhm), color='y', alpha=0.2)
pc ,= plt.plot(rm, Xmm, 'k-', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xmm+SXmm), abs(Xmm-SXmm), color='k', alpha=0.2)
pd ,= plt.plot(rm, Xvh, 'g-', ms=ms, mew=mew)
plt.plot(rm, -Xvh, 'g--', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xvh+SXvh), abs(Xvh-SXvh), color='g', alpha=0.2)
pe ,= plt.plot(rm, Xvm, 'm-', ms=ms, mew=mew)
plt.plot(rm, -Xvm, 'm--', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xvm+SXvm), abs(Xvm-SXvm), color='m', alpha=0.2)
pf ,= plt.plot(rm, Xvv, 'b-', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xvv+SXvv), abs(Xvv-SXvv), color='b', alpha=0.2)
plt.annotate(r'$\langle r_\mathrm{v}\rangle$', xy=(rv.mean(),1.01*10**np.floor(np.log10(abs(Xvh).min()))/margin), xytext=(10,300),textcoords='offset points',arrowprops=dict(fc='k',arrowstyle="->",connectionstyle="angle,angleA=0,angleB=90,rad=10"))
plt.xlabel(r'$r \;[h^{-1}\mathrm{Mpc}]$')
plt.ylabel(r'$\xi(r)$')
plt.title(r'Correlation functions')
plt.xscale('log')
plt.yscale('log')
plt.xlim(rmin,rmax)
plt.ylim(10**np.floor(np.log10(abs(Xvh).min()))/margin, max(10**np.ceil(np.log10(Xhh.max())),10**np.ceil(np.log10(Xvv.max())))*margin)
plt.legend([pa, pb, pc, pd, pe, pf],['hh', 'hm', 'mm', 'vh', 'vm', 'vv'],'lower left' )
plt.savefig(outputDir+'/correlation_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
pa, = plt.plot(km, b_hh, 'r-', ms=ms, mew=mew)
pb, = plt.plot(km, b_hm, 'r--', ms=ms, mew=mew)
pc, = plt.plot(km, b_vv, 'b-', ms=ms, mew=mew)
pd, = plt.plot(km, b_vm, 'b--', ms=ms, mew=mew)
pe, = plt.plot(km, b_vh/bm_vh, 'g-', ms=ms, mew=mew)
plt.plot(km, np.sin(km*rv.mean())/(km*rv.mean()), 'k:', ms=ms, mew=mew)
plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$')
plt.ylabel(r'$b(k)$')
plt.title(r'Bias')
plt.xscale('log')
plt.yscale('linear')
plt.xlim(kmin,kmax)
plt.ylim(np.floor(b_vm.min()),np.ceil(max(b_hh.max(),b_vv.max())))
plt.legend([pa,pb,pc,pd,pe],['hh', 'hm', 'vv', 'vm', r'$\bar{u}_\mathrm{v}(k)$'],'best' )
plt.savefig(outputDir+'/bias_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
pa, = plt.plot(km, sn_hh, 'r-', ms=ms, mew=mew)
pb, = plt.plot(km, sn_vh, 'g-', ms=ms, mew=mew)
pc, = plt.plot(km, sn_vv, 'b-', ms=ms, mew=mew)
plt.plot(km, abs(sn_vh), 'g:')
pd, = plt.plot(km, 1/nh, 'r--')
pe, = plt.plot(km, 1/nv, 'b-.')
plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$')
plt.ylabel(r'$\sigma^2(k)$')
plt.title(r'Shotnoise')
plt.xscale('log')
plt.yscale('log')
plt.xlim(kmin,kmax)
plt.ylim(10**np.floor(np.log10(abs(sn_vh).min())), 10**np.ceil(np.log10(sn_vv.max())))
plt.legend([pa,pb,pc,pd,pe],['hh', 'vh', 'vv', r'$\bar{n}_\mathrm{h}^{-1}$', r'$\bar{n}_\mathrm{v}^{-1}$'],'best' )
plt.savefig(outputDir+'/shotnoise_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()

7
c_tools/analysis/voidOverlap.cpp
View file

@ -271,7 +271,7 @@ int main(int argc, char **argv) {
filename = string(args.outfile_arg); filename = string(args.outfile_arg);
filename = filename.append("summary.out"); filename = filename.append("summary.out");
fp = fopen(filename.c_str(), "w"); fp = fopen(filename.c_str(), "w");
fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, common volume ratio 2, relative dist, num matches, num significant matches\n"); fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, common volume ratio 2, relative dist, num matches, num significant matches, match ID\n");
for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
int voidID = catalog1.voids[iVoid1].voidID; int voidID = catalog1.voids[iVoid1].voidID;
if (catalog1.voids[iVoid1].numMatches > 0) { if (catalog1.voids[iVoid1].numMatches > 0) {
@ -286,14 +286,15 @@ int main(int argc, char **argv) {
rdist = catalog1.voids[iVoid1].matches[0].dist; rdist = catalog1.voids[iVoid1].matches[0].dist;
rdist /= catalog1.voids[iVoid1].radius; rdist /= catalog1.voids[iVoid1].radius;
fprintf(fp, "%d %.4f %.4f %.4f %.4f %.4f %d %d\n", voidID, fprintf(fp, "%d %.4f %.4f %.4f %.4f %.4f %d %d %d\n", voidID,
catalog1.voids[iVoid1].radiusMpc, catalog1.voids[iVoid1].radiusMpc,
rRatio, rRatio,
commonVolRatio, commonVolRatio,
volRatio, volRatio,
rdist, rdist,
catalog1.voids[iVoid1].numMatches, catalog1.voids[iVoid1].numMatches,
catalog1.voids[iVoid1].numBigMatches); catalog1.voids[iVoid1].numBigMatches,
catalog2.voids[iVoid2].voidID);
} else { } else {
fprintf(fp, "%d %.2f 0.0 0.0 0.0 0.0 0 0\n", voidID, fprintf(fp, "%d %.2f 0.0 0.0 0.0 0.0 0 0\n", voidID,

21
c_tools/mock/generateMock.cpp
View file

@ -260,7 +260,7 @@ void selectBox(SimuData *simu, std::vector<long>& targets, generateMock_info& ar
acceptance = acceptance =
acceptance && acceptance &&
(simu->Pos[j][i] > ranges[j][0]) && (simu->Pos[j][i] > ranges[j][0]) &&
(simu->Pos[j][i] < ranges[j][1]); (simu->Pos[j][i] <= ranges[j][1]);
p.Pos[j] = simu->Pos[j][i]; p.Pos[j] = simu->Pos[j][i];
p.Vel[j] = (simu->Vel[j] != 0) ? simu->Vel[j][i] : 0; p.Vel[j] = (simu->Vel[j] != 0) ? simu->Vel[j][i] : 0;
} }
@ -274,7 +274,7 @@ void selectBox(SimuData *simu, std::vector<long>& targets, generateMock_info& ar
numAccepted++; numAccepted++;
} }
} }
cout << "Num accepted here = " << numAccepted << " / input = " << simu->NumPart << endl; cout << "SELECTBOX: Num accepted here = " << numAccepted << " / input = " << simu->NumPart << " (after resubsampling)" << endl;
} }
class PreselectParticles: public SimulationPreprocessor class PreselectParticles: public SimulationPreprocessor
@ -385,7 +385,7 @@ void buildBox(SimuData *simu, long num_targets, long loaded,
} }
} }
void saveBox(SimuData *&boxed, const std::string& outbox) void saveBox(SimuData *&boxed, const std::string& outbox, generateMock_info& args_info)
{ {
double *ranges = boxed->as<double>("ranges"); double *ranges = boxed->as<double>("ranges");
NcFile f(outbox.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4); NcFile f(outbox.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4);
@ -395,6 +395,11 @@ void saveBox(SimuData *&boxed, const std::string& outbox)
int num_snapshots = *boxed->as<int>("num_snapshots"); int num_snapshots = *boxed->as<int>("num_snapshots");
long *uniqueID = boxed->as<long>("uniqueID"); long *uniqueID = boxed->as<long>("uniqueID");
if (!f.is_valid())
{
cerr << "Could not create parameter file '" << outbox << "'. Aborting." << endl;
exit(1);
}
f.add_att("range_x_min", ranges[0]); f.add_att("range_x_min", ranges[0]);
f.add_att("range_x_max", ranges[1]); f.add_att("range_x_max", ranges[1]);
f.add_att("range_y_min", ranges[2]); f.add_att("range_y_min", ranges[2]);
@ -403,6 +408,7 @@ void saveBox(SimuData *&boxed, const std::string& outbox)
f.add_att("range_z_max", ranges[5]); f.add_att("range_z_max", ranges[5]);
f.add_att("mask_index", -1); f.add_att("mask_index", -1);
f.add_att("is_observation", 0); f.add_att("is_observation", 0);
f.add_att("data_subsampling", args_info.subsample_arg);
NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart); NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
NcDim *NumSnap_dim = f.add_dim("numsnap_dim", num_snapshots); NcDim *NumSnap_dim = f.add_dim("numsnap_dim", num_snapshots);
@ -446,6 +452,13 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a
boxed->time = simu->time; boxed->time = simu->time;
boxed->BoxSize = simu->BoxSize; boxed->BoxSize = simu->BoxSize;
NcAtt *d_sub = f.get_att("data_subsampling");
if (d_sub == 0 || d_sub->as_double(0)/args_info.subsample_arg-1. > 1.e-5)
{
cerr << "Parameter file was not generated with the same simulation subsampling argument. Particles will be different. Stop here." << d_sub->as_double(0) << " " << args_info.subsample_arg <<endl;
exit(1);
}
NcVar *v_id = f.get_var("particle_ids"); NcVar *v_id = f.get_var("particle_ids");
NcVar *v_snap = f.get_var("snapshot_split"); NcVar *v_snap = f.get_var("snapshot_split");
long *edges1; long *edges1;
@ -710,7 +723,7 @@ int main(int argc, char **argv)
delete[] efac; delete[] efac;
} }
saveBox(simuOut, args_info.outputParameter_arg); saveBox(simuOut, args_info.outputParameter_arg, args_info);
generateOutput(simuOut, args_info.axis_arg, generateOutput(simuOut, args_info.axis_arg,
args_info.output_arg); args_info.output_arg);
delete preselector; delete preselector;

29
c_tools/mock/loaders/multidark_loader.cpp
View file

@ -20,13 +20,16 @@
#include <cassert> #include <cassert>
#include <boost/format.hpp>
#include <string> #include <string>
#include <fstream> #include <fstream>
#include <iostream> #include <iostream>
#include <CosmoTool/loadRamses.hpp> #include <CosmoTool/loadRamses.hpp>
#include <CosmoTool/fortran.hpp> #include <CosmoTool/fortran.hpp>
#include "simulation_loader.hpp" #include "simulation_loader.hpp"
#include <libsdf/cosmo.h>
using boost::format;
using namespace std; using namespace std;
using namespace CosmoTool; using namespace CosmoTool;
@ -82,6 +85,12 @@ public:
double tempData; double tempData;
double shift = 0.5*simu->BoxSize;
double rescale_position = simu->Hubble*1.e-5*100./simu->time;
double rescale_velocity = one_kpc/one_Gyr;
#define INFINITY std::numeric_limits<float>::max()
float min_pos[3] = {INFINITY,INFINITY, INFINITY}, max_pos[3] = {-INFINITY,-INFINITY,-INFINITY};
cout << "loading multidark particles" << endl; cout << "loading multidark particles" << endl;
long actualNumPart = 0; long actualNumPart = 0;
@ -96,6 +105,17 @@ public:
p.Pos[2] == -99 && p.Vel[2] == -99) p.Pos[2] == -99 && p.Vel[2] == -99)
break; break;
//p.Pos[0] = p.Pos[0]*rescale_position + shift;
//p.Pos[1] = p.Pos[1]*rescale_position + shift;
//p.Pos[2] = p.Pos[2]*rescale_position + shift;
//p.Vel[2] = p.Vel[2]*rescale_velocity;
// enforce box size in case of roundoff error
for (int k = 0; k < 3; k++) {
if (p.Pos[k] < 0) p.Pos[k] += simu->BoxSize;
if (p.Pos[k] >= simu->BoxSize) p.Pos[k] -= simu->BoxSize;
}
if (preproc != 0 && !preproc->accept(p)) if (preproc != 0 && !preproc->accept(p))
continue; continue;
@ -110,7 +130,14 @@ public:
reallocArray(uniqueID, allocated, actualNumPart); reallocArray(uniqueID, allocated, actualNumPart);
reallocArray(index, allocated, actualNumPart); reallocArray(index, allocated, actualNumPart);
} }
}
for (int k = 0; k < 3; k++) {
min_pos[k] = std::min(min_pos[k], p.Pos[k]);
max_pos[k] = std::max(max_pos[k], p.Pos[k]);
}
}
for (int k = 0; k < 3; k++) cout << boost::format("min[%d] = %g, max[%d] = %g") % k % min_pos[k] % k %max_pos[k] << endl;
applyTransformations(simu); applyTransformations(simu);
simu->NumPart = actualNumPart; simu->NumPart = actualNumPart;
simu->TotalNumPart = actualNumPart; simu->TotalNumPart = actualNumPart;

2
c_tools/mock/loaders/sdf_loader.cpp
View file

@ -190,7 +190,7 @@ public:
if (load_flags & (NEED_POSITION | NEED_VELOCITY)) if (load_flags & (NEED_POSITION | NEED_VELOCITY))
rescaleParticles(d, d->Hubble*1e-5, one_kpc/one_Gyr); rescaleParticles(d, d->Hubble*1e-5, one_kpc/one_Gyr);
// enforceBoxSize(d); enforceBoxSize(d);
applyTransformations(d); applyTransformations(d);
basicPreprocessing(d, preproc); basicPreprocessing(d, preproc);

234
c_tools/stacking/pruneVoids.cpp
View file

@ -40,6 +40,9 @@
#include <netcdfcpp.h> #include <netcdfcpp.h>
#include "pruneVoids_conf.h" #include "pruneVoids_conf.h"
#include <vector> #include <vector>
#include "assert.h"
#include "voidTree.hpp"
#include "loadZobov.hpp"
#define LIGHT_SPEED 299792.458 #define LIGHT_SPEED 299792.458
#define MPC2Z 100./LIGHT_SPEED #define MPC2Z 100./LIGHT_SPEED
@ -48,6 +51,8 @@
#define CENTRAL_VOID 1 #define CENTRAL_VOID 1
#define EDGE_VOID 2 #define EDGE_VOID 2
using namespace std;
typedef struct partStruct { typedef struct partStruct {
float x, y, z, vol; float x, y, z, vol;
} PART; } PART;
@ -70,12 +75,26 @@ typedef struct voidStruct {
float center[3], barycenter[3]; float center[3], barycenter[3];
int accepted; int accepted;
int voidType; int voidType;
int parentID, numChildren;
bool isLeaf, hasHighCentralDen;
gsl_vector *eval; gsl_vector *eval;
gsl_matrix *evec; gsl_matrix *evec;
} VOID; } VOID;
void openFiles(string outputDir, string sampleName, string name,
int mockIndex, int numKept,
FILE** fpZobov, FILE** fpCenters,
FILE** fpCentersNoCut,
FILE** fpBarycenter, FILE** fpDistances, FILE** fpShapes,
FILE** fpSkyPositions);
void closeFiles(FILE* fpZobov, FILE* fpCenters,
FILE* fpCentersNoCut,
FILE* fpBarycenter, FILE* fpDistances, FILE* fpShapes,
FILE* fpSkyPositions);
void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters, void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters,
FILE* fpCenterNoCut, FILE* fpCentersNoCut,
FILE* fpSkyPositions, FILE* fpBarycenters, FILE* fpDistances, FILE* fpSkyPositions, FILE* fpBarycenters, FILE* fpDistances,
FILE* fpShapes, bool isObservation, double *boxLen); FILE* fpShapes, bool isObservation, double *boxLen);
@ -108,15 +127,12 @@ int main(int argc, char **argv) {
int i, p, p2, numPartTot, numZonesTot, dummy, iVoid, iZ; int i, p, p2, numPartTot, numZonesTot, dummy, iVoid, iZ;
int numVoids, mockIndex, numKept; int numVoids, mockIndex, numKept;
double tolerance; double tolerance;
FILE *fp, *fpZobovCentral, *fpZobovAll, *fpCentersCentral, *fpCentersAll, FILE *fp, *fpZobov, *fpCenters, *fpCentersNoCut, *fpBarycenter,
*fpCentersNoCutCentral, *fpCentersNoCutAll, *fpBarycenterCentral, *fpDistances, *fpShapes, *fpSkyPositions;
*fpBarycenterAll, *fpDistancesCentral, *fpDistancesAll,
*fpShapesCentral, *fpShapesAll, *fpSkyPositionsCentral,
*fpSkyPositionsAll;
PART *part, *voidPart; PART *part, *voidPart;
ZONE2PART *zones2Parts; ZONE2PART *zones2Parts;
VOID2ZONE *void2Zones; VOID2ZONE *void2Zones;
std::vector<VOID> voids; vector<VOID> voids;
float *temp, junk, voidVol; float *temp, junk, voidVol;
int junkInt, voidID, numPart, numZones, zoneID, partID, maxNumPart; int junkInt, voidID, numPart, numZones, zoneID, partID, maxNumPart;
int coreParticle, zoneNumPart; int coreParticle, zoneNumPart;
@ -124,6 +140,7 @@ int main(int argc, char **argv) {
float centralRad, centralDen; float centralRad, centralDen;
double nearestEdge, redshift; double nearestEdge, redshift;
char line[500], junkStr[10]; char line[500], junkStr[10];
string outputDir, sampleName, name;
int mask_index; int mask_index;
double ranges[3][2], boxLen[3], mul; double ranges[3][2], boxLen[3], mul;
double volNorm, radius; double volNorm, radius;
@ -250,6 +267,9 @@ int main(int argc, char **argv) {
voids[i-1].radius = pow(voidVol/volNorm*3./4./M_PI, 1./3.); voids[i-1].radius = pow(voidVol/volNorm*3./4./M_PI, 1./3.);
voids[i-1].accepted = 1; voids[i-1].accepted = 1;
voids[i-1].isLeaf = true;
voids[i-1].hasHighCentralDen = false;
voids[i-1].eval = gsl_vector_alloc(3); voids[i-1].eval = gsl_vector_alloc(3);
voids[i-1].evec = gsl_matrix_alloc(3,3); voids[i-1].evec = gsl_matrix_alloc(3,3);
} }
@ -306,6 +326,27 @@ int main(int argc, char **argv) {
fclose(fp); fclose(fp);
free(temp); free(temp);
// load voids *again* using Guilhem's code so we can get tree
printf(" Re-loading voids and building tree..\n");
ZobovRep zobovCat;
if (!loadZobov(args.voidDesc_arg, args.zone2Part_arg,
args.void2Zone_arg,
0, zobovCat)) {
printf("Error loading catalog!\n");
return -1;
}
VoidTree *tree;
tree = new VoidTree(zobovCat);
zobovCat.allZones.erase(zobovCat.allZones.begin(), zobovCat.allZones.end());
//zobovCat.allZones.erase(zobovCat.allVoids.begin(), zobovCat.allVoids.end());
// copy tree information to our own data structures
for (iVoid = 0; iVoid < numVoids; iVoid++) {
voidID = voids[iVoid].voidID;
voids[iVoid].parentID = tree->getParent(voidID);
voids[iVoid].numChildren = tree->getChildren(voidID).size();
}
// check boundaries // check boundaries
printf(" Computing void properties...\n"); printf(" Computing void properties...\n");
@ -541,6 +582,7 @@ int main(int argc, char **argv) {
int numCentral = 0; int numCentral = 0;
int numEdge = 0; int numEdge = 0;
int numNearZ = 0; int numNearZ = 0;
int numAreParents = 0;
int numTooSmall = 0; int numTooSmall = 0;
printf(" Picking winners and losers...\n"); printf(" Picking winners and losers...\n");
@ -607,12 +649,30 @@ int main(int argc, char **argv) {
} }
} }
voids.resize(iGood); voids.resize(iGood);
printf(" 4th filter: rejected %d too close to low redshift boundaries\n", numNearZ); printf(" 4th filter: rejected %d below redshift boundaries\n", numNearZ);
numAreParents = 0;
iGood = 0;
for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].numChildren > 0) {
numAreParents++;
voids[iVoid].isLeaf = false;
} else {
//voids[iGood++] = voids[iVoid];
voids[iVoid].isLeaf = true;
}
}
//voids.resize(iGood);
//printf(" 5th filter: rejected %d that are not leaf nodes\n", numAreParents);
for (iVoid = 0; iVoid < voids.size(); iVoid++) { for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].centralDen > args.maxCentralDen_arg) { if (voids[iVoid].centralDen > args.maxCentralDen_arg) {
voids[iVoid].accepted = -1; voids[iVoid].accepted = -1;
voids[iVoid].hasHighCentralDen = true;
numHighDen++; numHighDen++;
} else {
voids[iVoid].hasHighCentralDen = false;
} }
} }
@ -630,80 +690,46 @@ int main(int argc, char **argv) {
printf(" We have %d edge voids\n", numEdge); printf(" We have %d edge voids\n", numEdge);
printf(" We have %d central voids\n", numCentral); printf(" We have %d central voids\n", numCentral);
printf(" We have %d too high central density\n", numHighDen); printf(" We have %d too high central density\n", numHighDen);
printf(" We have %d that are not leaf nodes\n", numAreParents);
printf(" Output...\n"); printf(" Output...\n");
fpZobovCentral = fopen((std::string(args.outputDir_arg)+"/voidDesc_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpZobovCentral, "%d particles, %d voids.\n", mockIndex, numKept);
fprintf(fpZobovCentral, "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
fpZobovAll = fopen((std::string(args.outputDir_arg)+"/voidDesc_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w"); outputDir = string(args.outputDir_arg);
fprintf(fpZobovAll, "%d particles, %d voids.\n", mockIndex, numKept); sampleName = (string(args.sampleName_arg)+".out");
fprintf(fpZobovAll, "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
fpBarycenterCentral = fopen((std::string(args.outputDir_arg)+"/barycenters_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fpBarycenterAll = fopen((std::string(args.outputDir_arg)+"/barycenters_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fpCentersCentral = fopen((std::string(args.outputDir_arg)+"/centers_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpCentersCentral, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part\n");
fpCentersAll = fopen((std::string(args.outputDir_arg)+"/centers_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpCentersAll, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part\n");
fpCentersNoCutCentral = fopen((std::string(args.outputDir_arg)+"/centers_nocut_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpCentersNoCutCentral, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part\n");
fpCentersNoCutAll = fopen((std::string(args.outputDir_arg)+"/centers_nocut_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpCentersNoCutAll, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part\n");
fpDistancesCentral = fopen((std::string(args.outputDir_arg)+"boundaryDistancesCentral_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fpDistancesAll = fopen((std::string(args.outputDir_arg)+"boundaryDistancesAll_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fpSkyPositionsCentral = fopen((std::string(args.outputDir_arg)+"/sky_positions_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpSkyPositionsCentral, "# RA, dec, redshift, radius (Mpc/h), void ID\n");
fpSkyPositionsAll = fopen((std::string(args.outputDir_arg)+"/sky_positions_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpSkyPositionsAll, "# RA, dec, redshift, radius (Mpc/h), void ID\n");
fpShapesCentral = fopen((std::string(args.outputDir_arg)+"/shapes_central_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpShapesCentral, "# void ID, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n");
fpShapesAll = fopen((std::string(args.outputDir_arg)+"/shapes_all_"+std::string(args.sampleName_arg)+".out").c_str(), "w");
fprintf(fpShapesAll, "# void ID, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n");
name = "central";
openFiles(outputDir, sampleName, name,
mockIndex, numKept,
&fpZobov, &fpCenters, &fpCentersNoCut, &fpBarycenter,
&fpDistances, &fpShapes, &fpSkyPositions);
for (iVoid = 0; iVoid < voids.size(); iVoid++) { for (iVoid = 0; iVoid < voids.size(); iVoid++) {
if (voids[iVoid].voidType == CENTRAL_VOID) { if (voids[iVoid].voidType == CENTRAL_VOID) {
outputVoid(iVoid, voids[iVoid], fpZobovCentral, fpCentersCentral, outputVoid(iVoid, voids[iVoid], fpZobov, fpCenters,
fpCentersNoCutCentral, fpSkyPositionsCentral, fpCentersNoCut, fpSkyPositions,
fpBarycenterCentral, fpDistancesCentral, fpShapesCentral, fpBarycenter, fpDistances, fpShapes,
args.isObservation_flag, boxLen);
}
if (voids[iVoid].voidType == EDGE_VOID ||
voids[iVoid].voidType == CENTRAL_VOID) {
outputVoid(iVoid, voids[iVoid], fpZobovAll, fpCentersAll,
fpCentersNoCutAll, fpSkyPositionsAll,
fpBarycenterAll, fpDistancesAll, fpShapesAll,
args.isObservation_flag, boxLen); args.isObservation_flag, boxLen);
} }
} }
closeFiles(fpZobov, fpCenters, fpCentersNoCut, fpBarycenter,
fpDistances, fpShapes, fpSkyPositions);
fclose(fpZobovCentral); name = "all";
fclose(fpZobovAll); openFiles(outputDir, sampleName, name,
fclose(fpCentersCentral); mockIndex, numKept,
fclose(fpCentersAll); &fpZobov, &fpCenters, &fpCentersNoCut, &fpBarycenter,
fclose(fpCentersNoCutCentral); &fpDistances, &fpShapes, &fpSkyPositions);
fclose(fpCentersNoCutAll); for (iVoid = 0; iVoid < voids.size(); iVoid++) {
fclose(fpBarycenterCentral); if (voids[iVoid].voidType == EDGE_VOID ||
fclose(fpBarycenterAll); voids[iVoid].voidType == CENTRAL_VOID) {
fclose(fpDistancesCentral); outputVoid(iVoid, voids[iVoid], fpZobov, fpCenters,
fclose(fpDistancesAll); fpCentersNoCut, fpSkyPositions,
fclose(fpShapesCentral); fpBarycenter, fpDistances, fpShapes,
fclose(fpShapesAll); args.isObservation_flag, boxLen);
fclose(fpSkyPositionsCentral); }
fclose(fpSkyPositionsAll); }
closeFiles(fpZobov, fpCenters, fpCentersNoCut, fpBarycenter,
fpDistances, fpShapes, fpSkyPositions);
clock2 = clock(); clock2 = clock();
printf(" Time: %f sec (for %d voids)\n", printf(" Time: %f sec (for %d voids)\n",
@ -716,9 +742,57 @@ int main(int argc, char **argv) {
} // end main } // end main
// ----------------------------------------------------------------------------
void openFiles(string outputDir, string sampleName, string name,
int mockIndex, int numKept,
FILE** fpZobov, FILE** fpCenters,
FILE** fpCentersNoCut,
FILE** fpBarycenter, FILE** fpDistances, FILE** fpShapes,
FILE** fpSkyPositions) {
*fpZobov = fopen((outputDir+"/voidDesc_"+name+"_"+sampleName).c_str(), "w");
fprintf(*fpZobov, "%d particles, %d voids.\n", mockIndex, numKept);
fprintf(*fpZobov, "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
*fpBarycenter = fopen((outputDir+"/barycenters_"+name+"_"+sampleName).c_str(), "w");
*fpCenters = fopen((outputDir+"/centers_"+name+"_"+sampleName).c_str(), "w");
fprintf(*fpCenters, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part, parent ID, number of children\n");
*fpCentersNoCut = fopen((outputDir+"/centers_nocut_"+name+"_"+sampleName).c_str(), "w");
fprintf(*fpCentersNoCut, "# center x,y,z (Mpc/h), volume (normalized), radius (Mpc/h), redshift, volume (Mpc/h^3), void ID, density contrast, num part, parent ID, number of children\n");
*fpDistances = fopen((outputDir+"boundaryDistances_"+name+"_"+sampleName).c_str(), "w");
*fpSkyPositions = fopen((outputDir+"/sky_positions_"+name+"_"+sampleName).c_str(), "w");
fprintf(*fpSkyPositions, "# RA, dec, redshift, radius (Mpc/h), void ID\n");
*fpShapes = fopen((outputDir+"/shapes_"+name+"_"+sampleName).c_str(), "w");
fprintf(*fpShapes, "# void ID, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n");
} // end openFiles
// ----------------------------------------------------------------------------
void closeFiles(FILE* fpZobov, FILE* fpCenters,
FILE* fpCentersNoCut,
FILE* fpBarycenter, FILE* fpDistances, FILE* fpShapes,
FILE* fpSkyPositions) {
fclose(fpZobov);
fclose(fpCenters);
fclose(fpCentersNoCut);
fclose(fpBarycenter);
fclose(fpDistances);
fclose(fpShapes);
fclose(fpSkyPositions);
} // end closeFile
// ---------------------------------------------------------------------------- // ----------------------------------------------------------------------------
void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters, void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters,
FILE* fpCenterNoCut, FILE* fpSkyPositions, FILE* fpCentersNoCut, FILE* fpSkyPositions,
FILE* fpBarycenters, FILE* fpDistances, FILE* fpShapes, FILE* fpBarycenters, FILE* fpDistances, FILE* fpShapes,
bool isObservation, double *boxLen) { bool isObservation, double *boxLen) {
@ -756,8 +830,8 @@ void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters,
outCenter[2] = (outVoid.barycenter[2]-boxLen[2]/2.)*100.; outCenter[2] = (outVoid.barycenter[2]-boxLen[2]/2.)*100.;
} }
if (outVoid.accepted == 1) { if (!outVoid.hasHighCentralDen && outVoid.isLeaf) {
fprintf(fpCenters, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d\n", fprintf(fpCenters, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d %d %d\n",
outCenter[0], outCenter[0],
outCenter[1], outCenter[1],
outCenter[2], outCenter[2],
@ -767,10 +841,12 @@ void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters,
4./3.*M_PI*pow(outVoid.radius, 3), 4./3.*M_PI*pow(outVoid.radius, 3),
outVoid.voidID, outVoid.voidID,
outVoid.densCon, outVoid.densCon,
outVoid.numPart); outVoid.numPart,
outVoid.parentID,
outVoid.numChildren);
} }
fprintf(fpCenterNoCut, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d\n", fprintf(fpCentersNoCut, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d %d %d\n",
outCenter[0], outCenter[0],
outCenter[1], outCenter[1],
outCenter[2], outCenter[2],
@ -780,7 +856,9 @@ void outputVoid(int iVoid, VOID outVoid, FILE* fpZobov, FILE* fpCenters,
4./3.*M_PI*pow(outVoid.radius, 3), 4./3.*M_PI*pow(outVoid.radius, 3),
outVoid.voidID, outVoid.voidID,
outVoid.densCon, outVoid.densCon,
outVoid.numPart); outVoid.numPart,
outVoid.parentID,
outVoid.numChildren);
fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n", fprintf(fpSkyPositions, "%.2f %.2f %.5f %.2f %d\n",
atan((outVoid.barycenter[1]-boxLen[1]/2.) / atan((outVoid.barycenter[1]-boxLen[1]/2.) /

121
crossCompare/analysis/compareProfiles.py Executable file
View file

@ -0,0 +1,121 @@
#!/usr/bin/env python
#+
# VIDE -- Void IDEntification pipeline -- ./pipeline/apAnalysis.py
# Copyright (C) 2010-2013 Guilhem Lavaux
# Copyright (C) 2011-2013 P. M. Sutter
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; version 2 of the License.
#
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+
# takes a set of voids with given radii, and stacks and plots matches in
# other catalogs
from void_python_tools.backend import *
import imp
import pickle
import os
import numpy as np
import argparse
# ------------------------------------------------------------------------------
mergerNameBase = "mergerTree"
parser = argparse.ArgumentParser(description='Analyze.')
parser.add_argument('--parm', dest='parm', default='datasetsToAnalyze.py', help='path to parameter file')
args = parser.parse_args()
# ------------------------------------------------------------------------------
filename = args.parm
print " Loading parameters from", filename
if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename
exit(-1)
parms = imp.load_source("name", filename)
globals().update(vars(parms))
if not os.access(dataDir, os.F_OK):
os.makedirs(dataDir)
mergerFileBase = dataDir + "/" + mergerNameBase
baseIDList = []
# get list of base voids
with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input:
baseSample = pickle.load(input)
baseSampleName = baseSample.fullName
baseVoidList = np.loadtxt(workDir+baseSampleDir+"/centers_nocut_central_"+\
baseSampleName+".out")
for stack in baseSample.stacks:
print " Stack:", stack.rMin
accepted = (baseVoidList[:,4] > stack.rMin) & (baseVoidList[:,4] < stack.rMax)
baseIDList = baseVoidList[accepted][:,7]
for (iSample, sampleDir) in enumerate(sampleDirList):
with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input:
sample = pickle.load(input)
print " Working with", sample.fullName, "...",
sys.stdout.flush()
sampleName = sample.fullName
# get list of appropriate voids
if sample == baseSample:
idList = baseIDList
else:
idList = 2
matchList = np.loadtxt(mergerFileBase+"_"+baseSampleName+"_"+sampleName+\
"_summary.out")
accepted = (matchList[:,0] == baseIDList).any()
idList = matchList[accepted][:,8]
voidBaseDir = workDir+"/"+sampleDir+"stacks"
runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin,
stack.rMax, thisDataPortion,
customLine="selected")
stack.rMin = 0.
stack.rMax = 1000.
voidDir = voidBaseDir+"_"+runSuffix
if not os.access(voidDir,os.F_OK): os.makedirs(voidDir)
print " Stacking voids...",
sys.stdout.flush()
STACK_PATH = CTOOLS_PATH+"/stacking/stackVoidsZero"
launchStack(sample, stack, STACK_PATH,
thisDataPortion=thisDataPortion,
logDir=logDir, voidDir=voidDir,
zobovDir=workDir+"/"+sampleDir,
freshStack=freshStack, INCOHERENT=False,
ranSeed=101010, summaryFile=None,
continueRun=continueRun,
dataType=sample.dataType,
idList=idList)
print " Profiling stacks...",
sys.stdout.flush()
logFile = logDir+"/profile_"+sampleName+"_"+runSuffix+".out"
launchProfile(sample, stack, voidDir=voidDir,
logFile=logFile, continueRun=False)
print " Done!"

12
crossCompare/analysis/mergerTree.py Normal file → Executable file
View file

@ -1,3 +1,4 @@
#!/usr/bin/env python
#+ #+
# VIDE -- Void IDEntification pipeline -- ./crossCompare/analysis/mergerTree.py # VIDE -- Void IDEntification pipeline -- ./crossCompare/analysis/mergerTree.py
# Copyright (C) 2010-2013 Guilhem Lavaux # Copyright (C) 2010-2013 Guilhem Lavaux
@ -17,9 +18,6 @@
# with this program; if not, write to the Free Software Foundation, Inc., # with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+ #+
#!/usr/bin/env python
# plots cumulative distributions of number counts
from void_python_tools.backend import * from void_python_tools.backend import *
from void_python_tools.plotting import * from void_python_tools.plotting import *
@ -35,13 +33,13 @@ import argparse
dataNameBase = "mergerTree" dataNameBase = "mergerTree"
parser = argparse.ArgumentParser(description='Analyze.') parser = argparse.ArgumentParser(description='Analyze.')
parser.add_argument('--parmFile', dest='parmFile', default='datasetsToAnalyze.py', parser.add_argument('--parm', dest='parm', default='datasetsToAnalyze.py',
help='path to parameter file') help='path to parameter file')
args = parser.parse_args() args = parser.parse_args()
# ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------
filename = args.parmFile filename = args.parm
print " Loading parameters from", filename print " Loading parameters from", filename
if not os.access(filename, os.F_OK): if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename print " Cannot find parameter file %s!" % filename
@ -78,7 +76,7 @@ for (iSample, sampleDir) in enumerate(sampleDirList):
thisDataPortion="central", logFile=logFile, thisDataPortion="central", logFile=logFile,
continueRun=False, workDir=workDir, continueRun=False, workDir=workDir,
outputFile=stepOutputFileName, outputFile=stepOutputFileName,
matchMethod="proximity") matchMethod="useID")
# attach columns to summary file # attach columns to summary file
#if iSample == 1: #if iSample == 1:
@ -101,3 +99,5 @@ for (iSample, sampleDir) in enumerate(sampleDirList):
#if os.access("mergerTree.log", os.F_OK): os.unlink("mergerTree.log") #if os.access("mergerTree.log", os.F_OK): os.unlink("mergerTree.log")
#if os.access("temp.out", os.F_OK): os.unlink("temp.out") #if os.access("temp.out", os.F_OK): os.unlink("temp.out")
#if os.access("thisStep.out", os.F_OK): os.unlink("thisStep.out") #if os.access("thisStep.out", os.F_OK): os.unlink("thisStep.out")
print " Done!"

47
crossCompare/plotting/plotNumberFunc.py Normal file → Executable file
View file

@ -1,3 +1,4 @@
#!/usr/bin/env python
#+ #+
# VIDE -- Void IDEntification pipeline -- ./crossCompare/plotting/plotNumberFunc.py # VIDE -- Void IDEntification pipeline -- ./crossCompare/plotting/plotNumberFunc.py
# Copyright (C) 2010-2013 Guilhem Lavaux # Copyright (C) 2010-2013 Guilhem Lavaux
@ -17,7 +18,6 @@
# with this program; if not, write to the Free Software Foundation, Inc., # with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+ #+
#!/usr/bin/env python
# plots cumulative distributions of number counts # plots cumulative distributions of number counts
@ -36,8 +36,8 @@ from scipy.stats import ks_2samp
plotNameBase = "compdist" plotNameBase = "compdist"
obsFudgeFactor = 1.0 # what fraction of the volume are we *reall* capturing? #obsFudgeFactor = 1.0 # what fraction of the volume are we *reall* capturing?
#obsFudgeFactor = .66 # what fraction of the volume are we *reall* capturing? obsFudgeFactor = .15 # what fraction of the volume are we *reall* capturing?
linewidth = 1 linewidth = 1
@ -45,7 +45,7 @@ parser = argparse.ArgumentParser(description='Plot.')
parser.add_argument('--show', dest='showPlot', action='store_const', parser.add_argument('--show', dest='showPlot', action='store_const',
const=True, default=False, const=True, default=False,
help='display the plot (default: just write eps)') help='display the plot (default: just write eps)')
parser.add_argument('--parmFile', dest='parmFile', default='datasetsToPlot.py', parser.add_argument('--parm', dest='parm', default='datasetsToPlot.py',
help='path to parameter file') help='path to parameter file')
args = parser.parse_args() args = parser.parse_args()
@ -55,7 +55,7 @@ errorBarsX = np.linspace(0, plotMax, num=nErrorBars)
# ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------
filename = args.parmFile filename = args.parm
print " Loading parameters from", filename print " Loading parameters from", filename
if not os.access(filename, os.F_OK): if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename print " Cannot find parameter file %s!" % filename
@ -102,7 +102,7 @@ for dataPortion in dataPortions:
boxVol *= 1.e-9 # Mpc->Gpc boxVol *= 1.e-9 # Mpc->Gpc
filename = workDir+"/"+sampleDirList[iSample]+"/centers_"+dataPortion+"_"+\ filename = workDir+"/"+sampleDirList[iSample]+"/centers_nocut_"+dataPortion+"_"+\
sampleName+".out" sampleName+".out"
if not os.access(filename, os.F_OK): if not os.access(filename, os.F_OK):
print "File not found: ", filename print "File not found: ", filename
@ -137,8 +137,13 @@ for dataPortion in dataPortions:
ecolor=colorList[iSample], ecolor=colorList[iSample],
fmt=None, label='_nolegend_', capsize=0) fmt=None, label='_nolegend_', capsize=0)
hist, bin_edges = np.histogram(data, bins=100, range=(0,100)) hist, bin_edges = np.histogram(data, bins=40, range=(0,100))
allData.append(hist) #allData.append(hist)
binCenters = 0.5*(bin_edges[1:] + bin_edges[:-1])
#plt.plot(binCenters, hist, '-',
# label=lineTitle, color=colorList[iSample],
# linewidth=linewidth)
plt.legend(title = "Samples", loc = "upper right", prop={'size':8}) plt.legend(title = "Samples", loc = "upper right", prop={'size':8})
#plt.title(plotTitle) #plt.title(plotTitle)
@ -155,19 +160,19 @@ plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight") plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight") plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
dataFile = figDir+"/data_"+plotName+".dat" #dataFile = figDir+"/data_"+plotName+".dat"
fp = open(dataFile, 'w') #fp = open(dataFile, 'w')
fp.write("# R [Mpc/h], N [h^3 Gpc^-3]\n") #fp.write("# R [Mpc/h], N [h^3 Gpc^-3]\n")
fp.write("# ") #fp.write("# ")
for sample in dataSampleList: #for sample in dataSampleList:
fp.write(sample.fullName+" ") # fp.write(sample.fullName+" ")
fp.write("\n") #fp.write("\n")
for i in xrange(100): #for i in xrange(100):
fp.write(str(bin_edges[i]) + " ") # fp.write(str(bin_edges[i]) + " ")
for iSample in xrange(len(dataSampleList)): # for iSample in xrange(len(dataSampleList)):
fp.write(str(allData[iSample][i])+" ") # fp.write(str(allData[iSample][i])+" ")
fp.write("\n") # fp.write("\n")
fp.close() #fp.close()
if args.showPlot: if args.showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png") os.system("display %s" % figDir+"/fig_"+plotName+".png")

1
external/external_build.cmake vendored
View file

@ -234,6 +234,7 @@ ExternalProject_Add(cosmotool
-DGSLCBLAS_LIBRARY=${GSLCBLAS_LIBRARY} -DGSLCBLAS_LIBRARY=${GSLCBLAS_LIBRARY}
-DNETCDF_LIBRARY=${NETCDF_LIBRARY} -DNETCDF_LIBRARY=${NETCDF_LIBRARY}
-DNETCDFCPP_LIBRARY=${NETCDFCPP_LIBRARY} -DNETCDFCPP_LIBRARY=${NETCDFCPP_LIBRARY}
-DENABLE_SHARP=OFF
) )
SET(COSMOTOOL_LIBRARY ${CMAKE_BINARY_DIR}/ext_build/cosmotool/lib/libCosmoTool.a) SET(COSMOTOOL_LIBRARY ${CMAKE_BINARY_DIR}/ext_build/cosmotool/lib/libCosmoTool.a)
set(COSMOTOOL_INCLUDE_PATH ${CMAKE_BINARY_DIR}/ext_build/cosmotool/include) set(COSMOTOOL_INCLUDE_PATH ${CMAKE_BINARY_DIR}/ext_build/cosmotool/include)

9
pipeline/generateCatalog.py
View file

@ -29,7 +29,7 @@ import pickle
# ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------
if (len(sys.argv) == 1): if (len(sys.argv) == 1):
print "Usage: ./analyzeVoids.py parameter_file.py" print "Usage: ./generateCatalog.py parameter_file.py"
exit(-1) exit(-1)
if (len(sys.argv) > 1): if (len(sys.argv) > 1):
@ -88,7 +88,10 @@ for sample in dataSampleList:
if (sample.dataType == "simulation"): if (sample.dataType == "simulation"):
fp.write("Particles placed on lightcone: %g\n" % sample.useLightCone) fp.write("Particles placed on lightcone: %g\n" % sample.useLightCone)
fp.write("Peculiar velocities included: %g\n" % sample.usePecVel) fp.write("Peculiar velocities included: %g\n" % sample.usePecVel)
fp.write("Additional subsampling fraction: %s\n" % sample.subsample[-1]) if (len(sample.subsample) == 1):
fp.write("Additional subsampling fraction: %s\n" % sample.subsample)
else:
fp.write("Additional subsampling fraction: %s\n" % sample.subsample[-1])
fp.write("Simulation box length (Mpc/h): %g\n" % sample.boxLen) fp.write("Simulation box length (Mpc/h): %g\n" % sample.boxLen)
fp.write("Simulation Omega_M: %g\n" % sample.omegaM) fp.write("Simulation Omega_M: %g\n" % sample.omegaM)
fp.write("Number of simulation subvolumes: %s\n" % sample.numSubvolumes) fp.write("Number of simulation subvolumes: %s\n" % sample.numSubvolumes)
@ -147,5 +150,7 @@ if (startCatalogStage <= 4) and (endCatalogStage >= 4):
dataPortion=thisDataPortion, setName=setName) dataPortion=thisDataPortion, setName=setName)
plotNumberDistribution(workDir, dataSampleList, figDir, showPlot=False, plotNumberDistribution(workDir, dataSampleList, figDir, showPlot=False,
dataPortion=thisDataPortion, setName=setName) dataPortion=thisDataPortion, setName=setName)
plotVoidDistribution(workDir, dataSampleList, figDir, showPlot=False,
dataPortion=thisDataPortion, setName=setName)
print "\n Done!" print "\n Done!"

154
python_tools/pipeline_source/prepareCatalogs.in.py
View file

@ -50,13 +50,13 @@ parser.add_argument('--hod', dest='hod', action='store_const',
parser.add_argument('--all', dest='all', action='store_const', parser.add_argument('--all', dest='all', action='store_const',
const=True, default=False, const=True, default=False,
help='write everything') help='write everything')
parser.add_argument('--parmFile', dest='parmFile', parser.add_argument('--parm', dest='parm',
default="", default="",
help='path to parameter file') help='path to parameter file')
args = parser.parse_args() args = parser.parse_args()
filename = args.parmFile filename = args.parm
print " Loading parameters from", filename print " Loading parameters from", filename
if not os.access(filename, os.F_OK): if not os.access(filename, os.F_OK):
print " Cannot find parameter file %s!" % filename print " Cannot find parameter file %s!" % filename
@ -76,6 +76,26 @@ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1):
return sampleName return sampleName
#------------------------------------------------------------------------------
def getNickName(sampleName):
splitName = sampleName.split('_')
if "ss" in splitName[1]:
nickName = "Subsample = " + splitName[1].replace("ss","")
nickName += ", z = " + splitName[2].replace("z","")
elif "hod" in splitName[1]:
nickName = "HOD = " + splitName[2]
nickName += ", z = " + splitName[3].replace("z","")
elif "halos" in splitName[1]:
nickName = "Halos, min mass = " + splitName[2].replace("min","")
nickName += ", z = " + splitName[3].replace("z","")
else:
nickName = sampleName
return nickName
#------------------------------------------------------------------------------ #------------------------------------------------------------------------------
# for given dataset parameters, outputs a script for use with analyzeVoids # for given dataset parameters, outputs a script for use with analyzeVoids
def writeScript(setName, dataFileNameBase, dataFormat, def writeScript(setName, dataFileNameBase, dataFormat,
@ -137,7 +157,7 @@ newSample = Sample(dataFile = "{dataFile}",
dataFormat = "{dataFormat}", dataFormat = "{dataFormat}",
dataUnit = {dataUnit}, dataUnit = {dataUnit},
fullName = "{sampleName}", fullName = "{sampleName}",
nickName = "{sampleName}", nickName = "{nickName}",
dataType = "simulation", dataType = "simulation",
zBoundary = ({zMin}, {zMax}), zBoundary = ({zMin}, {zMax}),
zRange = ({zMin}, {zMax}), zRange = ({zMin}, {zMax}),
@ -218,11 +238,14 @@ newSample.addStack(0.0, 5.0, 90, 95, False, False)
sampleName = getSampleName(setName, sliceMin, useVel, sampleName = getSampleName(setName, sliceMin, useVel,
iSlice=iSlice, iVol=mySubvolume) iSlice=iSlice, iVol=mySubvolume)
nickName = getNickName(sampleName)
scriptFile.write(sampleInfo.format(dataFile=dataFileName, scriptFile.write(sampleInfo.format(dataFile=dataFileName,
dataFormat=dataFormat, dataFormat=dataFormat,
dataUnit=dataUnit, dataUnit=dataUnit,
sampleName=sampleName, sampleName=sampleName,
nickName=nickName,
zMin=sliceMin, zMin=sliceMin,
zMax=sliceMax, zMax=sliceMax,
zMinMpc=sliceMinMpc, zMinMpc=sliceMinMpc,
@ -285,24 +308,39 @@ for iSubSample in xrange(len(subSamples)):
scriptDir, catalogDir, fileNums, redshifts, scriptDir, catalogDir, fileNums, redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM, numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
subsample=subSampleToUse, suffix=suffix) subsample=subSampleToUse, suffix=suffix)
writeScript(setName, fileToUse, dataFormat, if doPecVel:
scriptDir, catalogDir, fileNums, redshifts, writeScript(setName, fileToUse, dataFormat,
numSubvolumes, numSlices, True, lbox, minRadius, omegaM, scriptDir, catalogDir, fileNums, redshifts,
subsample=subSampleToUse, suffix=suffix) numSubvolumes, numSlices, True, lbox, minRadius, omegaM,
subsample=subSampleToUse, suffix=suffix)
else: else:
subSampleToUse = keepFractionList if doSubSampling:
suffix = "" # prepare scripts using our own format
fileToUse = partFileList[0] dataFormatToUse = "multidark"
writeScript(setName, fileToUse, dataFormat, subSampleToUse = 1.0
fileToUse = prefix+"ss"+str(thisSubSample)+"_z"
partFileList = []
suffix = ".dat"
for (iRedshift, redshift) in enumerate(redshifts):
sampleName = getSampleName(setName, redshift, False)
outFileName = sampleName+".dat"
partFileList.append(outFileName)
else:
dataFormatToUse = dataFormat
subSampleToUse = keepFractionList
suffix = ""
fileToUse = partFileList[0]
writeScript(setName, fileToUse, dataFormatToUse,
scriptDir, catalogDir, fileNums, redshifts, scriptDir, catalogDir, fileNums, redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM, numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
subsample=subSampleToUse, suffix=suffix, subsample=subSampleToUse, suffix=suffix,
dataFileNameList=partFileList) dataFileNameList=partFileList)
writeScript(setName, fileToUse, dataFormat, if doPecVel:
scriptDir, catalogDir, fileNums, redshifts, writeScript(setName, fileToUse, dataFormatToUse,
numSubvolumes, numSlices, True, lbox, minRadius, omegaM, scriptDir, catalogDir, fileNums, redshifts,
subsample=subSampleToUse, suffix=suffix, numSubvolumes, numSlices, True, lbox, minRadius, omegaM,
dataFileNameList=partFileList) subsample=subSampleToUse, suffix=suffix,
dataFileNameList=partFileList)
if args.subsample or args.all: if args.subsample or args.all:
@ -313,17 +351,46 @@ for iSubSample in xrange(len(subSamples)):
print " redshift", redshift print " redshift", redshift
sys.stdout.flush() sys.stdout.flush()
if dataFormat == "multidark": if dataFormat == "multidark" or dataFormat == "sdf":
# reuse previous subamples in order to: # reuse previous subamples in order to:
# - preserve unique IDs across multiple subsamples # - preserve unique IDs across multiple subsamples
# - reuse smaller files for faster processing # - reuse smaller files for faster processing
if prevSubSample == -1: if prevSubSample == -1:
dataFile = catalogDir+"/"+particleFileBase+fileNums[iRedshift] if particleFileDummy == '':
dataFile = catalogDir+"/"+particleFileBase+fileNums[iRedshift]
else:
dataFile = particleFileBase.replace(particleFileDummy,
fileNums[iRedshift])
dataFile = catalogDir+"/"+dataFile
keepFraction = float(thisSubSample) / baseResolution
else: else:
sampleName = prefix+"ss"+str(prevSubSample)+"_z"+redshift sampleName = prefix+"ss"+str(prevSubSample)+"_z"+redshift
dataFile = catalogDir+"/"+sampleName+".dat" dataFile = catalogDir+"/"+sampleName+".dat"
keepFraction = float(thisSubSample) / float(prevSubSample) keepFraction = float(thisSubSample) / float(prevSubSample)
if prevSubSample == -1 and dataFormat == "sdf":
convertedFile = dataFile + "_temp"
SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64"
scale = 1./(1.+float(redshift))
rescale_position = hubble/1000./scale
shift = lbox/2.
rescale_velocity = 3.08567802e16/3.1558149984e16
command = "%s %s x y z vz vy vx | awk '{print $1*%g+%g, $2*%g+%g, $3*%g+%g, $4*%g, $5*%g, $6*%g}' > %s" % (SDFcvt_PATH, dataFile,
rescale_position,
shift,
rescale_position,
shift,
rescale_position,
shift,
rescale_velocity,
rescale_velocity,
rescale_velocity,
convertedFile )
os.system(command)
os.system(command)
dataFile = convertedFile
inFile = open(dataFile, 'r') inFile = open(dataFile, 'r')
sampleName = prefix+"ss"+str(thisSubSample)+"_z"+redshift sampleName = prefix+"ss"+str(thisSubSample)+"_z"+redshift
@ -335,6 +402,11 @@ for iSubSample in xrange(len(subSamples)):
outFile.write("%s\n" %(redshift)) outFile.write("%s\n" %(redshift))
outFile.write("%d\n" %(maxKeep)) outFile.write("%d\n" %(maxKeep))
if dataFormat == "sdf":
splitter = ' '
if dataFormat == "multidark":
splitter = ','
numKept = 0 numKept = 0
for (i,line) in enumerate(inFile): for (i,line) in enumerate(inFile):
if (prevSubSample != -1 and i < 5): continue # skip header if (prevSubSample != -1 and i < 5): continue # skip header
@ -344,7 +416,7 @@ for iSubSample in xrange(len(subSamples)):
#if numKept > maxKeep: break #if numKept > maxKeep: break
if (prevSubSample == -1): if (prevSubSample == -1):
line = line.split(',') line = line.split(splitter)
x = float(line[0]) x = float(line[0])
y = float(line[1]) y = float(line[1])
z = float(line[2]) z = float(line[2])
@ -360,6 +432,9 @@ for iSubSample in xrange(len(subSamples)):
inFile.close() inFile.close()
outFile.close() outFile.close()
if prevSubSample == -1 and dataFormat == "sdf":
os.unlink(dataFile)
elif dataFormat == "random": elif dataFormat == "random":
sampleName = "ran.ss"+str(thisSubSample)+"_z"+redshift sampleName = "ran.ss"+str(thisSubSample)+"_z"+redshift
outFile = open(catalogDir+"/"+sampleName+".dat", 'w') outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
@ -409,14 +484,25 @@ if (args.script or args.all) and haloFileBase != "":
minRadies = 10 minRadies = 10
setName = prefix+"halos_min"+str(minHaloMass) setName = prefix+"halos_min"+str(minHaloMass)
fileList = []
for (iRedshift, redshift) in enumerate(redshifts):
sampleName = getSampleName(setName, redshift, False)
outFileName = sampleName+".dat"
fileList.append(outFileName)
writeScript(setName, prefix+"halos_min"+str(minHaloMass)+"_z", "multidark", writeScript(setName, prefix+"halos_min"+str(minHaloMass)+"_z", "multidark",
scriptDir, catalogDir, fileNums, scriptDir, catalogDir, fileNums,
redshifts, redshifts,
numSubvolumes, numSlices, False, lbox, minRadius, omegaM) numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
writeScript(setName, prefix+"halos_min"+str(minHaloMass)+"_z", "multidark", dataFileNameList = fileList)
scriptDir, catalogDir, fileNums,
redshifts, if doPecVel:
numSubvolumes, numSlices, True, lbox, minRadius, omegaM) writeScript(setName, prefix+"halos_min"+str(minHaloMass)+"_z",
"multidark",
scriptDir, catalogDir, fileNums,
redshifts,
numSubvolumes, numSlices, True, lbox, minRadius, omegaM,
dataFileNameList = fileList)
if (args.halos or args.all) and haloFileBase != "": if (args.halos or args.all) and haloFileBase != "":
print " Doing halos" print " Doing halos"
@ -451,7 +537,7 @@ if (args.halos or args.all) and haloFileBase != "":
numPart += 1 numPart += 1
inFile.close() inFile.close()
sampleName = prefix+"halos_min"+str(minHaloMass)+"_z"+fileNums[iRedshift] sampleName = prefix+"halos_min"+str(minHaloMass)+"_z"+redshifts[iRedshift]
outFileName = catalogDir+"/"+sampleName+".dat" outFileName = catalogDir+"/"+sampleName+".dat"
outFile = open(outFileName, 'w') outFile = open(outFileName, 'w')
outFile.write("%f\n" %(lbox)) outFile.write("%f\n" %(lbox))
@ -536,15 +622,25 @@ root_filename {workDir}/hod
if (args.script or args.all) and haloFileBase != "": if (args.script or args.all) and haloFileBase != "":
print " Doing HOD scripts" print " Doing HOD scripts"
sys.stdout.flush() sys.stdout.flush()
for thisHod in hodParmList: for thisHod in hodParmList:
fileList = []
for (iRedshift, redshift) in enumerate(redshifts):
sampleName = getSampleName(prefix+"hod_"+thisHod['name'], redshift, False)
outFileName = sampleName+".dat"
fileList.append(outFileName)
print " ", thisHod['name'] print " ", thisHod['name']
setName = prefix+"hod_"+thisHod['name'] setName = prefix+"hod_"+thisHod['name']
writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark", writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark",
scriptDir, catalogDir, fileNums, redshifts, scriptDir, catalogDir, fileNums, redshifts,
numSubvolumes, numSlices, False, lbox, 15, omegaM) numSubvolumes, numSlices, False, lbox, 15, omegaM,
writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark", dataFileNameList = fileList)
scriptDir, catalogDir, fileNums, redshifts, if doPecVel:
numSubvolumes, numSlices, True, lbox, 15, omegaM) writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark",
scriptDir, catalogDir, fileNums, redshifts,
numSubvolumes, numSlices, True, lbox, 15, omegaM,
dataFileNameList = fileList)
if (args.hod or args.all) and haloFileBase != "": if (args.hod or args.all) and haloFileBase != "":
print " Doing HOD" print " Doing HOD"

2
python_tools/setup.py
View file

@ -31,7 +31,7 @@ setup(
cmdclass = {'build_ext': build_ext}, cmdclass = {'build_ext': build_ext},
include_dirs = [np.get_include()], include_dirs = [np.get_include()],
packages= packages=
['void_python_tools','void_python_tools.backend','void_python_tools.apTools', ['void_python_tools','void_python_tools.backend','void_python_tools.apTools', 'void_python_tools.xcor',
'void_python_tools.apTools.profiles','void_python_tools.apTools.chi2', 'void_python_tools.plotting'], 'void_python_tools.apTools.profiles','void_python_tools.apTools.chi2', 'void_python_tools.plotting'],
#ext_modules = [Extension("void_python_tools.chi2.velocityProfileFitNative", ["void_python_tools/chi2/velocityProfileFitNative.pyx"], libraries=["gsl", "gslcblas"]), Extension("void_python_tools.chi2.likelihoo", ["void_python_tools/chi2/likelihood.pyx"], libraries=["gsl", "gslcblas"])] #ext_modules = [Extension("void_python_tools.chi2.velocityProfileFitNative", ["void_python_tools/chi2/velocityProfileFitNative.pyx"], libraries=["gsl", "gslcblas"]), Extension("void_python_tools.chi2.likelihoo", ["void_python_tools/chi2/likelihood.pyx"], libraries=["gsl", "gslcblas"])]
ext_modules = [ ext_modules = [

5
python_tools/void_python_tools/__init__.py
View file

@ -17,6 +17,7 @@
# with this program; if not, write to the Free Software Foundation, Inc., # with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+ #+
from void_python_tools.backend import * from void_python_tools.backend import *
from void_python_tools.apTools import * from void_python_tools.apTools import *
from void_python_tools.plotting import * from void_python_tools.plotting import *
from void_python_tools.xcor import *

4
python_tools/void_python_tools/backend/classes.py
View file

@ -184,9 +184,9 @@ def jobSuccessful(logFile, doneString):
if doneString in line: jobDone = True if doneString in line: jobDone = True
return jobDone return jobDone
def getStackSuffix(zMin, zMax, rMin, rMax, dataPortion): def getStackSuffix(zMin, zMax, rMin, rMax, dataPortion, customLine=""):
return "z"+str(zMin)+"-"+str(zMax)+"_"+str(rMin)+"-"+str(rMax)+\ return "z"+str(zMin)+"-"+str(zMax)+"_"+str(rMin)+"-"+str(rMax)+\
"Mpc"+"_"+dataPortion "Mpc"+customLine+"_"+dataPortion
def my_import(name): def my_import(name):
mod = __import__(name) mod = __import__(name)

55
python_tools/void_python_tools/backend/launchers.py
View file

@ -132,6 +132,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
return return
prevSubSample = -1 prevSubSample = -1
firstSubSample = -1
for thisSubSample in sample.subsample.split(', '): for thisSubSample in sample.subsample.split(', '):
if prevSubSample == -1: if prevSubSample == -1:
@ -141,13 +142,14 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
keepFraction = float(thisSubSample) keepFraction = float(thisSubSample)
subSampleLine = "subsample %g" % keepFraction subSampleLine = "subsample %g" % keepFraction
resubSampleLine = "" resubSampleLine = ""
firstSubSample = keepFraction
else: else:
inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+\ inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+\
sampleName+"_ss"+prevSubSample+".par" sampleName+"_ss"+prevSubSample+".par"
outputFile = zobovDir+"/zobov_slice_" + sampleName + "_ss" + \ outputFile = zobovDir+"/zobov_slice_" + sampleName + "_ss" + \
thisSubSample thisSubSample
keepFraction = float(thisSubSample)/float(prevSubSample) keepFraction = float(thisSubSample)/float(prevSubSample)
subSampleLine = "subsample %g" % keepFraction subSampleLine = "subsample %s" % firstSubSample
resubSampleLine = "resubsample %g" % keepFraction resubSampleLine = "resubsample %g" % keepFraction
includePecVelString = "" includePecVelString = ""
@ -162,6 +164,8 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
dataFileLine = "gadget " + datafile dataFileLine = "gadget " + datafile
elif sample.dataFormat == "sdf": elif sample.dataFormat == "sdf":
dataFileLine = "sdf " + datafile dataFileLine = "sdf " + datafile
else:
raise ValueError("unknown dataFormat '%s'" % sample.dataFormat)
iX = float(sample.mySubvolume[0]) iX = float(sample.mySubvolume[0])
iY = float(sample.mySubvolume[1]) iY = float(sample.mySubvolume[1])
@ -320,8 +324,8 @@ def launchPrune(sample, binPath,
periodicLine = " --periodic='" periodicLine = " --periodic='"
if sample.numSubvolumes == 1: periodicLine += "xy" if sample.numSubvolumes == 1: periodicLine += "xy"
if sample.zBoundaryMpc[1] - sample.zBoundaryMpc[0] - \ if np.abs(sample.zBoundaryMpc[1] - sample.zBoundaryMpc[0] - \
sample.boxLen <= 1.e-1: sample.boxLen) <= 1.e-1:
periodicLine += "z" periodicLine += "z"
periodicLine += "' " periodicLine += "' "
@ -376,7 +380,7 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir,
periodicLine = " --periodic='" periodicLine = " --periodic='"
if sample1.dataType != "observation": if sample1.dataType != "observation":
if sample1.numSubvolumes == 1: periodicLine += "xy" if sample1.numSubvolumes == 1: periodicLine += "xy"
if sample1.zBoundaryMpc[1] - sample1.zBoundaryMpc[0] - sample1.boxLen <= 1.e-1: if np.abs(sample1.zBoundaryMpc[1] - sample1.zBoundaryMpc[0] - sample1.boxLen) <= 1.e-1:
periodicLine += "z" periodicLine += "z"
periodicLine += "' " periodicLine += "' "
@ -416,6 +420,7 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir,
cmd += periodicLine cmd += periodicLine
cmd += " --outfile=" + outputFile cmd += " --outfile=" + outputFile
cmd += " &> " + logFile cmd += " &> " + logFile
open("temp.par",'w').write(cmd)
os.system(cmd) os.system(cmd)
if jobSuccessful(logFile, "Done!\n"): if jobSuccessful(logFile, "Done!\n"):
@ -433,12 +438,19 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
voidDir=None, freshStack=True, runSuffix=None, voidDir=None, freshStack=True, runSuffix=None,
zobovDir=None, zobovDir=None,
INCOHERENT=False, ranSeed=None, summaryFile=None, INCOHERENT=False, ranSeed=None, summaryFile=None,
continueRun=None, dataType=None, prefixRun=""): continueRun=None, dataType=None, prefixRun="",
idList=None, rescaleOverride=None):
sampleName = sample.fullName sampleName = sample.fullName
if idList != None:
customLine = "selected"
else:
customLine = ""
runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin,
stack.rMax, thisDataPortion) stack.rMax, thisDataPortion,
customLine=customLine)
logFile = logDir+("/%sstack_"%prefixRun)+sampleName+"_"+runSuffix+".out" logFile = logDir+("/%sstack_"%prefixRun)+sampleName+"_"+runSuffix+".out"
@ -474,11 +486,26 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
else: else:
nullTestFlag = "" nullTestFlag = ""
if stack.rescaleMode == "rmax": if rescaleOverride == None: rescaleOverride = stack.rescaleMode
if rescaleOverride == "rmax":
rescaleFlag = "rescale" rescaleFlag = "rescale"
else: else:
rescaleFlag = "" rescaleFlag = ""
idListFile = "idlist.temp"
if idList != None:
idListFlag = "idList " + idListFile
idFile = open(idListFile, 'w')
for id in idList: idFile.write(str(id)+"\n")
idFile.close()
rMinToUse = 0.
rMaxToUse = 100000.
centralRadius = rMaxToUse
else:
idListFlag = ""
rMinToUse = stack.rMin
rMaxToUse = stack.rMax
conf=""" conf="""
desc %s desc %s
partzone %s partzone %s
@ -501,13 +528,14 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
barycenters %s barycenters %s
boundaryDistances %s boundaryDistances %s
%s %s
%s
""" % \ """ % \
(zobovDir+"/voidDesc_"+thisDataPortion+"_"+sampleName+".out", (zobovDir+"/voidDesc_"+thisDataPortion+"_"+sampleName+".out",
zobovDir+"/voidPart_"+sampleName+".dat", zobovDir+"/voidPart_"+sampleName+".dat",
zobovDir+"/voidZone_"+sampleName+".dat", zobovDir+"/voidZone_"+sampleName+".dat",
zobovDir+"/vol_"+sampleName+".dat", zobovDir+"/vol_"+sampleName+".dat",
stack.rMin, rMinToUse,
stack.rMax, rMaxToUse,
zobovDir+("/%szobov_slice_"%prefixRun)+sampleName, zobovDir+("/%szobov_slice_"%prefixRun)+sampleName,
zobovDir+"/zobov_slice_"+sampleName+".par", zobovDir+"/zobov_slice_"+sampleName+".par",
maxDen, maxDen,
@ -522,7 +550,9 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
thisDataPortion, thisDataPortion,
zobovDir+"/barycenters_"+thisDataPortion+"_"+sampleName+".out", zobovDir+"/barycenters_"+thisDataPortion+"_"+sampleName+".out",
zobovDir+"/boundaryDistances_"+thisDataPortion+"_"+sampleName+".out", zobovDir+"/boundaryDistances_"+thisDataPortion+"_"+sampleName+".out",
rescaleFlag) rescaleFlag,
idListFlag
)
parmFile = os.getcwd()+("/%sstack_"%prefixRun)+sample.fullName+".par" parmFile = os.getcwd()+("/%sstack_"%prefixRun)+sample.fullName+".par"
@ -672,8 +702,9 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
os.system("mv %s %s" % ("normalizations.txt", voidDir+"/")) os.system("mv %s %s" % ("normalizations.txt", voidDir+"/"))
os.system("mv %s %s" % ("boundaryDistances.txt", voidDir+"/")) os.system("mv %s %s" % ("boundaryDistances.txt", voidDir+"/"))
if os.access(parmFile, os.F_OK): if os.access(idListFile, os.F_OK): os.unlink(idListFile)
os.unlink(parmFile)
if os.access(parmFile, os.F_OK): os.unlink(parmFile)
return return

50
python_tools/void_python_tools/plotting/plotTools.py
View file

@ -18,7 +18,7 @@
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+ #+
__all__=['plotRedshiftDistribution', 'plotSizeDistribution', 'plot1dProfiles', __all__=['plotRedshiftDistribution', 'plotSizeDistribution', 'plot1dProfiles',
'plotMarg1d', 'plotNumberDistribution'] 'plotMarg1d', 'plotNumberDistribution', 'plotVoidDistribution']
from void_python_tools.backend.classes import * from void_python_tools.backend.classes import *
from plotDefs import * from plotDefs import *
@ -26,6 +26,7 @@ import numpy as np
import os import os
import pylab as plt import pylab as plt
import void_python_tools.apTools as vp import void_python_tools.apTools as vp
import void_python_tools.xcor as xcor
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
def plotRedshiftDistribution(workDir=None, sampleList=None, figDir=None, def plotRedshiftDistribution(workDir=None, sampleList=None, figDir=None,
@ -266,7 +267,7 @@ def plotNumberDistribution(workDir=None, sampleList=None, figDir=None,
boxVol *= 1.e-9 boxVol *= 1.e-9
filename = workDir+"/sample_"+sampleName+"/centers_"+dataPortion+"_"+\ filename = workDir+"/sample_"+sampleName+"/centers_nocut_"+dataPortion+"_"+\
sampleName+".out" sampleName+".out"
if not os.access(filename, os.F_OK): if not os.access(filename, os.F_OK):
print "File not found: ", filename print "File not found: ", filename
@ -294,5 +295,50 @@ def plotNumberDistribution(workDir=None, sampleList=None, figDir=None,
if showPlot: if showPlot:
os.system("display %s" % figDir+"/fig_"+plotName+".png") os.system("display %s" % figDir+"/fig_"+plotName+".png")
# -----------------------------------------------------------------------------
def plotVoidDistribution(workDir=None, sampleList=None, figDir=None,
plotNameBase="dv",
showPlot=False, dataPortion=None, setName=None):
Nmesh = 256
for (iSample,sample) in enumerate(sampleList):
plt.clf()
plt.xlabel("x (Mpc/h)")
plt.ylabel("y (Mpc/h)")
sampleName = sample.fullName
plotName = plotNameBase+"_"+sampleName
filename = workDir+"/sample_"+sampleName+"/centers_nocut_"+dataPortion+"_"+\
sampleName+".out"
if not os.access(filename, os.F_OK):
print "File not found: ", filename
continue
void_file = open(filename,'r')
void_header1 = void_file.readline().split(",")
void_data1 = np.reshape(void_file.read().split(),
(-1,len(void_header1))).astype(np.float32)
void_file.close()
xv = void_data1[:,0:3]
rv = void_data1[:,4]
vv = void_data1[:,6]
dv, wm, ws = xcor.cic(xv, sample.boxLen, Lboxcut = 0., Nmesh = Nmesh)
plt.imshow(np.sum(dv+1,2)/Nmesh,extent=[0,sample.boxLen,0,sample.boxLen],
aspect='equal',
cmap='YlGnBu_r',interpolation='gaussian')
plt.colorbar()
plt.title("Voids: "+sampleName)
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")

20
python_tools/void_python_tools/xcor/__init__.py Normal file
View file

@ -0,0 +1,20 @@
#+
# VIDE -- Void IDEntification pipeline -- ./python_tools/void_python_tools/plotting/__init__.py
# Copyright (C) 2010-2013 Guilhem Lavaux
# Copyright (C) 2011-2013 P. M. Sutter
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; version 2 of the License.
#
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+
from xcorlib import *

73
python_tools/void_python_tools/xcor/xcorlib.py Normal file
View file

@ -0,0 +1,73 @@
import numpy as np
def cic( x, Lbox, Lboxcut = 0, Nmesh = 128, weights = None ):
if weights == None: weights = 1
wm = np.mean(weights)
ws = np.mean(weights**2)
d = np.mod(x/(Lbox+2*Lboxcut)*Nmesh,1)
box = ([Lboxcut,Lbox+Lboxcut],[Lboxcut,Lbox+Lboxcut],[Lboxcut,Lbox+Lboxcut])
rho = np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*(1-d[:,1])*(1-d[:,2]))[0] \
+ np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*(1-d[:,1])*(1-d[:,2]))[0],1,0) \
+ np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*d[:,1]*(1-d[:,2]))[0],1,1) \
+ np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*(1-d[:,1])*d[:,2])[0],1,2) \
+ np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*d[:,1]*(1-d[:,2]))[0],1,0),1,1) \
+ np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*(1-d[:,1])*d[:,2])[0],1,0),1,2) \
+ np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*d[:,1]*d[:,2])[0],1,1),1,2) \
+ np.roll(np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*d[:,1]*d[:,2])[0],1,0),1,1),1,2)
rho /= wm
rho = rho/rho.mean() - 1.
return (rho, wm, ws)
def powcor( d1, d2, Lbox, Nmesh = 128, Nbin = 100, scale = 'lin', cor = False ):
# Fourier transform
d1 = np.fft.fftn(d1)
d2 = np.fft.fftn(d2)
# CIC correction
wid = np.indices(np.shape(d1)) - Nmesh/2
#wid[np.where(wid >= Nmesh/2)] -= Nmesh
wid = wid*np.pi/Nmesh + 1e-100
wcic = np.prod(np.sin(wid)/wid,0)**2
# Shell average power spectrum
dk = 2*np.pi/Lbox
Pk = np.conj(d1)*d2*(Lbox/Nmesh**2)**3
Pk = np.fft.fftshift(Pk)/wcic**2
(Nm, km, Pkm, SPkm) = shellavg(np.real(Pk), dk, Nmesh, Nbin = Nbin, xmin = 0., xmax = Nmesh*dk/2, scale = scale)
# Inverse Fourier transform and shell average correlation function
if cor == True:
dx = Lbox/Nmesh
Xr = np.fft.ifftshift(np.fft.ifftn(np.fft.ifftshift(Pk)))*(Nmesh/Lbox)**3
(Nmx, rm, Xrm, SXrm) = shellavg(np.real(Xr), dx, Nmesh, Nbin = Nbin, xmin = dx, xmax = 140., scale = scale)
return ((Nm, km, Pkm, SPkm),(Nmx, rm, Xrm, SXrm))
else: return (Nm, km, Pkm, SPkm)
def shellavg( f, dx, Nmesh, Nbin = 100, xmin = 0., xmax = 1., scale = 'lin' ):
x = np.indices(np.shape(f)) - Nmesh/2
#x[np.where(x >= Nmesh/2)] -= Nmesh
x = dx*np.sqrt(np.sum(x**2,0))
if scale == 'lin': bins = xmin+(xmax-xmin)* (np.arange(Nbin+1)/float(Nbin))
if scale == 'log': bins = xmin*(xmax/xmin)**(np.arange(Nbin+1)/float(Nbin))
Nm = np.histogram(x, bins = bins)[0]
xm = np.histogram(x, bins = bins, weights = x)[0]/Nm
fm = np.histogram(x, bins = bins, weights = f)[0]/Nm
fs = np.sqrt((np.histogram(x, bins = bins, weights = f**2)[0]/Nm - fm**2)/(Nm-1))
return (Nm, xm, fm, fs)

75
zobov/jozov.c
View file

@ -85,7 +85,7 @@ int main(int argc,char **argv) {
printf("Bad density threshold.\n"); printf("Bad density threshold.\n");
exit(0); exit(0);
} }
if (sscanf(argv[7],"%f",&mockIndex) == 0) { if (sscanf(argv[7],"%d",&mockIndex) == 0) {
printf("Bad mock galaxy index.\n"); printf("Bad mock galaxy index.\n");
exit(0); exit(0);
} }
@ -101,6 +101,8 @@ int main(int argc,char **argv) {
exit(0); exit(0);
} }
fread(&np,1, sizeof(int),adj); fread(&np,1, sizeof(int),adj);
if (mockIndex < 0)
mockIndex = np;
printf("adj: %d particles\n", np); printf("adj: %d particles\n", np);
FF; FF;
@ -108,33 +110,40 @@ int main(int argc,char **argv) {
p = (PARTICLE *)malloc(np*sizeof(PARTICLE)); p = (PARTICLE *)malloc(np*sizeof(PARTICLE));
/* Adjacencies*/ /* Adjacencies*/
for (i=0;i<np;i++) { for (i=0;i<np;i++) {
fread(&p[i].nadj,1,sizeof(int),adj); fread(&p[i].nadj,1,sizeof(pid_t),adj);
/* The number of adjacencies per particle */ /* The number of adjacencies per particle */
if (p[i].nadj > 0) if (p[i].nadj > 0)
p[i].adj = (int *)malloc(p[i].nadj*sizeof(int)); p[i].adj = (pid_t *)malloc(p[i].nadj*sizeof(pid_t));
else p[i].adj = 0; else p[i].adj = 0;
p[i].ncnt = 0; /* Temporarily, it's an adj counter */ p[i].ncnt = 0; /* Temporarily, it's an adj counter */
} }
for (i=0;i<np;i++) { for (i=0;i<np;i++) {
fread(&nin,1,sizeof(int),adj); fread(&nin,1,sizeof(pid_t),adj);
if (nin > 0) if (nin > 0)
for (k=0;k<nin;k++) { for (k=0;k<nin;k++) {
fread(&j,1,sizeof(int),adj); fread(&j,1,sizeof(pid_t),adj);
if (j < np) { if (j < np) {
/* Set both halves of the pair */ /* Set both halves of the pair */
assert(i < j); assert(i < j);
if (p[i].ncnt == p[i].nadj) if (p[i].ncnt == p[i].nadj)
{ {
p[i].adj = (int *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(int)); int q;
p[i].nadj++; printf("OVERFLOW for particle %d (pending %d). List of accepted:\n", i, j);
for (q=0;q<p[i].nadj;q++)
printf(" %d\n", p[i].adj[q]);
abort();
} }
if (p[j].ncnt == p[j].nadj) if (p[j].ncnt == p[j].nadj)
{ {
p[j].adj = (int *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(int)); int q;
p[j].nadj++; printf("OVERFLOW for particle %d (pending %d). List of accepted:\n", j, i);
} for (q=0;q<p[j].nadj;q++)
p[i].adj[p[i].ncnt] = j; printf(" %d\n", p[j].adj[q]);
p[j].adj[p[j].ncnt] = i; abort();
}
p[i].adj[p[i].ncnt] = j;
p[j].adj[p[j].ncnt] = i;
p[i].ncnt++; p[j].ncnt++; p[i].ncnt++; p[j].ncnt++;
} else { } else {
printf("%d: adj = %d\n",i,j); printf("%d: adj = %d\n",i,j);
@ -144,13 +153,13 @@ int main(int argc,char **argv) {
fclose(adj); fclose(adj);
/* Check that we got all the pairs */ /* Check that we got all the pairs */
// adj = fopen(adjfile, "r"); /* adj = fopen(adjfile, "r");
// fread(&np,1, sizeof(int),adj); fread(&np,1, sizeof(int),adj);*/
for (i=0;i<np;i++) { for (i=0;i<np;i++) {
// fread(&nin,1,sizeof(int),adj); /* actually nadj */ /* fread(&nin,1,sizeof(int),adj); /* actually nadj */
// PMS // PMS
if (p[i].ncnt != p[i].nadj && i < mockIndex) { if (p[i].ncnt != p[i].nadj && i < mockIndex) {
//if (p[i].ncnt != p[i].nadj) { /*if (p[i].ncnt != p[i].nadj) {*/
// END PMS // END PMS
p[i].nadj = p[i].ncnt; p[i].nadj = p[i].ncnt;
printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj); printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
@ -185,9 +194,9 @@ int main(int argc,char **argv) {
} }
fclose(vol); fclose(vol);
jumped = (int *)malloc(np*sizeof(int)); jumped = (pid_t *)malloc(np*sizeof(pid_t));
jumper = (int *)malloc(np*sizeof(int)); jumper = (pid_t *)malloc(np*sizeof(pid_t));
numinh = (int *)malloc(np*sizeof(int)); numinh = (pid_t *)malloc(np*sizeof(pid_t));
/* find jumper */ /* find jumper */
for (i = 0; i < np; i++) { for (i = 0; i < np; i++) {
@ -252,7 +261,7 @@ int main(int argc,char **argv) {
} }
for (h=0;h<nzones;h++) { for (h=0;h<nzones;h++) {
zt[h].adj = (int *)malloc(zt[h].nadj*sizeof(int)); zt[h].adj = (pid_t *)malloc(zt[h].nadj*sizeof(pid_t));
if (zt[h].adj == NULL) { if (zt[h].adj == NULL) {
printf("Unable to allocate %d adj's of zone %d\n",zt[h].nadj,h); printf("Unable to allocate %d adj's of zone %d\n",zt[h].nadj,h);
exit(0); exit(0);
@ -313,7 +322,7 @@ int main(int argc,char **argv) {
for (h=0;h<nzones;h++) { for (h=0;h<nzones;h++) {
/*printf("%d ",zt[h].nadj);*/ /*printf("%d ",zt[h].nadj);*/
z[h].nadj = zt[h].nadj; z[h].nadj = zt[h].nadj;
z[h].adj = (int *)malloc(zt[h].nadj*sizeof(int)); z[h].adj = (pid_t *)malloc(zt[h].nadj*sizeof(pid_t));
z[h].slv = (float *)malloc(zt[h].nadj*sizeof(float)); z[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
for (za = 0; za<zt[h].nadj; za++) { for (za = 0; za<zt[h].nadj; za++) {
z[h].adj[za] = zt[h].adj[za]; z[h].adj[za] = zt[h].adj[za];
@ -326,12 +335,12 @@ int main(int argc,char **argv) {
free(zt); free(zt);
free(numinh); free(numinh);
m = (int **)malloc(nzones*sizeof(int *)); m = (pid_t **)malloc(nzones*sizeof(pid_t *));
/* Not in the zone struct since it'll be freed up (contiguously, we hope) /* Not in the zone struct since it'll be freed up (contiguously, we hope)
soon */ soon */
nm = (int *)malloc(nzones*sizeof(int)); nm = (pid_t *)malloc(nzones*sizeof(pid_t));
for (h=0; h<nzones; h++) { for (h=0; h<nzones; h++) {
m[h] = (int *)malloc(z[h].np*sizeof(int)); m[h] = (pid_t *)malloc(z[h].np*sizeof(pid_t));
nm[h] = 0; nm[h] = 0;
z[h].vol = 0.; z[h].vol = 0.;
} }
@ -354,14 +363,14 @@ int main(int argc,char **argv) {
printf("Problem opening zonefile %s.\n\n",zonfile); printf("Problem opening zonefile %s.\n\n",zonfile);
exit(0); exit(0);
} }
fwrite(&np,1,4,zon); fwrite(&np,1,sizeof(pid_t),zon);
fwrite(&nzones,1,4,zon); fwrite(&nzones,1,sizeof(int),zon);
for (h=0; h<nzones; h++) { for (h=0; h<nzones; h++) {
// PMS // PMS
//printf("%d %d %d", &(z[h].np), m[h], z[h].np); //printf("%d %d %d", &(z[h].np), m[h], z[h].np);
// END PMS // END PMS
fwrite(&(z[h].np),1,4,zon); fwrite(&(z[h].np),1,sizeof(pid_t),zon);
fwrite(m[h],z[h].np,4,zon); fwrite(m[h],z[h].np,sizeof(pid_t),zon);
free(m[h]); free(m[h]);
} }
free(m); free(m);
@ -394,7 +403,7 @@ int main(int argc,char **argv) {
printf("Problem opening zonefile %s.\n\n",zonfile2); printf("Problem opening zonefile %s.\n\n",zonfile2);
exit(0); exit(0);
} }
fwrite(&nzones,1,4,zon2); fwrite(&nzones,1,sizeof(int),zon2);
for (h = 0; h<nzones; h++) { for (h = 0; h<nzones; h++) {
nhlcount = 0; nhlcount = 0;
@ -532,8 +541,8 @@ int main(int argc,char **argv) {
z[h].nhl = nhl; z[h].nhl = nhl;
fwrite(&nhl,1,4,zon2); fwrite(&nhl,1,sizeof(int),zon2);
fwrite(zonelist,nhl,4,zon2); fwrite(zonelist,nhl,sizeof(int),zon2);
} }
fclose(zon2); fclose(zon2);

6
zobov/voz.h
View file

@ -3,7 +3,9 @@
##define NGUARD 42 /*Actually, the number of SPACES between guard points ##define NGUARD 42 /*Actually, the number of SPACES between guard points
in each dim */ in each dim */
typedef int pid_t;
typedef struct Partadj { typedef struct Partadj {
int nadj; pid_t nadj;
int *adj; pid_t *adj;
} PARTADJ; } PARTADJ;

63
zobov/voz1b1.c
View file

@ -12,7 +12,7 @@ int posread(char *posfile, float ***p, float fact);
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
int exitcode; int exitcode;
int i, j, np; pid_t i, j, np;
float **r; float **r;
coordT rtemp[3], *parts; coordT rtemp[3], *parts;
coordT deladjs[3*MAXVERVER], points[3*MAXVERVER]; coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
@ -21,15 +21,16 @@ int main(int argc, char *argv[]) {
char *posfile, outfile[200], *suffix, *outDir; char *posfile, outfile[200], *suffix, *outDir;
PARTADJ *adjs; PARTADJ *adjs;
float *vols; float *vols;
float predict, xmin,xmax,ymin,ymax,zmin,zmax; realT predict, xmin,xmax,ymin,ymax,zmin,zmax;
int *orig; pid_t *orig;
int isitinbuf; int isitinbuf;
char isitinmain, d; char isitinmain, d;
int numdiv, nvp, nvpall, nvpbuf; int numdiv;
float width, width2, totwidth, totwidth2, bf, s, g; pid_t nvp, nvpall, nvpbuf;
realT width, width2, totwidth, totwidth2, bf, s, g;
float border, boxsize; float border, boxsize;
float c[3]; realT c[3];
int b[3]; int b[3];
double totalvol; double totalvol;
@ -117,7 +118,7 @@ int main(int argc, char *argv[]) {
exit(0); exit(0);
} }
DL c[d] = ((float)b[d]+0.5)*width; DL c[d] = ((float)b[d])*width;
printf("c: %f,%f,%f\n",c[0],c[1],c[2]); printf("c: %f,%f,%f\n",c[0],c[1],c[2]);
/* Assign temporary array*/ /* Assign temporary array*/
nvpbuf = 0; /* Number of particles to tesselate, including nvpbuf = 0; /* Number of particles to tesselate, including
@ -142,7 +143,7 @@ int main(int argc, char *argv[]) {
points */ points */
parts = (coordT *)malloc(3*nvpbuf*sizeof(coordT)); parts = (coordT *)malloc(3*nvpbuf*sizeof(coordT));
orig = (int *)malloc(nvpbuf*sizeof(int)); orig = (pid_t *)malloc(nvpbuf*sizeof(pid_t));
if (parts == NULL) { if (parts == NULL) {
printf("Unable to allocate parts\n"); printf("Unable to allocate parts\n");
@ -158,7 +159,7 @@ int main(int argc, char *argv[]) {
for (i=0; i<np; i++) { for (i=0; i<np; i++) {
isitinmain = 1; isitinmain = 1;
DL { DL {
rtemp[d] = r[i][d] - c[d]; rtemp[d] = (realT)r[i][d] - (realT)c[d];
if (rtemp[d] > 0.5) rtemp[d] --; if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++; if (rtemp[d] < -0.5) rtemp[d] ++;
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
@ -182,15 +183,16 @@ int main(int argc, char *argv[]) {
nvpbuf = nvp; nvpbuf = nvp;
for (i=0; i<np; i++) { for (i=0; i<np; i++) {
isitinbuf = 1; isitinbuf = 1;
isitinmain = 1;
DL { DL {
rtemp[d] = r[i][d] - c[d]; rtemp[d] = (realT)r[i][d] - (realT)c[d];
if (rtemp[d] > 0.5) rtemp[d] --; if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++; if (rtemp[d] < -0.5) rtemp[d] ++;
isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2); isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2);
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
} }
if ((isitinbuf > 0) && if (isitinbuf && !isitinmain) {
((fabs(rtemp[0])>width2)||(fabs(rtemp[1])>width2)||(fabs(rtemp[2])>width2))) {
/*printf("%3.3f ",sqrt(rtemp[0]*rtemp[0] + rtemp[1]*rtemp[1] + /*printf("%3.3f ",sqrt(rtemp[0]*rtemp[0] + rtemp[1]*rtemp[1] +
rtemp[2]*rtemp[2])); rtemp[2]*rtemp[2]));
printf("|%2.2f,%2.2f,%2.2f,%f,%f",r[i][0],r[i][1],r[i][2],width2,totwidth2);*/ printf("|%2.2f,%2.2f,%2.2f,%f,%f",r[i][0],r[i][1],r[i][2],width2,totwidth2);*/
@ -221,40 +223,40 @@ int main(int argc, char *argv[]) {
for (i=0; i<NGUARD+1; i++) { for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) { for (j=0; j<NGUARD+1; j++) {
/* Bottom */ /* Bottom */
parts[3*nvpall] = -width2 + (float)i * s; parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (float)j * s; parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = -width2 - g; parts[3*nvpall+2] = -width2 - g;
nvpall++; nvpall++;
/* Top */ /* Top */
parts[3*nvpall] = -width2 + (float)i * s; parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (float)j * s; parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = width2 + g; parts[3*nvpall+2] = width2 + g;
nvpall++; nvpall++;
} }
} }
for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/ for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
for (j=0; j<NGUARD+1; j++) { for (j=0; j<NGUARD+1; j++) {
parts[3*nvpall] = -width2 + (float)i * s; parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 - g; parts[3*nvpall+1] = -width2 - g;
parts[3*nvpall+2] = -width2 + (float)j * s; parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++; nvpall++;
parts[3*nvpall] = -width2 + (float)i * s; parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = width2 + g; parts[3*nvpall+1] = width2 + g;
parts[3*nvpall+2] = -width2 + (float)j * s; parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++; nvpall++;
} }
} }
for (i=0; i<NGUARD+1; i++) { for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) { for (j=0; j<NGUARD+1; j++) {
parts[3*nvpall] = -width2 - g; parts[3*nvpall] = -width2 - g;
parts[3*nvpall+1] = -width2 + (float)i * s; parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (float)j * s; parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++; nvpall++;
parts[3*nvpall] = width2 + g; parts[3*nvpall] = width2 + g;
parts[3*nvpall+1] = -width2 + (float)i * s; parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (float)j * s; parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++; nvpall++;
} }
} }
@ -275,6 +277,11 @@ int main(int argc, char *argv[]) {
/* Do tesselation*/ /* Do tesselation*/
printf("File read. Tessellating ...\n"); fflush(stdout); printf("File read. Tessellating ...\n"); fflush(stdout);
exitcode = delaunadj(parts, nvp, nvpbuf, nvpall, &adjs); exitcode = delaunadj(parts, nvp, nvpbuf, nvpall, &adjs);
if (exitcode != 0)
{
printf("Error while tesselating. Stopping here."); fflush(stdout);
exit(1);
}
/* Now calculate volumes*/ /* Now calculate volumes*/
printf("Now finding volumes ...\n"); fflush(stdout); printf("Now finding volumes ...\n"); fflush(stdout);
@ -323,14 +330,14 @@ int main(int argc, char *argv[]) {
printf("nvp = %d\n",nvp); printf("nvp = %d\n",nvp);
/* Tell us where the original particles were */ /* Tell us where the original particles were */
fwrite(orig,sizeof(int),nvp,out); fwrite(orig,sizeof(pid_t),nvp,out);
/* Volumes*/ /* Volumes*/
fwrite(vols,sizeof(float),nvp,out); fwrite(vols,sizeof(float),nvp,out);
/* Adjacencies */ /* Adjacencies */
for (i=0;i<nvp;i++) { for (i=0;i<nvp;i++) {
fwrite(&(adjs[i].nadj),1,sizeof(int),out); fwrite(&(adjs[i].nadj),1,sizeof(pid_t),out);
if (adjs[i].nadj > 0) if (adjs[i].nadj > 0)
fwrite(adjs[i].adj,adjs[i].nadj,sizeof(int),out); fwrite(adjs[i].adj,adjs[i].nadj,sizeof(pid_t),out);
else printf("0"); else printf("0");
} }
fclose(out); fclose(out);

14
zobov/vozinit.c
View file

@ -5,6 +5,8 @@
#define DL for (d=0;d<3;d++) #define DL for (d=0;d<3;d++)
#define BF 1e30 #define BF 1e30
#define QHULL_MAX_PARTICLES ((1L<<24)-1)
int posread(char *posfile, float ***p, float fact); int posread(char *posfile, float ***p, float fact);
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
@ -22,6 +24,7 @@ int main(int argc, char *argv[]) {
float width, width2, totwidth, totwidth2, bf, s, g; float width, width2, totwidth, totwidth2, bf, s, g;
float border, boxsize; float border, boxsize;
float c[3]; float c[3];
int numGuards;
int b[3]; int b[3];
int numThreads; int numThreads;
int mockIndex; int mockIndex;
@ -145,6 +148,15 @@ int main(int argc, char *argv[]) {
printf("Nvp range: %d,%d\n",nvpmin,nvpmax); printf("Nvp range: %d,%d\n",nvpmin,nvpmax);
printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax); printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax);
numGuards = 6*(NGUARD+1)*(NGUARD+1);
printf("Total max particles: %d\n" , nvpbufmax+numGuards);
if (nvpbufmax+numGuards >= QHULL_MAX_PARTICLES)
{
printf("Too many particles to tesselate per division (Qhull will overflow). Increase divisions or reduce buffer size.\n");
fflush(stdout);
exit(1);
}
/* Output script file */ /* Output script file */
sprintf(scrfile,"scr%s",suffix); sprintf(scrfile,"scr%s",suffix);
printf("Writing script file to %s.\n",scrfile);fflush(stdout); printf("Writing script file to %s.\n",scrfile);fflush(stdout);
@ -152,7 +164,7 @@ int main(int argc, char *argv[]) {
if (scr == NULL) { if (scr == NULL) {
printf("Problem opening script file.\n"); printf("Problem opening script file.\n");
fflush(stdout); fflush(stdout);
exit(0); exit(1);
} }
fprintf(scr,"#!/bin/bash -f\n"); fprintf(scr,"#!/bin/bash -f\n");
p = 0; p = 0;

127
zobov/voztie.c
View file

@ -1,3 +1,4 @@
#include <assert.h>
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
@ -14,10 +15,12 @@ int main(int argc, char *argv[]) {
int numdiv,np,np2,na; int numdiv,np,np2,na;
int i,j,k,p,nout; pid_t i,j,k,p,nout;
int nvp,npnotdone,nvpmax, nvpsum, *orig; pid_t nvp,npnotdone,nvpmax, nvpsum, *orig;
double avgnadj, avgvol; double avgnadj, avgvol;
int numRemoved = 0;
// PMS // PMS
int mockIndex; int mockIndex;
// END PMS // END PMS
@ -87,15 +90,20 @@ int main(int argc, char *argv[]) {
if (adjs == NULL) printf("Couldn't allocate adjs.\n"); if (adjs == NULL) printf("Couldn't allocate adjs.\n");
vols = (float *)malloc(np*sizeof(float)); vols = (float *)malloc(np*sizeof(float));
if (vols == NULL) printf("Couldn't allocate vols.\n"); if (vols == NULL) printf("Couldn't allocate vols.\n");
orig = (int *)malloc(nvpmax*sizeof(int)); orig = (pid_t *)malloc(nvpmax*sizeof(pid_t));
if (orig == NULL) printf("Couldn't allocate orig.\n"); if (orig == NULL) printf("Couldn't allocate orig.\n");
if ((cnt_adj==NULL) || (vols == NULL) || (orig == NULL) || (adjs == NULL)) { if ((cnt_adj==NULL) || (vols == NULL) || (orig == NULL) || (adjs == NULL)) {
printf("Not enough memory to allocate. Exiting.\n"); printf("Not enough memory to allocate. Exiting.\n");
exit(0); exit(0);
} }
for (p=0;p<np;p++) for (p=0;p<np;p++)
vols[p] = -1.; {
vols[p] = -1.;
adjs[p].nadj = 0;
adjs[p].adj = 0;
}
nvpsum = 0;
for (i = 0; i < numdiv; i++) { for (i = 0; i < numdiv; i++) {
for (j = 0; j < numdiv; j++) { for (j = 0; j < numdiv; j++) {
for (k = 0; k < numdiv; k++) { for (k = 0; k < numdiv; k++) {
@ -105,23 +113,19 @@ int main(int argc, char *argv[]) {
printf("Unable to open file %s.\n\n",partfile); printf("Unable to open file %s.\n\n",partfile);
exit(0); exit(0);
} }
fread(&np2,1,sizeof(int),part); fread(&np2,1,sizeof(pid_t),part);
fread(&nvp,1,sizeof(int),part); fread(&nvp,1,sizeof(pid_t),part);
/*printf("nvp = %d\n",nvp);fflush(stdout);*/ /*printf("nvp = %d\n",nvp);fflush(stdout);*/
nvpsum += nvp; nvpsum += nvp;
fread(orig,nvp,sizeof(int),part); fread(orig,nvp,sizeof(pid_t),part);
for (p=0;p<nvp;p++) { for (p=0;p<nvp;p++) {
fread(&volstemp,1,sizeof(float),part); fread(&volstemp,1,sizeof(float),part);
if (vols[orig[p]] > -1.) if (vols[orig[p]] > -1.)
// PMS
if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3 && orig[p] < mockIndex) { if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3 && orig[p] < mockIndex) {
//if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) {
// END PMS
printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n",
orig[p],vols[orig[p]],volstemp); orig[p],vols[orig[p]],volstemp);
// TEST
//exit(0); //exit(0);
} }
vols[orig[p]] = volstemp; vols[orig[p]] = volstemp;
@ -129,14 +133,17 @@ int main(int argc, char *argv[]) {
for (p=0;p<nvp;p++) { for (p=0;p<nvp;p++) {
fread(&na,1,sizeof(int),part); fread(&na,1,sizeof(int),part);
pid_t pid = orig[p];
if (na > 0) { if (na > 0) {
assert(adjs[pid].nadj == 0);// || adjs[pid].nadj == na);
adjs[orig[p]].nadj = na; adjs[orig[p]].nadj = na;
adjs[orig[p]].adj = (int *)malloc(na*sizeof(int)); if (adjs[pid].adj == 0)
if (adjs[orig[p]].adj == NULL) { adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t));
if (adjs[orig[p]].adj == 0) {
printf("Couldn't allocate adjs[orig[%d]].adj.\n",p); printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
exit(0); exit(0);
} }
fread(adjs[orig[p]].adj,na,sizeof(int),part); fread(adjs[orig[p]].adj,na,sizeof(pid_t),part);
} else { } else {
printf("0"); fflush(stdout); printf("0"); fflush(stdout);
} }
@ -157,53 +164,37 @@ int main(int argc, char *argv[]) {
adjs[i].nadj = 0; adjs[i].nadj = 0;
} }
int numRemoved = 0; // unlink particles adjacent to mock galaxies
// unlink particles adjacent to mock galaxies
for (i = 0; i < mockIndex; i++) {
for (j = 0; j < adjs[i].nadj; j++) {
if (adjs[i].adj[j] > mockIndex) {
//printf("KILLING %d\n", i);
vols[i] = 1.e-29;
adjs[i].nadj = 0;
numRemoved++;
break;
}
}
}
// update all other adjacencies
for (i = 0; i < mockIndex; i++) {
int numAdjSaved = 0;
for (j = 0; j < adjs[i].nadj; j++) {
//if ( vols[adjs[i].adj[j]] != -99) {
if ( adjs[adjs[i].adj[j]].nadj != 0) {
adjs[i].adj[numAdjSaved] = adjs[i].adj[j];
numAdjSaved++;
}
}
adjs[i].nadj = numAdjSaved;
}
/*
for (i = 0; i < mockIndex; i++) { for (i = 0; i < mockIndex; i++) {
printf("ADJ: %d %d : ", i, adjs[i].nadj);
for (j = 0; j < adjs[i].nadj; j++) { for (j = 0; j < adjs[i].nadj; j++) {
printf(" %d", adjs[i].adj[j]); if (adjs[i].adj[j] > mockIndex) {
//printf("KILLING %d\n", i);
vols[i] = 1.e-29;
adjs[i].nadj = 0;
numRemoved++;
break;
}
} }
printf("\n");
} }
*/
// update all other adjacencies
for (i = 0; i < mockIndex; i++) {
int numAdjSaved = 0;
for (j = 0; j < adjs[i].nadj; j++) {
printf("Removed %d mock galaxies and %d adjacent galaxies.\n", np-mockIndex, if ( adjs[adjs[i].adj[j]].nadj != 0) {
adjs[i].adj[numAdjSaved] = adjs[i].adj[j];
numAdjSaved++;
}
}
adjs[i].nadj = numAdjSaved;
}
printf("Removed %d mock galaxies and %d adjacent galaxies.\n", np-mockIndex,
numRemoved); numRemoved);
printf("There are %d galaxies remaining.\n", mockIndex-numRemoved); printf("There are %d galaxies remaining.\n", mockIndex-numRemoved);
// END PMS // END PMS
@ -217,7 +208,7 @@ printf("\n");
avgvol += (double)(vols[p]); avgvol += (double)(vols[p]);
} }
if (npnotdone > 0) if (npnotdone > 0)
printf("%d particles not done!\n"); printf("%d particles not done!\n", npnotdone);
printf("%d particles done more than once.\n",nvpsum-np); printf("%d particles done more than once.\n",nvpsum-np);
avgnadj /= (double)np; avgnadj /= (double)np;
avgvol /= (double)np; avgvol /= (double)np;
@ -237,44 +228,28 @@ printf("\n");
printf("Unable to open %s\n",adjfile); printf("Unable to open %s\n",adjfile);
exit(0); exit(0);
} }
// PMS
fwrite(&mockIndex,1, sizeof(int),adj); fwrite(&mockIndex,1, sizeof(int),adj);
//fwrite(&np,1, sizeof(int),adj);
// END OMS
/* Adjacencies: first the numbers of adjacencies, /* Adjacencies: first the numbers of adjacencies,
and the number we're actually going to write per particle */ and the number we're actually going to write per particle */
// Recount the number of adjacencies after merge // Recount the number of adjacencies after merge
for(i=0;i<mockIndex;i++) for(i=0;i<mockIndex;i++)
cnt_adj[i] = 0; cnt_adj[i] = adjs[i].nadj;
for(i=0;i<mockIndex;i++)
{
for (j=0;j<adjs[i].nadj;j++)
{
cnt_adj[adjs[i].adj[j]]++;
cnt_adj[i]++;
}
}
// PMS
for (i=0;i<mockIndex;i++) for (i=0;i<mockIndex;i++)
//for (i=0;i<np;i++)
// END PMS
fwrite(&cnt_adj[i],1,sizeof(int),adj); fwrite(&cnt_adj[i],1,sizeof(int),adj);
/* Now the lists of adjacencies (without double counting) */ /* Now the lists of adjacencies (without double counting) */
// PMS
for (i=0;i<mockIndex;i++) { for (i=0;i<mockIndex;i++) {
//for (i=0;i<np;i++) {
//if (adjs[i].nadj > 0) {
// END PMS
nout = 0; nout = 0;
for (j=0;j<adjs[i].nadj; j++) if (adjs[i].adj[j] > i) nout++; for (j=0;j<adjs[i].nadj; j++)
if (adjs[i].adj[j] > i)
nout++;
fwrite(&nout,1,sizeof(int),adj); fwrite(&nout,1,sizeof(int),adj);
for (j=0;j<adjs[i].nadj; j++) { for (j=0;j<adjs[i].nadj; j++) {
int id = adjs[i].adj[j]; pid_t id = adjs[i].adj[j];
if (id > i) if (id > i)
fwrite(&(adjs[i].adj[j]),1,sizeof(int),adj); fwrite(&id,1,sizeof(pid_t),adj);
} }
} }

2
zobov/vozutil.c
View file

@ -46,7 +46,7 @@ int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs)
} }
/* Delaunay triangulation*/ /* Delaunay triangulation*/
sprintf (flags, "qhull s d QJ"); sprintf (flags, "qhull s d Qt");
exitcode= qh_new_qhull (dim, nvpall, points, ismalloc, exitcode= qh_new_qhull (dim, nvpall, points, ismalloc,
flags, outfile, errfile); flags, outfile, errfile);