[Yt-svn] yt-commit r1516 - in branches/yt-1.5/yt/extensions: . lightcone
britton at wrangler.dreamhost.com
britton at wrangler.dreamhost.com
Tue Nov 3 15:45:19 PST 2009
Author: britton
Date: Tue Nov 3 15:45:17 2009
New Revision: 1516
URL: http://yt.enzotools.org/changeset/1516
Log:
Added new light cone generator, EnzoSimulation, and HaloProfiler, all with
docstrings.
Added:
branches/yt-1.5/yt/extensions/EnzoSimulation.py
Removed:
branches/yt-1.5/yt/extensions/lightcone/README
Modified:
branches/yt-1.5/yt/extensions/HaloProfiler.py
branches/yt-1.5/yt/extensions/lightcone/HaloMask.py
branches/yt-1.5/yt/extensions/lightcone/LightCone.py
branches/yt-1.5/yt/extensions/lightcone/LightConeProjection.py
branches/yt-1.5/yt/extensions/lightcone/UniqueSolution.py
Added: branches/yt-1.5/yt/extensions/EnzoSimulation.py
==============================================================================
--- (empty file)
+++ branches/yt-1.5/yt/extensions/EnzoSimulation.py Tue Nov 3 15:45:17 2009
@@ -0,0 +1,463 @@
+import yt.lagos as lagos
+from yt.logger import lagosLogger as mylog
+import numpy as na
+
+dt_Tolerance = 1e-3
+
+class EnzoSimulation(object):
+ """
+ Super class for performing the same operation over all data dumps in
+ a simulation from one redshift to another.
+ """
+ def __init__(self, EnzoParameterFile, initial_time=None, final_time=None, initial_redshift=None, final_redshift=None,
+ links=False, enzo_parameters=None, get_time_outputs=True, get_redshift_outputs=True):
+ """
+ Initialize an EnzoSimulation object.
+ :param initial_time (float): the initial time in code units for the dataset list. Default: None.
+ :param final_time (float): the final time in code units for the dataset list. Default: None.
+ :param initial_redshift (float): the initial (highest) redshift for the dataset list. Only for
+ cosmological simulations. Default: None.
+ :param final_redshift (float): the final (lowest) redshift for the dataset list. Only for cosmological
+ simulations. Default: None.
+ :param links (bool): if True, each entry in the dataset list will contain entries, previous and next, that
+ point to the previous and next entries on the dataset list. Default: False.
+ :param enzo_parameters (dict): a dictionary specify additional parameters to be retrieved from the
+ parameter file. The format should be the name of the parameter as the key and the variable type as
+ the value. For example, {'CosmologyComovingBoxSize':float}. All parameter values will be stored in
+ the dictionary attribute, enzoParameters. Default: None.
+ :param get_time_outputs (bool): if False, the time datasets, specified in Enzo with the dtDataDump, will not
+ be added to the dataset list. Default: True.
+ :param get_redshift_outputs (bool): if False, the redshift datasets will not be added to the dataset list. Default: True.
+ """
+ self.EnzoParameterFile = EnzoParameterFile
+ self.enzoParameters = {}
+ self.redshiftOutputs = []
+ self.timeOutputs = []
+ self.allOutputs = []
+ self.InitialTime = initial_time
+ self.FinalTime = final_time
+ self.InitialRedshift = initial_redshift
+ self.FinalRedshift = final_redshift
+ self.links = links
+ self.get_time_outputs = get_time_outputs
+ self.get_redshift_outputs = get_redshift_outputs
+
+ # Add any extra parameters to parameter dict.
+ if enzo_parameters is None: enzo_parameters = {}
+ EnzoParameterDict.update(enzo_parameters)
+
+ # Set some parameter defaults.
+ self._SetParameterDefaults()
+
+ # Read parameters.
+ self._ReadEnzoParameterFile()
+
+ # Check for sufficient starting/ending parameters.
+ if self.InitialTime is None and self.InitialRedshift is None:
+ if self.enzoParameters['ComovingCoordinates'] and \
+ self.enzoParameters.has_key('CosmologyInitialRedshift'):
+ self.InitialRedshift = self.enzoParameters['CosmologyInitialRedshift']
+ elif self.enzoParameters.has_key('InitialTime'):
+ self.InitialTime = self.enzoParameters['InitialTime']
+ else:
+ mylog.error("Couldn't find parameter for initial time or redshift from parameter file.")
+ return None
+ if self.FinalTime is None and self.FinalRedshift is None:
+ if self.enzoParameters['ComovingCoordinates'] and \
+ self.enzoParameters.has_key('CosmologyFinalRedshift'):
+ self.FinalRedshift = self.enzoParameters['CosmologyFinalRedshift']
+ elif self.enzoParameters.has_key('StopTime'):
+ self.FinalTime = self.enzoParameters['StopTime']
+ else:
+ mylog.error("Couldn't find parameter for final time or redshift from parameter file.")
+ return None
+
+ # Convert initial/final redshifts to times.
+ if self.enzoParameters['ComovingCoordinates']:
+ # Instantiate a cosmology calculator.
+ self.cosmology = lagos.Cosmology(HubbleConstantNow =
+ (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
+ OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
+ OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'])
+
+ # Instantiate EnzoCosmology object for units and time conversions.
+ self.enzo_cosmology = lagos.EnzoCosmology(HubbleConstantNow =
+ (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
+ OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
+ OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'],
+ InitialRedshift = self.enzoParameters['CosmologyInitialRedshift'])
+ if self.InitialRedshift is not None:
+ self.InitialTime = self.enzo_cosmology.ComputeTimeFromRedshift(self.InitialRedshift) / \
+ self.enzo_cosmology.TimeUnits
+ if self.FinalRedshift is not None:
+ self.FinalTime = self.enzo_cosmology.ComputeTimeFromRedshift(self.FinalRedshift) / \
+ self.enzo_cosmology.TimeUnits
+
+ # Get initial time of simulation.
+ if self.enzoParameters['ComovingCoordinates'] and \
+ self.enzoParameters.has_key('CosmologyInitialRedshift'):
+ self.SimulationInitialTime = self.enzo_cosmology.InitialTime / self.enzo_cosmology.TimeUnits
+ elif self.enzoParameters.has_key('InitialTime'):
+ self.SimulationInitialTime = self.enzoParameters['InitialTime']
+ else:
+ self.SimulationInitialTime = 0.0
+
+ # Combine all data dumps.
+ self._CombineDataOutputs()
+
+ def _CalculateRedshiftDumpTimes(self):
+ "Calculates time from redshift of redshift dumps."
+
+ for output in self.redshiftOutputs:
+ output['time'] = self.enzo_cosmology.ComputeTimeFromRedshift(output['redshift']) / \
+ self.enzo_cosmology.TimeUnits
+
+ def _CalculateTimeDumps(self):
+ "Calculates time dumps and their redshifts if cosmological."
+
+ if self.enzoParameters['dtDataDump'] <= 0.0: return
+
+ index = 0
+ current_time = self.SimulationInitialTime
+ while (current_time <= self.FinalTime) or \
+ (abs(self.FinalTime - current_time) / self.FinalTime) < dt_Tolerance:
+ filename = "%s/%s%04d/%s%04d" % (self.enzoParameters['GlobalDir'],
+ self.enzoParameters['DataDumpDir'],index,
+ self.enzoParameters['DataDumpName'],index)
+
+ self.timeOutputs.append({'index':index,'filename':filename,'time':current_time})
+ if self.enzoParameters['ComovingCoordinates']:
+ t = self.enzo_cosmology.InitialTime + (self.enzoParameters['dtDataDump'] * self.enzo_cosmology.TimeUnits * index)
+ self.timeOutputs[-1]['redshift'] = self.enzo_cosmology.ComputeRedshiftFromTime(t)
+
+ current_time += self.enzoParameters['dtDataDump']
+ index += 1
+
+ def _CombineDataOutputs(self):
+ "Combines redshift and time data into one sorted list."
+
+ # Calculate redshifts for dt data dumps.
+ if self.enzoParameters.has_key('dtDataDump') and self.get_time_outputs:
+ self._CalculateTimeDumps()
+
+ # Calculate times for redshift dumps.
+ if self.enzoParameters['ComovingCoordinates'] and self.get_redshift_outputs:
+ self._CalculateRedshiftDumpTimes()
+
+ self.allOutputs = self.redshiftOutputs + self.timeOutputs
+ self.allOutputs.sort(key=lambda obj:obj['time'])
+
+ start_index = None
+ end_index = None
+
+ for q in range(len(self.allOutputs)):
+ del self.allOutputs[q]['index']
+
+ if start_index is None:
+ if self.allOutputs[q]['time'] >= self.InitialTime or \
+ abs(self.allOutputs[q]['time'] - self.InitialTime) < \
+ self.allOutputs[q]['time'] * dt_Tolerance:
+ start_index = q
+
+ if end_index is None:
+ if self.allOutputs[q]['time'] > self.FinalTime:
+ end_index = q - 1
+ if abs(self.allOutputs[q]['time'] - self.FinalTime) < \
+ self.allOutputs[q]['time'] * dt_Tolerance:
+ end_index = q
+ if q == len(self.allOutputs) - 1:
+ end_index = q
+
+ if self.links and start_index is not None:
+ if q == start_index:
+ self.allOutputs[q]['previous'] = None
+ else:
+ self.allOutputs[q]['previous'] = self.allOutputs[q-1]
+
+ if q == end_index:
+ self.allOutputs[q]['next'] = None
+ else:
+ self.allOutputs[q]['next'] = self.allOutputs[q+1]
+
+ del self.redshiftOutputs
+ del self.timeOutputs
+
+ if end_index is None:
+ end_index = len(self.allOutputs)-1
+
+ self.allOutputs = self.allOutputs[start_index:end_index+1]
+
+ def _ReadEnzoParameterFile(self):
+ "Reads an Enzo parameter file looking for cosmology and output parameters."
+ lines = open(self.EnzoParameterFile).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().split("#", 1)[0].split("//", 1)[0]
+ except ValueError:
+ continue
+ if EnzoParameterDict.has_key(param):
+ t = map(EnzoParameterDict[param], vals.split())
+ if len(t) == 1:
+ self.enzoParameters[param] = t[0]
+ else:
+ self.enzoParameters[param] = t
+ elif param.startswith("CosmologyOutputRedshift["):
+ index = param[param.find("[")+1:param.find("]")]
+ self.redshiftOutputs.append({'index':int(index),
+ 'redshift':float(vals)})
+
+ # Add filenames to redshift outputs.
+ for output in self.redshiftOutputs:
+ output["filename"] = "%s/%s%04d/%s%04d" % (self.enzoParameters['GlobalDir'],
+ self.enzoParameters['RedshiftDumpDir'],output['index'],
+ self.enzoParameters['RedshiftDumpName'],output['index'])
+
+ def _SetParameterDefaults(self):
+ "Set some default parameters to avoid problems if they are not in the parameter file."
+ self.enzoParameters['GlobalDir'] = "."
+ self.enzoParameters['RedshiftDumpName'] = "RD"
+ self.enzoParameters['RedshiftDumpDir'] = "RD"
+ self.enzoParameters['DataDumpName'] = "DD"
+ self.enzoParameters['DataDumpDir'] = "DD"
+ self.enzoParameters['ComovingCoordinates'] = 0
+
+ def imagine_minimal_splice(self, initial_redshift, final_redshift, decimals=3, filename=None,
+ redshift_output_string='CosmologyOutputRedshift', start_index=0):
+ """
+ Create imaginary list of redshift outputs to maximally span a redshift interval.
+ :param decimals (int): The decimal place to which the output redshift will be rounded.
+ If the decimal place in question is nonzero, the redshift will be rounded up to
+ ensure continuity of the splice. Default: 3.
+ :param filename (str): If provided, a file will be written with the redshift outputs in
+ the form in which they should be given in the enzo parameter file. Default: None.
+ :param redshift_output_string (str): The parameter accompanying the redshift outputs in the
+ enzo parameter file. Default: "CosmologyOutputRedshift".
+ :param start_index (int): The index of the first redshift output. Default: 0.
+ """
+
+ z = initial_redshift
+ outputs = []
+
+ while z > final_redshift:
+ rounded = na.round(z, decimals=decimals)
+ if rounded - z < 0:
+ rounded += na.power(10.0,(-1.0*decimals))
+ z = rounded
+ deltaz_max = deltaz_forward(self.cosmology, z, self.enzoParameters['CosmologyComovingBoxSize'])
+ outputs.append({'redshift': z, 'deltazMax': deltaz_max})
+ z -= deltaz_max
+
+ mylog.info("imagine_maximal_splice: Needed %d data dumps to get from z = %f to %f." %
+ (len(outputs), initial_redshift, final_redshift))
+
+ if filename is not None:
+ mylog.info("Writing redshift dump list to %s." % filename)
+ f = open(filename,'w')
+ for q, output in enumerate(outputs):
+ z_string = "%%s[%%d] = %%.%df" % decimals
+ f.write(("%s[%d] = %." + str(decimals) + "f\n") % (redshift_output_string, (q+start_index), output['redshift']))
+ f.close()
+
+ return outputs
+
+ def create_cosmology_splice(self, minimal=True, deltaz_min=0.0, initial_redshift=None, final_redshift=None):
+ """
+ Create list of datasets to be used for LightCones or LightRays.
+ :param minimal (bool): if True, the minimum number of datasets is used to connect the initial and final
+ redshift. If false, the list will contain as many entries as possible within the redshift
+ interval. Default: True.
+ :param deltaz_min (float): specifies the minimum delta z between consecutive datasets in the returned
+ list. Default: 0.0.
+ :param initial_redshift (float): the initial (highest) redshift in the cosmology splice list. If none
+ given, the highest redshift dataset present will be used. Default: None.
+ :param final_redshift (float): the final (lowest) redshift in the cosmology splice list. If none given,
+ the lowest redshift dataset present will be used. Default: None.
+ """
+
+ if initial_redshift is None: initial_redshift = self.InitialRedshift
+ if final_redshift is None: final_redshift = self.FinalRedshift
+
+ # Calculate maximum delta z for each data dump.
+ self._calculate_deltaz_max()
+
+ # Calculate minimum delta z for each data dump.
+ self._calculate_deltaz_min(deltaz_min=deltaz_min)
+
+ cosmology_splice = []
+
+ # Use minimum number of datasets to go from z_i to z_f.
+ if minimal:
+
+ z_Tolerance = 1e-3
+ z = initial_redshift
+
+ # fill redshift space with datasets
+ while ((z > final_redshift) and
+ (na.fabs(z - final_redshift) > z_Tolerance)):
+
+ # For first data dump, choose closest to desired redshift.
+ if (len(cosmology_splice) == 0):
+ # Sort data outputs by proximity to current redsfhit.
+ self.allOutputs.sort(key=lambda obj:na.fabs(z - obj['redshift']))
+ cosmology_splice.append(self.allOutputs[0])
+
+ # Move forward from last slice in stack until z > z_max.
+ else:
+ current_slice = cosmology_splice[-1]
+ while current_slice['next'] is not None and \
+ (z < current_slice['next']['redshift'] or \
+ na.abs(z - current_slice['next']['redshift']) < z_Tolerance):
+ current_slice = current_slice['next']
+
+ if current_slice is cosmology_splice[-1]:
+ final_redshift = cosmology_splice[-1]['redshift'] - cosmology_splice[-1]['deltazMax']
+ mylog.error("Cosmology splice incomplete due to insufficient data outputs.")
+ break
+ else:
+ cosmology_splice.append(current_slice)
+
+ z = cosmology_splice[-1]['redshift'] - cosmology_splice[-1]['deltazMax']
+
+ # Make light ray using maximum number of datasets (minimum spacing).
+ else:
+ # Sort data outputs by proximity to current redsfhit.
+ self.allOutputs.sort(key=lambda obj:na.fabs(initial_redshift - obj['redshift']))
+ # For first data dump, choose closest to desired redshift.
+ cosmology_splice.append(self.allOutputs[0])
+
+ nextOutput = cosmology_splice[-1]['next']
+ while (nextOutput is not None):
+ if (nextOutput['redshift'] <= final_redshift):
+ break
+ if ((cosmology_splice[-1]['redshift'] - nextOutput['redshift']) > cosmology_splice[-1]['deltazMin']):
+ cosmology_splice.append(nextOutput)
+ nextOutput = nextOutput['next']
+ if (cosmology_splice[-1]['redshift'] - cosmology_splice[-1]['deltazMax']) > final_redshift:
+ mylog.error("Cosmology splice incomplete due to insufficient data outputs.")
+ final_redshift = cosmology_splice[-1]['redshift'] - cosmology_splice[-1]['deltazMax']
+
+ mylog.info("create_cosmology_splice: Used %d data dumps to get from z = %f to %f." %
+ (len(cosmology_splice),initial_redshift,final_redshift))
+
+ return cosmology_splice
+
+ def _calculate_deltaz_max(self):
+ "Calculate delta z that corresponds to full box length going from z to (z - delta z)."
+
+ d_Tolerance = 1e-4
+ max_Iterations = 100
+
+ targetDistance = self.enzoParameters['CosmologyComovingBoxSize']
+
+ for output in self.allOutputs:
+ z = output['redshift']
+
+ # Calculate delta z that corresponds to the length of the box at a given redshift.
+ # Use Newton's method to calculate solution.
+ z1 = z
+ z2 = z1 - 0.1 # just an initial guess
+ distance1 = 0.0
+ iteration = 1
+
+ # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
+ distance2 = self.cosmology.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+
+ while ((na.fabs(distance2-targetDistance)/distance2) > d_Tolerance):
+ m = (distance2 - distance1) / (z2 - z1)
+ z1 = z2
+ distance1 = distance2
+ z2 = ((targetDistance - distance2) / m) + z2
+ distance2 = self.cosmology.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+ iteration += 1
+ if (iteration > max_Iterations):
+ mylog.error("calculate_deltaz_max: Warning - max iterations exceeded for z = %f (delta z = %f)." %
+ (z,na.fabs(z2-z)))
+ break
+ output['deltazMax'] = na.fabs(z2-z)
+
+ def _calculate_deltaz_min(self, deltaz_min=0.0):
+ "Calculate delta z that corresponds to a single top grid pixel going from z to (z - delta z)."
+
+ d_Tolerance = 1e-4
+ max_Iterations = 100
+
+ targetDistance = self.enzoParameters['CosmologyComovingBoxSize'] / self.enzoParameters['TopGridDimensions'][0]
+
+ for output in self.allOutputs:
+ z = output['redshift']
+
+ # Calculate delta z that corresponds to the length of a top grid pixel at a given redshift.
+ # Use Newton's method to calculate solution.
+ z1 = z
+ z2 = z1 - 0.01 # just an initial guess
+ distance1 = 0.0
+ iteration = 1
+
+ # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
+ distance2 = self.cosmology.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+
+ while ((na.fabs(distance2 - targetDistance) / distance2) > d_Tolerance):
+ m = (distance2 - distance1) / (z2 - z1)
+ z1 = z2
+ distance1 = distance2
+ z2 = ((targetDistance - distance2) / m) + z2
+ distance2 = self.cosmology.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+ iteration += 1
+ if (iteration > max_Iterations):
+ mylog.error("calculate_deltaz_max: Warning - max iterations exceeded for z = %f (delta z = %f)." %
+ (z,na.fabs(z2-z)))
+ break
+ # Use this calculation or the absolute minimum specified by the user.
+ output['deltazMin'] = max(na.fabs(z2-z),deltaz_min)
+
+EnzoParameterDict = {"CosmologyOmegaMatterNow": float,
+ "CosmologyOmegaLambdaNow": float,
+ "CosmologyHubbleConstantNow": float,
+ "CosmologyInitialRedshift": float,
+ "CosmologyFinalRedshift": float,
+ "ComovingCoordinates": int,
+ "InitialTime": float,
+ "StopTime": float,
+ "TopGridDimensions": float,
+ "dtDataDump": float,
+ "RedshiftDumpName": str,
+ "RedshiftDumpDir": str,
+ "DataDumpName": str,
+ "DataDumpDir": str,
+ "GlobalDir" : str}
+
+def deltaz_forward(cosmology, z, target_distance):
+ "Calculate deltaz corresponding to moving a comoving distance starting from some redshift."
+
+ d_Tolerance = 1e-4
+ max_Iterations = 100
+
+ # Calculate delta z that corresponds to the length of the box at a given redshift.
+ # Use Newton's method to calculate solution.
+ z1 = z
+ z2 = z1 - 0.1 # just an initial guess
+ distance1 = 0.0
+ iteration = 1
+
+ # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
+ distance2 = cosmology.ComovingRadialDistance(z2,z) * cosmology.HubbleConstantNow / 100.0
+
+ while ((na.fabs(distance2-target_distance)/distance2) > d_Tolerance):
+ m = (distance2 - distance1) / (z2 - z1)
+ z1 = z2
+ distance1 = distance2
+ z2 = ((target_distance - distance2) / m) + z2
+ distance2 = cosmology.ComovingRadialDistance(z2,z) * cosmology.HubbleConstantNow / 100.0
+ iteration += 1
+ if (iteration > max_Iterations):
+ mylog.error("deltaz_forward: Warning - max iterations exceeded for z = %f (delta z = %f)." %
+ (z,na.fabs(z2-z)))
+ break
+ return na.fabs(z2-z)
Modified: branches/yt-1.5/yt/extensions/HaloProfiler.py
==============================================================================
--- branches/yt-1.5/yt/extensions/HaloProfiler.py (original)
+++ branches/yt-1.5/yt/extensions/HaloProfiler.py Tue Nov 3 15:45:17 2009
@@ -27,21 +27,99 @@
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'):
+ "Radial profiling, filtering, and projections for halos in cosmological simulations."
+ 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'], filter_quantities=['id','center']):
+ """
+ Initialize a HaloProfiler object.
+ :param halos (str): "multiple" for profiling more than one halo. In this mode halos are read in
+ from a list or identified with a halo finder. In "single" mode, the one and only halo
+ center is identified automatically as the location of the peak in the density field.
+ Default: "multiple".
+ :param halo_list_file (str): name of file containing the list of halos. The HaloProfiler will
+ look for this file in the data directory. Default: "HopAnalysis.out".
+ :param halo_list_format (str or dict): the format of the halo list file. "yt_hop" for the format
+ given by yt's halo finders. "enzo_hop" for the format written by enzo_hop. This keyword
+ can also be given in the form of a dictionary specifying the column in which various
+ properties can be found. For example, {"id": 0, "center": [1, 2, 3], "mass": 4, "radius": 5}.
+ Default: "yt_hop".
+ :param halo_finder_function (function): If halos is set to multiple and the file given by
+ halo_list_file does not exit, the halo finding function specified here will be called.
+ Default: HaloFinder (yt_hop).
+ :param halo_finder_args (tuple): args given with call to halo finder function. Default: None.
+ :param halo_finder_kwargs (dict): kwargs given with call to halo finder function. Default: None.
+ :param use_density_center (bool): re-center halos before performing profiles with an center of mass
+ weighted by overdensity. This is generally not needed. Default: False.
+ :param density_center_exponent (float): when use_density_center set to True, this specifies the
+ exponent, alpha, such that the halo center calculation is weighted by overdensity^alpha.
+ Default: 1.0.
+ :param use_field_max_center (str): another alternative for halo re-centering by selecting the
+ location of the maximum of the field given by this keyword. This is generally not needed.
+ Default: None.
+ :param halo_radius (float): if no halo radii are provided in the halo list file, this parameter is
+ used to specify the radius out to which radial profiles will be made. This keyword is also
+ used when halos is set to single. Default: 0.1.
+ :param radius_units (str): the units of halo_radius. Default: "1" (code units).
+ :param n_profile_bins (int): the number of bins in the radial profiles. Default: 50.
+ :param profile_output_dir (str): the subdirectory, inside the data directory, in which radial profile
+ output files will be created. The directory will be created if it does not exist.
+ Default: "radial_profiles".
+ :param projection_output_dir (str): the subdirectory, inside the data directory, in which projection
+ output files will be created. The directory will be created if it does not exist.
+ Default: "projections".
+ :param projection_width (float): the width of halo projections. Default: 8.0.
+ :param projection_width_units (str): the units of projection_width. Default: "mpc".
+ :param project_at_level (int or "max"): the maximum refinement level to be included in projections.
+ Default: "max" (maximum level within the dataset).
+ :param velocity_center (list): the method in which the halo bulk velocity is calculated (used for
+ calculation of radial and tangential velocities. Valid options are:
+ - ["bulk", "halo"] (Default): the velocity provided in the halo list
+ - ["bulk", "sphere"]: the bulk velocity of the sphere centered on the halo center.
+ - ["max", field]: the velocity of the cell that is the location of the maximum of the field
+ specified (used only when halos set to single).
+ :param filter_quantities (list): quantities from the original halo list file to be written out in the
+ filtered list file. Default: ['id','center'].
+ """
+
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.filter_quantities = filter_quantities
+ if self.filter_quantities is None: self.filter_quantities = []
+
+ 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,36 +129,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._need_bulk_velocity = False
+ for field in [hp['field'] for hp in self.profile_fields]:
+ if 'Velocity' in field or 'Mach' in field:
+ 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
@@ -91,144 +195,258 @@
# 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."
+
+ # Reset filtered halo list.
+ self.filtered_halos = []
+
+ # 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_prefilter = None
+ virial_prefilter_safety_factor = 0.5
+ all_filter_functions = [hf['function'] for hf in self._halo_filters]
+ virial_filter = VirialFilter in all_filter_functions
+ 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'])
+ # Profile all halos.
+ for halo in self._get_objs('all_halos', round_robin=True):
+
+ # Apply prefilters to avoid profiling unwanted halos.
+ filter_result = True
+ haloQuantities = {}
+ if prefilters is not None:
+ for prefilter in prefilters:
+ if not eval(prefilter):
+ filter_result = False
+ break
+
+ if filter_result and len(self.profile_fields) > 0:
+
+ profile_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))
+ profiledHalo = self._get_halo_profile(halo, profile_filename, virial_filter=virial_filter)
+
+ if profiledHalo is None:
continue
- sphere = self.pf.h.sphere(halo['center'],halo['r_max']/self.pf.units['mpc'])
- if len(sphere._grids) == 0: continue
+ # Apply filter and keep track of the quantities that are returned.
+ for hFilter in self._halo_filters:
+ filter_result, filterQuantities = hFilter['function'](profiledHalo, *hFilter['args'], **hFilter['kwargs'])
- # Set 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 not filter_result: break
+
+ if filterQuantities is not None:
+ haloQuantities.update(filterQuantities)
+
+ if filter_result:
+ for quantity in self.filter_quantities:
+ if halo.has_key(quantity): haloQuantities[quantity] = halo[quantity]
+
+ self.filtered_halos.append(haloQuantities)
+
+ self.filtered_halos = self._mpi_catlist(self.filtered_halos)
+ self.filtered_halos.sort(key = lambda a:a['id'])
- if newProfile:
+ if filename is not None:
+ self._write_filtered_halo_list(filename, **kwargs)
+
+ def _get_halo_profile(self, halo, filename, virial_filter=True, 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.
+ """
+
+ # 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
+
+ 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" % halo['id'])
+ profile.write_out(filename, format='%0.6e')
+ elif force_write:
+ mylog.info("Re-writing halo %d" % halo['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=True,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'])):
@@ -238,191 +456,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()
- center = [float(onLine[4]),float(onLine[5]),float(onLine[6])]
- r_max = self.haloRadius * self.pf.units['mpc']
- halo = {'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()
@@ -434,12 +619,14 @@
fields.pop(0)
profile = {}
+ profile_obj = FakeProfile(self.pf)
for field in fields:
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.
@@ -452,138 +639,94 @@
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])
+
+ profile_obj._data = profile
+
if len(profile[fields[0]]) > 1:
- return profile
+ return profile_obj
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])]
+ halo_fields = []
+ for halo_field in self.filter_quantities:
+ if halo_field in fields:
+ fields.remove(halo_field)
+ halo_fields.append(halo_field)
+ # Make it so number of fields in header is same as number of data columns.
+ header_fields = []
+ for halo_field in halo_fields:
+ if isinstance(self.filtered_halos[0][halo_field], types.ListType):
+ header_fields.extend(["%s[%d]" % (halo_field, q)
+ for q in range(len(self.filtered_halos[0][halo_field]))])
+ else:
+ header_fields.append(halo_field)
+ file.write("# ")
+ file.write("\t".join(header_fields + fields + ["\n"]))
+
+ for halo in self.filtered_halos:
+ for halo_field in halo_fields:
+ if isinstance(halo[halo_field], types.ListType):
+ field_data = na.array(halo[halo_field])
+ field_data.tofile(file, sep="\t", format=format)
+ else:
+ if halo_field == 'id':
+ file.write("%04d" % halo[halo_field])
+ else:
+ file.write("%s" % halo[halo_field])
+ file.write("\t")
+ 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.
@@ -596,7 +739,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
@@ -652,18 +795,32 @@
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
+
+class FakeProfile(lagos.ParallelAnalysisInterface):
+ """
+ This is used to mimic a profile object when reading profile data from disk.
+ """
+ def __init__(self, pf):
+ self.pf = pf
+ self._data = {}
+
+ def __getitem__(self, key):
+ return self._data[key]
+
+ def keys(self):
+ return self._data.keys()
Modified: branches/yt-1.5/yt/extensions/lightcone/HaloMask.py
==============================================================================
--- branches/yt-1.5/yt/extensions/lightcone/HaloMask.py (original)
+++ branches/yt-1.5/yt/extensions/lightcone/HaloMask.py Tue Nov 3 15:45:17 2009
@@ -25,81 +25,100 @@
from yt.extensions.HaloProfiler import *
from yt.logger import lagosLogger as mylog
+from yt.config import ytcfg
import yt.lagos as lagos
import copy
import numpy as na
import h5py
-#### Note: assumption of box width 1. I'll fix it someday.
-
-def MakeLightConeHaloMask(lightCone,HaloMaskParameterFile,cube_file=None,mask_file=None,**kwargs):
+def light_cone_halo_mask(lightCone, cube_file=None, mask_file=None, **kwargs):
"Make a boolean mask to cut clusters out of light cone projections."
- pixels = int(lightCone.lightConeParameters['FieldOfViewInArcMinutes'] * 60.0 / \
- lightCone.lightConeParameters['ImageResolutionInArcSeconds'])
+ pixels = int(lightCone.field_of_view_in_arcminutes * 60.0 /
+ lightCone.image_resolution_in_arcseconds)
- lightConeMask = []
+ light_cone_mask = []
# Loop through files in light cone solution and get virial quantities.
- for slice in lightCone.lightConeSolution:
- hp = HaloProfiler(slice['filename'],HaloMaskParameterFile,**kwargs)
- hp._LoadVirialData()
-
- lightConeMask.append(_MakeSliceMask(slice,hp.virialQuantities,pixels))
+ for slice in lightCone.light_cone_solution:
+ halo_list = _get_halo_list(slice['filename'], **kwargs)
+ light_cone_mask.append(_make_slice_mask(slice, halo_list, pixels))
# Write out cube of masks from each slice.
- if cube_file is not None:
+ if cube_file is not None and ytcfg.getint("yt", "__parallel_rank") == 0:
mylog.info("Saving halo mask cube to %s." % cube_file)
- output = h5py.File(cube_file,'a')
- output.create_dataset('haloMaskCube',data=na.array(lightConeMask))
+ output = h5py.File(cube_file, 'a')
+ output.create_dataset('haloMaskCube', data=na.array(light_cone_mask))
output.close()
# Write out final mask.
- if mask_file is not None:
+ if mask_file is not None and ytcfg.getint("yt", "__parallel_rank") == 0:
# Final mask is simply the product of the mask from each slice.
mylog.info("Saving halo mask to %s." % mask_file)
- finalMask = na.ones(shape=(pixels,pixels))
- for mask in lightConeMask:
+ finalMask = na.ones(shape=(pixels, pixels))
+ for mask in light_cone_mask:
finalMask *= mask
- output = h5py.File(mask_file,'a')
- output.create_dataset('haloMask',data=na.array(finalMask))
+ output = h5py.File(mask_file, 'a')
+ output.create_dataset('HaloMask', data=na.array(finalMask))
output.close()
- return lightConeMask
+ return light_cone_mask
-def MakeLightConeHaloMap(lightCone,HaloMaskParameterFile,map_file='halo_map.dat',**kwargs):
+def light_cone_halo_map(lightCone, map_file='halo_map.out', **kwargs):
"Make a text list of location of halos in a light cone image with virial quantities."
haloMap = []
# Loop through files in light cone solution and get virial quantities.
- for slice in lightCone.lightConeSolution:
- hp = HaloProfiler(slice['filename'],HaloMaskParameterFile,**kwargs)
- hp._LoadVirialData()
-
- haloMap.extend(_MakeSliceHaloMap(slice,hp.virialQuantities))
+ for slice in lightCone.light_cone_solution:
+ halo_list = _get_halo_list(slice['filename'], **kwargs)
+ haloMap.extend(_make_slice_halo_map(slice, halo_list))
# Write out file.
- mylog.info("Saving halo map to %s." % map_file)
- f = open(map_file,'w')
- f.write("#z x y M [Msun] R [Mpc] R [image]\n")
- for halo in haloMap:
- f.write("%7.4f %9.6f %9.6f %9.3e %9.3e %9.3e\n" % \
- (halo['redshift'],halo['x'],halo['y'],
- halo['mass'],halo['radiusMpc'],halo['radiusImage']))
- f.close()
+ if ytcfg.getint("yt", "__parallel_rank") == 0:
+ mylog.info("Saving halo map to %s." % map_file)
+ f = open(map_file, 'w')
+ f.write("#z x y M [Msun] R [Mpc] R [image]\n")
+ for halo in haloMap:
+ f.write("%7.4f %9.6f %9.6f %9.3e %9.3e %9.3e\n" % \
+ (halo['redshift'], halo['x'], halo['y'],
+ halo['mass'], halo['radiusMpc'], halo['radiusImage']))
+ f.close()
+
+def _get_halo_list(dataset, halo_profiler_kwargs=None, halo_profiler_actions=None, halo_list='all'):
+ "Load a list of halos for the dataset."
+
+ if halo_profiler_kwargs is None: halo_profiler_kwargs = {}
+ if halo_profiler_actions is None: halo_profiler_actions = []
+
+ hp = HaloProfiler(dataset, **halo_profiler_kwargs)
+ for action in halo_profiler_actions:
+ if not action.has_key('args'): action['args'] = ()
+ if not action.has_key('kwargs'): action['kwargs'] = {}
+ action['function'](hp, *action['args'], **action['kwargs'])
+
+ if halo_list == 'all':
+ return_list = copy.deepcopy(hp.all_halos)
+ elif halo_list == 'filtered':
+ return_list = copy.deepcopy(hp.filtered_halos)
+ else:
+ mylog.error("Keyword, halo_list, must be either 'all' or 'filtered'.")
+ return_list = None
+
+ del hp
+ return return_list
-def _MakeSliceMask(slice,virialQuantities,pixels):
+def _make_slice_mask(slice, halo_list, pixels):
"Make halo mask for one slice in light cone solution."
# Get shifted, tiled halo list.
- all_halo_x,all_halo_y,all_halo_radius,all_halo_mass = _MakeSliceHaloList(slice,virialQuantities)
+ all_halo_x, all_halo_y, all_halo_radius, all_halo_mass = _make_slice_halo_list(slice, halo_list)
# Make boolean mask and cut out halos.
dx = slice['WidthBoxFraction'] / pixels
x = [(q + 0.5) * dx for q in range(pixels)]
- haloMask = na.ones(shape=(pixels,pixels),dtype=bool)
+ haloMask = na.ones(shape=(pixels, pixels), dtype=bool)
# Cut out any pixel that has any part at all in the circle.
for q in range(len(all_halo_radius)):
@@ -114,7 +133,7 @@
return haloMask
-def _MakeSliceHaloMap(slice,virialQuantities):
+def _make_slice_halo_map(slice, halo_list):
"Make list of halos for one slice in light cone solution."
# Get units to convert virial radii back to physical units.
@@ -123,7 +142,7 @@
del dataset_object
# Get shifted, tiled halo list.
- all_halo_x,all_halo_y,all_halo_radius,all_halo_mass = _MakeSliceHaloList(slice,virialQuantities)
+ all_halo_x, all_halo_y, all_halo_radius, all_halo_mass = _make_slice_halo_list(slice, halo_list)
# Construct list of halos
haloMap = []
@@ -142,7 +161,7 @@
return haloMap
-def _MakeSliceHaloList(slice,virialQuantities):
+def _make_slice_halo_list(slice, halo_list):
"Make shifted, tiled list of halos for halo mask and halo map."
# Make numpy arrays for halo centers and virial radii.
@@ -157,7 +176,7 @@
Mpc_units = dataset_object.units['mpc']
del dataset_object
- for halo in virialQuantities:
+ for halo in halo_list:
if halo is not None:
center = copy.deepcopy(halo['center'])
halo_depth.append(center.pop(slice['ProjectionAxis']))
@@ -182,13 +201,13 @@
add_left = (halo_depth + halo_radius) > 1 # should be box width
add_right = (halo_depth - halo_radius) < 0
- halo_depth = na.concatenate([halo_depth,(halo_depth[add_left]-1),(halo_depth[add_right]+1)])
- halo_x = na.concatenate([halo_x,halo_x[add_left],halo_x[add_right]])
- halo_y = na.concatenate([halo_y,halo_y[add_left],halo_y[add_right]])
- halo_radius = na.concatenate([halo_radius,halo_radius[add_left],halo_radius[add_right]])
- halo_mass = na.concatenate([halo_mass,halo_mass[add_left],halo_mass[add_right]])
+ halo_depth = na.concatenate([halo_depth, (halo_depth[add_left]-1), (halo_depth[add_right]+1)])
+ halo_x = na.concatenate([halo_x, halo_x[add_left], halo_x[add_right]])
+ halo_y = na.concatenate([halo_y, halo_y[add_left], halo_y[add_right]])
+ halo_radius = na.concatenate([halo_radius, halo_radius[add_left], halo_radius[add_right]])
+ halo_mass = na.concatenate([halo_mass, halo_mass[add_left], halo_mass[add_right]])
- del add_left,add_right
+ del add_left, add_right
# Cut out the halos outside the region of interest.
if (slice['DepthBoxFraction'] < 1):
@@ -217,12 +236,12 @@
# Copy original into offset positions to make tiles.
for x in range(int(na.ceil(slice['WidthBoxFraction']))):
for y in range(int(na.ceil(slice['WidthBoxFraction']))):
- all_halo_x = na.concatenate([all_halo_x,halo_x+x])
- all_halo_y = na.concatenate([all_halo_y,halo_y+y])
- all_halo_radius = na.concatenate([all_halo_radius,halo_radius])
- all_halo_mass = na.concatenate([all_halo_mass,halo_mass])
+ all_halo_x = na.concatenate([all_halo_x, halo_x+x])
+ all_halo_y = na.concatenate([all_halo_y, halo_y+y])
+ all_halo_radius = na.concatenate([all_halo_radius, halo_radius])
+ all_halo_mass = na.concatenate([all_halo_mass, halo_mass])
- del halo_x,halo_y,halo_radius,halo_mass
+ del halo_x, halo_y, halo_radius, halo_mass
# Shift centers laterally.
offset = copy.deepcopy(slice['ProjectionCenter'])
@@ -276,17 +295,17 @@
del add_y_left
# Add the hanging centers back to the projection data.
- all_halo_x = na.concatenate([all_halo_x,add_x_halo_x,add2_x_halo_x,add_y_halo_x,add2_y_halo_x])
- all_halo_y = na.concatenate([all_halo_y,add_x_halo_y,add2_x_halo_y,add_y_halo_y,add2_y_halo_y])
- all_halo_radius = na.concatenate([all_halo_radius,add_x_halo_radius,add2_x_halo_radius,
- add_y_halo_radius,add2_y_halo_radius])
- all_halo_mass = na.concatenate([all_halo_mass,add_x_halo_mass,add2_x_halo_mass,
- add_y_halo_mass,add2_y_halo_mass])
-
- del add_x_halo_x,add_x_halo_y,add_x_halo_radius
- del add2_x_halo_x,add2_x_halo_y,add2_x_halo_radius
- del add_y_halo_x,add_y_halo_y,add_y_halo_radius
- del add2_y_halo_x,add2_y_halo_y,add2_y_halo_radius
+ all_halo_x = na.concatenate([all_halo_x, add_x_halo_x, add2_x_halo_x, add_y_halo_x, add2_y_halo_x])
+ all_halo_y = na.concatenate([all_halo_y, add_x_halo_y, add2_x_halo_y, add_y_halo_y, add2_y_halo_y])
+ all_halo_radius = na.concatenate([all_halo_radius, add_x_halo_radius, add2_x_halo_radius,
+ add_y_halo_radius, add2_y_halo_radius])
+ all_halo_mass = na.concatenate([all_halo_mass, add_x_halo_mass, add2_x_halo_mass,
+ add_y_halo_mass, add2_y_halo_mass])
+
+ del add_x_halo_x, add_x_halo_y, add_x_halo_radius
+ del add2_x_halo_x, add2_x_halo_y, add2_x_halo_radius
+ del add_y_halo_x, add_y_halo_y, add_y_halo_radius
+ del add2_y_halo_x, add2_y_halo_y, add2_y_halo_radius
# Cut edges to proper width.
cut_mask = (all_halo_x - all_halo_radius < slice['WidthBoxFraction']) & \
@@ -297,4 +316,4 @@
all_halo_mass = all_halo_mass[cut_mask]
del cut_mask
- return (all_halo_x,all_halo_y,all_halo_radius,all_halo_mass)
+ return (all_halo_x, all_halo_y, all_halo_radius, all_halo_mass)
Modified: branches/yt-1.5/yt/extensions/lightcone/LightCone.py
==============================================================================
--- branches/yt-1.5/yt/extensions/lightcone/LightCone.py (original)
+++ branches/yt-1.5/yt/extensions/lightcone/LightCone.py Tue Nov 3 15:45:17 2009
@@ -24,6 +24,7 @@
"""
from yt.extensions.lightcone import *
+from yt.extensions.EnzoSimulation import *
from yt.logger import lagosLogger as mylog
from yt.config import ytcfg
from yt.funcs import *
@@ -32,24 +33,50 @@
import copy
import os
import numpy as na
-import random as rand
-import yt.lagos.Cosmology as cosmo
-class LightCone(object):
- def __init__(self,EnzoParameterFile,LightConeParameterFile,verbose=True):
- self.verbose = verbose
- self.EnzoParameterFile = EnzoParameterFile
- self.LightConeParameterFile = LightConeParameterFile
- self.enzoParameters = {}
- self.lightConeParameters = {}
- self.redshiftOutputs = []
- self.timeOutputs = []
- self.allOutputs = []
- self.lightConeSolution = []
- self.masterSolution = [] # kept to compare with recycled solutions
- self.projectionStack = []
- self.projectionWeightFieldStack = []
- self.haloMask = []
+class LightCone(EnzoSimulation):
+ def __init__(self, EnzoParameterFile, initial_redshift=1.0, final_redshift=0.0, observer_redshift=0.0,
+ field_of_view_in_arcminutes=600.0, image_resolution_in_arcseconds=60.0,
+ use_minimum_datasets=True, deltaz_min=0.0, minimum_coherent_box_fraction=0.0,
+ output_dir='LC', output_prefix='LightCone'):
+ """
+ Initialize a LightCone object.
+ :param initial_redshift (float): the initial (highest) redshift for the light cone. Default: 1.0.
+ :param final_redshift (float): the final (lowest) redshift for the light cone. Default: 0.0.
+ :param observer_redshift (float): the redshift of the observer. Default: 0.0.
+ :param field_of_view_in_arcminutes (float): the field of view of the image in units of arcminutes.
+ Default: 600.0.
+ :param image_resolution_in_arcseconds (float): the size of each image pixel in units of arcseconds.
+ Default: 60.0.
+ :param use_minimum_datasets (bool): if True, the minimum number of datasets is used to connect the
+ initial and final redshift. If false, the light cone solution will contain as many entries
+ as possible within the redshift interval. Default: True.
+ :param deltaz_min (float): specifies the minimum :math:`\Delta z` between consecutive datasets in
+ the returned list. Default: 0.0.
+ :param minimum_coherent_box_fraction (float): used with use_minimum_datasets set to False, this
+ parameter specifies the fraction of the total box size to be traversed before rerandomizing
+ the projection axis and center. This was invented to allow light cones with thin slices to
+ sample coherent large scale structure, but in practice does not work so well. Try setting
+ this parameter to 1 and see what happens. Default: 0.0.
+ :param output_dir (str): the directory in which images and data files will be written. Default: 'LC'.
+ :param output_prefix (str): the prefix of all images and data files. Default: 'LightCone'.
+ """
+
+ self.initial_redshift = initial_redshift
+ self.final_redshift = final_redshift
+ self.observer_redshift = observer_redshift
+ self.field_of_view_in_arcminutes = field_of_view_in_arcminutes
+ self.image_resolution_in_arcseconds = image_resolution_in_arcseconds
+ self.use_minimum_datasets = use_minimum_datasets
+ self.deltaz_min = deltaz_min
+ self.minimum_coherent_box_fraction = minimum_coherent_box_fraction
+ self.output_dir = output_dir
+ self.output_prefix = output_prefix
+
+ self.master_solution = [] # kept to compare with recycled solutions
+ self.projection_stack = []
+ self.projection_weight_field_stack = []
+ self.halo_mask = []
# Original random seed of the first solution.
self.originalRandomSeed = 0
@@ -58,231 +85,189 @@
self.recycleSolution = False
self.recycleRandomSeed = 0
- # Set some parameter defaults.
- self._SetParameterDefaults()
-
- # Read parameters.
- self._ReadLightConeParameterFile()
- self._ReadEnzoParameterFile()
+ # Initialize EnzoSimulation machinery for getting dataset list.
+ EnzoSimulation.__init__(self, EnzoParameterFile, initial_redshift=self.initial_redshift,
+ final_redshift=self.final_redshift, links=True,
+ enzo_parameters={'CosmologyComovingBoxSize':float})
# Calculate number of pixels.
- self.pixels = int(self.lightConeParameters['FieldOfViewInArcMinutes'] * 60.0 / \
- self.lightConeParameters['ImageResolutionInArcSeconds'])
+ self.pixels = int(self.field_of_view_in_arcminutes * 60.0 / \
+ self.image_resolution_in_arcseconds)
- if ytcfg.getint("yt","__parallel_rank") == 0:
+ if ytcfg.getint("yt", "__parallel_rank") == 0:
# Create output directory.
- if (os.path.exists(self.lightConeParameters['OutputDir'])):
- if not(os.path.isdir(self.lightConeParameters['OutputDir'])):
- mylog.error("Output directory exists, but is not a directory: %s." % self.lightConeParameters['OutputDir'])
- self.lightConeParameters['OutputDir'] = './'
- else:
- os.mkdir(self.lightConeParameters['OutputDir'])
-
- # Calculate redshifts for dt data dumps.
- if (self.enzoParameters.has_key('dtDataDump')):
- if self.lightConeParameters['dtSplit'] is None:
- self._CalculateTimeDumpRedshifts()
- else:
- self._CalculateSplitTimeDumpRedshifts(self.lightConeParameters['dtSplit'],self.lightConeParameters['dt_new'])
- # Combine all data dumps.
- self._CombineDataOutputs()
-
- # Calculate maximum delta z for each data dump.
- self._CalculateDeltaZMax()
-
- if not self.lightConeParameters['UseMinimumNumberOfProjections']:
- # Calculate minimum delta z for each data dump.
- self._CalculateDeltaZMin()
+ if (os.path.exists(self.output_dir)):
+ if not(os.path.isdir(self.output_dir)):
+ mylog.error("Output directory exists, but is not a directory: %s." % self.output_dir)
+ self.output_dir = './'
+ else:
+ os.mkdir(self.output_dir)
+
+ # Get list of datasets for light cone solution.
+ self.light_cone_solution = self.create_cosmology_splice(minimal=self.use_minimum_datasets,
+ deltaz_min=self.deltaz_min)
- def CalculateLightConeSolution(self,seed=None):
- "Create list of projections to be added together to make the light cone."
+ def calculate_light_cone_solution(self, seed=None, filename=None):
+ """
+ Create list of projections to be added together to make the light cone.
+ :param seed (int): the seed for the random number generator. Any light cone solution
+ can be reproduced by giving the same random seed. Default: None (each solution
+ will be distinct).
+ :param filename (str): if given, a text file detailing the solution will be written out. Default: None.
+ """
# Don't use box coherence with maximum projection depths.
- if self.lightConeParameters['UseMinimumNumberOfProjections'] and \
- self.lightConeParameters['MinimumCoherentBoxFraction'] > 0:
- mylog.info("Setting MinimumCoherentBoxFraction to 0 with minimal light cone.")
- self.lightConeParameters['MinimumCoherentBoxFraction'] = 0
+ if self.use_minimum_datasets and \
+ self.minimum_coherent_box_fraction > 0:
+ mylog.info("Setting minimum_coherent_box_fraction to 0 with minimal light cone.")
+ self.minimum_coherent_box_fraction = 0
# Make sure recycling flag is off.
self.recycleSolution = False
# Get rid of old halo mask, if one was there.
- self.haloMask = []
+ self.halo_mask = []
if seed is not None:
self.originalRandomSeed = int(seed)
- # Make light cone using minimum number of projections.
- if (self.lightConeParameters['UseMinimumNumberOfProjections']):
-
- z_Tolerance = 1e-4
- z = self.lightConeParameters['InitialRedshift']
-
- # fill redshift space with projections
- while ((z > self.lightConeParameters['FinalRedshift']) and
- (na.fabs(z - self.lightConeParameters['FinalRedshift']) > z_Tolerance)):
- # Sort data outputs by proximity to current redsfhit.
- self.allOutputs.sort(key=lambda obj:na.fabs(z - obj['redshift']))
- # For first data dump, choose closest to desired redshift.
- if (len(self.lightConeSolution) == 0):
- self.lightConeSolution.append(self.allOutputs[0])
- # Start with data dump closest to desired redshift and move backward
- # until one is within max delta z of last output in solution list.
- else:
- output = self.allOutputs[0]
- while (z > output['redshift']):
- output = output['previous']
- if (output is None):
- if self.verbose: mylog.error("CalculateLightConeSolution: search for data output went off the end the stack.")
- if self.verbose: mylog.error("Could not calculate light cone solution.")
- return
- if (output['redshift'] == self.lightConeSolution[-1]['redshift']):
- if self.verbose: mylog.error("CalculateLightConeSolution: No data dump between z = %f and %f." % \
- ((self.lightConeSolution[-1]['redshift'] - self.lightConeSolution[-1]['deltazMax']),
- self.lightConeSolution[-1]['redshift']))
- if self.verbose: mylog.error("Could not calculate light cone solution.")
- return
- self.lightConeSolution.append(output)
- z = self.lightConeSolution[-1]['redshift'] - self.lightConeSolution[-1]['deltazMax']
-
- # Make light cone using maximum number of projections (minimum spacing).
- else:
- # Sort data outputs by proximity to current redsfhit.
- self.allOutputs.sort(key=lambda obj:na.fabs(self.lightConeParameters['InitialRedshift'] - obj['redshift']))
- # For first data dump, choose closest to desired redshift.
- self.lightConeSolution.append(self.allOutputs[0])
-
- nextOutput = self.lightConeSolution[-1]['next']
- while (nextOutput is not None):
- if (nextOutput['redshift'] <= self.lightConeParameters['FinalRedshift']):
- break
- if ((self.lightConeSolution[-1]['redshift'] - nextOutput['redshift']) > self.lightConeSolution[-1]['deltazMin']):
- self.lightConeSolution.append(nextOutput)
- nextOutput = nextOutput['next']
-
- if self.verbose: mylog.info("CalculateLightConeSolution: Used %d data dumps to get from z = %f to %f." %
- (len(self.lightConeSolution),
- self.lightConeParameters['InitialRedshift'],
- self.lightConeParameters['FinalRedshift']))
-
# Calculate projection sizes, and get random projection axes and centers.
- rand.seed(self.originalRandomSeed)
- co = lagos.Cosmology(HubbleConstantNow = (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
- OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
- OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'])
+ na.random.seed(self.originalRandomSeed)
# For box coherence, keep track of effective depth travelled.
boxFractionUsed = 0.0
- for q in range(len(self.lightConeSolution)):
- del self.lightConeSolution[q]['previous']
- del self.lightConeSolution[q]['next']
- if (q == len(self.lightConeSolution) - 1):
- z_next = self.lightConeParameters['FinalRedshift']
+ for q in range(len(self.light_cone_solution)):
+ del self.light_cone_solution[q]['previous']
+ del self.light_cone_solution[q]['next']
+ if (q == len(self.light_cone_solution) - 1):
+ z_next = self.final_redshift
else:
- z_next = self.lightConeSolution[q+1]['redshift']
+ z_next = self.light_cone_solution[q+1]['redshift']
# Calculate fraction of box required for a depth of delta z
- self.lightConeSolution[q]['DepthBoxFraction'] = co.ComovingRadialDistance(z_next,self.lightConeSolution[q]['redshift']) * \
+ self.light_cone_solution[q]['DepthBoxFraction'] = self.cosmology.ComovingRadialDistance(z_next, self.light_cone_solution[q]['redshift']) * \
self.enzoParameters['CosmologyHubbleConstantNow'] / self.enzoParameters['CosmologyComovingBoxSize']
# Simple error check to make sure more than 100% of box depth is never required.
- if (self.lightConeSolution[q]['DepthBoxFraction'] > 1.0):
- if self.verbose: mylog.error("Warning: box fraction required to go from z = %f to %f is %f" %
- (self.lightConeSolution[q]['redshift'],z_next,
- self.lightConeSolution[q]['DepthBoxFraction']))
- if self.verbose: mylog.error("Full box delta z is %f, but it is %f to the next data dump." %
- (self.lightConeSolution[q]['deltazMax'],
- self.lightConeSolution[q]['redshift']-z_next))
+ if (self.light_cone_solution[q]['DepthBoxFraction'] > 1.0):
+ mylog.debug("Warning: box fraction required to go from z = %f to %f is %f" %
+ (self.light_cone_solution[q]['redshift'], z_next,
+ self.light_cone_solution[q]['DepthBoxFraction']))
+ mylog.debug("Full box delta z is %f, but it is %f to the next data dump." %
+ (self.light_cone_solution[q]['deltazMax'],
+ self.light_cone_solution[q]['redshift']-z_next))
# Calculate fraction of box required for width corresponding to requested image size.
- scale = co.AngularScale_1arcsec_kpc(self.lightConeParameters['ObserverRedshift'],self.lightConeSolution[q]['redshift'])
- size = self.lightConeParameters['FieldOfViewInArcMinutes'] * 60.0 * scale / 1000.0
+ scale = self.cosmology.AngularScale_1arcsec_kpc(self.observer_redshift, self.light_cone_solution[q]['redshift'])
+ size = self.field_of_view_in_arcminutes * 60.0 * scale / 1000.0
boxSizeProper = self.enzoParameters['CosmologyComovingBoxSize'] / (self.enzoParameters['CosmologyHubbleConstantNow'] *
- (1.0 + self.lightConeSolution[q]['redshift']))
- self.lightConeSolution[q]['WidthBoxFraction'] = size / boxSizeProper
+ (1.0 + self.light_cone_solution[q]['redshift']))
+ self.light_cone_solution[q]['WidthBoxFraction'] = size / boxSizeProper
# Get projection axis and center.
# If using box coherence, only get random axis and center if enough of the box has been used,
# or if boxFractionUsed will be greater than 1 after this slice.
- if (q == 0) or (self.lightConeParameters['MinimumCoherentBoxFraction'] == 0) or \
- (boxFractionUsed > self.lightConeParameters['MinimumCoherentBoxFraction']) or \
- (boxFractionUsed + self.lightConeSolution[q]['DepthBoxFraction'] > 1.0):
+ if (q == 0) or (self.minimum_coherent_box_fraction == 0) or \
+ (boxFractionUsed > self.minimum_coherent_box_fraction) or \
+ (boxFractionUsed + self.light_cone_solution[q]['DepthBoxFraction'] > 1.0):
# Random axis and center.
- self.lightConeSolution[q]['ProjectionAxis'] = rand.randint(0,2)
- self.lightConeSolution[q]['ProjectionCenter'] = [rand.random(),rand.random(),rand.random()]
+ self.light_cone_solution[q]['ProjectionAxis'] = na.random.randint(0, 3)
+ self.light_cone_solution[q]['ProjectionCenter'] = [na.random.random(), na.random.random(), na.random.random()]
boxFractionUsed = 0.0
else:
# Same axis and center as previous slice, but with depth center shifted.
- self.lightConeSolution[q]['ProjectionAxis'] = self.lightConeSolution[q-1]['ProjectionAxis']
- self.lightConeSolution[q]['ProjectionCenter'] = copy.deepcopy(self.lightConeSolution[q-1]['ProjectionCenter'])
- self.lightConeSolution[q]['ProjectionCenter'][self.lightConeSolution[q]['ProjectionAxis']] += \
- 0.5 * (self.lightConeSolution[q]['DepthBoxFraction'] + self.lightConeSolution[q-1]['DepthBoxFraction'])
- if self.lightConeSolution[q]['ProjectionCenter'][self.lightConeSolution[q]['ProjectionAxis']] >= 1.0:
- self.lightConeSolution[q]['ProjectionCenter'][self.lightConeSolution[q]['ProjectionAxis']] -= 1.0
+ self.light_cone_solution[q]['ProjectionAxis'] = self.light_cone_solution[q-1]['ProjectionAxis']
+ self.light_cone_solution[q]['ProjectionCenter'] = copy.deepcopy(self.light_cone_solution[q-1]['ProjectionCenter'])
+ self.light_cone_solution[q]['ProjectionCenter'][self.light_cone_solution[q]['ProjectionAxis']] += \
+ 0.5 * (self.light_cone_solution[q]['DepthBoxFraction'] + self.light_cone_solution[q-1]['DepthBoxFraction'])
+ if self.light_cone_solution[q]['ProjectionCenter'][self.light_cone_solution[q]['ProjectionAxis']] >= 1.0:
+ self.light_cone_solution[q]['ProjectionCenter'][self.light_cone_solution[q]['ProjectionAxis']] -= 1.0
- boxFractionUsed += self.lightConeSolution[q]['DepthBoxFraction']
+ boxFractionUsed += self.light_cone_solution[q]['DepthBoxFraction']
# Store this as the master solution.
- self.masterSolution = [copy.deepcopy(q) for q in self.lightConeSolution]
+ self.master_solution = [copy.deepcopy(q) for q in self.light_cone_solution]
- # Clear out some stuff.
- del co
+ # Write solution to a file.
+ if filename is not None:
+ self._save_light_cone_solution(filename=filename)
- @rootonly
- def GetHaloMask(self,HaloMaskParameterFile,mask_file=None,map_file=None,**kwargs):
- "Gets a halo mask from a file or makes a new one."
+ def get_halo_mask(self, mask_file=None, map_file=None, **kwargs):
+ """
+ Gets a halo mask from a file or makes a new one.
+ :param mask_file (str): specify an hdf5 file to output the halo mask.
+ :param map_file (str): specify a text file to output the halo map
+ (locations in image of halos).
+ """
# Get halo map if map_file given.
- if (map_file is not None):
- MakeLightConeHaloMap(self,HaloMaskParameterFile,map_file=map_file,**kwargs)
+ if map_file is not None and not os.path.exists(map_file):
+ light_cone_halo_map(self, map_file=map_file, **kwargs)
# Check if file already exists.
- if (mask_file is not None) and os.path.exists(mask_file):
- input = h5.openFile(mask_file,'r')
- self.haloMask = input.root.haloMask.read()
+ if mask_file is not None and os.path.exists(mask_file):
+ mylog.info('Reading halo mask from %s.' % mask_file)
+ input = h5py.File(mask_file, 'r')
+ self.halo_mask = input['HaloMask'].value
input.close()
# Otherwise, make a halo mask.
else:
- haloMaskCube = MakeLightConeHaloMask(self,HaloMaskParameterFile,mask_file=mask_file,**kwargs)
+ halo_mask_cube = light_cone_halo_mask(self, mask_file=mask_file, **kwargs)
# Collapse cube into final mask.
- self.haloMask = na.ones(shape=(self.pixels,self.pixels),dtype=bool)
- for mask in haloMaskCube:
- self.haloMask *= mask
- del haloMaskCube
-
- def ProjectLightCone(self,field,weight_field=None,apply_halo_mask=False,node=None,
- save_stack=True,save_slice_images=False,flatten_stack=False,photon_field=False,**kwargs):
- "Create projections for light cone, then add them together."
+ if ytcfg.getint("yt", "__parallel_rank") == 0:
+ self.halo_mask = na.ones(shape=(self.pixels, self.pixels), dtype=bool)
+ for mask in halo_mask_cube:
+ self.halo_mask *= mask
+ del halo_mask_cube
+
+ def project_light_cone(self, field, weight_field=None, apply_halo_mask=False, node=None,
+ save_stack=True, save_slice_images=False, flatten_stack=False, photon_field=False, **kwargs):
+ """
+ Create projections for light cone, then add them together.
+ :param weight_field (str): the weight field of the projection. This has the same meaning as in standard
+ projections. Default: None.
+ :param apply_halo_mask (bool): if True, a boolean mask is apply to the light cone projection. See below for a
+ description of halo masks. Default: False.
+ :param node (str): a prefix to be prepended to the node name under which the projection data is serialized.
+ Default: None.
+ :param save_stack (bool): if True, the unflatted light cone data including each individual slice is written to
+ an hdf5 file. Default: True.
+ :param save_slice_images (bool): save images for each individual projection slice. Default: False.
+ :param flatten_stack (bool): if True, the light cone stack is continually flattened each time a slice is added
+ in order to save memory. This is generally not necessary. Default: False.
+ :param photon_field (bool): if True, the projection data for each slice is decremented by 4 Pi R^2`, where R
+ is the luminosity distance between the observer and the slice redshift. Default: False.
+ """
# Clear projection stack.
- self.projectionStack = []
- self.projectionWeightFieldStack = []
- if (self.lightConeSolution[-1].has_key('object')):
- del self.lightConeSolution[-1]['object']
+ self.projection_stack = []
+ self.projection_weight_field_stack = []
+ if (self.light_cone_solution[-1].has_key('object')):
+ del self.light_cone_solution[-1]['object']
- if not(self.lightConeParameters['OutputDir'].endswith("/")):
- self.lightConeParameters['OutputDir'] += "/"
+ if not(self.output_dir.endswith("/")):
+ self.output_dir += "/"
- for q,output in enumerate(self.lightConeSolution):
+ for q, output in enumerate(self.light_cone_solution):
if node is None:
- name = "%s%s_%04d_%04d" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'],
- q,len(self.lightConeSolution))
+ name = "%s%s_%04d_%04d" % (self.output_dir, self.output_prefix,
+ q, len(self.light_cone_solution))
else:
- name = "%s%s_%s_%04d_%04d" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'],
- node,q,len(self.lightConeSolution))
+ name = "%s%s_%s_%04d_%04d" % (self.output_dir, self.output_prefix,
+ node, q, len(self.light_cone_solution))
output['object'] = lagos.EnzoStaticOutput(output['filename'])
- frb = LightConeProjection(output,field,self.pixels,weight_field=weight_field,
+ frb = LightConeProjection(output, field, self.pixels, weight_field=weight_field,
save_image=save_slice_images,
- name=name,node=node,**kwargs)
- if ytcfg.getint("yt","__parallel_rank") == 0:
+ name=name, node=node, **kwargs)
+ if ytcfg.getint("yt", "__parallel_rank") == 0:
if photon_field:
- #need to decrement the flux by the luminosity distance. Assume field in frb is in erg/s/cm^2/Hz
+ # Decrement the flux by the luminosity distance. Assume field in frb is in erg/s/cm^2/Hz
co = lagos.Cosmology(HubbleConstantNow = (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'])
- dL = co.LuminosityDistance(self.lightConeParameters['ObserverRedshift'],output['redshift']) #in Mpc
+ dL = self.cosmology.LuminosityDistance(self.observer_redshift, output['redshift']) #in Mpc
boxSizeProper = self.enzoParameters['CosmologyComovingBoxSize'] / (self.enzoParameters['CosmologyHubbleConstantNow'] *
(1.0 + output['redshift']))
pixelarea = (boxSizeProper/self.pixels)**2 #in proper cm^2
@@ -290,44 +275,44 @@
mylog.info("Distance to slice = %e" % dL)
frb[field] *= factor #in erg/s/cm^2/Hz on observer's image plane.
- if ytcfg.getint("yt","__parallel_rank") == 0:
- if (weight_field is not None):
+ if ytcfg.getint("yt", "__parallel_rank") == 0:
+ if weight_field is not None:
# Data come back normalized by the weight field.
# Undo that so it can be added up for the light cone.
- self.projectionStack.append(frb[field]*frb['weight_field'])
- self.projectionWeightFieldStack.append(frb['weight_field'])
+ self.projection_stack.append(frb[field]*frb['weight_field'])
+ self.projection_weight_field_stack.append(frb['weight_field'])
else:
- self.projectionStack.append(frb[field])
+ self.projection_stack.append(frb[field])
# Delete the frb. This saves a decent amount of ram.
- if (q < len(self.lightConeSolution) - 1):
+ if (q < len(self.light_cone_solution) - 1):
del frb
# Flatten stack to save memory.
- if flatten_stack and (len(self.projectionStack) > 1):
- self.projectionStack = [sum(self.projectionStack)]
+ if flatten_stack and (len(self.projection_stack) > 1):
+ self.projection_stack = [sum(self.projection_stack)]
if weight_field is not None:
- self.projectionWeightFieldStack = [sum(self.projectionWeightFieldStack)]
+ self.projection_weight_field_stack = [sum(self.projection_weight_field_stack)]
# Delete the plot collection now that the frb is deleted.
del output['pc']
# Unless this is the last slice, delete the dataset object.
# The last one will be saved to make the plot collection.
- if (q < len(self.lightConeSolution) - 1):
+ if (q < len(self.light_cone_solution) - 1):
del output['object']
- if ytcfg.getint("yt","__parallel_rank") == 0:
+ if ytcfg.getint("yt", "__parallel_rank") == 0:
# Add up slices to make light cone projection.
if (weight_field is None):
- lightConeProjection = sum(self.projectionStack)
+ lightConeProjection = sum(self.projection_stack)
else:
- lightConeProjection = sum(self.projectionStack) / sum(self.projectionWeightFieldStack)
+ lightConeProjection = sum(self.projection_stack) / sum(self.projection_weight_field_stack)
if node is None:
- filename = "%s%s" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
+ filename = "%s%s" % (self.output_dir, self.output_prefix)
else:
- filename = "%s%s_%s" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'],node)
+ filename = "%s%s_%s" % (self.output_dir, self.output_prefix, node)
# Save the last fixed resolution buffer for the plot collection,
# but replace the data with the full light cone projection data.
@@ -335,28 +320,28 @@
# Write stack to hdf5 file.
if save_stack:
- self._SaveLightConeStack(field=field,weight_field=weight_field,filename=filename)
+ self._save_light_cone_stack(field=field, weight_field=weight_field, filename=filename)
# Apply halo mask.
if apply_halo_mask:
- if len(self.haloMask) > 0:
+ if len(self.halo_mask) > 0:
mylog.info("Applying halo mask.")
- frb.data[field] *= self.haloMask
+ frb.data[field] *= self.halo_mask
else:
- mylog.error("No halo mask loaded, call GetHaloMask.")
+ mylog.error("No halo mask loaded, call get_halo_mask.")
# Make a plot collection for the light cone projection.
- center = [0.5 * (self.lightConeSolution[-1]['object'].parameters['DomainLeftEdge'][w] +
- self.lightConeSolution[-1]['object'].parameters['DomainRightEdge'][w])
- for w in range(self.lightConeSolution[-1]['object'].parameters['TopGridRank'])]
- pc = raven.PlotCollection(self.lightConeSolution[-1]['object'],center=center)
- pc.add_fixed_resolution_plot(frb,field,**kwargs)
+ center = [0.5 * (self.light_cone_solution[-1]['object'].parameters['DomainLeftEdge'][w] +
+ self.light_cone_solution[-1]['object'].parameters['DomainRightEdge'][w])
+ for w in range(self.light_cone_solution[-1]['object'].parameters['TopGridRank'])]
+ pc = raven.PlotCollection(self.light_cone_solution[-1]['object'], center=center)
+ pc.add_fixed_resolution_plot(frb, field, **kwargs)
pc.save(filename)
# Return the plot collection so the user can remake the plot if they want.
return pc
- def RerandomizeLightConeSolution(self,newSeed,recycle=True):
+ def rerandomize_light_cone_solution(self, newSeed, recycle=True, filename=None):
"""
When making a projection for a light cone, only randomizations along the line of sight make any
given projection unique, since the lateral shifting and tiling is done after the projection is made.
@@ -368,23 +353,29 @@
This routine has now been updated to be a general solution rescrambler. If the keyword recycle is set to
True, then it will recycle. Otherwise, it will create a completely new solution.
+
+ :param recycle (bool): if True, the new solution will have the same shift in the line of sight as the original
+ solution. Since the projections of each slice are serialized and stored for the entire width of the
+ box (even if the width used is left than the total box), the projection data can be deserialized
+ instead of being remade from scratch. This can greatly speed up the creation of a large number of light
+ cone projections. Default: True.
+ :param filename (str): if given, a text file detailing the solution will be written out. Default: None.
"""
# Get rid of old halo mask, if one was there.
- self.haloMask = []
+ self.halo_mask = []
# Clean pf objects out of light cone solution.
- for slice in self.lightConeSolution:
+ for slice in self.light_cone_solution:
if slice.has_key('object'):
del slice['object']
if recycle:
- if self.verbose: mylog.info("Recycling solution made with %s with new seed %s." %
- (self.originalRandomSeed,
- newSeed))
+ mylog.debug("Recycling solution made with %s with new seed %s." %
+ (self.originalRandomSeed, newSeed))
self.recycleRandomSeed = int(newSeed)
else:
- if self.verbose: mylog.info("Creating new solution with random seed %s." % newSeed)
+ mylog.debug("Creating new solution with random seed %s." % newSeed)
self.originalRandomSeed = int(newSeed)
self.recycleRandomSeed = 0
@@ -398,86 +389,90 @@
boxFractionUsed = 0.0
# Seed random number generator with new seed.
- rand.seed(int(newSeed))
+ na.random.seed(int(newSeed))
- for q,output in enumerate(self.lightConeSolution):
+ for q, output in enumerate(self.light_cone_solution):
# It is necessary to make the same number of calls to the random number generator
# so the original solution willbe produced if the same seed is given.
# Get projection axis and center.
# If using box coherence, only get random axis and center if enough of the box has been used,
# or if boxFractionUsed will be greater than 1 after this slice.
- if (q == 0) or (self.lightConeParameters['MinimumCoherentBoxFraction'] == 0) or \
- (boxFractionUsed > self.lightConeParameters['MinimumCoherentBoxFraction']) or \
- (boxFractionUsed + self.lightConeSolution[q]['DepthBoxFraction'] > 1.0):
+ if (q == 0) or (self.minimum_coherent_box_fraction == 0) or \
+ (boxFractionUsed > self.minimum_coherent_box_fraction) or \
+ (boxFractionUsed + self.light_cone_solution[q]['DepthBoxFraction'] > 1.0):
# Get random projection axis and center.
# If recycling, axis will get thrown away since it is used in creating a unique projection object.
- newAxis = rand.randint(0,2)
+ newAxis = na.random.randint(0, 3)
- newCenter = [rand.random(),rand.random(),rand.random()]
+ newCenter = [na.random.random(), na.random.random(), na.random.random()]
boxFractionUsed = 0.0
else:
# Same axis and center as previous slice, but with depth center shifted.
- newAxis = self.lightConeSolution[q-1]['ProjectionAxis']
- newCenter = copy.deepcopy(self.lightConeSolution[q-1]['ProjectionCenter'])
+ newAxis = self.light_cone_solution[q-1]['ProjectionAxis']
+ newCenter = copy.deepcopy(self.light_cone_solution[q-1]['ProjectionCenter'])
newCenter[newAxis] += \
- 0.5 * (self.lightConeSolution[q]['DepthBoxFraction'] + self.lightConeSolution[q-1]['DepthBoxFraction'])
+ 0.5 * (self.light_cone_solution[q]['DepthBoxFraction'] + self.light_cone_solution[q-1]['DepthBoxFraction'])
if newCenter[newAxis] >= 1.0:
newCenter[newAxis] -= 1.0
if recycle:
- output['ProjectionAxis'] = self.masterSolution[q]['ProjectionAxis']
+ output['ProjectionAxis'] = self.master_solution[q]['ProjectionAxis']
else:
output['ProjectionAxis'] = newAxis
- boxFractionUsed += self.lightConeSolution[q]['DepthBoxFraction']
+ boxFractionUsed += self.light_cone_solution[q]['DepthBoxFraction']
# Make list of rectangle corners to calculate common volume.
- newCube = na.zeros(shape=(len(newCenter),2))
- oldCube = na.zeros(shape=(len(newCenter),2))
+ newCube = na.zeros(shape=(len(newCenter), 2))
+ oldCube = na.zeros(shape=(len(newCenter), 2))
for w in range(len(newCenter)):
- if (w == self.masterSolution[q]['ProjectionAxis']):
- oldCube[w] = [self.masterSolution[q]['ProjectionCenter'][w] - 0.5 * self.masterSolution[q]['DepthBoxFraction'],
- self.masterSolution[q]['ProjectionCenter'][w] + 0.5 * self.masterSolution[q]['DepthBoxFraction']]
+ if (w == self.master_solution[q]['ProjectionAxis']):
+ oldCube[w] = [self.master_solution[q]['ProjectionCenter'][w] - 0.5 * self.master_solution[q]['DepthBoxFraction'],
+ self.master_solution[q]['ProjectionCenter'][w] + 0.5 * self.master_solution[q]['DepthBoxFraction']]
else:
- oldCube[w] = [self.masterSolution[q]['ProjectionCenter'][w] - 0.5 * self.masterSolution[q]['WidthBoxFraction'],
- self.masterSolution[q]['ProjectionCenter'][w] + 0.5 * self.masterSolution[q]['WidthBoxFraction']]
+ oldCube[w] = [self.master_solution[q]['ProjectionCenter'][w] - 0.5 * self.master_solution[q]['WidthBoxFraction'],
+ self.master_solution[q]['ProjectionCenter'][w] + 0.5 * self.master_solution[q]['WidthBoxFraction']]
if (w == output['ProjectionAxis']):
if recycle:
newCube[w] = oldCube[w]
else:
- newCube[w] = [newCenter[w] - 0.5 * self.masterSolution[q]['DepthBoxFraction'],
- newCenter[w] + 0.5 * self.masterSolution[q]['DepthBoxFraction']]
+ newCube[w] = [newCenter[w] - 0.5 * self.master_solution[q]['DepthBoxFraction'],
+ newCenter[w] + 0.5 * self.master_solution[q]['DepthBoxFraction']]
else:
- newCube[w] = [newCenter[w] - 0.5 * self.masterSolution[q]['WidthBoxFraction'],
- newCenter[w] + 0.5 * self.masterSolution[q]['WidthBoxFraction']]
+ newCube[w] = [newCenter[w] - 0.5 * self.master_solution[q]['WidthBoxFraction'],
+ newCenter[w] + 0.5 * self.master_solution[q]['WidthBoxFraction']]
- commonVolume += commonNVolume(oldCube,newCube,periodic=na.array([[0,1],[0,1],[0,1]]))
+ commonVolume += commonNVolume(oldCube, newCube, periodic=na.array([[0, 1], [0, 1], [0, 1]]))
totalVolume += output['DepthBoxFraction'] * output['WidthBoxFraction']**2
# Replace centers for every axis except the line of sight axis.
for w in range(len(newCenter)):
- if not(recycle and (w == self.lightConeSolution[q]['ProjectionAxis'])):
- self.lightConeSolution[q]['ProjectionCenter'][w] = newCenter[w]
+ if not(recycle and (w == self.light_cone_solution[q]['ProjectionAxis'])):
+ self.light_cone_solution[q]['ProjectionCenter'][w] = newCenter[w]
if recycle:
- if self.verbose: mylog.info("Fractional common volume between master and recycled solution is %.2e" % (commonVolume/totalVolume))
+ mylog.debug("Fractional common volume between master and recycled solution is %.2e" % (commonVolume/totalVolume))
else:
- if self.verbose: mylog.info("Fraction of total volume in common with old solution is %.2e." % (commonVolume/totalVolume))
- self.masterSolution = [copy.deepcopy(q) for q in self.lightConeSolution]
+ mylog.debug("Fraction of total volume in common with old solution is %.2e." % (commonVolume/totalVolume))
+ self.master_solution = [copy.deepcopy(q) for q in self.light_cone_solution]
+
+ # Write solution to a file.
+ if filename is not None:
+ self._save_light_cone_solution(filename=filename)
- def RestoreMasterSolution(self):
+ def restore_master_solution(self):
"Reset the active light cone solution to the master solution."
- self.lightConeSolution = [copy.deepcopy(q) for q in self.masterSolution]
+ self.light_cone_solution = [copy.deepcopy(q) for q in self.master_solution]
@rootonly
- def SaveLightConeSolution(self,file="light_cone.dat"):
+ def _save_light_cone_solution(self, filename="light_cone.dat"):
"Write out a text file with information on light cone solution."
- if self.verbose: mylog.info("Saving light cone solution to %s." % file)
+ mylog.info("Saving light cone solution to %s." % filename)
- f = open(file,'w')
+ f = open(filename, 'w')
if self.recycleSolution:
f.write("Recycled Solution\n")
f.write("OriginalRandomSeed = %s\n" % self.originalRandomSeed)
@@ -486,243 +481,33 @@
f.write("Original Solution\n")
f.write("OriginalRandomSeed = %s\n" % self.originalRandomSeed)
f.write("EnzoParameterFile = %s\n" % self.EnzoParameterFile)
- f.write("LightConeParameterFile = %s\n" % self.LightConeParameterFile)
f.write("\n")
- for parameter in self.lightConeParameters:
- f.write("%s = %s\n" % (parameter,str(self.lightConeParameters[parameter])))
- f.write("\n")
- for q,output in enumerate(self.lightConeSolution):
+ for q, output in enumerate(self.light_cone_solution):
f.write("Proj %04d, %s, z = %f, depth/box = %f, width/box = %f, axis = %d, center = %f, %f, %f\n" %
- (q,output['filename'],output['redshift'],output['DepthBoxFraction'],output['WidthBoxFraction'],
- output['ProjectionAxis'],output['ProjectionCenter'][0],output['ProjectionCenter'][1],output['ProjectionCenter'][2]))
+ (q, output['filename'], output['redshift'], output['DepthBoxFraction'], output['WidthBoxFraction'],
+ output['ProjectionAxis'], output['ProjectionCenter'][0], output['ProjectionCenter'][1], output['ProjectionCenter'][2]))
f.close()
- def _CalculateDeltaZMax(self):
- "Calculate delta z that corresponds to full box length going from z to (z - delta z)."
- co = lagos.Cosmology(HubbleConstantNow = (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
- OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
- OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'])
-
- d_Tolerance = 1e-4
- max_Iterations = 100
-
- targetDistance = self.enzoParameters['CosmologyComovingBoxSize']
-
- for output in self.allOutputs:
- z = output['redshift']
-
- # Calculate delta z that corresponds to the length of the box at a given redshift.
- # Use Newton's method to calculate solution.
- z1 = z
- z2 = z1 - 0.1 # just an initial guess
- distance1 = 0.0
- iteration = 1
-
- # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
- distance2 = co.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
-
- while ((na.fabs(distance2-targetDistance)/distance2) > d_Tolerance):
- m = (distance2 - distance1) / (z2 - z1)
- z1 = z2
- distance1 = distance2
- z2 = ((targetDistance - distance2) / m) + z2
- distance2 = co.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
- iteration += 1
- if (iteration > max_Iterations):
- if self.verbose: mylog.error("CalculateDeltaZMax: Warning - max iterations exceeded for z = %f (delta z = %f)." % (z,na.fabs(z2-z)))
- break
- output['deltazMax'] = na.fabs(z2-z)
-
- del co
-
- def _CalculateDeltaZMin(self):
- "Calculate delta z that corresponds to a single top grid pixel going from z to (z - delta z)."
- co = lagos.Cosmology(HubbleConstantNow = (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
- OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
- OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'])
-
- d_Tolerance = 1e-4
- max_Iterations = 100
-
- targetDistance = self.enzoParameters['CosmologyComovingBoxSize'] / self.enzoParameters['TopGridDimensions'][0]
-
- for output in self.allOutputs:
- z = output['redshift']
-
- # Calculate delta z that corresponds to the length of a top grid pixel at a given redshift.
- # Use Newton's method to calculate solution.
- z1 = z
- z2 = z1 - 0.01 # just an initial guess
- distance1 = 0.0
- iteration = 1
-
- # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
- distance2 = co.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
-
- while ((na.fabs(distance2 - targetDistance) / distance2) > d_Tolerance):
- m = (distance2 - distance1) / (z2 - z1)
- z1 = z2
- distance1 = distance2
- z2 = ((targetDistance - distance2) / m) + z2
- distance2 = co.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
- iteration += 1
- if (iteration > max_Iterations):
- if self.verbose: mylog.error("CalculateDeltaZMax: Warning - max iterations exceeded for z = %f (delta z = %f)." %
- (z,na.fabs(z2-z)))
- break
- # Use this calculation or the absolute minimum specified by the user.
- output['deltazMin'] = max(na.fabs(z2-z),self.lightConeParameters['MinimumSliceDeltaz'])
-
- del co
-
- def _CalculateTimeDumpRedshifts(self):
- "Calculates redshift values for the dt data dumps."
- ec = lagos.EnzoCosmology(HubbleConstantNow = \
- (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
- OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
- OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'],
- InitialRedshift = self.enzoParameters['CosmologyInitialRedshift'])
-
- z_Tolerance = 1e-3
- z = self.enzoParameters['CosmologyInitialRedshift']
- index = 0
- while ((z > self.enzoParameters['CosmologyFinalRedshift']) and
- (na.fabs(z-self.enzoParameters['CosmologyFinalRedshift']) > z_Tolerance)):
- t = ec.InitialTime + (self.enzoParameters['dtDataDump'] * ec.TimeUnits * index)
- z = ec.ComputeRedshiftFromTime(t)
- filename = "%s/%s%04d/%s%04d" % (self.enzoParameters['GlobalDir'],
- self.enzoParameters['DataDumpDir'],
- index,
- self.enzoParameters['DataDumpDir'],
- index)
- self.timeOutputs.append({'index':index,'redshift':z,'filename':filename})
- index += 1
-
- def _CalculateSplitTimeDumpRedshifts(self,break_index,new_dt):
- "Calculates redshift values for the dt data dumps."
- ec = lagos.EnzoCosmology(HubbleConstantNow = \
- (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
- OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
- OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'],
- InitialRedshift = self.enzoParameters['CosmologyInitialRedshift'])
-
- z_Tolerance = 1e-3
- z = self.enzoParameters['CosmologyInitialRedshift']
- index = 0
- while ((z > self.enzoParameters['CosmologyFinalRedshift']) and
- (na.fabs(z-self.enzoParameters['CosmologyFinalRedshift']) > z_Tolerance)):
- if index <= break_index:
- t = ec.InitialTime + (self.enzoParameters['dtDataDump'] * ec.TimeUnits * index)
- z = ec.ComputeRedshiftFromTime(t)
- else:
- t = ec.InitialTime + (self.enzoParameters['dtDataDump'] * ec.TimeUnits * (break_index))
- t += (new_dt * ec.TimeUnits * (index-break_index))
- z = ec.ComputeRedshiftFromTime(t)
- filename = "%s/%s%04d/%s%04d" % (self.enzoParameters['GlobalDir'],
- self.enzoParameters['DataDumpDir'],
- index,
- self.enzoParameters['DataDumpDir'],
- index)
- self.timeOutputs.append({'index':index,'redshift':z,'filename':filename})
- index += 1
-
-
- def _CombineDataOutputs(self):
- "Combines redshift and time data into one sorted list."
- self.allOutputs = self.redshiftOutputs + self.timeOutputs
- self.allOutputs.sort(reverse=True,key=lambda obj:obj['redshift'])
- for q in range(len(self.allOutputs)):
- del self.allOutputs[q]['index']
- if (q == 0):
- self.allOutputs[q]['previous'] = None
- self.allOutputs[q]['next'] = self.allOutputs[q+1]
- elif (q == len(self.allOutputs) - 1):
- self.allOutputs[q]['previous'] = self.allOutputs[q-1]
- self.allOutputs[q]['next'] = None
- else:
- self.allOutputs[q]['previous'] = self.allOutputs[q-1]
- self.allOutputs[q]['next'] = self.allOutputs[q+1]
- del self.redshiftOutputs
- del self.timeOutputs
-
- def _ReadEnzoParameterFile(self):
- "Reads an Enzo parameter file looking for cosmology and output parameters."
- lines = open(self.EnzoParameterFile).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:
- if self.verbose: mylog.error("Skipping line: %s" % line)
- continue
- if EnzoParameterDict.has_key(param):
- t = map(EnzoParameterDict[param], vals.split())
- if len(t) == 1:
- self.enzoParameters[param] = t[0]
- else:
- self.enzoParameters[param] = t
- elif param.startswith("CosmologyOutputRedshift["):
- index = param[param.find("[")+1:param.find("]")]
- self.redshiftOutputs.append({'index':int(index),
- 'redshift':float(vals)})
-
- # Add filenames to redshift outputs.
- for output in self.redshiftOutputs:
- output["filename"] = "%s/%s%04d/%s%04d" % (self.enzoParameters['GlobalDir'],
- self.enzoParameters['RedshiftDumpDir'],
- output['index'],
- self.enzoParameters['RedshiftDumpDir'],
- output['index'])
-
- def _ReadLightConeParameterFile(self):
- "Reads a light cone parameter file."
- lines = open(self.LightConeParameterFile).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:
- if self.verbose: mylog.error("Skipping line: %s" % line)
- continue
- if LightConeParameterDict.has_key(param):
- t = map(LightConeParameterDict[param], vals.split())
- if len(t) == 1:
- self.lightConeParameters[param] = t[0]
- else:
- self.lightConeParameters[param] = t
-
- def _SaveLightConeStack(self,field=None,weight_field=None,filename=None,over_write=True):
+ def _save_light_cone_stack(self, field=None, weight_field=None, filename=None, over_write=True):
"Save the light cone projection stack as a 3d array in and hdf5 file."
# Make list of redshifts to include as a dataset attribute.
- redshiftList = na.array([slice['redshift'] for slice in self.lightConeSolution])
+ redshiftList = na.array([slice['redshift'] for slice in self.light_cone_solution])
- field_node = "%s_%s" % (field,weight_field)
+ field_node = "%s_%s" % (field, weight_field)
weight_field_node = "weight_field_%s" % weight_field
import h5py
if (filename is None):
- filename = "%s/%s_data" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
+ filename = "%s/%s_data" % (self.output_dir, self.output_prefix)
if not(filename.endswith('.h5')):
filename += ".h5"
- if (len(self.projectionStack) == 0):
- if self.verbose: mylog.error("SaveLightConeStack: no projection data loaded.")
+ if (len(self.projection_stack) == 0):
+ mylog.debug("save_light_cone_stack: no projection data loaded.")
return
- if self.verbose: mylog.info("Writing light cone data to %s." % filename)
+ mylog.info("Writing light cone data to %s." % filename)
output = h5py.File(filename, "a")
@@ -734,82 +519,38 @@
write_data = True
del output[field_node]
else:
- mylog.info("Dataset, %s, already exists in %s, not saving." % (field_node,filename))
+ mylog.info("Dataset, %s, already exists in %s, not saving." % (field_node, filename))
write_data = False
else:
write_data = True
if write_data:
mylog.info("Saving %s to %s." % (field_node, filename))
- self.projectionStack = na.array(self.projectionStack)
- field_dataset = output.create_dataset(field_node,data=self.projectionStack)
+ self.projection_stack = na.array(self.projection_stack)
+ field_dataset = output.create_dataset(field_node, data=self.projection_stack)
field_dataset.attrs['redshifts'] = redshiftList
- field_dataset.attrs['ObserverRedshift'] = na.float(self.lightConeParameters['ObserverRedshift'])
- field_dataset.attrs['FieldOfViewInArcMinutes'] = na.float(self.lightConeParameters['FieldOfViewInArcMinutes'])
- field_dataset.attrs['ImageResolutionInArcSeconds'] = na.float(self.lightConeParameters['ImageResolutionInArcSeconds'])
+ field_dataset.attrs['observer_redshift'] = na.float(self.observer_redshift)
+ field_dataset.attrs['field_of_view_in_arcminutes'] = na.float(self.field_of_view_in_arcminutes)
+ field_dataset.attrs['image_resolution_in_arcseconds'] = na.float(self.image_resolution_in_arcseconds)
- if (len(self.projectionWeightFieldStack) > 0):
+ if (len(self.projection_weight_field_stack) > 0):
if node_exists:
if over_write:
mylog.info("Dataset, %s, already exists, overwriting." % weight_field_node)
del output[field_node]
else:
- mylog.info("Dataset, %s, already exists in %s, not saving." % (weight_field_node,filename))
+ mylog.info("Dataset, %s, already exists in %s, not saving." % (weight_field_node, filename))
write_data = False
else:
write_data = True
if write_data:
mylog.info("Saving %s to %s." % (weight_field_node, filename))
- self.projectionWeightFieldStack = na.array(self.projectionWeightFieldStack)
- weight_field_dataset = output.create_dataset(weight_field_node,data=self.projectionWeightFieldStack)
+ self.projection_weight_field_stack = na.array(self.projection_weight_field_stack)
+ weight_field_dataset = output.create_dataset(weight_field_node, data=self.projection_weight_field_stack)
weight_field_dataset.attrs['redshifts'] = redshiftList
- weight_field_dataset.attrs['ObserverRedshift'] = na.float(self.lightConeParameters['ObserverRedshift'])
- weight_field_dataset.attrs['FieldOfViewInArcMinutes'] = na.float(self.lightConeParameters['FieldOfViewInArcMinutes'])
- weight_field_dataset.attrs['ImageResolutionInArcSeconds'] = na.float(self.lightConeParameters['ImageResolutionInArcSeconds'])
+ weight_field_dataset.attrs['observer_redshift'] = na.float(self.observer_redshift)
+ weight_field_dataset.attrs['field_of_view_in_arcminutes'] = na.float(self.field_of_view_in_arcminutes)
+ weight_field_dataset.attrs['image_resolution_in_arcseconds'] = na.float(self.image_resolution_in_arcseconds)
output.close()
-
- def _SetParameterDefaults(self):
- "Set some default parameters to avoid problems if they are not in the parameter file."
- self.enzoParameters['GlobalDir'] = "./"
- self.enzoParameters['RedshiftDumpName'] = "RD"
- self.enzoParameters['RedshiftDumpDir'] = "RD"
- self.enzoParameters['DataDumpName'] = "DD"
- self.enzoParameters['DataDumpDir'] = "DD"
- self.lightConeParameters['UseMinimumNumberOfProjections'] = 1
- self.lightConeParameters['MinimumCoherentBoxFraction'] = 0.0
- self.lightConeParameters['MinimumSliceDeltaz'] = 0.0
- self.lightConeParameters['ObserverRedshift'] = 0.0
- self.lightConeParameters['OutputDir'] = "./"
- self.lightConeParameters['OutputPrefix'] = "LightCone"
- self.lightConeParameters['dtSplit'] = None
- self.lightConeParameters['dt_new'] = 0.0
-
-EnzoParameterDict = {"CosmologyCurrentRedshift": float,
- "CosmologyComovingBoxSize": float,
- "CosmologyOmegaMatterNow": float,
- "CosmologyOmegaLambdaNow": float,
- "CosmologyHubbleConstantNow": float,
- "CosmologyInitialRedshift": float,
- "CosmologyFinalRedshift": float,
- "TopGridDimensions": float,
- "dtDataDump": float,
- "RedshiftDumpName": str,
- "RedshiftDumpDir": str,
- "DataDumpName": str,
- "DataDumpDir": str,
- "GlobalDir" : str}
-
-LightConeParameterDict = {"InitialRedshift": float,
- "FinalRedshift": float,
- "ObserverRedshift": float,
- "FieldOfViewInArcMinutes": float,
- "ImageResolutionInArcSeconds": float,
- "UseMinimumNumberOfProjections": int,
- "MinimumCoherentBoxFraction": float,
- "MinimumSliceDeltaz": float,
- "OutputDir": str,
- "OutputPrefix": str,
- "dtSplit": float,
- "dt_new": float}
Modified: branches/yt-1.5/yt/extensions/lightcone/LightConeProjection.py
==============================================================================
--- branches/yt-1.5/yt/extensions/lightcone/LightConeProjection.py (original)
+++ branches/yt-1.5/yt/extensions/lightcone/LightConeProjection.py Tue Nov 3 15:45:17 2009
@@ -31,13 +31,9 @@
import numpy as na
import copy
-#### Note: There is an assumption that the box width is 1 in every direction. This doesn't have to be, but
-#### I don't have time to fix it now. All that is required is a multiplication by the box width in
-#### the direction in question wherever DepthBoxFraction or WidthBoxFraction appears.
-#### To be fixed before turning 30. - Britton
-
@lagos.parallel_blocking_call
-def LightConeProjection(lightConeSlice,field,pixels,weight_field=None,save_image=False,name="",node=None,field_cuts=None,**kwargs):
+def LightConeProjection(lightConeSlice, field, pixels, weight_field=None, save_image=False, name="", node=None, field_cuts=None,
+ add_redshift_label=False, **kwargs):
"Create a single projection to be added into the light cone stack."
# Use some projection parameters to seed random number generator to make unique node name.
@@ -50,11 +46,11 @@
# DepthBoxFraction
# Name node with user specified keyword if given with 'node' keyword.
- node_name = "LightCone_%s_%d_%f_%f" % (node,lightConeSlice['ProjectionAxis'],
+ node_name = "LightCone_%s_%d_%f_%f" % (node, lightConeSlice['ProjectionAxis'],
lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']],
lightConeSlice['DepthBoxFraction'])
- mylog.info("Making projection at z = %f from %s." % (lightConeSlice['redshift'],lightConeSlice['filename']))
+ mylog.info("Making projection at z = %f from %s." % (lightConeSlice['redshift'], lightConeSlice['filename']))
region_center = [0.5 * (lightConeSlice['object'].parameters['DomainRightEdge'][q] +
lightConeSlice['object'].parameters['DomainLeftEdge'][q]) \
@@ -62,7 +58,7 @@
# Make the plot collection and put it in the slice so we can delete it cleanly in the same scope
# as where the frb will be deleted.
- lightConeSlice['pc'] = raven.PlotCollection(lightConeSlice['object'],center=region_center)
+ lightConeSlice['pc'] = raven.PlotCollection(lightConeSlice['object'], center=region_center)
# 1. The Depth Problem
# Use coordinate field cut in line of sight to cut projection to proper depth.
@@ -72,26 +68,26 @@
these_field_cuts = copy.deepcopy(field_cuts)
if (lightConeSlice['DepthBoxFraction'] < 1):
- axis = ('x','y','z')[lightConeSlice['ProjectionAxis']]
+ axis = ('x', 'y', 'z')[lightConeSlice['ProjectionAxis']]
depthLeft = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] - 0.5 * lightConeSlice['DepthBoxFraction']
depthRight = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] + 0.5 * lightConeSlice['DepthBoxFraction']
if (depthLeft < 0):
cut_mask = "((grid[\"%s\"] + 0.5*grid[\"d%s\"] >= 0) & (grid[\"%s\"] - 0.5*grid[\"d%s\"] <= %f)) | ((grid[\"%s\"] + 0.5*grid[\"d%s\"] >= %f) & (grid[\"%s\"] - 0.5*grid[\"d%s\"] <= 1))" % \
- (axis,axis,axis,axis,depthRight,axis,axis,(depthLeft+1),axis,axis)
+ (axis, axis, axis, axis, depthRight, axis, axis, (depthLeft+1), axis, axis)
elif (depthRight > 1):
cut_mask = "((grid[\"%s\"] + 0.5*grid[\"d%s\"] >= 0) & (grid[\"%s\"] - 0.5*grid[\"d%s\"] <= %f)) | ((grid[\"%s\"] + 0.5*grid[\"d%s\"] >= %f) & (grid[\"%s\"] - 0.5*grid[\"d%s\"] <= 1))" % \
- (axis,axis,axis,axis,(depthRight-1),axis,axis,depthLeft,axis,axis)
+ (axis, axis, axis, axis, (depthRight-1), axis, axis, depthLeft, axis, axis)
else:
- cut_mask = "(grid[\"%s\"] + 0.5*grid[\"d%s\"] >= %f) & (grid[\"%s\"] - 0.5*grid[\"%s\"] <= %f)" % (axis,axis,depthLeft,axis,axis,depthRight)
+ cut_mask = "(grid[\"%s\"] + 0.5*grid[\"d%s\"] >= %f) & (grid[\"%s\"] - 0.5*grid[\"%s\"] <= %f)" % (axis, axis, depthLeft, axis, axis, depthRight)
these_field_cuts.append(cut_mask)
# Make projection.
- lightConeSlice['pc'].add_projection(field,lightConeSlice['ProjectionAxis'],weight_field=weight_field,field_cuts=these_field_cuts,
- node_name=node_name,**kwargs)
+ lightConeSlice['pc'].add_projection(field, lightConeSlice['ProjectionAxis'], weight_field=weight_field, field_cuts=these_field_cuts,
+ node_name=node_name, **kwargs)
# If parallel: all the processes have the whole projection object, but we only need to do the tiling, shifting, and cutting once.
- if ytcfg.getint("yt","__parallel_rank") == 0:
+ if ytcfg.getint("yt", "__parallel_rank") == 0:
# 2. The Tile Problem
# Tile projection to specified width.
@@ -108,12 +104,12 @@
for x in range(int(na.ceil(lightConeSlice['WidthBoxFraction']))):
for y in range(int(na.ceil(lightConeSlice['WidthBoxFraction']))):
if ((x + y) > 0):
- lightConeSlice['pc'].plots[0].data['px'] = na.concatenate([lightConeSlice['pc'].plots[0].data['px'],original_px+x])
- lightConeSlice['pc'].plots[0].data['py'] = na.concatenate([lightConeSlice['pc'].plots[0].data['py'],original_py+y])
- lightConeSlice['pc'].plots[0].data['pdx'] = na.concatenate([lightConeSlice['pc'].plots[0].data['pdx'],original_pdx])
- lightConeSlice['pc'].plots[0].data['pdy'] = na.concatenate([lightConeSlice['pc'].plots[0].data['pdy'],original_pdy])
- lightConeSlice['pc'].plots[0].data[field] = na.concatenate([lightConeSlice['pc'].plots[0].data[field],original_field])
- lightConeSlice['pc'].plots[0].data['weight_field'] = na.concatenate([lightConeSlice['pc'].plots[0].data['weight_field'],original_weight_field])
+ lightConeSlice['pc'].plots[0].data['px'] = na.concatenate([lightConeSlice['pc'].plots[0].data['px'], original_px+x])
+ lightConeSlice['pc'].plots[0].data['py'] = na.concatenate([lightConeSlice['pc'].plots[0].data['py'], original_py+y])
+ lightConeSlice['pc'].plots[0].data['pdx'] = na.concatenate([lightConeSlice['pc'].plots[0].data['pdx'], original_pdx])
+ lightConeSlice['pc'].plots[0].data['pdy'] = na.concatenate([lightConeSlice['pc'].plots[0].data['pdy'], original_pdy])
+ lightConeSlice['pc'].plots[0].data[field] = na.concatenate([lightConeSlice['pc'].plots[0].data[field], original_field])
+ lightConeSlice['pc'].plots[0].data['weight_field'] = na.concatenate([lightConeSlice['pc'].plots[0].data['weight_field'], original_weight_field])
# Delete originals.
del original_px
@@ -186,21 +182,21 @@
del add_y_left
# Add the hanging cells back to the projection data.
- lightConeSlice['pc'].plots[0].data['px'] = na.concatenate([lightConeSlice['pc'].plots[0]['px'],add_x_px,add_y_px,add2_x_px,add2_y_px])
- lightConeSlice['pc'].plots[0].data['py'] = na.concatenate([lightConeSlice['pc'].plots[0]['py'],add_x_py,add_y_py,add2_x_py,add2_y_py])
- lightConeSlice['pc'].plots[0].data['pdx'] = na.concatenate([lightConeSlice['pc'].plots[0]['pdx'],add_x_pdx,add_y_pdx,add2_x_pdx,add2_y_pdx])
- lightConeSlice['pc'].plots[0].data['pdy'] = na.concatenate([lightConeSlice['pc'].plots[0]['pdy'],add_x_pdy,add_y_pdy,add2_x_pdy,add2_y_pdy])
- lightConeSlice['pc'].plots[0].data[field] = na.concatenate([lightConeSlice['pc'].plots[0][field],add_x_field,add_y_field,add2_x_field,add2_y_field])
- lightConeSlice['pc'].plots[0].data['weight_field'] = na.concatenate([lightConeSlice['pc'].plots[0]['weight_field'],add_x_weight_field,add_y_weight_field,
- add2_x_weight_field,add2_y_weight_field])
+ lightConeSlice['pc'].plots[0].data['px'] = na.concatenate([lightConeSlice['pc'].plots[0]['px'], add_x_px, add_y_px, add2_x_px, add2_y_px])
+ lightConeSlice['pc'].plots[0].data['py'] = na.concatenate([lightConeSlice['pc'].plots[0]['py'], add_x_py, add_y_py, add2_x_py, add2_y_py])
+ lightConeSlice['pc'].plots[0].data['pdx'] = na.concatenate([lightConeSlice['pc'].plots[0]['pdx'], add_x_pdx, add_y_pdx, add2_x_pdx, add2_y_pdx])
+ lightConeSlice['pc'].plots[0].data['pdy'] = na.concatenate([lightConeSlice['pc'].plots[0]['pdy'], add_x_pdy, add_y_pdy, add2_x_pdy, add2_y_pdy])
+ lightConeSlice['pc'].plots[0].data[field] = na.concatenate([lightConeSlice['pc'].plots[0][field], add_x_field, add_y_field, add2_x_field, add2_y_field])
+ lightConeSlice['pc'].plots[0].data['weight_field'] = na.concatenate([lightConeSlice['pc'].plots[0]['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
# Tiles were made rounding up the width to the nearest integer.
# Cut off the edges to get the specified width.
@@ -226,14 +222,16 @@
# Save an image if requested.
if save_image:
- lightConeSlice['pc'].set_xlim(0,lightConeSlice['WidthBoxFraction'])
- lightConeSlice['pc'].set_ylim(0,lightConeSlice['WidthBoxFraction'])
+ lightConeSlice['pc'].set_xlim(0, lightConeSlice['WidthBoxFraction'])
+ lightConeSlice['pc'].set_ylim(0, lightConeSlice['WidthBoxFraction'])
+ if add_redshift_label:
+ lightConeSlice['pc'].plots[-1].modify['text']((0.5, 0.03), "z = %.3f" % lightConeSlice['redshift'], dict(color='black',size=50))
lightConeSlice['pc'].save(name)
# Create fixed resolution buffer to return back to the light cone object.
# These buffers will be stacked together to make the light cone.
- frb = raven.FixedResolutionBuffer(lightConeSlice['pc'].plots[0].data,(0,lightConeSlice['WidthBoxFraction'],0,lightConeSlice['WidthBoxFraction']),
- (pixels,pixels),antialias=False)
+ frb = raven.FixedResolutionBuffer(lightConeSlice['pc'].plots[0].data, (0, lightConeSlice['WidthBoxFraction'], 0, lightConeSlice['WidthBoxFraction']),
+ (pixels, pixels), antialias=False)
return frb
else:
Modified: branches/yt-1.5/yt/extensions/lightcone/UniqueSolution.py
==============================================================================
--- branches/yt-1.5/yt/extensions/lightcone/UniqueSolution.py (original)
+++ branches/yt-1.5/yt/extensions/lightcone/UniqueSolution.py Tue Nov 3 15:45:17 2009
@@ -28,36 +28,35 @@
from yt.logger import lagosLogger as mylog
import numpy as na
import random as rand
+import copy
import sys
-def ProjectUniqueLightCones(EnzoParameterFile,LightConeParameterFile,seedFile,field,**kwargs):
+def project_unique_light_cones(lightcone, seed_file, field, **kwargs):
"Make light cone projections using list of random seeds in a file."
- seedList = _ReadSeedFile(seedFile)
+ seedList = _read_seed_file(seed_file)
- lc = LightCone(EnzoParameterFile,LightConeParameterFile)
- prefix = lc.lightConeParameters['OutputPrefix']
- lc.CalculateLightConeSolution(seed=0)
+ prefix = lightcone.output_prefix
+ lightcone.calculate_light_cone_solution(seed=0)
lastSeed = None
for seed in seedList:
if (seed['master'] != lastSeed):
- lc.RerandomizeLightConeSolution(seed['master'],recycle=False)
+ lightcone.rerandomize_light_cone_solution(seed['master'], recycle=False)
lastSeed = seed['master']
if (seed['recycle'] is not None):
- lc.RerandomizeLightConeSolution(seed['recycle'],recycle=True)
+ lightcone.rerandomize_light_cone_solution(seed['recycle'], recycle=True)
- lc.lightConeParameters['OutputPrefix'] = "%s_%s_%s" % (prefix,seed['master'],seed['recycle'])
- lc.ProjectLightCone(field,**kwargs)
+ lightcone.output_prefix = "%s_%s_%s" % (prefix, seed['master'], seed['recycle'])
+ lightcone.project_light_cone(field, **kwargs)
-def FindUniqueSolutions(EnzoParameterFile,LightConeParameterFile,solutions=100,seed=None,
- max_overlap=0.25,failures=10,recycle=True,filename='unique.dat'):
+def find_unique_solutions(lightcone1, solutions=100, seed=None, max_overlap=0.25, failures=10,
+ recycle=True, filename='unique.dat'):
"Find a set of random seeds that will give light cones will minimal volume overlap."
- solution1 = LightCone(EnzoParameterFile,LightConeParameterFile,verbose=False)
- solution2 = LightCone(EnzoParameterFile,LightConeParameterFile,verbose=False)
- solution1.CalculateLightConeSolution(seed=0)
- solution2.CalculateLightConeSolution(seed=0)
+ lightcone2 = copy.deepcopy(lightcone1)
+ lightcone1.calculate_light_cone_solution(seed=0)
+ lightcone2.calculate_light_cone_solution(seed=0)
uniqueSeeds = []
if recycle:
@@ -83,40 +82,41 @@
if (recycle and master is not None):
newSeed = master
if state is not None: rand.setstate(state)
- newRecycleSeed = rand.randint(1,1e9)
+ newRecycleSeed = rand.randint(1, 1e9)
state = rand.getstate()
else:
if state is not None: rand.setstate(state)
- newSeed = rand.randint(1,1e9)
+ newSeed = rand.randint(1, 1e9)
state = rand.getstate()
if recycle:
master = newSeed
recycleFails = 0
newRecycleSeed = None
- sys.stderr.write(("Unique solutions: %d, consecutive failures: %"+failDigits+"d, %"+failDigits+"d.\r") % (len(uniqueSeeds),fails,recycleFails))
+ sys.stderr.write(("Unique solutions: %d, consecutive failures: %"+failDigits+"d, %"+failDigits+"d.\r") % \
+ (len(uniqueSeeds), fails, recycleFails))
- solution1.RerandomizeLightConeSolution(newSeed,recycle=False)
+ lightcone1.rerandomize_light_cone_solution(newSeed, recycle=False)
if newRecycleSeed is not None:
- solution1.RerandomizeLightConeSolution(newRecycleSeed,recycle=True)
+ lightcone1.rerandomize_light_cone_solution(newRecycleSeed, recycle=True)
# Compare with all other seeds.
testPass = True
for uniqueSeed in uniqueSeeds:
- solution2.RerandomizeLightConeSolution(uniqueSeed['master'],recycle=False)
+ lightcone2.rerandomize_light_cone_solution(uniqueSeed['master'], recycle=False)
if uniqueSeed['recycle'] is not None:
- solution2.RerandomizeLightConeSolution(uniqueSeed['recycle'],recycle=True)
+ lightcone2.rerandomize_light_cone_solution(uniqueSeed['recycle'], recycle=True)
- common = CompareSolutions(solution1.lightConeSolution,solution2.lightConeSolution)
+ common = _compare_solutions(lightcone1.light_cone_solution, lightcone2.light_cone_solution)
if (common > max_overlap):
testPass = False
break
else:
- maxCommon = max(maxCommon,common)
+ maxCommon = max(maxCommon, common)
if testPass:
- uniqueSeeds.append({'master':newSeed,'recycle':newRecycleSeed})
+ uniqueSeeds.append({'master':newSeed, 'recycle':newRecycleSeed})
fails = 0
recycleFails = 0
@@ -127,21 +127,23 @@
fails += 1
if (recycleFails >= failures):
- sys.stderr.write(("Unique solutions: %d, consecutive failures: %"+failDigits+"d, %"+failDigits+"d.\n") % (len(uniqueSeeds),fails,recycleFails))
+ sys.stderr.write(("Unique solutions: %d, consecutive failures: %"+failDigits+"d, %"+failDigits+"d.\n") % \
+ (len(uniqueSeeds), fails, recycleFails))
fails += 1
mylog.info("Max recycled failures reached with master seed %d." % newSeed)
master = None
if (fails >= failures):
- sys.stderr.write(("Unique solutions: %d, consecutive failures: %"+failDigits+"d, %"+failDigits+"d.\n") % (len(uniqueSeeds),fails,recycleFails))
+ sys.stderr.write(("Unique solutions: %d, consecutive failures: %"+failDigits+"d, %"+failDigits+"d.\n") % \
+ (len(uniqueSeeds), fails, recycleFails))
mylog.error("Max consecutive failures reached.")
break
mylog.info("Created %d unique solutions." % len(uniqueSeeds))
mylog.info("Maximum common volume is %.2e." % maxCommon)
- _WriteSeedFile(uniqueSeeds,filename)
+ _write_seed_file(uniqueSeeds, filename)
return uniqueSeeds
-def CompareSolutions(solution1,solution2):
+def _compare_solutions(solution1, solution2):
"Calculate common volume between two light cone solutions."
if (len(solution1) != len(solution2)):
@@ -157,7 +159,7 @@
mylog.error("Light cone solutions do not have equal volumes, will use the smaller one.")
for q in range(len(solution1)):
- cube1 = na.zeros(shape=(len(solution1[q]['ProjectionCenter']),2))
+ cube1 = na.zeros(shape=(len(solution1[q]['ProjectionCenter']), 2))
volume1 = 1.0
for w in range(len(cube1)):
if (w == solution1[q]['ProjectionAxis']):
@@ -168,7 +170,7 @@
cube1[w] = [solution1[q]['ProjectionCenter'][w] - 0.5 * width,
solution1[q]['ProjectionCenter'][w] + 0.5 * width]
- cube2 = na.zeros(shape=(len(solution2[q]['ProjectionCenter']),2))
+ cube2 = na.zeros(shape=(len(solution2[q]['ProjectionCenter']), 2))
volume2 = 1.0
for w in range(len(cube2)):
if (w == solution2[q]['ProjectionAxis']):
@@ -179,12 +181,12 @@
cube2[w] = [solution2[q]['ProjectionCenter'][w] - 0.5 * width,
solution2[q]['ProjectionCenter'][w] + 0.5 * width]
- totalVolume += min(volume1,volume2)
- commonVolume += commonNVolume(cube1,cube2,periodic=na.array([[0,1],[0,1],[0,1]]))
+ totalVolume += min(volume1, volume2)
+ commonVolume += commonNVolume(cube1, cube2, periodic=na.array([[0, 1], [0, 1], [0, 1]]))
return (commonVolume/totalVolume)
-def _ReadSeedFile(filename):
+def _read_seed_file(filename):
"Read list of random seeds from a file."
mylog.info("Reading random seed list from %s." % filename)
@@ -203,15 +205,16 @@
return seedList
-def _WriteSeedFile(seedList,filename):
+ at rootonly
+def _write_seed_file(seedList, filename):
"Write list of random seeds to a file."
mylog.info("Writing random seed list to %s." % filename)
- f = open(filename,'w')
+ f = open(filename, 'w')
for seed in seedList:
if seed['recycle'] is None:
f.write("%s\n" % seed['master'])
else:
- f.write("%s,%s\n" % (seed['master'],seed['recycle']))
+ f.write("%s, %s\n" % (seed['master'], seed['recycle']))
f.close()
More information about the yt-svn
mailing list