[Yt-svn] yt-commit r1416 - trunk/yt/extensions
britton at wrangler.dreamhost.com
britton at wrangler.dreamhost.com
Mon Aug 31 10:40:57 PDT 2009
Author: britton
Date: Mon Aug 31 10:40:56 2009
New Revision: 1416
URL: http://yt.spacepope.org/changeset/1416
Log:
Committing new halo profiler and virial filter function.
Added:
trunk/yt/extensions/HaloFilters.py
Modified:
trunk/yt/extensions/HaloProfiler.py
Added: trunk/yt/extensions/HaloFilters.py
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/HaloFilters.py Mon Aug 31 10:40:56 2009
@@ -0,0 +1,115 @@
+"""
+Halo filters to be used with the HaloProfiler.
+
+Author: Britton Smith <brittons at origins.colorado.edu>
+Affiliation: CASA/University of CO, Boulder
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008-2009 Britton Smith. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt 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; either version 3 of the License, or
+ (at your option) any later version.
+
+ 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, see <http://www.gnu.org/licenses/>.
+"""
+
+from yt.logger import lagosLogger as mylog
+from copy import deepcopy
+import numpy as na
+
+def VirialFilter(profile,overdensity_field='ActualOverdensity',
+ virial_overdensity=200.,must_be_virialized=True,
+ virial_filters=[['TotalMassMsun','>=','1e14']],
+ virial_quantities=['TotalMassMsun','RadiusMpc'],
+ virial_index=None):
+ """
+ Filter halos by virial quantities.
+ Return values are a True or False whether the halo passed the filter,
+ along with a dictionary of virial quantities for the fields specified in
+ the virial_quantities keyword. Thresholds for virial quantities are
+ given with the virial_filters keyword in the following way:
+ [field, condition, value].
+ """
+
+ fields = deepcopy(virial_quantities)
+ if virial_filters is None: virial_filters = []
+ for vfilter in virial_filters:
+ if not vfilter[0] in fields:
+ fields.append(vfilter[0])
+
+ overDensity = []
+ temp_profile = {}
+ for field in fields:
+ temp_profile[field] = []
+
+ for q in range(len(profile[overdensity_field])):
+ good = True
+ if (profile[overdensity_field][q] != profile[overdensity_field][q]):
+ good = False
+ continue
+ for field in fields:
+ if (profile[field][q] != profile[field][q]):
+ good = False
+ break
+ if good:
+ overDensity.append(profile[overdensity_field][q])
+ for field in fields:
+ temp_profile[field].append(profile[field][q])
+
+ virial = {}
+ for field in fields:
+ virial[field] = 0.0
+
+ if (not (na.array(overDensity) >= virial_overdensity).any()) and \
+ must_be_virialized:
+ mylog.error("This halo is not virialized!")
+ return [False, {}]
+
+ if (len(overDensity) < 2):
+ mylog.error("Skipping halo with no valid points in profile.")
+ return [False, {}]
+
+ if (overDensity[1] <= virial_overdensity):
+ index = 0
+ elif (overDensity[-1] >= virial_overdensity):
+ index = -2
+ else:
+ for q in (na.arange(len(overDensity)-2))+2:
+ if (overDensity[q] < virial_overdensity):
+ index = q - 1
+ break
+
+ if type(virial_index) is list:
+ virial_index.append(index)
+
+ for field in fields:
+ if (overDensity[index+1] - overDensity[index]) == 0:
+ mylog.error("Overdensity profile has slope of zero.")
+ return [False, {}]
+ else:
+ slope = (temp_profile[field][index+1] - temp_profile[field][index]) / \
+ (overDensity[index+1] - overDensity[index])
+ value = slope * (virial_overdensity - overDensity[index]) + \
+ temp_profile[field][index]
+ virial[field] = value
+
+ for vfilter in virial_filters:
+ if eval("%s %s %s" % (virial[vfilter[0]],vfilter[1],vfilter[2])):
+ mylog.debug("(%s %s %s) returned True for %s." % (vfilter[0],vfilter[1],vfilter[2],virial[vfilter[0]]))
+ continue
+ else:
+ mylog.debug("(%s %s %s) returned False for %s." % (vfilter[0],vfilter[1],vfilter[2],virial[vfilter[0]]))
+ return [False, {}]
+
+ return [True, dict((q,virial[q]) for q in virial_quantities)]
+
Modified: trunk/yt/extensions/HaloProfiler.py
==============================================================================
--- trunk/yt/extensions/HaloProfiler.py (original)
+++ trunk/yt/extensions/HaloProfiler.py Mon Aug 31 10:40:56 2009
@@ -27,21 +27,46 @@
from yt.lagos.HaloFinding import HaloFinder
from yt.logger import lagosLogger as mylog
import yt.raven as raven
+from HaloFilters import *
import numpy as na
import os
import h5py
+import types
PROFILE_RADIUS_THRESHOLD = 2
class HaloProfiler(lagos.ParallelAnalysisInterface):
- def __init__(self,dataset,HaloProfilerParameterFile,halos='multiple',radius=0.1,radius_units='1',hop_style='new'):
+ def __init__(self, dataset, halos='multiple', halo_list_file='HopAnalysis.out', halo_list_format='yt_hop',
+ halo_finder_function=HaloFinder, halo_finder_args=None, halo_finder_kwargs=None,
+ use_density_center=False, density_center_exponent=1.0, use_field_max_center=None,
+ halo_radius=0.1, radius_units='1', n_profile_bins=50,
+ profile_output_dir='radial_profiles', projection_output_dir='projections',
+ projection_width=8.0, projection_width_units='mpc', project_at_level='max',
+ velocity_center=['bulk', 'halo']):
+
self.dataset = dataset
- self.HaloProfilerParameterFile = HaloProfilerParameterFile
- self.haloProfilerParameters = {}
- self.profileFields = {}
- self.projectionFields = {}
- self.hopHalos = []
- self.virialQuantities = []
+
+ self.profile_output_dir = profile_output_dir
+ self.projection_output_dir = projection_output_dir
+ self.n_profile_bins = n_profile_bins
+ self.projection_width = projection_width
+ self.projection_width_units = projection_width_units
+ self.project_at_level = project_at_level
+
+ self.profile_fields = []
+ self.projection_fields = []
+
+ self._halo_filters = []
+ self.all_halos = []
+ self.filtered_halos = []
+ self._projection_halo_list = []
+
+ # Set halo finder function and parameters, if needed.
+ self.halo_finder_function = halo_finder_function
+ self.halo_finder_args = halo_finder_args
+ if self.halo_finder_args is None: self.halo_finder_args = ()
+ self.halo_finder_kwargs = halo_finder_kwargs
+ if self.halo_finder_kwargs is None: self.halo_finder_kwargs = {}
# Set option to get halos from hop or single halo at density maximum.
# multiple: get halos from hop
@@ -51,43 +76,62 @@
mylog.error("Keyword, halos, must be either 'single' or 'multiple'.")
return None
- # Set hop file style.
- # old: enzo_hop output.
- # new: yt hop output.
- self.hop_style = hop_style
- if not(self.hop_style is 'old' or self.hop_style is 'new'):
- mylog.error("Keyword, hop_style, must be either 'old' or 'new'.")
+ # Set halo list format.
+ # 'yt_hop': yt hop output.
+ # 'enzo_hop': enzo_hop output.
+ # dictionary: a dictionary containing fields and their corresponding columns.
+ self.halo_list_file = halo_list_file
+ if halo_list_format == 'yt_hop':
+ self.halo_list_format = {'id':0, 'mass':1, 'center':[7, 8, 9], 'velocity':[10, 11, 12], 'r_max':13}
+ elif halo_list_format == 'enzo_hop':
+ self.halo_list_format = {'id':0, 'center':[4, 5, 6]}
+ elif isinstance(halo_list_format, types.DictType):
+ self.halo_list_format = halo_list_format
+ else:
+ mylog.error("Keyword, halo_list_format, must be 'yt_hop', 'enzo_hop', or a dictionary of custom settings.")
return None
- # Set some parameter defaults.
- self._SetParameterDefaults()
-
- # Read parameter file.
- self._ReadHaloProfilerParameterFile()
+ # Option to recenter sphere on density center.
+ self.use_density_center = use_density_center
+ self.density_center_exponent = density_center_exponent
+ if self.use_density_center:
+ def _MatterDensityXTotalMass(field, data):
+ return na.power((data['Matter_Density'] * data['TotalMassMsun']), self.density_center_exponent)
+ def _Convert_MatterDensityXTotalMass(data):
+ return 1
+ lagos.add_field("MatterDensityXTotalMass", units=r"",
+ function=_MatterDensityXTotalMass,
+ convert_function=_Convert_MatterDensityXTotalMass)
+
+ # Option to recenter sphere on the location of a field max.
+ self.use_field_max_center = use_field_max_center
+ if self.use_field_max_center is not None:
+ self.use_density_center = False
# Look for any field that might need to have the bulk velocity set.
- self.needBulkVelocity = False
- for field in self.profileFields:
+ self._need_bulk_velocity = False
+ for field in [hp['field'] for hp in self.profile_fields]:
if 'Velocity' in field or 'Mach' in field:
- self.needBulkVelocity = True
+ self._need_bulk_velocity = True
break
# Check validity for VelocityCenter parameter which toggles how the
# velocity is zeroed out for radial velocity profiles.
- if self.haloProfilerParameters['VelocityCenter'][0] == 'bulk':
- if self.haloProfilerParameters['VelocityCenter'][1] == 'halo' and \
+ self.velocity_center = velocity_center[:]
+ if self.velocity_center[0] == 'bulk':
+ if self.velocity_center[1] == 'halo' and \
self.halos is 'single':
mylog.error("Parameter, VelocityCenter, must be set to 'bulk sphere' or 'max <field>' with halos flag set to 'single'.")
return None
- if self.haloProfilerParameters['VelocityCenter'][1] == 'halo' and \
- self.hop_style is 'old':
+ if self.velocity_center[1] == 'halo' and \
+ self.halo_list_format is 'enzo_hop':
mylog.error("Parameter, VelocityCenter, must be 'bulk sphere' for old style hop output files.")
return None
- if not(self.haloProfilerParameters['VelocityCenter'][1] == 'halo' or
- self.haloProfilerParameters['VelocityCenter'][1] == 'sphere'):
+ if not(self.velocity_center[1] == 'halo' or
+ self.velocity_center[1] == 'sphere'):
mylog.error("Second value of VelocityCenter must be either 'halo' or 'sphere' if first value is 'bulk'.")
return None
- elif self.haloProfilerParameters['VelocityCenter'][0] == 'max':
+ elif self.velocity_center[0] == 'max':
if self.halos is 'multiple':
mylog.error("Getting velocity center from a max field value only works with halos='single'.")
return None
@@ -98,145 +142,250 @@
# Create dataset object.
self.pf = lagos.EnzoStaticOutput(self.dataset)
self.pf.h
- if self.halos is 'single' or hop_style is 'old':
- self.haloRadius = radius / self.pf[radius_units]
-
- def makeProfiles(self):
- "Make radial profiles for all halos on the list."
+ if self.halos is 'single' or not 'r_max' in self.halo_list_format:
+ self.halo_radius = halo_radius / self.pf[radius_units]
# Get halo(s).
if self.halos is 'single':
v, center = self.pf.h.find_max('Density')
singleHalo = {}
singleHalo['center'] = center
- singleHalo['r_max'] = self.haloRadius * self.pf.units['mpc']
+ singleHalo['r_max'] = self.halo_radius * self.pf.units['mpc']
singleHalo['id'] = 0
- self.hopHalos.append(singleHalo)
+ self.all_halos.append(singleHalo)
elif self.halos is 'multiple':
# Get hop data.
- self._LoadHopData()
+ self._load_halo_data()
+ if len(self.all_halos) == 0:
+ mylog.error("No halos loaded, there will be nothing to do.")
+ return None
else:
mylog.error("I don't know whether to get halos from hop or from density maximum. This should not have happened.")
- return
+ return None
+
+ def add_halo_filter(self, function, *args, **kwargs):
+ "Add a halo filter to the filter list."
+
+ self._halo_filters.append({'function':function, 'args':args, 'kwargs':kwargs})
+
+ def add_profile(self, field, weight_field=None, accumulation=False):
+ "Add a field for profiling."
+
+ self.profile_fields.append({'field':field, 'weight_field':weight_field, 'accumulation':accumulation})
+
+ def add_projection(self, field, weight_field=None):
+ "Add a field for projection."
+
+ self.projection_fields.append({'field':field, 'weight_field':weight_field})
+
+ def make_profiles(self, filename=None, prefilters=None, **kwargs):
+ "Make radial profiles for all halos on the list."
+
+ # Check to see if the VirialFilter has been added to the filter list.
+ # If a lower mass cutoff is being used, use it to make a pre-filter.
+ if prefilters is None: prefilters = []
+ virial_filter = True
+ virial_prefilter = None
+ virial_prefilter_safety_factor = 0.1
+ all_filter_functions = [hf['function'] for hf in self._halo_filters]
+ if 'mass' in self.halo_list_format and VirialFilter in all_filter_functions:
+ vFilter = self._halo_filters[all_filter_functions.index(VirialFilter)]
+ if vFilter['kwargs'].has_key('virial_filters') and \
+ vFilter['kwargs']['virial_filters'] is not None:
+ all_vqFilters = [vqf[0] for vqf in vFilter['kwargs']['virial_filters']]
+ if 'TotalMassMsun' in all_vqFilters:
+ mass_filter = vFilter['kwargs']['virial_filters'][all_vqFilters.index('TotalMassMsun')]
+ if '>' in mass_filter[1]:
+ virial_prefilter = "halo['mass'] %s %f * %s" % (mass_filter[1], virial_prefilter_safety_factor, mass_filter[2])
+ prefilters.append(virial_prefilter)
+ elif '<' in mass_filter[1]:
+ virial_prefilter = "halo['mass'] %s %f * %s" % (mass_filter[1], (1./virial_prefilter_safety_factor), mass_filter[2])
+ prefilters.append(virial_prefilter)
# Add profile fields necessary for calculating virial quantities.
- self._CheckForNeededProfileFields()
+ if virial_filter: self._check_for_needed_profile_fields()
- outputDir = "%s/%s" % (self.pf.fullpath,self.haloProfilerParameters['ProfileOutputDir'])
+ outputDir = "%s/%s" % (self.pf.fullpath, self.profile_output_dir)
self.__check_directory(outputDir)
- for q,halo in enumerate(self._get_objs('hopHalos', round_robin=True)):
- filename = "%s/Halo_%04d_profile.dat" % (outputDir,halo['id'])
-
- # Read profile from file if it already exists.
- # If not, profile will be None.
- profile = self._ReadProfile(filename)
-
- # Make profile if necessary.
- newProfile = profile is None
- if profile is None:
-
- r_min = 2*self.pf.h.get_smallest_dx() * self.pf['mpc']
- if (halo['r_max'] / r_min < PROFILE_RADIUS_THRESHOLD):
- mylog.error("Skipping halo with r_max / r_min = %f." % (halo['r_max']/r_min))
- continue
-
- sphere = self.pf.h.sphere(halo['center'],halo['r_max']/self.pf.units['mpc'])
- if len(sphere._grids) == 0: continue
+ # Profile all halos.
+ for halo in self._get_objs('all_halos', round_robin=True):
+ profile_filename = "%s/Halo_%04d_profile.dat" % (outputDir, halo['id'])
+
+ profiledHalo = self._get_halo_profile(halo, profile_filename, virial_filter=virial_filter, prefilters=prefilters)
- if self.needBulkVelocity:
- # Set bulk velocity to zero out radial velocity profiles.
- if self.haloProfilerParameters['VelocityCenter'][0] == 'bulk':
- if self.haloProfilerParameters['VelocityCenter'][1] == 'halo':
- sphere.set_field_parameter('bulk_velocity',halo['velocity'])
- elif self.haloProfilerParameters['VelocityCenter'][1] == 'sphere':
- sphere.set_field_parameter('bulk_velocity',sphere.quantities['BulkVelocity']())
- else:
- mylog.error("Invalid parameter: VelocityCenter.")
- elif self.haloProfilerParameters['VelocityCenter'][0] == 'max':
- max_grid,max_cell,max_value,max_location = self.pf.h.find_max_cell_location(self.haloProfilerParameters['VelocityCenter'][1])
- sphere.set_field_parameter('bulk_velocity',[max_grid['x-velocity'][max_cell],
- max_grid['y-velocity'][max_cell],
- max_grid['z-velocity'][max_cell]])
-
- profile = lagos.BinnedProfile1D(sphere,self.haloProfilerParameters['n_bins'],"RadiusMpc",
- r_min,halo['r_max'],
- log_space=True, lazy_reader=False)
- for field in self.profileFields.keys():
- profile.add_fields(field,weight=self.profileFields[field][0],
- accumulation=self.profileFields[field][1])
-
- self._AddActualOverdensity(profile)
-
- virial = self._CalculateVirialQuantities(profile)
- virial['center'] = halo['center']
- virial['id'] = halo['id']
-
- if (virial['TotalMassMsun'] >= self.haloProfilerParameters['VirialMassCutoff']):
- self.virialQuantities.append(virial)
- if newProfile:
- mylog.info("Writing halo %d" % virial['id'])
- profile.write_out(filename, format='%0.6e')
- del profile
+ if profiledHalo is None:
+ continue
+
+ # Apply filter and keep track of the quantities that are returned.
+ haloQuantities = {}
+ filterResult = True
+ for hFilter in self._halo_filters:
+ filterResult, filterQuantities = hFilter['function'](profiledHalo, *hFilter['args'], **hFilter['kwargs'])
+
+ if not filterResult: break
+
+ if filterQuantities is not None:
+ haloQuantities.update(filterQuantities)
+
+ if filterResult:
+ haloQuantities['id'] = halo['id']
+ if halo.has_key('center'): haloQuantities['center'] = halo['center']
+ self.filtered_halos.append(haloQuantities)
+
+ self.filtered_halos = self._mpi_catlist(self.filtered_halos)
+ self.filtered_halos.sort(key = lambda a:a['id'])
+
+ if filename is not None:
+ self._write_filtered_halo_list(filename, **kwargs)
+
+ def _get_halo_profile(self, halo, filename, virial_filter=True, prefilters=None, force_write=False):
+ """
+ Profile a single halo and write profile data to a file.
+ If file already exists, read profile data from file.
+ Return a dictionary of id, center, and virial quantities if virial_filter is True.
+ """
+
+ # Apply prefilters to avoid profiling unwanted halos.
+ if prefilters is not None:
+ for prefilter in prefilters:
+ if not eval(prefilter):
+ return None
+
+ # Read profile from file if it already exists.
+ # If not, profile will be None.
+ profile = self._read_profile(filename)
+
+ # Make profile if necessary.
+ newProfile = profile is None
+ if newProfile:
+
+ r_min = 2 * self.pf.h.get_smallest_dx() * self.pf['mpc']
+ if (halo['r_max'] / r_min < PROFILE_RADIUS_THRESHOLD):
+ mylog.error("Skipping halo with r_max / r_min = %f." % (halo['r_max']/r_min))
+ return None
- if newProfile:
+ sphere = self.pf.h.sphere(halo['center'], halo['r_max']/self.pf.units['mpc'])
+ if len(sphere._grids) == 0: return None
+ new_sphere = False
+
+ if self.use_density_center:
+ dc_x = sphere.quantities['WeightedAverageQuantity']('x', 'MatterDensityXTotalMass')
+ dc_y = sphere.quantities['WeightedAverageQuantity']('y', 'MatterDensityXTotalMass')
+ dc_z = sphere.quantities['WeightedAverageQuantity']('z', 'MatterDensityXTotalMass')
+ mylog.info("Moving halo center from %s to %s." % (halo['center'], [dc_x, dc_y, dc_z]))
+ halo['center'] = [dc_x, dc_y, dc_z]
+ new_sphere = True
+
+ if self.use_field_max_center is not None:
+ ma, maxi, mx, my, mz, mg = sphere.quantities['MaxLocation'](self.use_field_max_center)
+ mylog.info("Moving halo center from %s to %s." % (halo['center'], [mx, my, mz]))
+ halo['center'] = [mx, my, mz]
+ new_sphere = True
+
+ if new_sphere:
# Temporary solution to memory leak.
for g in self.pf.h.grids:
g.clear_data()
sphere.clear_data()
del sphere
+ sphere = self.pf.h.sphere(halo['center'], halo['r_max']/self.pf.units['mpc'])
- self._WriteVirialQuantities()
+ if self._need_bulk_velocity:
+ # Set bulk velocity to zero out radial velocity profiles.
+ if self.velocity_center[0] == 'bulk':
+ if self.velocity_center[1] == 'halo':
+ sphere.set_field_parameter('bulk_velocity', halo['velocity'])
+ elif self.velocity_center[1] == 'sphere':
+ sphere.set_field_parameter('bulk_velocity', sphere.quantities['BulkVelocity']())
+ else:
+ mylog.error("Invalid parameter: VelocityCenter.")
+ elif self.velocity_center[0] == 'max':
+ max_grid, max_cell, max_value, max_location = self.pf.h.find_max_cell_location(self.velocity_center[1])
+ sphere.set_field_parameter('bulk_velocity', [max_grid['x-velocity'][max_cell],
+ max_grid['y-velocity'][max_cell],
+ max_grid['z-velocity'][max_cell]])
+
+ profile = lagos.BinnedProfile1D(sphere, self.n_profile_bins, "RadiusMpc",
+ r_min, halo['r_max'],
+ log_space=True, lazy_reader=False)
+ for hp in self.profile_fields:
+ profile.add_fields(hp['field'], weight=hp['weight_field'], accumulation=hp['accumulation'])
+
+ profiledHalo = {}
+ if virial_filter:
+ self._add_actual_overdensity(profile)
+
+ profiledHalo['center'] = halo['center']
+ profiledHalo['id'] = halo['id']
+
+ if newProfile:
+ mylog.info("Writing halo %d" % profiledHalo['id'])
+ profile.write_out(filename, format='%0.6e')
+ elif force_write:
+ mylog.info("Re-writing halo %d" % profiledHalo['id'])
+ self._write_profile(profile, filename, format='%0.6e')
+
+ if newProfile:
+ # Temporary solution to memory leak.
+ for g in self.pf.h.grids:
+ g.clear_data()
+ sphere.clear_data()
+ del sphere
- def _finalize_parallel(self):
- self.virialQuantities = self._mpi_catlist(self.virialQuantities)
- self.virialQuantities.sort(key = lambda a:a['id'])
+ return profile
- @lagos.parallel_root_only
- def __check_directory(self, outputDir):
- if (os.path.exists(outputDir)):
- if not(os.path.isdir(outputDir)):
- mylog.error("Output directory exists, but is not a directory: %s." % outputDir)
- raise IOError(outputDir)
+ def make_projections(self, axes=[0, 1, 2], halo_list='filtered', save_images=False, save_cube=True, **kwargs):
+ "Make projections of all halos using specified fields."
+
+ # Get list of halos for projecting.
+ if halo_list == 'filtered':
+ self._halo_projection_list = self.filtered_halos
+ elif halo_list == 'all':
+ self._halo_projection_list = self.all_halos
+ elif isinstance(halo_list, types.StringType):
+ self._halo_projection_list = self._read_halo_list(halo_list)
+ elif isinstance(halo_list, types.ListType):
+ self._halo_projection_list = halo_list
else:
- os.mkdir(outputDir)
+ mylog.error("Keyword, halo_list', must be 'filtered', 'all', a filename, or an actual list.")
+ return
- def makeProjections(self,save_images=False,save_cube=True,**kwargs):
- "Make projections of all halos using specified fields."
- # Get virial quantities.
- self._LoadVirialData()
+ if len(self._halo_projection_list) == 0:
+ mylog.error("Halo list for projections is empty.")
+ return
# Set resolution for fixed resolution output.
if save_cube:
- if (str(self.haloProfilerParameters['ProjectAtLevel']).find('max') >= 0):
+ if self.project_at_level == 'max':
proj_level = self.pf.h.maxLevel
else:
- proj_level = int(self.haloProfilerParameters['ProjectAtLevel'])
- proj_dx = self.pf.units['mpc'] / self.pf.parameters['TopGridDimensions'][0] / \
+ proj_level = int(self.project_at_level)
+ proj_dx = self.pf.units[self.projection_width_units] / self.pf.parameters['TopGridDimensions'][0] / \
(self.pf.parameters['RefineBy']**proj_level)
- projectionResolution = int(self.haloProfilerParameters['ProjectionWidth'] / proj_dx)
+ projectionResolution = int(self.projection_width / proj_dx)
- outputDir = "%s/%s" % (self.pf.fullpath,self.haloProfilerParameters['ProjectionOutputDir'])
+ outputDir = "%s/%s" % (self.pf.fullpath, self.projection_output_dir)
self.__check_directory(outputDir)
center = [0.5 * (self.pf.parameters['DomainLeftEdge'][w] + self.pf.parameters['DomainRightEdge'][w])
for w in range(self.pf.parameters['TopGridRank'])]
# Create a plot collection.
- pc = raven.PlotCollection(self.pf,center=center)
+ pc = raven.PlotCollection(self.pf, center=center)
- for halo in self._get_objs('virialQuantities', round_robin=True):
+ for halo in self._get_objs('_halo_projection_list', round_robin=True):
if halo is None:
continue
# Check if region will overlap domain edge.
# Using non-periodic regions is faster than using periodic ones.
- leftEdge = [(halo['center'][w] - 0.5 * self.haloProfilerParameters['ProjectionWidth']/self.pf.units['mpc'])
+ leftEdge = [(halo['center'][w] - 0.5 * self.projection_width/self.pf.units[self.projection_width_units])
for w in range(len(halo['center']))]
- rightEdge = [(halo['center'][w] + 0.5 * self.haloProfilerParameters['ProjectionWidth']/self.pf.units['mpc'])
+ rightEdge = [(halo['center'][w] + 0.5 * self.projection_width/self.pf.units[self.projection_width_units])
for w in range(len(halo['center']))]
mylog.info("Projecting halo %04d in region: [%f, %f, %f] to [%f, %f, %f]." %
- (halo['id'],leftEdge[0],leftEdge[1],leftEdge[2],rightEdge[0],rightEdge[1],rightEdge[2]))
+ (halo['id'], leftEdge[0], leftEdge[1], leftEdge[2], rightEdge[0], rightEdge[1], rightEdge[2]))
need_per = False
for w in range(len(halo['center'])):
@@ -246,192 +395,158 @@
break
if need_per:
- region = self.pf.h.periodic_region(halo['center'],leftEdge,rightEdge)
+ region = self.pf.h.periodic_region(halo['center'], leftEdge, rightEdge)
else:
- region = self.pf.h.region(halo['center'],leftEdge,rightEdge)
+ region = self.pf.h.region(halo['center'], leftEdge, rightEdge)
# Make projections.
- for w in range(self.pf.parameters['TopGridRank']):
+ if not isinstance(axes, types.ListType): axes = list([axes])
+ for w in axes:
# YT projections do not follow the right-hand rule.
coords = range(3)
del coords[w]
x_axis = coords[0]
y_axis = coords[1]
- for field in self.projectionFields.keys():
- pc.add_projection(field,w,weight_field=self.projectionFields[field],source=region,lazy_reader=False,
- serialize=False,**kwargs)
+ for hp in self.projection_fields:
+ pc.add_projection(hp['field'], w, weight_field=hp['weight_field'], source=region, lazy_reader=False,
+ serialize=False, **kwargs)
# Set x and y limits, shift image if it overlaps domain boundary.
if need_per:
- pw = self.haloProfilerParameters['ProjectionWidth']/self.pf.units['mpc']
- ShiftProjections(self.pf,pc,halo['center'],center,w)
+ pw = self.projection_width/self.pf.units[self.projection_width_units]
+ shift_projections(self.pf, pc, halo['center'], center, w)
# Projection has now been shifted to center of box.
proj_left = [center[x_axis]-0.5*pw, center[y_axis]-0.5*pw]
proj_right = [center[x_axis]+0.5*pw, center[y_axis]+0.5*pw]
else:
- proj_left = [leftEdge[x_axis],leftEdge[y_axis]]
- proj_right = [rightEdge[x_axis],rightEdge[y_axis]]
+ proj_left = [leftEdge[x_axis], leftEdge[y_axis]]
+ proj_right = [rightEdge[x_axis], rightEdge[y_axis]]
- pc.set_xlim(proj_left[0],proj_right[0])
- pc.set_ylim(proj_left[1],proj_right[1])
+ pc.set_xlim(proj_left[0], proj_right[0])
+ pc.set_ylim(proj_left[1], proj_right[1])
# Save projection data to hdf5 file.
if save_cube:
- axes = ['x','y','z']
+ axis_labels = ['x', 'y', 'z']
dataFilename = "%s/Halo_%04d_%s_data.h5" % \
- (outputDir,halo['id'],axes[w])
+ (outputDir, halo['id'], axis_labels[w])
mylog.info("Saving projection data to %s." % dataFilename)
output = h5py.File(dataFilename, "a")
# Create fixed resolution buffer for each projection and write them out.
- for e,field in enumerate(self.projectionFields.keys()):
- frb = raven.FixedResolutionBuffer(pc.plots[e].data,(proj_left[0],proj_right[0],proj_left[1],proj_right[1]),
- (projectionResolution,projectionResolution),
+ for e, hp in enumerate(self.projection_fields):
+ frb = raven.FixedResolutionBuffer(pc.plots[e].data, (proj_left[0], proj_right[0], proj_left[1], proj_right[1]),
+ (projectionResolution, projectionResolution),
antialias=False)
- dataset_name = "%s_%s" % (field,self.projectionFields[field])
+ dataset_name = "%s_%s" % (hp['field'], hp['weight_field'])
if dataset_name in output.listnames(): del output[dataset_name]
- output.create_dataset(dataset_name,data=frb[field])
+ output.create_dataset(dataset_name, data=frb[hp['field']])
output.close()
if save_images:
- pc.save("%s/Halo_%04d" % (outputDir,halo['id']),force_save=True)
+ pc.save("%s/Halo_%04d" % (outputDir, halo['id']), force_save=True)
pc.clear_plots()
del region
del pc
- @lagos.parallel_root_only
- def _WriteVirialQuantities(self):
- "Write out file with halo centers and virial masses and radii."
- filename = "%s/%s" % (self.pf.fullpath,self.haloProfilerParameters['VirialQuantitiesOutputFile'])
- mylog.info("Writing virial quantities to %s." % filename)
- file = open(filename,'w')
- file.write("#Index\tx\ty\tz\tMass [Msolar]\tRadius [Mpc]\n")
- for vq in self.virialQuantities:
- if vq is not None:
- file.write("%04d %.10f %.10f %.10f %.6e %.6e\n" % (vq['id'],vq['center'][0],vq['center'][1],vq['center'][2],
- vq['TotalMassMsun'], vq['RadiusMpc']))
- file.close()
-
- def _AddActualOverdensity(self,profile):
+ def _add_actual_overdensity(self, profile):
"Calculate overdensity from TotalMassMsun and CellVolume fields."
- if (profile.keys()).count('ActualOverdensity') > 0:
+ if 'ActualOverdensity' in profile.keys():
return
- rho_crit_now = 1.8788e-29 * self.pf['CosmologyHubbleConstantNow']**2.0 # g cm^-3
+ rho_crit_now = 1.8788e-29 * self.pf['CosmologyHubbleConstantNow']**2.0 * \
+ self.pf['CosmologyOmegaMatterNow'] # g cm^-3
Msun2g = 1.989e33
rho_crit = rho_crit_now * ((1.0 + self.pf['CosmologyCurrentRedshift'])**3.0)
profile['ActualOverdensity'] = (Msun2g * profile['TotalMassMsun']) / \
profile['CellVolume'] / rho_crit
- def _CalculateVirialQuantities(self,profile,overdensity_field='ActualOverdensity'):
- "Calculate virial radius and virial mass from radial profiles."
+ def _check_for_needed_profile_fields(self):
+ "Make sure CellVolume and TotalMass fields are added so virial quantities can be calculated."
+ all_profile_fields = [hp['field'] for hp in self.profile_fields]
+ if not 'CellVolume' in all_profile_fields:
+ mylog.info("Adding CellVolume field to so virial quantities can be calculated")
+ self.add_profile('CellVolume', weight_field=None, accumulation=True)
+ if not 'TotalMassMsun' in all_profile_fields:
+ mylog.info("Adding TotalMassMsun field to so virial quantities can be calculated")
+ self.add_profile('TotalMassMsun', weight_field=None, accumulation=True)
- fields = ['TotalMassMsun','RadiusMpc']
- overDensity = []
- temp_profile = {}
- for field in fields:
- temp_profile[field] = []
+ def _load_halo_data(self, filename=None):
+ "Read hop output file or run hop if it doesn't exist."
- for q in range(len(profile[overdensity_field])):
- good = True
- if (profile[overdensity_field][q] != profile[overdensity_field][q]):
- good = False
- continue
- for field in fields:
- if (profile[field][q] != profile[field][q]):
- good = False
- break
- if good:
- overDensity.append(profile[overdensity_field][q])
- for field in fields:
- temp_profile[field].append(profile[field][q])
+ # Don't run if hop data already loaded.
+ if self.all_halos:
+ return
- virial = {}
- if (len(overDensity) < 2):
- mylog.error("Skipping halo with no valid points in profile.")
- for field in fields:
- virial[field] = 0.0
- return virial
-
- if (overDensity[1] <= self.haloProfilerParameters['VirialOverdensity']):
- index = 0
- elif (overDensity[-1] >= self.haloProfilerParameters['VirialOverdensity']):
- index = -2
- else:
- for q in (na.array(range(len(overDensity)-2))+2):
- if (overDensity[q] < self.haloProfilerParameters['VirialOverdensity']):
- index = q - 1
- break
+ if filename is None:
+ filename = self.halo_list_file
- for field in fields:
- if (overDensity[index+1] - overDensity[index]) == 0:
- virial[field] = 0.0
- else:
- slope = (temp_profile[field][index+1] - temp_profile[field][index]) / \
- (overDensity[index+1] - overDensity[index])
- value = slope * (self.haloProfilerParameters['VirialOverdensity'] - overDensity[index]) + \
- temp_profile[field][index]
- virial[field] = value
+ hopFile = "%s/%s" % (self.pf.fullpath, filename)
- return virial
+ if not(os.path.exists(hopFile)):
+ mylog.info("Hop file not found, running hop to get halos.")
+ self._run_hop(hopFile)
- def _CheckForNeededProfileFields(self):
- "Make sure CellVolume and TotalMass fields are added so virial quantities can be calculated."
- if not(self.profileFields.has_key('CellVolume')):
- mylog.info("Adding CellVolume field to so virial quantities can be calculated")
- self.profileFields['CellVolume'] = [None,True]
- if not(self.profileFields.has_key('TotalMassMsun')):
- mylog.info("Adding TotalMassMsun field to so virial quantities can be calculated")
- self.profileFields['TotalMassMsun'] = [None,True]
+ self.all_halos = self._read_halo_list(hopFile)
- def _ReadHopFile(self,hopFile):
- "Read hop file to get halo information."
- mylog.info("Reading halo information from %s." % hopFile)
- hopLines = file(hopFile)
+ def _read_halo_list(self, listFile):
+ """
+ Read halo list from aue file.
+ Allow for columnar data in varying formats.
+ """
+
+ def __isE(arg):
+ parts = arg.lower().split('e')
+ if len(parts) != 2: return False
+ return not (True in [q.isalpha() for q in ''.join(parts)])
+
+ def __get_num(arg):
+ if __isE(arg):
+ return float(arg)
+ if arg != arg.swapcase():
+ return arg
+ return float(arg)
+
+ mylog.info("Reading halo information from %s." % listFile)
+ haloList = []
+ listLines = file(listFile)
+
+ fields = self.halo_list_format.keys()
+ getID = not 'id' in fields
+ getR_max = not 'r_max' in fields
- for line in hopLines:
+ for line in listLines:
line = line.strip()
if not(line.startswith('#')):
+ halo = {}
onLine = line.split()
- id = int(onLine[0])
- mass = float(onLine[1])
- if (mass >= self.haloProfilerParameters['VirialMassCutoff']):
- center = [float(onLine[7]),float(onLine[8]),float(onLine[9])]
- velocity = [float(onLine[10]),float(onLine[11]),float(onLine[12])]
- r_max = float(onLine[13]) * self.pf.units['mpc']
- halo = {'id': id, 'center': center, 'r_max': r_max, 'velocity': velocity}
- self.hopHalos.append(halo)
-
- mylog.info("Loaded %d halos with total dark matter mass af at least %e Msolar." %
- (len(self.hopHalos),self.haloProfilerParameters['VirialMassCutoff']))
-
- def _ReadOldHopFile(self,hopFile):
- "Read old style hop file made by enzo_hop."
- mylog.info("Reading halo information from old style hop file %s." % hopFile)
- hopLines = file(hopFile)
+ for field in fields:
+ if isinstance(self.halo_list_format[field], types.ListType):
+ halo[field] = [__get_num(onLine[q]) for q in self.halo_list_format[field]]
+ else:
+ halo[field] = __get_num(onLine[self.halo_list_format[field]])
+ if getID: halo['id'] = len(haloList)
+ if getR_max:
+ halo['r_max'] = self.halo_radius * self.pf.units['mpc']
+ else:
+ halo['r_max'] *= self.pf.units['mpc']
+ haloList.append(halo)
- for line in hopLines:
- line = line.strip()
- if not(line.startswith('#')):
- onLine = line.split()
- id = int(onLine[0])
- center = [float(onLine[4]),float(onLine[5]),float(onLine[6])]
- r_max = self.haloRadius * self.pf.units['mpc']
- halo = {'id': id, 'center': center, 'r_max': r_max}
- self.hopHalos.append(halo)
+ mylog.info("Loaded %d halos." % (len(haloList)))
+ return haloList
- def _ReadProfile(self,profileFile):
+ def _read_profile(self, profileFile):
"Read radial profile from file. Return None if it doesn't have all the fields requested."
# Check to see if file exists.
if not os.path.exists(profileFile):
return None
- f = open(profileFile,'r')
+ f = open(profileFile, 'r')
lines = f.readlines()
f.close()
@@ -447,8 +562,9 @@
profile[field] = []
# Check if all fields needed are present.
- for field in self.profileFields:
- if not profile.has_key(field):
+ all_profile_fields = [hp['field'] for hp in self.profile_fields]
+ for field in all_profile_fields:
+ if not field in profile:
return None
# Fill profile fields, skip bad values.
@@ -461,138 +577,72 @@
lineOK = False
break
if lineOK:
- for q,field in enumerate(fields):
+ for q, field in enumerate(fields):
profile[field].append(float(onLine[q]))
+ for field in fields:
+ profile[field] = na.array(profile[field])
+
if len(profile[fields[0]]) > 1:
return profile
else:
return None
- def _RunHop(self,hopFile):
+ def _run_hop(self, hopFile):
"Run hop to get halos."
- hop_results = HaloFinder(self.pf)
+ hop_results = self.halo_finder_function(self.pf, *self.halo_finder_args, **self.halo_finder_kwargs)
hop_results.write_out(hopFile)
del hop_results
- def _LoadHopData(self):
- "Read hop output file or run hop if it doesn't exist."
+ @lagos.parallel_root_only
+ def _write_filtered_halo_list(self, filename, format="%s"):
+ """
+ Write out list of filtered halos along with any quantities
+ picked up during the filtering process.
+ """
- # Don't run if hop data already loaded.
- if self.hopHalos:
+ if len(self.filtered_halos) == 0:
+ mylog.error("No halos in filtered list.")
return
- hopFile = "%s/%s" % (self.pf.fullpath,
- self.haloProfilerParameters['HopOutputFile'])
+ filename = "%s/%s" % (self.pf.fullpath, filename)
+ mylog.info("Writing filtered halo list to %s." % filename)
+ file = open(filename, "w")
+ fields = [field for field in sorted(self.filtered_halos[0])]
+ fields.pop(fields.index('id'))
+ fields.pop(fields.index('center'))
+ file.write("\t".join(["#Index", 'x', 'y', 'z'] + fields + ["\n"]))
+
+ for halo in self.filtered_halos:
+ file.write("%04d\t" % halo['id'])
+ file.write("%.12f\t%.12f\t%.12f\t" % (halo['center'][0], halo['center'][1], halo['center'][2]))
+ field_data = na.array([halo[field] for field in fields])
+ field_data.tofile(file, sep="\t", format=format)
+ file.write("\n")
+ file.close()
- if not(os.path.exists(hopFile)):
- mylog.info("Hop file not found, running hop to get halos.")
- self._RunHop(hopFile)
+ def _write_profile(self, profile, filename, format="%0.16e"):
+ fid = open(filename, "w")
+ fields = [field for field in sorted(profile.keys()) if field != "UsedBins"]
+ fid.write("\t".join(["#"] + fields + ["\n"]))
+ field_data = na.array([profile[field] for field in fields])
+ for line in range(field_data.shape[1]):
+ field_data[:, line].tofile(fid, sep="\t", format=format)
+ fid.write("\n")
+ fid.close()
- if self.hop_style is 'new':
- self._ReadHopFile(hopFile)
+ @lagos.parallel_root_only
+ def __check_directory(self, outputDir):
+ if (os.path.exists(outputDir)):
+ if not(os.path.isdir(outputDir)):
+ mylog.error("Output directory exists, but is not a directory: %s." % outputDir)
+ raise IOError(outputDir)
else:
- self._ReadOldHopFile(hopFile)
-
- def _LoadVirialData(self):
- "Read virial quantities data or ask for profiles to be run if it doesn't exist."
-
- if self.virialQuantities:
- return
-
- virialFile = "%s/%s" % (self.pf.fullpath,
- self.haloProfilerParameters['VirialQuantitiesOutputFile'])
- if not(os.path.exists(virialFile)):
- mylog.info("Virial quantities file not found. Making profiles to calculate virial quantities.")
- self.makeProfiles()
-
- mylog.info("Reading virial quantities from %s." % virialFile)
- virialLines = file(virialFile)
-
- halos = 0
- for line in virialLines:
- line = line.strip()
- if not(line.startswith('#')):
- onLine = line.split()
- index = int(onLine[0])
- virial = {}
- virial['id'] = int(onLine[0])
- virial['center'] = [float(onLine[1]),float(onLine[2]),float(onLine[3])]
- virial['TotalMassMsun'] = float(onLine[4])
- virial['RadiusMpc'] = float(onLine[5])
- if (virial['TotalMassMsun'] >= self.haloProfilerParameters['VirialMassCutoff']):
- self.virialQuantities.append(virial)
- halos += 1
-
- mylog.info("Loaded virial quantities for %d halos." % halos)
-
- def _ReadHaloProfilerParameterFile(self):
- "Read a halo profiler parameter file."
- lines = open(self.HaloProfilerParameterFile).readlines()
- for line in lines:
- if line.find("#") >= 0: # Keep the commented lines
- line=line[:line.find("#")]
- line=line.strip()
- if len(line) < 2:
- continue
- try:
- param, vals = map(str,line.split("="))
- param = param.strip()
- vals = vals.strip()
- except ValueError:
- mylog.error("Skipping line: %s" % line)
- continue
- if haloProfilerParameterDict.has_key(param):
- t = map(haloProfilerParameterDict[param], vals.split())
- if len(t) == 1:
- self.haloProfilerParameters[param] = t[0]
- else:
- self.haloProfilerParameters[param] = t
- elif param.startswith("Profile["):
- field = param[param.find("[")+1:param.find("]")]
- field_vals = vals.split(',')
- for val in field_vals:
- val.strip()
- if (field_vals[0].find('None') >= 0):
- field_vals[0] = None
- if (field_vals[1].find('True') >= 0):
- field_vals[1] = True
- else:
- field_vals[1] = False
- self.profileFields[field] = field_vals
- elif param.startswith("Projection["):
- field = param[param.find("[")+1:param.find("]")]
- if (vals.find('None') >= 0):
- vals = None
- self.projectionFields[field] = vals
-
- def _SetParameterDefaults(self):
- "Set some default parameters."
- self.haloProfilerParameters['VirialMassCutoff'] = 1e13 # M_solar
- self.haloProfilerParameters['VirialOverdensity'] = 200
- self.haloProfilerParameters['VelocityCenter'] = ['bulk','halo']
- self.haloProfilerParameters['n_bins'] = 50
- self.haloProfilerParameters['ProfileOutputDir'] = 'radial_profiles'
- self.haloProfilerParameters['ProjectionOutputDir'] = 'projections'
- self.haloProfilerParameters['HopOutputFile'] = "HopAnalysis.out"
- self.haloProfilerParameters['VirialQuantitiesOutputFile'] = "VirialQuantities.out"
- self.haloProfilerParameters['ProjectionWidth'] = 4.0 # Mpc
- self.haloProfilerParameters['ProjectAtLevel'] = 'max'
-
-haloProfilerParameterDict = {"ProfileOutputDir": str,
- "VirialQuantitiesOutputFile": str,
- "HopOutputFile": str,
- "VirialMassCutoff": float,
- "VirialOverdensity": float,
- "VelocityCenter": str,
- "n_bins": int,
- "ProjectionOutputDir": str,
- "ProjectionWidth": float,
- "ProjectAtLevel": str}
+ os.mkdir(outputDir)
-def ShiftProjections(pf,pc,oldCenter,newCenter,axis):
+def shift_projections(pf, pc, oldCenter, newCenter, axis):
"""
Shift projection data around.
This is necessary when projecting a preiodic region.
@@ -605,7 +655,7 @@
for plot in pc.plots:
# Get name of data field.
- other_fields = {'px':True,'py':True,'pdx':True,'pdy':True,'weight_field':True}
+ other_fields = {'px':True, 'py':True, 'pdx':True, 'pdy':True, 'weight_field':True}
for pfield in plot.data.data.keys():
if not(other_fields.has_key(pfield)):
field = pfield
@@ -661,18 +711,18 @@
add2_y_weight_field = plot['weight_field'][plot['py'] - 0.5 * plot['pdy'] < 0]
# Add the hanging cells back to the projection data.
- plot.data['px'] = na.concatenate([plot['px'],add_x_px,add_y_px,add2_x_px,add2_y_px])
- plot.data['py'] = na.concatenate([plot['py'],add_x_py,add_y_py,add2_x_py,add2_y_py])
- plot.data['pdx'] = na.concatenate([plot['pdx'],add_x_pdx,add_y_pdx,add2_x_pdx,add2_y_pdx])
- plot.data['pdy'] = na.concatenate([plot['pdy'],add_x_pdy,add_y_pdy,add2_x_pdy,add2_y_pdy])
- plot.data[field] = na.concatenate([plot[field],add_x_field,add_y_field,add2_x_field,add2_y_field])
+ plot.data['px'] = na.concatenate([plot['px'], add_x_px, add_y_px, add2_x_px, add2_y_px])
+ plot.data['py'] = na.concatenate([plot['py'], add_x_py, add_y_py, add2_x_py, add2_y_py])
+ plot.data['pdx'] = na.concatenate([plot['pdx'], add_x_pdx, add_y_pdx, add2_x_pdx, add2_y_pdx])
+ plot.data['pdy'] = na.concatenate([plot['pdy'], add_x_pdy, add_y_pdy, add2_x_pdy, add2_y_pdy])
+ plot.data[field] = na.concatenate([plot[field], add_x_field, add_y_field, add2_x_field, add2_y_field])
plot.data['weight_field'] = na.concatenate([plot['weight_field'],
- add_x_weight_field,add_y_weight_field,add2_x_weight_field,add2_y_weight_field])
+ add_x_weight_field, add_y_weight_field, add2_x_weight_field, add2_y_weight_field])
# Delete original copies of hanging cells.
- del add_x_px,add_y_px,add2_x_px,add2_y_px
- del add_x_py,add_y_py,add2_x_py,add2_y_py
- del add_x_pdx,add_y_pdx,add2_x_pdx,add2_y_pdx
- del add_x_pdy,add_y_pdy,add2_x_pdy,add2_y_pdy
- del add_x_field,add_y_field,add2_x_field,add2_y_field
- del add_x_weight_field,add_y_weight_field,add2_x_weight_field,add2_y_weight_field
+ del add_x_px, add_y_px, add2_x_px, add2_y_px
+ del add_x_py, add_y_py, add2_x_py, add2_y_py
+ del add_x_pdx, add_y_pdx, add2_x_pdx, add2_y_pdx
+ del add_x_pdy, add_y_pdy, add2_x_pdy, add2_y_pdy
+ del add_x_field, add_y_field, add2_x_field, add2_y_field
+ del add_x_weight_field, add_y_weight_field, add2_x_weight_field, add2_y_weight_field
More information about the yt-svn
mailing list