renamed python_tools to python_src to maintain consistency

python_src/vide/voidUtil/__init__.py

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

python_src/vide/voidUtil/catalogUtil.py

@@ -0,0 +1,586 @@
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/voidUtil/catalogUtil.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+#
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+
+# Various utilities for loading and modifying particle datasets
+
+import numpy as np
+from netCDF4 import Dataset
+import sys
+from vide.backend import *
+import vide.apTools as vp
+import pickle
+from .periodic_kdtree import PeriodicCKDTree
+import os
+
+NetCDFFile = Dataset
+ncFloat = 'f8'
+
+# -----------------------------------------------------------------------------
+def loadPart(sampleDir):
+    print("    Loading particle data...")
+    sys.stdout.flush()
+
+    with open(sampleDir+"/sample_info.dat", 'rb') as input:
+      sample = pickle.load(input)
+
+    infoFile = sampleDir+"/zobov_slice_"+sample.fullName+".par"
+    File = NetCDFFile(infoFile, 'r')
+    ranges = np.zeros((3,2))
+    ranges[0][0] = getattr(File, 'range_x_min')
+    ranges[0][1] = getattr(File, 'range_x_max')
+    ranges[1][0] = getattr(File, 'range_y_min')
+    ranges[1][1] = getattr(File, 'range_y_max')
+    ranges[2][0] = getattr(File, 'range_z_min')
+    ranges[2][1] = getattr(File, 'range_z_max')
+    isObservation = getattr(File, 'is_observation')
+    maskIndex = getattr(File, 'mask_index')
+    File.close()
+    mul = np.zeros((3))
+    mul[:] = ranges[:,1] - ranges[:,0]
+
+    partFile = sampleDir+"/zobov_slice_"+sample.fullName
+    iLine = 0
+    partData = []
+    part = np.zeros((3))
+    with open(partFile, mode="rb") as File:
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      Np = np.fromfile(File, dtype=np.int32,count=1)[0]
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      x = np.fromfile(File, dtype=np.float32,count=Np)
+      x *= mul[0]
+      if isObservation != 1:
+        x += ranges[0][0]
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      y = np.fromfile(File, dtype=np.float32,count=Np)
+      y *= mul[1]
+      if isObservation != 1:
+        y += ranges[1][0]
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      z = np.fromfile(File, dtype=np.float32,count=Np)
+      z *= mul[2]
+      if isObservation != 1:
+        z += ranges[2][0]
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      RA = np.fromfile(File, dtype=np.float32,count=Np)
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      Dec = np.fromfile(File, dtype=np.float32,count=Np)
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      redshift = np.fromfile(File, dtype=np.float32,count=Np)
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      uniqueID = np.fromfile(File, dtype=np.int64,count=Np)
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+
+    if isObservation == 1:
+      x = x[0:maskIndex]# * 100/300000
+      y = y[0:maskIndex]# * 100/300000
+      z = z[0:maskIndex]# * 100/300000
+      RA = RA[0:maskIndex]
+      Dec = Dec[0:maskIndex]
+      redshift = redshift[0:maskIndex]
+      uniqueID = uniqueID[0:maskIndex]
+
+    partData = np.column_stack((x,y,z))
+
+    extraData = np.column_stack((RA,Dec,redshift,uniqueID))
+
+    boxLen = mul
+
+    #if isObservation == 1:
+    #  # look for the mask file
+    #  if os.access(sample.maskFile, os.F_OK):
+    #    maskFile = sample.maskFile
+    #  else:
+    #    maskFile = sampleDir+"/"+os.path.basename(sample.maskFile)
+    #    print "Using maskfile found in:", maskFile
+    #  props = vp.getSurveyProps(maskFile, sample.zBoundary[0],
+    #                            sample.zBoundary[1],
+    #                            sample.zBoundary[0],
+    #                            sample.zBoundary[1], "all",
+    #                            selectionFuncFile=sample.selFunFile,
+    #                            useComoving=sample.useComoving)
+    #  boxVol = props[0]
+    #  volNorm = maskIndex/boxVol
+    #else:
+    boxVol = np.prod(boxLen)
+    volNorm = len(x)/boxVol
+
+    isObservationData = isObservation == 1
+
+    return partData, boxLen, volNorm, isObservationData, ranges, extraData
+
+# -----------------------------------------------------------------------------
+def getVolNorm(sampleDir):
+    with open(sampleDir+"/sample_info.dat", 'rb') as input:
+      sample = pickle.load(input)
+
+    infoFile = sampleDir+"/zobov_slice_"+sample.fullName+".par"
+    File = NetCDFFile(infoFile, 'r')
+    ranges = np.zeros((3,2))
+    ranges[0][0] = getattr(File, 'range_x_min')
+    ranges[0][1] = getattr(File, 'range_x_max')
+    ranges[1][0] = getattr(File, 'range_y_min')
+    ranges[1][1] = getattr(File, 'range_y_max')
+    ranges[2][0] = getattr(File, 'range_z_min')
+    ranges[2][1] = getattr(File, 'range_z_max')
+    isObservation = getattr(File, 'is_observation')
+    maskIndex = getattr(File, 'mask_index')
+    File.close()
+    mul = np.zeros((3))
+    mul[:] = ranges[:,1] - ranges[:,0]
+
+    partFile = sampleDir+"/zobov_slice_"+sample.fullName
+    with open(partFile, mode="rb") as File:
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      Np = np.fromfile(File, dtype=np.int32,count=1)[0]
+
+    boxLen = mul
+
+    #if isObservation == 1:
+    #  # look for the mask file
+    #  if os.access(sample.maskFile, os.F_OK):
+    #    maskFile = sample.maskFile
+    #  else:
+    #    maskFile = sampleDir+"/"+os.path.basename(sample.maskFile)
+    #    print "Using maskfile found in:", maskFile
+    #  props = vp.getSurveyProps(maskFile, sample.zBoundary[0],
+    #                            sample.zBoundary[1],
+    #                            sample.zBoundary[0],
+    #                            sample.zBoundary[1], "all",
+    #                            selectionFuncFile=sample.selFunFile,
+    #                            useComoving=sample.useComoving)
+    #  boxVol = props[0]
+    #  volNorm = maskIndex/boxVol
+    #else:
+    boxVol = np.prod(boxLen)
+    volNorm = Np/boxVol
+
+    return volNorm
+
+
+# -----------------------------------------------------------------------------
+def loadPartVel(sampleDir):
+    #print "    Loading particle velocities..."
+    sys.stdout.flush()
+
+    with open(sampleDir+"/sample_info.dat", 'rb') as input:
+      sample = pickle.load(input)
+
+    infoFile = sampleDir+"/zobov_slice_"+sample.fullName+".par"
+    File = NetCDFFile(infoFile, 'r')
+    isObservation = getattr(File, 'is_observation')
+
+    if isObservation:
+      print("No velocities for observations!")
+      return -1
+
+    vx = File.variables['vel_x'][0:]
+    vy = File.variables['vel_y'][0:]
+    vz = File.variables['vel_z'][0:]
+
+    File.close()
+
+    partVel = np.column_stack((vx,vy,vz))
+
+    return partVel
+
+# -----------------------------------------------------------------------------
+def getPartTree(catalog):
+
+  sample = catalog.sampleInfo
+  partData = catalog.partPos
+  boxLen = catalog.boxLen
+
+  periodicLine = getPeriodic(sample)
+
+  periodic = 1.*boxLen
+
+  if not "x" in periodicLine: periodic[0] = -1
+  if not "y" in periodicLine: periodic[1] = -1
+  if not "z" in periodicLine: periodic[2] = -1
+
+  return PeriodicCKDTree(periodic, partData)
+
+# -----------------------------------------------------------------------------
+def getBall(partTree, center, radius):
+
+  return partTree.query_ball_point(center, r=radius)
+
+
+# -----------------------------------------------------------------------------
+def shiftPart(inPart, center, periodicLine, ranges):
+
+  part = inPart.copy()
+  newCenter = 1.*center;
+
+  boxLen = np.zeros((3))
+  boxLen[0] = ranges[0][1] - ranges[0][0]
+  boxLen[1] = ranges[1][1] - ranges[1][0]
+  boxLen[2] = ranges[2][1] - ranges[2][0]
+
+  # shift to box coordinates
+  part[:,0] -= ranges[0][0]
+  part[:,1] -= ranges[1][0]
+  part[:,2] -= ranges[2][0]
+
+  newCenter[:] -= ranges[:,0]
+
+  part[:,0] -= newCenter[0]
+  part[:,1] -= newCenter[1]
+  part[:,2] -= newCenter[2]
+
+  shiftUs = np.abs(part[:,0]) > boxLen[0]/2.
+  if ("x" in periodicLine): part[shiftUs,0] -= \
+                            np.copysign(boxLen[0], part[shiftUs,0])
+  shiftUs = np.abs(part[:,1]) > boxLen[1]/2.
+  if ("y" in periodicLine): part[shiftUs,1] -= \
+                            np.copysign(boxLen[1], part[shiftUs,1])
+  shiftUs = np.abs(part[:,2]) > boxLen[2]/2.
+  if ("z" in periodicLine): part[shiftUs,2] -= \
+                            np.copysign(boxLen[2], part[shiftUs,2])
+
+
+  #part[:,0] += ranges[0][0]
+  #part[:,1] += ranges[1][0]
+  #part[:,2] += ranges[2][0]
+
+  return part
+
+# -----------------------------------------------------------------------------
+
+class Bunch:
+    def __init__(self, **kwds):
+        self.__dict__.update(kwds)
+
+class Catalog:
+  numVoids = 0
+  numPartTot = 0
+  numZonesTot = 0
+  volNorm = 0
+  boxLen = np.zeros((3))
+  ranges = np.zeros((3,2))
+  part = None
+  zones2Parts = None
+  void2Zones = None
+  voids = None
+  sampleInfo = None
+
+# -----------------------------------------------------------------------------
+def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True,
+                    untrimmed=False):
+
+# loads a void catalog
+# by default, loads parent-level voids with central densities greater than 0.2*mean
+#   sampleDir: path to VIDE output directory
+#   dataPortion: "central" or "all"
+#   loadParticles: if True, also load particle information
+#   untrimmed: if True, catalog contains all voids, regardless of density or hierarchy level
+
+  sys.stdout.flush()
+
+  catalog = Catalog()
+  with open(sampleDir+"/sample_info.dat", 'rb') as input:
+    sample = pickle.load(input)
+  catalog.sampleInfo = sample
+
+  print("Loading info...")
+  infoFile = sampleDir+"/zobov_slice_"+sample.fullName+".par"
+  File = NetCDFFile(infoFile, 'r')
+  ranges = np.zeros((3,2))
+  ranges[0][0] = getattr(File, 'range_x_min')
+  ranges[0][1] = getattr(File, 'range_x_max')
+  ranges[1][0] = getattr(File, 'range_y_min')
+  ranges[1][1] = getattr(File, 'range_y_max')
+  ranges[2][0] = getattr(File, 'range_z_min')
+  ranges[2][1] = getattr(File, 'range_z_max')
+  catalog.boxLen[0] = ranges[0][1] - ranges[0][0]
+  catalog.boxLen[1] = ranges[1][1] - ranges[1][0]
+  catalog.boxLen[2] = ranges[2][1] - ranges[2][0]
+  catalog.ranges = ranges
+  File.close()
+
+  volNorm = getVolNorm(sampleDir)
+  catalog.volNorm = volNorm
+
+  if untrimmed:
+    prefix = "untrimmed_"
+  else:
+    prefix = ""
+
+  print("Loading voids...")
+  fileName = sampleDir+"/"+prefix+"voidDesc_"+dataPortion+"_"+sample.fullName+".out"
+  catData = np.loadtxt(fileName, comments="#", skiprows=2)
+  catalog.voids = []
+  for line in catData:
+    catalog.voids.append(Bunch(iVoid = int(line[0]),
+                               voidID = int(line[1]),
+                               coreParticle = line[2],
+                               coreDens = line[3],
+                               zoneVol = line[4],
+                               zoneNumPart = line[5],
+                               numZones = int(line[6]),
+                               voidVol = line[7],
+                               numPart = int(line[8]),
+                               densCon = line[9],
+                               voidProb = line[10],
+                               radius = 0., # this is read in later
+                               macrocenter = np.zeros((3)),
+                               redshift = 0,
+                               RA = 0,
+                               Dec = 0,
+                               parentID = 0,
+                               treeLevel = 0,
+                               numChildren = 0,
+                               centralDen = 0.,
+                               ellipticity = 0.,
+                               eigenVals = np.zeros((3)),
+                               eigenVecs = np.zeros((3,3)),
+                               ))
+
+  catalog.numVoids = len(catalog.voids)
+  print("Read %d voids" % catalog.numVoids)
+
+  print("Loading macrocenters...")
+  iLine = 0
+  for line in open(sampleDir+"/"+prefix+"macrocenters_"+dataPortion+"_"+sample.fullName+".out"):
+    line = line.split()
+    catalog.voids[iLine].macrocenter[0] = float(line[1])
+    catalog.voids[iLine].macrocenter[1] = float(line[2])
+    catalog.voids[iLine].macrocenter[2] = float(line[3])
+    iLine += 1
+
+  iLine = 0
+  fileName = sampleDir+"/"+prefix+"sky_positions_"+dataPortion+"_"+sample.fullName+".out"
+  catData = np.loadtxt(fileName, comments="#")
+  for line in catData:
+    catalog.voids[iLine].RA = float(line[0])
+    catalog.voids[iLine].Dec = float(line[1])
+    iLine += 1
+
+
+  print("Loading derived void information...")
+  fileName = sampleDir+"/"+prefix+"centers_"+dataPortion+"_"+sample.fullName+".out"
+  catData = np.loadtxt(fileName, comments="#")
+  for (iLine,line) in enumerate(catData):
+    catalog.voids[iLine].volume = float(line[6])
+    catalog.voids[iLine].radius = float(line[4])
+    catalog.voids[iLine].redshift = float(line[5])
+    catalog.voids[iLine].parentID = float(line[10])
+    catalog.voids[iLine].treeLevel = float(line[11])
+    catalog.voids[iLine].numChildren = float(line[12])
+    catalog.voids[iLine].centralDen = float(line[13])
+    iLine += 1
+
+  fileName = sampleDir+"/"+prefix+"shapes_"+dataPortion+"_"+sample.fullName+".out"
+  catData = np.loadtxt(fileName, comments="#")
+  for (iLine,line) in enumerate(catData):
+    catalog.voids[iLine].ellipticity = float(line[1])
+
+    catalog.voids[iLine].eigenVals[0] = float(line[2])
+    catalog.voids[iLine].eigenVals[1] = float(line[3])
+    catalog.voids[iLine].eigenVals[2] = float(line[4])
+
+    catalog.voids[iLine].eigenVecs[0][0] = float(line[5])
+    catalog.voids[iLine].eigenVecs[0][1] = float(line[6])
+    catalog.voids[iLine].eigenVecs[0][2] = float(line[7])
+
+    catalog.voids[iLine].eigenVecs[1][0] = float(line[8])
+    catalog.voids[iLine].eigenVecs[1][1] = float(line[9])
+    catalog.voids[iLine].eigenVecs[1][2] = float(line[10])
+
+    catalog.voids[iLine].eigenVecs[2][0] = float(line[11])
+    catalog.voids[iLine].eigenVecs[2][1] = float(line[12])
+    catalog.voids[iLine].eigenVecs[2][2] = float(line[13])
+
+    iLine += 1
+
+  if loadParticles:
+    print("Loading all particles...")
+    partData, boxLen, volNorm, isObservationData, ranges, extraData = loadPart(sampleDir)
+    numPartTot = len(partData)
+    catalog.numPartTot = numPartTot
+    catalog.partPos = partData
+    catalog.part = []
+    for i in range(len(partData)):
+      catalog.part.append(Bunch(x = partData[i][0],
+                                y = partData[i][1],
+                                z = partData[i][2],
+                                volume = 0,
+                                ra = extraData[i][0],
+                                dec = extraData[i][1],
+                                redshift = extraData[i][2],
+                                uniqueID = extraData[i][3]))
+
+
+    print("Loading volumes...")
+    volFile = sampleDir+"/vol_"+sample.fullName+".dat"
+    with open(volFile, mode="rb") as File:
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      vols = np.fromfile(File, dtype=np.float32,count=numPartTot)
+    for ivol in range(len(vols)):
+      catalog.part[ivol].volume = vols[ivol] / volNorm
+
+    print("Loading zone-void membership info...")
+    zoneFile = sampleDir+"/voidZone_"+sample.fullName+".dat"
+    catalog.void2Zones = []
+    with open(zoneFile, mode="rb") as File:
+      numZonesTot = np.fromfile(File, dtype=np.int32,count=1)[0]
+      catalog.numZonesTot = numZonesTot
+      for iZ in range(numZonesTot):
+        numZones = np.fromfile(File, dtype=np.int32,count=1)[0]
+        catalog.void2Zones.append(Bunch(numZones = numZones,
+                                        zoneIDs = []))
+
+        for p in range(numZones):
+          zoneID = np.fromfile(File, dtype=np.int32,count=1)[0]
+          catalog.void2Zones[iZ].zoneIDs.append(zoneID)
+
+    print("Loading particle-zone membership info...")
+    zonePartFile = sampleDir+"/voidPart_"+sample.fullName+".dat"
+    catalog.zones2Parts = []
+    with open(zonePartFile) as File:
+      chk = np.fromfile(File, dtype=np.int32,count=1)
+      numZonesTot = np.fromfile(File, dtype=np.int32,count=1)[0]
+      for iZ in range(numZonesTot):
+        numPart = np.fromfile(File, dtype=np.int32,count=1)[0]
+        catalog.zones2Parts.append(Bunch(numPart = numPart,
+                                         partIDs = []))
+
+        for p in range(numPart):
+          partID = np.fromfile(File, dtype=np.int32,count=1)[0]
+          catalog.zones2Parts[iZ].partIDs.append(partID)
+
+  return catalog
+
+
+# -----------------------------------------------------------------------------
+def getArray(objectList, attr):
+
+  if hasattr(objectList[0], attr):
+    ndim = np.shape( np.atleast_1d( getattr(objectList[0], attr) ) )[0]
+    attrArr = np.zeros(( len(objectList), ndim ))
+
+    for idim in range(ndim):
+      attrArr[:,idim] = np.fromiter((np.atleast_1d(getattr(v, attr))[idim] \
+                                    for v in objectList), float )
+
+    if ndim == 1: attrArr = attrArr[:,0]
+
+    return attrArr
+  else:
+    print(" Attribute", attr, "not found!")
+    return -1
+
+
+# -----------------------------------------------------------------------------
+def getVoidPart(catalog, voidID):
+
+  partOut = []
+  for iZ in range(catalog.void2Zones[voidID].numZones):
+    zoneID = catalog.void2Zones[voidID].zoneIDs[iZ]
+    for p in range(catalog.zones2Parts[zoneID].numPart):
+      partID = catalog.zones2Parts[zoneID].partIDs[p]
+      partOut.append(catalog.part[partID])
+
+  return partOut
+
+# -----------------------------------------------------------------------------
+def filterVoidsOnSize(catalog, rMin):
+  catalog.voids = [v for v in catalog.voids if v.radius >= rMin]
+  catalog.numVoids = len(catalog.voids)
+  return catalog
+
+# -----------------------------------------------------------------------------
+def filterVoidsOnTreeLevel(catalog, level):
+  catalog.voids = [v for v in catalog.voids if v.treeLevel == level]
+  if level == -1:
+    catalog.voids = [v for v in catalog.voids if v.numChildren == 0]
+  catalog.numVoids = len(catalog.voids)
+  return catalog
+
+# -----------------------------------------------------------------------------
+def filterVoidsOnCentralDen(catalog, maxCentralDen):
+  catalog.voids = [v for v in catalog.voids if v.centralDen <= maxCentralDen]
+  catalog.numVoids = len(catalog.voids)
+  return catalog
+
+# -----------------------------------------------------------------------------
+def filterVoidsOnCoreDen(catalog, maxCoreDen):
+  catalog.voids = [v for v in catalog.voids if v.coreDens <= maxCoreDen]
+  catalog.numVoids = len(catalog.voids)
+  return catalog
+
+# -----------------------------------------------------------------------------
+def filterVoidsOnDenCon(catalog, minDenCon):
+  catalog.voids = [v for v in catalog.voids if v.densCon >= minDenCon]
+
+  return catalog
+
+
+
+# -----------------------------------------------------------------------------
+def stackVoids(catalog, stackMode = "ball"):
+
+# builds a stack of voids from the given catalog
+#   catalog: void catalog
+#   stackMode:
+#     "ball": spherical cut
+#     "voronoi": only void member particles
+#
+# returns:
+#   stackedPart: array of relative particle positions in the stack
+
+  rMax = 100.
+  periodicLine = getPeriodic(catalog.sampleInfo)
+
+  if stackMode == "ball":
+    partTree = getPartTree(catalog)
+
+  stackedPart = []
+  for void in catalog.voids:
+    center = void.macrocenter
+
+    if stackMode == "ball":
+      localPart = catalog.partPos[ getBall(partTree, center, rMax) ]
+    else:
+      voidPart = getVoidPart(catalog, void.voidID)
+      localPart = np.zeros((len(voidPart),3))
+      localPart[:,0] = getArray(voidPart, 'x')
+      localPart[:,1] = getArray(voidPart, 'y')
+      localPart[:,2] = getArray(voidPart, 'z')
+
+    shiftedPart = shiftPart(localPart, center, periodicLine, catalog.ranges)
+
+    stackedPart.extend(shiftedPart)
+
+  return stackedPart

python_src/vide/voidUtil/matchUtil.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/voidUtil/matchUtil.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+# 
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+__all__=['compareCatalogs',]
+
+from vide.backend import *
+import imp
+import pickle
+import os
+import matplotlib.pyplot as plt
+import numpy as np
+import argparse
+
+def compareCatalogs(baseCatalogDir, compareCatalogDir, 
+                    outputDir="./", logDir="./", 
+                    matchMethod="useID", dataPortion="central",
+                    strictMatch=True, 
+                    overlapFrac=0.25,
+                    pathToCTools="../../../c_tools"):
+
+# reports the overlap between two void catalogs
+#  baseCatalogDir: directory of catalog 1
+#  compareCatagDir: directory of catalog 2
+#  outputDir: directory to place outputs
+#  logDir: directory to place log files
+#  matchMethod: "useID" to use unique IDs, "prox" to use overlap of Voronoi cells
+#  dataPortion: "central" or "all"
+#  strictMatch: if True, only attempt to match to trimmed catalog
+#  overlapFrac: threshold fraction of Voronoi radius to count as matched
+#  pathToCTools: path to location of VIDE c_tools directory
+
+  if not os.access(outputDir, os.F_OK):
+    os.makedirs(outputDir)
+  
+  if not os.access(logDir, os.F_OK):
+    os.makedirs(logDir)
+  
+  outFileName = outputDir + "/" + "voidOverlap" #+ ".dat"
+  
+  with open(baseCatalogDir+"/sample_info.dat", 'rb') as input:
+    baseSample = pickle.load(input)
+  
+  with open(compareCatalogDir+"/sample_info.dat", 'rb') as input:
+    sample = pickle.load(input)
+  
+  print(" Comparing", baseSample.fullName, "to", sample.fullName, "...", end=' ')
+  sys.stdout.flush()
+  
+  sampleName = sample.fullName
+  
+  binPath = pathToCTools+"/analysis/voidOverlap"
+  logFile = logDir+"/compare_"+baseSample.fullName+"_"+sampleName+".out"
+  stepOutputFileName = outFileName + "_" + baseSample.fullName + "_" + \
+                       sampleName + "_"
+  
+  launchVoidOverlap(baseSample, sample, baseCatalogDir, 
+                    compareCatalogDir, binPath, 
+                    thisDataPortion=dataPortion, logFile=logFile,
+                    continueRun=False, 
+                    outputFile=stepOutputFileName,
+                    matchMethod=matchMethod,
+                    overlapFrac=overlapFrac,
+                    strictMatch=strictMatch)
+  
+  print(" Done!")
+  return

python_src/vide/voidUtil/periodic_kdtree.py

@@ -0,0 +1,394 @@
+# periodic_kdtree.py
+#
+# A wrapper around scipy.spatial.kdtree to implement periodic boundary
+# conditions
+#
+# Written by Patrick Varilly, 6 Jul 2012
+# Released under the scipy license
+
+import numpy as np
+from scipy.spatial import KDTree, cKDTree
+import itertools
+import heapq
+
+def _gen_relevant_images(x, bounds, distance_upper_bound):
+    # Map x onto the canonical unit cell, then produce the relevant
+    # mirror images
+    real_x = x - np.where(bounds > 0.0,
+                          np.floor(x / bounds) * bounds, 0.0)
+    m = len(x)
+
+    xs_to_try = [real_x]
+    for i in range(m):
+        if bounds[i] > 0.0:
+            disp = np.zeros(m)
+            disp[i] = bounds[i]
+
+            if distance_upper_bound == np.inf:
+                xs_to_try = list(
+                    itertools.chain.from_iterable(
+                        (_ + disp, _, _ - disp) for _ in xs_to_try))
+            else:
+                extra_xs = []
+
+                # Point near lower boundary, include image on upper side
+                if abs(real_x[i]) < distance_upper_bound:
+                    extra_xs.extend(_ + disp for _ in xs_to_try)
+
+                # Point near upper boundary, include image on lower side
+                if abs(bounds[i] - real_x[i]) < distance_upper_bound:
+                    extra_xs.extend(_ - disp for _ in xs_to_try)
+
+                xs_to_try.extend(extra_xs)
+
+    return xs_to_try
+
+
+class PeriodicKDTree(KDTree):
+    """
+    kd-tree for quick nearest-neighbor lookup with periodic boundaries
+
+    See scipy.spatial.kdtree for details on kd-trees.
+
+    Searches with periodic boundaries are implemented by mapping all
+    initial data points to one canonical periodic image, building an
+    ordinary kd-tree with these points, then querying this kd-tree multiple
+    times, if necessary, with all the relevant periodic images of the
+    query point.
+
+    Note that to ensure that no two distinct images of the same point
+    appear in the results, it is essential to restrict the maximum
+    distance between a query point and a data point to half the smallest
+    box dimension.
+    """
+
+    def __init__(self, bounds, data, leafsize=10):
+        """Construct a kd-tree.
+
+        Parameters
+        ----------
+        bounds : array_like, shape (k,)
+            Size of the periodic box along each spatial dimension.  A
+            negative or zero size for dimension k means that space is not
+            periodic along k.
+        data : array_like, shape (n,k)
+            The data points to be indexed. This array is not copied, and
+            so modifying this data will result in bogus results.
+        leafsize : positive int
+            The number of points at which the algorithm switches over to
+            brute-force.
+        """
+
+        # Map all points to canonical periodic image
+        self.bounds = np.array(bounds)
+        self.real_data = np.asarray(data)
+        wrapped_data = (
+            self.real_data - np.where(bounds > 0.0,
+                (np.floor(self.real_data / bounds) * bounds), 0.0))
+
+        # Calculate maximum distance_upper_bound
+        self.max_distance_upper_bound = np.min(
+            np.where(self.bounds > 0, 0.5 * self.bounds, np.inf))
+
+        # Set up underlying kd-tree
+        super(PeriodicKDTree, self).__init__(wrapped_data, leafsize)
+
+    # The following name is a kludge to override KDTree's private method
+    def _KDTree__query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):
+        # This is the internal query method, which guarantees that x
+        # is a single point, not an array of points
+        #
+        # A slight complication: k could be "None", which means "return
+        # all neighbors within the given distance_upper_bound".
+
+        # Cap distance_upper_bound
+        distance_upper_bound = min([distance_upper_bound,
+                                      self.max_distance_upper_bound])
+
+        # Run queries over all relevant images of x
+        hits_list = []
+        for real_x in _gen_relevant_images(x, self.bounds, distance_upper_bound):
+            hits_list.append(
+                super(PeriodicKDTree, self)._KDTree__query(
+                    real_x, k, eps, p, distance_upper_bound))
+
+        # Now merge results
+        if k is None:
+            return list(heapq.merge(*hits_list))
+        elif k > 1:
+            return heapq.nsmallest(k, itertools.chain(*hits_list))
+        elif k == 1:
+            return [min(itertools.chain(*hits_list))]
+        else:
+            raise ValueError("Invalid k in periodic_kdtree._KDTree__query")
+
+    # The following name is a kludge to override KDTree's private method
+    def _KDTree__query_ball_point(self, x, r, p=2., eps=0):
+        # This is the internal query method, which guarantees that x
+        # is a single point, not an array of points
+
+        # Cap r
+        r = min(r, self.max_distance_upper_bound)
+
+        # Run queries over all relevant images of x
+        results = []
+        for real_x in _gen_relevant_images(x, self.bounds, r):
+            results.extend(
+                super(PeriodicKDTree, self)._KDTree__query_ball_point(
+                    real_x, r, p, eps))
+        return results
+
+    def query_ball_tree(self, other, r, p=2., eps=0):
+        raise NotImplementedError()
+
+    def query_pairs(self, r, p=2., eps=0):
+        raise NotImplementedError()
+
+    def count_neighbors(self, other, r, p=2.):
+        raise NotImplementedError()
+
+    def sparse_distance_matrix(self, other, max_distance, p=2.):
+        raise NotImplementedError()
+
+
+class PeriodicCKDTree(cKDTree):
+    """
+    Cython kd-tree for quick nearest-neighbor lookup with periodic boundaries
+
+    See scipy.spatial.ckdtree for details on kd-trees.
+
+    Searches with periodic boundaries are implemented by mapping all
+    initial data points to one canonical periodic image, building an
+    ordinary kd-tree with these points, then querying this kd-tree multiple
+    times, if necessary, with all the relevant periodic images of the
+    query point.
+
+    Note that to ensure that no two distinct images of the same point
+    appear in the results, it is essential to restrict the maximum
+    distance between a query point and a data point to half the smallest
+    box dimension.
+    """
+
+    def __init__(self, bounds, data, leafsize=10):
+        """Construct a kd-tree.
+
+        Parameters
+        ----------
+        bounds : array_like, shape (k,)
+            Size of the periodic box along each spatial dimension.  A
+            negative or zero size for dimension k means that space is not
+            periodic along k.
+        data : array-like, shape (n,m)
+            The n data points of dimension mto be indexed. This array is
+            not copied unless this is necessary to produce a contiguous
+            array of doubles, and so modifying this data will result in
+            bogus results.
+        leafsize : positive integer
+            The number of points at which the algorithm switches over to
+            brute-force.
+        """
+
+        # Map all points to canonical periodic image
+        self.bounds = np.array(bounds)
+        self.real_data = np.asarray(data)
+        wrapped_data = (
+            self.real_data - np.where(bounds > 0.0,
+                (np.floor(self.real_data / bounds) * bounds), 0.0))
+
+        # Calculate maximum distance_upper_bound
+        self.max_distance_upper_bound = np.min(
+            np.where(self.bounds > 0, 0.5 * self.bounds, np.inf))
+
+        # Set up underlying kd-tree
+        super(PeriodicCKDTree, self).__init__(wrapped_data, leafsize)
+
+    # Ideally, KDTree and cKDTree would expose identical query and __query
+    # interfaces.  But they don't, and cKDTree.__query is also inaccessible
+    # from Python.  We do our best here to cope.
+    def __query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):
+        # This is the internal query method, which guarantees that x
+        # is a single point, not an array of points
+        #
+        # A slight complication: k could be "None", which means "return
+        # all neighbors within the given distance_upper_bound".
+
+        # Cap distance_upper_bound
+        distance_upper_bound = np.min([distance_upper_bound,
+                                      self.max_distance_upper_bound])
+
+        # Run queries over all relevant images of x
+        hits_list = []
+        for real_x in _gen_relevant_images(x, self.bounds, distance_upper_bound):
+            d, i = super(PeriodicCKDTree, self).query(
+                    real_x, k, eps, p, distance_upper_bound)
+            if k > 1:
+                hits_list.append(list(zip(d, i)))
+            else:
+                hits_list.append([(d, i)])
+
+        # Now merge results
+        if k > 1:
+            return heapq.nsmallest(k, itertools.chain(*hits_list))
+        elif k == 1:
+            return [min(itertools.chain(*hits_list))]
+        else:
+            raise ValueError("Invalid k in periodic_kdtree._KDTree__query")
+
+    def query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):
+        """
+        Query the kd-tree for nearest neighbors
+
+        Parameters
+        ----------
+        x : array_like, last dimension self.m
+            An array of points to query.
+        k : integer
+            The number of nearest neighbors to return.
+        eps : non-negative float
+            Return approximate nearest neighbors; the kth returned value
+            is guaranteed to be no further than (1+eps) times the
+            distance to the real k-th nearest neighbor.
+        p : float, 1<=p<=infinity
+            Which Minkowski p-norm to use.
+            1 is the sum-of-absolute-values "Manhattan" distance
+            2 is the usual Euclidean distance
+            infinity is the maximum-coordinate-difference distance
+        distance_upper_bound : nonnegative float
+            Return only neighbors within this distance.  This is used to prune
+            tree searches, so if you are doing a series of nearest-neighbor
+            queries, it may help to supply the distance to the nearest neighbor
+            of the most recent point.
+
+        Returns
+        -------
+        d : array of floats
+            The distances to the nearest neighbors.
+            If x has shape tuple+(self.m,), then d has shape tuple+(k,).
+            Missing neighbors are indicated with infinite distances.
+        i : ndarray of ints
+            The locations of the neighbors in self.data.
+            If `x` has shape tuple+(self.m,), then `i` has shape tuple+(k,).
+            Missing neighbors are indicated with self.n.
+
+        """
+        x = np.asarray(x)
+        if np.shape(x)[-1] != self.m:
+            raise ValueError("x must consist of vectors of length %d but has shape %s" % (self.m, np.shape(x)))
+        if p<1:
+            raise ValueError("Only p-norms with 1<=p<=infinity permitted")
+        retshape = np.shape(x)[:-1]
+        if retshape!=():
+            if k>1:
+                dd = np.empty(retshape+(k,),dtype=np.float)
+                dd.fill(np.inf)
+                ii = np.empty(retshape+(k,),dtype=np.int)
+                ii.fill(self.n)
+            elif k==1:
+                dd = np.empty(retshape,dtype=np.float)
+                dd.fill(np.inf)
+                ii = np.empty(retshape,dtype=np.int)
+                ii.fill(self.n)
+            else:
+                raise ValueError("Requested %s nearest neighbors; acceptable numbers are integers greater than or equal to one, or None")
+            for c in np.ndindex(retshape):
+                hits = self.__query(x[c], k=k, eps=eps, p=p, distance_upper_bound=distance_upper_bound)
+                if k>1:
+                    for j in range(len(hits)):
+                        dd[c+(j,)], ii[c+(j,)] = hits[j]
+                elif k==1:
+                    if len(hits)>0:
+                        dd[c], ii[c] = hits[0]
+                    else:
+                        dd[c] = np.inf
+                        ii[c] = self.n
+            return dd, ii
+        else:
+            hits = self.__query(x, k=k, eps=eps, p=p, distance_upper_bound=distance_upper_bound)
+            if k==1:
+                if len(hits)>0:
+                    return hits[0]
+                else:
+                    return np.inf, self.n
+            elif k>1:
+                dd = np.empty(k,dtype=np.float)
+                dd.fill(np.inf)
+                ii = np.empty(k,dtype=np.int)
+                ii.fill(self.n)
+                for j in range(len(hits)):
+                    dd[j], ii[j] = hits[j]
+                return dd, ii
+            else:
+                raise ValueError("Requested %s nearest neighbors; acceptable numbers are integers greater than or equal to one, or None")
+
+    # Ideally, KDTree and cKDTree would expose identical __query_ball_point
+    # interfaces.  But they don't, and cKDTree.__query_ball_point is also
+    # inaccessible from Python.  We do our best here to cope.
+    def __query_ball_point(self, x, r, p=2., eps=0):
+        # This is the internal query method, which guarantees that x
+        # is a single point, not an array of points
+
+        # Cap r
+        r = min(r, self.max_distance_upper_bound)
+
+        # Run queries over all relevant images of x
+        results = []
+        for real_x in _gen_relevant_images(x, self.bounds, r):
+            results.extend(super(PeriodicCKDTree, self).query_ball_point(
+                real_x, r, p, eps))
+        return results
+
+    def query_ball_point(self, x, r, p=2., eps=0):
+        """
+        Find all points within distance r of point(s) x.
+
+        Parameters
+        ----------
+        x : array_like, shape tuple + (self.m,)
+            The point or points to search for neighbors of.
+        r : positive float
+            The radius of points to return.
+        p : float, optional
+            Which Minkowski p-norm to use.  Should be in the range [1, inf].
+        eps : nonnegative float, optional
+            Approximate search. Branches of the tree are not explored if their
+            nearest points are further than ``r / (1 + eps)``, and branches are
+            added in bulk if their furthest points are nearer than
+            ``r * (1 + eps)``.
+
+        Returns
+        -------
+        results : list or array of lists
+            If `x` is a single point, returns a list of the indices of the
+            neighbors of `x`. If `x` is an array of points, returns an object
+            array of shape tuple containing lists of neighbors.
+
+        Notes
+        -----
+        If you have many points whose neighbors you want to find, you may
+        save substantial amounts of time by putting them in a
+        PeriodicCKDTree and using query_ball_tree.
+        """
+        x = np.asarray(x).astype(np.float)
+        if x.shape[-1] != self.m:
+            raise ValueError("Searching for a %d-dimensional point in a " \
+                             "%d-dimensional KDTree" % (x.shape[-1], self.m))
+        if len(x.shape) == 1:
+            return self.__query_ball_point(x, r, p, eps)
+        else:
+            retshape = x.shape[:-1]
+            result = np.empty(retshape, dtype=np.object)
+            for c in np.ndindex(retshape):
+                result[c] = self.__query_ball_point(x[c], r, p, eps)
+            return result
+
+    def query_ball_tree(self, other, r, p=2., eps=0):
+        raise NotImplementedError()
+
+    def query_pairs(self, r, p=2., eps=0):
+        raise NotImplementedError()
+
+    def count_neighbors(self, other, r, p=2.):
+        raise NotImplementedError()
+
+    def sparse_distance_matrix(self, other, max_distance, p=2.):
+        raise NotImplementedError()

python_src/vide/voidUtil/plotDefs.py

@@ -0,0 +1,33 @@
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/voidUtil/plotDefs.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+# 
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+LIGHT_SPEED = 299792.458 
+
+colorList = ['r', 'b', 'g', 'y', 'c', 'm', 
+             'darkred', 'grey',
+             'orange', 'darkblue',
+             'indigo', 'lightseagreen', 'maroon', 'olive',
+             'royalblue', 'palevioletred', 'seagreen', 'tomato',
+             'aquamarine', 'darkslateblue',
+             'khaki', 'lawngreen', 'mediumorchid',
+             'orangered', 'thistle'
+             'yellowgreen']
+
+linewidth = 4
+fontsize = 12

python_src/vide/voidUtil/plotUtil.py

@@ -0,0 +1,298 @@
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/voidUtil/plotUtil.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+# 
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+__all__=['plotNumberFunction','plotEllipDist','plotVoidCells']
+
+from vide.backend.classes import *
+from .plotDefs import *
+import numpy as np
+import os
+import pylab as plt
+import vide.backend.cosmologyTools as vp
+from vide.voidUtil import getArray, shiftPart, getVoidPart
+
+def fill_between(x, y1, y2=0, ax=None, **kwargs):
+    """Plot filled region between `y1` and `y2`.
+
+    This function works exactly the same as matplotlib's fill_between, except
+    that it also plots a proxy artist (specifically, a rectangle of 0 size)
+    so that it can be added it appears on a legend.
+    """
+    ax = ax if ax is not None else plt.gca()
+    ax.fill_between(x, y1, y2, interpolate=True, **kwargs)
+    p = plt.Rectangle((0, 0), 0, 0, **kwargs)
+    ax.add_patch(p)
+
+# -----------------------------------------------------------------------------
+def plotNumberFunction(catalogList, 
+                       figDir="./", 
+                       plotName="numberfunc",
+                       cumulative=True,
+                       binWidth=1):
+
+# plots a cumulative number function
+#   catalogList: list of void catalogs to plot
+#   figDir: output directory for figures
+#   plotName: name to prefix to all outputs
+#   cumulative: if True, plots cumulative number function
+#   binWidth: width of histogram bins in Mpc/h
+# returns:
+#   ellipDistList: array of len(catalogList), 
+#                   each element has array of size bins, number, +/- 1 sigma
+
+  print("Plotting number function")
+  
+  catalogList = np.atleast_1d(catalogList)
+
+  plt.clf()
+  plt.xlabel("$R_{eff}$ [$h^{-1}Mpc$]", fontsize=14)
+  
+  if cumulative:
+    plt.ylabel(r"log ($n$ (> R) [$h^3$ Gpc$^{-3}$])", fontsize=14)
+  else:
+    plt.ylabel(r"log ($dn/dR$ [$h^3$ Gpc$^{-3}$])", fontsize=14)
+ 
+  ellipDistList = [] 
+
+  for (iSample,catalog) in enumerate(catalogList):
+    sample = catalog.sampleInfo
+    data = getArray(catalog.voids, 'radius')
+  
+    if sample.dataType == "observation":
+      maskFile = sample.maskFile
+  
+      boxVol = vp.getSurveyProps(maskFile,
+                                 sample.zBoundary[0], sample.zBoundary[1],
+                                 sample.zRange[0], sample.zRange[1], "all",
+                                 selectionFuncFile=sample.selFunFile)[0]
+    else:
+      boxVol = sample.boxLen*sample.boxLen*(sample.zBoundaryMpc[1] -
+                                            sample.zBoundaryMpc[0])
+  
+    boxVol *= 1.e-9 # Mpc->Gpc
+  
+    bins = 100./binWidth
+    hist, binEdges = np.histogram(data, bins=bins, range=(0., 100.))
+    binCenters = 0.5*(binEdges[1:] + binEdges[:-1])
+
+    if cumulative:
+      foundStart = False
+      for iBin in range(len(hist)):
+        if not foundStart and hist[iBin] == 0:
+          continue
+        foundStart = True
+        hist[iBin] = np.sum(hist[iBin:])
+  
+    nvoids = len(data)
+    var = hist * (1. - hist/nvoids)
+    sig = np.sqrt(var)
+  
+    lowerbound = hist - sig
+    upperbound = hist + sig
+  
+    mean = np.log10(hist/boxVol)
+    lowerbound = np.log10(lowerbound/boxVol)
+    upperbound = np.log10(upperbound/boxVol)
+  
+    lineColor = colorList[iSample]
+    lineTitle = sample.fullName
+    
+    trim = (lowerbound > .01)
+    mean = mean[trim]
+    binCentersToUse = binCenters[trim]
+    lower = lowerbound[trim]
+    upper = upperbound[trim]
+  
+    alpha = 0.55
+    fill_between(binCentersToUse, lower, upper,
+               label=lineTitle, color=lineColor,
+               alpha=alpha,
+               )
+  
+    lineStyle = '-'
+    plt.plot(binCentersToUse, mean, lineStyle,
+                color=lineColor,
+                linewidth=3)
+    
+    ellipDistList.append((binCentersToUse, mean, lower, upper))
+  
+  plt.legend(loc = "upper right", fancybox=True, prop={'size':14})
+  
+  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")
+
+  return ellipDistList 
+
+
+# -----------------------------------------------------------------------------
+def plotEllipDist(catalogList, 
+                  figDir="./", 
+                  plotName="ellipdist"):
+
+# plots ellipticity distributions
+#   catalogList: list of void catalogs to plot
+#   figDir: output directory for figures
+#   plotName: name to prefix to all outputs
+# returns:
+#   ellipDistList: array of len(catalogList), 
+#                  each element has array of ellipticity distributions
+
+  print("Plotting ellipticity distributions")
+  
+  plt.clf()
+  plt.xlabel(r"Ellipticity $\epsilon$", fontsize=14)
+  plt.ylabel(r"P($\epsilon$)", fontsize=14)
+  
+  ellipDistList = [] 
+
+  catalogList = np.atleast_1d(catalogList)
+
+  for (iSample,catalog) in enumerate(catalogList):
+    sample = catalog.sampleInfo
+    data = getArray(catalog.voids, 'ellipticity')
+  
+    dataWeights = np.ones_like(data)/len(data)
+    dataHist, dataBins = np.histogram(data, bins=10, weights=dataWeights,
+                                         range=(0.0,0.35)) 
+    binCenters = 0.5*(dataBins[1:] + dataBins[:-1])
+  
+    plt.plot(binCenters, dataHist, label=sample.fullName, 
+             color=colorList[iSample])
+  
+    ellipDistList.append((dataBins, dataHist,))
+    
+  plt.legend(loc = "upper right", fancybox=True, prop={'size':14})
+  
+  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")
+  
+  return
+
+
+# -----------------------------------------------------------------------------
+def plotVoidCells(catalog,
+              voidID,
+              figDir="./",
+              plotName="cells",
+              plotDensity=True,
+              sliceWidth=250):
+
+# plots the particles belonging to a void
+#   catalog: void catalog
+#   voidID: ID of void to plot
+#   figDir: output directory for figures
+#   plotName: name to prefix to all outputs
+#   sliceWidth: width of slice in Mpc/h
+
+  sample = catalog.sampleInfo
+  periodicLine = getPeriodic(sample)
+
+  plt.clf()
+
+  iVoid = -1
+  for i in range(len(catalog.voids)):
+    if catalog.voids[i].voidID == voidID:
+      iVoid = i
+      break
+
+  if iVoid == -1:
+    print("Void ID %d not found!" % voidID)
+    return
+  
+  sliceCenter = catalog.voids[iVoid].macrocenter
+
+  xwidth = sliceWidth
+  ywidth = sliceWidth
+  zwidth = max(sliceWidth/4., 50)
+
+  xmin = -xwidth/2.
+  xmax = xwidth/2.
+  ymin = -ywidth/2.
+  ymax = ywidth/2.
+  zmin = -zwidth/2.
+  zmax = zwidth/2.
+
+  #Slice Particles
+  if plotDensity:
+    part = catalog.partPos
+    part = shiftPart(part, sliceCenter, periodicLine, catalog.ranges) 
+
+    filter = (part[:,0] > xmin) & (part[:,0] < xmax) & \
+             (part[:,1] > ymin) & (part[:,1] < ymax) & \
+             (part[:,2] > zmin) & (part[:,2] < zmax)
+    part = part[filter]
+
+    extent = [xmin, xmax, ymin, ymax]
+    hist, xedges, yedges = np.histogram2d(part[:,0], part[:,1], normed=False,
+                                          bins=64)
+
+    hist = np.log10(hist+1)
+    plt.imshow(hist,
+               aspect='equal',
+               extent=extent,
+               interpolation='gaussian',
+               cmap='YlGnBu_r')
+
+
+  # overlay voids as circles
+  fig = plt.gcf()
+
+  voidPart = getVoidPart(catalog, voidID)
+
+  newpart = np.zeros((len(voidPart),3))
+  volume=np.zeros(len(voidPart))
+  radius=np.zeros(len(voidPart))
+
+  newpart[:,0] = getArray(voidPart, 'x')
+  newpart[:,1] = getArray(voidPart, 'y')
+  newpart[:,2] = getArray(voidPart, 'z')
+
+  volume = getArray(voidPart, 'volume')
+  radius = (3.*volume/(4.*np.pi))**(1./3.)
+
+  shiftedPartVoid =shiftPart(newpart,sliceCenter, periodicLine, catalog.ranges)
+
+  #Limiting plotted cells to cells into the slice
+  #Possibility to only plot bigger cells (through cellsradiuslim)
+  cellsMinlimz = zmin
+  cellsMaxlimz = zmax
+  cellsradiuslim = 0.0
+
+  for p in range(len(volume)):
+    if (shiftedPartVoid[p,2]>(cellsMinlimz) and \
+        shiftedPartVoid[p,2]<(cellsMaxlimz) and \
+        radius[p]>cellsradiuslim):
+      color = 'blue'
+      circle = plt.Circle((shiftedPartVoid[p,0], \
+                           shiftedPartVoid[p,1]), \
+                           radius[p],
+                           alpha =.2, fc=color,edgecolor=None,linewidth=1)
+      fig.gca().add_artist(circle)
+
+  title="cells"+str(voidID)
+  plt.title(title, fontsize=20)
+  plotName="cells"+str(voidID)
+
+  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")
+
+  return

python_src/vide/voidUtil/profileUtil.py

@@ -0,0 +1,213 @@
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/voidUtil/profileUtil.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+# 
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+__all__=['buildProfile','fitHSWProfile','getHSWProfile',]
+
+from vide.backend.classes import *
+from vide.voidUtil import *
+from .plotDefs import *
+import numpy as np
+import os
+import vide.apTools as vp
+from scipy.optimize import curve_fit
+from scipy.interpolate import interp1d
+
+def HSWProfile(r, rs, dc):
+  alpha = -2.0*rs + 4.0
+  if rs < 0.91:
+    beta = 17.5*rs - 6.5
+  else:
+    beta = -9.8*rs + 18.4
+  return dc * (1 - (r/rs)**alpha) / (1+ (r)**beta) + 1
+
+# -----------------------------------------------------------------------------
+def buildProfile(catalog, rMin, rMax, nBins=10):
+
+# builds a stacked radial density profile from the given catalog
+#   catalog: void catalog
+#   rMin: minimum void radius, in Mpc/h
+#   rMax: maximum void radius, in Mpc/h
+#   nBins: number of bins in profile (detaulf 10)
+#
+# returns:
+#   binCenters: array of radii in binned profile
+#   stackedProfile: the stacked density profile
+#   sigmas: 1-sigma uncertainty in each bin
+
+  rMaxProfile = rMin*3 + 2
+  periodicLine = getPeriodic(catalog.sampleInfo)
+
+  print("  Building particle tree...")
+  partTree = getPartTree(catalog)
+
+  print("  Selecting voids to stack...")
+  voidsToStack = [v for v in catalog.voids if (v.radius > rMin and v.radius < rMax)]
+
+  if len(voidsToStack) == 0:
+    print("  No voids to stack!")
+    return -1, -1, -1
+    
+  print("  Stacking voids...")
+  allProfiles = []
+  for void in voidsToStack:
+    center = void.macrocenter
+    
+    localPart = catalog.partPos[ getBall(partTree, center, rMaxProfile) ]
+    shiftedPart = shiftPart(localPart, center, periodicLine, catalog.ranges)
+
+    dist = np.sqrt(np.sum(shiftedPart[:,:]**2, axis=1))
+    thisProfile, radii = np.histogram(dist, bins=nBins, range=(0,rMaxProfile))
+    deltaV = 4*np.pi/3*(radii[1:]**3-radii[0:(radii.size-1)]**3)
+    thisProfile = np.float32(thisProfile)
+    thisProfile /= deltaV
+    thisProfile /= catalog.volNorm
+    allProfiles.append(thisProfile)
+    binCenters = 0.5*(radii[1:] + radii[:-1])
+
+  nVoids = len(voidsToStack)
+  stackedProfile = np.mean(allProfiles, axis=0) 
+  sigmas = np.std(allProfiles, axis=0) / np.sqrt(nVoids)
+  
+  return binCenters, stackedProfile, sigmas
+
+
+# -----------------------------------------------------------------------------
+def fitHSWProfile(radii, densities, sigmas, rV):
+
+# fits the given density profile to the HSW function
+#   radii: array of radii
+#   densities: array of densities in units of mean density
+#   sigmas: array of uncertainties
+#   rV: mean effective void radius
+# 
+# returns:
+#   popt: best-fit values of rs and dc
+#   pcov: covariance matrix
+#   rVals: array of radii for best-fit curve
+#   hswProfile: array of densities for best-fit curve in units of mean density
+
+  popt, pcov = curve_fit(HSWProfile, radii/rV, densities,
+                         sigma=sigmas,
+                         maxfev=10000, xtol=5.e-3,
+                         p0=[1.0,-1.0])
+
+  # return best-fits
+  rVals = np.linspace(0.0, radii[-1], 100) / rV
+  return popt, pcov, rVals*rV, HSWProfile(rVals,popt[0],popt[1])
+
+# -----------------------------------------------------------------------------
+def getHSWProfile(density, radius):
+
+# returns the HSW profile for the given sample density and void size
+#   will interpolate/extrapole the radius
+
+#  density: choice of sample (see arXiv:1309.5087):
+#           maxDM: DM at 1 particles per cubic Mpc/h
+#           fullDM: DM at 0.01 particles per cubic Mpc/h
+#           denseDM: DM at 4.e-3 particles per cubic Mpc/h
+#           sparseDM: DM at 3.e-4 particles per cubic Mpc/h
+#           
+#           denseHalos: halos at 4.e-3 particles per cubic Mpc/h
+#           sparseHalos: halos at 3.e-4 particles per cubic Mpc/h
+#          
+#           denseGal: galaxies at 4.e-3 particles per cubic Mpc/h
+#           sparseGal: galaxies at 3.e-4 particles per cubic Mpc/h
+#
+#  radius: void size in Mpc/h
+
+# returns:
+#   (rs, dc): best-fit values
+#   binCenters: array of radii in binned profile
+#   stackedProfile: the density profile
+
+  samples = [
+    {'name': 'maxDM',
+     'rv': [8.74041095,  11.65557095,  15.54746657,  20.94913774, 28.41063131,  38.61523696,  51.85944898,  69.42297033],
+     'rs': [0.74958738,  0.79650829,  0.86245251,  0.93960051, 1.01595177,  1.13159483,  1.31457096,  1.65611709],
+     'dc': [-0.95353184, -0.94861939, -0.91888455, -0.84480086, -0.73431544, -0.62614422, -0.54908132, -0.4912146],
+    },
+    {'name': 'fullDM',
+     'rv': [10, 15, 20, 25, 30, 35],
+     'rs': [ 0.76986779,  0.80980775,  0.86590177,  0.93732629,  1.02643542,
+        1.12875503],
+     'dc': [-0.81021   , -0.78115474, -0.78326026, -0.78670109, -0.76626508,
+       -0.72531084],
+    },
+    {'name': 'denseDM',
+     'rv': [10, 15, 20, 25, 30, 35, 40],
+     'rs': [0.7477462 ,  0.79932797,  0.84369297,  0.90309363,  0.92990401,
+        1.06970842,  1.16393474],
+     'dc': [-0.78333107, -0.75780925, -0.71586397, -0.74669512, -0.74902649,
+       -0.75342964, -0.76598043],
+    },
+    {'name': 'sparseDM',
+     'rv': [25, 30, 35, 40, 45, 50, 55, 60],
+     'rs': [0.78351117,  0.79751047,  0.8225573 ,  0.83751894,  0.85443167,
+        0.86346031,  0.85692501,  0.91470448],
+     'dc': [-0.67407548, -0.62389586, -0.59125414, -0.55341724, -0.54659457,
+       -0.5297821 , -0.534683  , -0.51055946],
+    },
+    {'name': 'denseHalos',
+     'rv': [15, 20, 25, 30, 35, 40, 45, 50, 55],
+     'rs': [0.75393528,  0.7758442 ,  0.79720886,  0.81560553,  0.83797299,
+        0.84377082,  0.84900783,  0.8709897 ,  0.86886687],
+     'dc': [-0.81348968, -0.77362777, -0.74336192, -0.72571135, -0.67928029,
+       -0.6279349 , -0.6313316 , -0.55188564, -0.48096026],
+    },
+    {'name': 'sparseHalos',
+     'rv': [25, 30, 35, 40, 45, 50, 55, 60, 65, 70],
+     'rs': [0.76420703,  0.78939067,  0.80480265,  0.82315275,  0.8158607 ,
+        0.82553517,  0.83843323,  0.85759226,  0.8471268 ,  0.89286939],
+     'dc': [-0.79317192, -0.81661415, -0.76770778, -0.7151494 , -0.718561  ,
+       -0.70858856, -0.68995608, -0.67415305, -0.63706798, -0.5303759],
+    },
+    {'name': 'denseGal',
+     'rv': [10, 15, 20, 25, 30, 35, 40, 45],
+     'rs': [0.70048139,  0.73717884,  0.75338516,  0.76782043,  0.79292536,
+        0.80157122,  0.8207239 ,  0.8091386],
+     'dc': [-0.83166549, -0.86505329, -0.81066899, -0.7662453 , -0.72840363,
+       -0.65163607, -0.57937656, -0.57125164],
+    },
+    {'name': 'sparseGal',
+     'rv': [25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75],
+     'rs': [0.75928922,  0.76622648,  0.77695425,  0.79963152,  0.8045125 ,
+        0.81892965,  0.83439691,  0.86600085,  0.83166875,  0.85283258,
+        0.83971344],
+     'dc': [-0.82413212, -0.87483536, -0.8221596 , -0.78459706, -0.75290061,
+       -0.77513988, -0.70012913, -0.67487994, -0.69903308, -0.65811992,
+       -0.57929526],
+    },
+  ]
+
+  mySample = next((item for item in samples if item['name'] == density), None)
+  if mySample == None:
+    print("Sample", density," not found! Use one of ", [item['name'] for item in samples])
+    return 
+
+  # interpolate the radii
+  rsFunc = interp1d( mySample['rv'], mySample['rs'], kind='cubic' )
+  dcFunc = interp1d( mySample['rv'], mySample['dc'], kind='cubic' )
+    
+  rs = rsFunc(radius)
+  dc = dcFunc(radius)
+
+  # return best-fits
+  rVals = np.linspace(0.0, 3*radius, 100) / radius
+  return (rs,dc), rVals*radius, HSWProfile(rVals,rs,dc)
+  

python_src/vide/voidUtil/xcorUtil.py

@@ -0,0 +1,147 @@
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/voidUtil/xcorUtil.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+# 
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+from matplotlib import rc
+from . import xcorlib
+from vide.voidUtil import getArray
+
+def computeXcor(catalog,
+                figDir="./",
+                Nmesh = 256,
+                Nbin = 100
+                ):
+
+# Computes and plots void-void and void-matter(galaxy) correlations
+#   catalog: catalog to analyze
+#   figDir: where to place plots
+#   Nmesh: number of grid cells in cic mesh-interpolation
+#   Nbin:  number of bins in final plots
+
+  # Parameters
+  Lbox = catalog.boxLen[0] # Boxlength
+  Lboxcut = 0.
+  Lbox -= 2*Lboxcut
+  
+  # Input particle arrays of shape (N,3)
+  xm = catalog.partPos # Halos / Galaxies / Dark matter
+  xv = getArray(catalog.voids, 'macrocenter')
+  
+  
+  # Interpolate to mesh
+  dm, wm, ws = xcorlib.cic(xm, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh, weights = None)
+  dv, wm, ws = xcorlib.cic(xv, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh, weights = None)
+  
+  # Fourier transform
+  dmk = np.fft.rfftn(dm)
+  dvk = np.fft.rfftn(dv)
+  
+  # 1D Power spectra & correlation functions
+  ((Nm, km, Pmm, SPmm),(Nmx, rm, Xmm, SXmm)) = xcorlib.powcor(dmk, dmk, Lbox, Nbin, 'lin', True, True, 1)
+  ((Nm, km, Pvm, SPvm),(Nmx, rm, Xvm, SXvm)) = xcorlib.powcor(dvk, dmk, Lbox, Nbin, 'lin', True, True, 1)
+  ((Nm, km, Pvv, SPvv),(Nmx, rm, Xvv, SXvv)) = xcorlib.powcor(dvk, dvk, Lbox, Nbin, 'lin', True, True, 1)
+  
+  # Number densities
+  nm = np.empty(len(km))
+  nv = np.empty(len(km))
+  nm[:] = len(xm)/Lbox**3
+  nv[:] = len(xv)/Lbox**3
+  
+  
+  # Plots
+  mpl.rc('font', family='serif')
+  ms = 2.5
+  fs = 16
+  mew = 0.1
+  margin = 1.2
+  kmin = km.min()/margin
+  kmax = km.max()*margin
+  rmin = rm.min()/margin
+  rmax = rm.max()*margin
+  
+  
+  # Density fields (projected)
+  plt.imshow(np.sum(dm[:,:,:]+1,2),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.savefig(figDir+'/dm.eps', bbox_inches="tight")
+  plt.savefig(figDir+'/dm.pdf', bbox_inches="tight")
+  plt.savefig(figDir+'/dm.png', 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.savefig(figDir+'/dv.eps', bbox_inches="tight") #, dpi=300
+  plt.savefig(figDir+'/dv.pdf', bbox_inches="tight") #, dpi=300
+  plt.savefig(figDir+'/dv.png', bbox_inches="tight") #, dpi=300
+  plt.clf()
+  
+  
+  # Power spectra & correlation functions
+  pa ,= plt.plot(km, Pmm, 'k-o', ms=0.8*ms, mew=mew, mec='k')
+  #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)
+  pb ,= plt.plot(km, Pvm, 'm-D', ms=ms, mew=mew, mec='k')
+  plt.plot(km, -Pvm, 'mD', ms=ms, mew=mew, mec='k')
+  plt.fill_between(km, abs(Pvm+SPvm), abs(Pvm-SPvm), color='m', alpha=0.2)
+  pc ,= plt.plot(km, Pvv, 'b-p', ms=1.3*ms, mew=mew, mec='k')
+  #plt.plot(km, Pvv-1./nv, 'b--', ms=ms, mew=mew)
+  plt.fill_between(km, Pvv+SPvv, abs(Pvv-SPvv), color='b', alpha=0.2)
+  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(Pvm[1:]).min()))/margin, max(10**np.ceil(np.log10(Pmm.max())),10**np.ceil(np.log10(Pvv.max())))*margin)
+  plt.legend([pa, pb, pc],['tt', 'vt', 'vv'],'best',prop={'size':12})
+  plt.savefig(figDir+'/power.eps', bbox_inches="tight") 
+  plt.savefig(figDir+'/power.pdf', bbox_inches="tight") 
+  plt.savefig(figDir+'/power.png', bbox_inches="tight") 
+  plt.clf()
+  
+  pa ,= plt.plot(rm, Xmm, 'k-o', ms=0.8*ms, mew=mew)
+  plt.fill_between(rm, abs(Xmm+SXmm), abs(Xmm-SXmm), color='k', alpha=0.2)
+  pb ,= plt.plot(rm, Xvm, 'm-D', ms=ms, mew=mew)
+  plt.plot(rm, -Xvm, 'mD', ms=ms, mew=mew)
+  plt.fill_between(rm, abs(Xvm+SXvm), abs(Xvm-SXvm), color='m', alpha=0.2)
+  pc ,= plt.plot(rm, Xvv, 'b-p', ms=1.3*ms, mew=mew)
+  plt.plot(rm, -Xvv, 'bp', ms=ms, mew=1.3*mew)
+  plt.fill_between(rm, abs(Xvv+SXvv), abs(Xvv-SXvv), color='b', alpha=0.2)
+  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(min(10**np.floor(np.log10(abs(Xvm).min())),10**np.floor(np.log10(abs(Xmm).min())))/margin, max(10**np.ceil(np.log10(Xmm.max())),10**np.ceil(np.log10(Xvv.max())))*margin)
+  plt.legend([pa, pb, pc],['tt', 'vt', 'vv'],'best',prop={'size':12})
+  plt.savefig(figDir+'/correlation.eps', bbox_inches="tight") 
+  plt.savefig(figDir+'/correlation.pdf', bbox_inches="tight")
+  plt.savefig(figDir+'/correlation.png', bbox_inches="tight")
+  plt.clf()
+  
+  
+  return 

python_src/vide/voidUtil/xcorlib.py

@@ -0,0 +1,107 @@
+#+
+#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/voidUtil/xcorlib.py
+#   Copyright (C) 2010-2014 Guilhem Lavaux
+#   Copyright (C) 2011-2014 P. M. Sutter
+#
+#   This program is free software; you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation; version 2 of the License.
+# 
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License along
+#   with this program; if not, write to the Free Software Foundation, Inc.,
+#   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#+
+import numpy as np
+
+# CIC interpolation
+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)
+
+
+# Power spectra & correlation functions
+def powcor(d1, d2, Lbox, Nbin = 10, scale = 'lin', cor = False, cic = True, dim = 1):
+
+ Nmesh = len(d1)
+
+ # CIC correction
+ if cic:
+  wid = np.indices(np.shape(d1))
+  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
+ if cic: Pk /= wcic**2
+ 
+ (Nm, km, Pkm, SPkm) = shellavg(np.real(Pk), dk, Nmesh, Nbin = Nbin, xmin = 0., xmax = Nmesh*dk/2, scale = scale, dim = dim)
+ 
+ # Inverse Fourier transform and shell average correlation function
+ if cor:
+  if cic: Pk *= wcic**2 # Undo cic-correction in correlation function
+  dx = Lbox/Nmesh
+  Xr = np.fft.irfftn(Pk)*(Nmesh/Lbox)**3
+
+  (Nmx, rm, Xrm, SXrm) = shellavg(np.real(Xr), dx, Nmesh, Nbin = Nbin/2, xmin = dx, xmax = 140., scale = scale, dim = dim)
+  
+  return ((Nm, km, Pkm, SPkm),(Nmx, rm, Xrm, SXrm))
+ 
+ else: return (Nm, km, Pkm, SPkm)
+
+
+# Shell averaging
+def shellavg(f, dx, Nmesh, Nbin = 10, xmin = 0., xmax = 1., scale = 'lin', dim = 1):
+ 
+ x = np.indices(np.shape(f))
+ x[np.where(x >= Nmesh/2)] -= Nmesh
+ f = f.flatten()
+ 
+ if scale == 'lin': bins = xmin+(xmax-xmin)* np.linspace(0,1,Nbin+1)
+ if scale == 'log': bins = xmin*(xmax/xmin)**np.linspace(0,1,Nbin+1)
+
+ if dim == 1: # 1D
+  x = dx*np.sqrt(np.sum(x**2,0)).flatten()
+  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)
+
+ elif dim == 2: # 2D
+  xper = dx*np.sqrt(x[0,:,:,:]**2 + x[1,:,:,:]**2 + 1e-100).flatten()
+  xpar = dx*np.abs(x[2,:,:,:]).flatten()
+  x = dx*np.sqrt(np.sum(x**2,0)).flatten()
+  Nm = np.histogram2d(xper, xpar, bins = [bins,bins])[0]
+  xmper = np.histogram2d(xper, xpar, bins = [bins,bins], weights = xper)[0]/Nm
+  xmpar = np.histogram2d(xper, xpar, bins = [bins,bins], weights = xpar)[0]/Nm
+  fm = np.histogram2d(xper, xpar, bins = [bins,bins], weights = f)[0]/Nm
+  return (Nm, xmper, xmpar, fm)