[Yt-svn] yt-commit r1017 - trunk/yt/extensions/lightcone
britton at wrangler.dreamhost.com
britton at wrangler.dreamhost.com
Tue Dec 16 16:41:22 PST 2008
Author: britton
Date: Tue Dec 16 16:41:22 2008
New Revision: 1017
URL: http://yt.spacepope.org/changeset/1017
Log:
Changed serialization of projection objects for the light cone. Additionally,
SaveLightConeStack is no longer called directly by the user. Instead,
the save_stack=True keyword is given to ProjectLightCone and it will be done
from there. The datasets in the hdf5 file are now named for the fields they
represent, instead of just being called 'data'. This allows multiple light
cone fields to be written to the same file. Also, the option to save images of
each slice is now handled in a similar manner, with the save_slice_images=True
keyword, instead of the paramter file parameter.
Modified:
trunk/yt/extensions/lightcone/LightCone.py
trunk/yt/extensions/lightcone/LightConeProjection.py
Modified: trunk/yt/extensions/lightcone/LightCone.py
==============================================================================
--- trunk/yt/extensions/lightcone/LightCone.py (original)
+++ trunk/yt/extensions/lightcone/LightCone.py Tue Dec 16 16:41:22 2008
@@ -25,13 +25,13 @@
from yt.extensions.lightcone import *
from yt.logger import lagosLogger as mylog
+from yt.config import ytcfg
from Common_nVolume import *
from HaloMask import *
import copy
import os
import numpy as na
import random as rand
-import tables as h5
class LightCone(object):
def __init__(self,EnzoParameterFile,LightConeParameterFile,verbose=True):
@@ -79,8 +79,12 @@
# Combine all data dumps.
self._CombineDataOutputs()
- # Calculate maximum delta z for each data dump.
- self._CalculateDeltaZMax()
+ if self.lightConeParameters['UseMinimumNumberOfProjections']:
+ # Calculate maximum delta z for each data dump.
+ self._CalculateDeltaZMax()
+ else:
+ # Calculate minimum delta z for each data dump.
+ self._CalculateDeltaZMin()
def CalculateLightConeSolution(self,seed=None):
"Create list of projections to be added together to make the light cone."
@@ -109,7 +113,7 @@
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 delta z of last output in solution list.
+ # until one is within max delta z of last output in solution list.
else:
output = self.allOutputs[0]
while (z > output['redshift']):
@@ -120,17 +124,15 @@
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]['deltaz']),
+ ((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]['deltaz']
+ z = self.lightConeSolution[-1]['redshift'] - self.lightConeSolution[-1]['deltazMax']
# Make light cone using maximum number of projections (minimum spacing).
else:
- deltazMin = 0.01
-
# 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.
@@ -140,7 +142,7 @@
while (nextOutput is not None):
if (nextOutput['redshift'] <= self.lightConeParameters['FinalRedshift']):
break
- if ((self.lightConeSolution[-1]['redshift'] - nextOutput['redshift']) > deltazMin):
+ if ((self.lightConeSolution[-1]['redshift'] - nextOutput['redshift']) > self.lightConeSolution[-1]['deltazMin']):
self.lightConeSolution.append(nextOutput)
nextOutput = nextOutput['next']
@@ -168,7 +170,7 @@
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]['deltaz'],
+ 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))
# Calculate fraction of box required for width corresponding to requested image size.
@@ -186,9 +188,6 @@
self.masterSolution = [copy.deepcopy(q) for q in self.lightConeSolution]
# Clear out some stuff.
- del self.allOutputs
- del self.redshiftOutputs
- del self.timeOutputs
del co
def GetHaloMask(self,HaloMaskParameterFile,mask_file=None,**kwargs):
@@ -209,7 +208,7 @@
self.haloMask *= mask
del haloMaskCube
- def ProjectLightCone(self,field,weight_field=None,apply_halo_mask=False,**kwargs):
+ def ProjectLightCone(self,field,weight_field=None,apply_halo_mask=False,save_stack=True,save_slice_images=False,**kwargs):
"Create projections for light cone, then add them together."
# Clear projection stack.
@@ -226,51 +225,59 @@
q,len(self.lightConeSolution))
output['object'] = lagos.EnzoStaticOutput(output['filename'])
frb = LightConeProjection(output,field,self.pixels,weight_field=weight_field,
- save_image=self.lightConeParameters['SaveLightConeSlices'],
+ save_image=save_slice_images,
name=name,**kwargs)
- 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'])
- else:
- self.projectionStack.append(frb[field])
+ 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'])
+ else:
+ self.projectionStack.append(frb[field])
# 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):
del output['object']
- # Add up slices to make light cone projection.
- if (weight_field is None):
- lightConeProjection = sum(self.projectionStack)
- else:
- lightConeProjection = sum(self.projectionStack) / sum(self.projectionWeightFieldStack)
+ if ytcfg.getint("yt","__parallel_rank") == 0:
+ # Add up slices to make light cone projection.
+ if (weight_field is None):
+ lightConeProjection = sum(self.projectionStack)
+ else:
+ lightConeProjection = sum(self.projectionStack) / sum(self.projectionWeightFieldStack)
- filename = "%s%s" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
+ filename = "%s%s" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
- # Save the last fixed resolution buffer for the plot collection,
- # but replace the data with the full light cone projection data.
- frb.data[field] = lightConeProjection
-
- # Apply halo mask.
- if apply_halo_mask:
- if len(self.haloMask) > 0:
- mylog.info("Applying halo mask.")
- frb.data[field] *= self.haloMask
- else:
- mylog.error("No halo mask loaded, call GetHaloMask.")
+ # Save the last fixed resolution buffer for the plot collection,
+ # but replace the data with the full light cone projection data.
+ frb.data[field] = lightConeProjection
+
+ # Write stack to hdf5 file.
+ if save_stack:
+ self._SaveLightConeStack(field=field,weight_field=weight_field,filename=filename)
+
+ # Apply halo mask.
+ if apply_halo_mask:
+ if len(self.haloMask) > 0:
+ mylog.info("Applying halo mask.")
+ frb.data[field] *= self.haloMask
+ else:
+ mylog.error("No halo mask loaded, call GetHaloMask.")
- # 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)
- pc.save(filename)
+ # 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)
+ pc.save(filename)
- # Return the plot collection so the user can remake the plot if they want.
- return pc
+ # Return the plot collection so the user can remake the plot if they want.
+ return pc
+ else:
+ mylog.info("I'm not the root process so I'm just going to chill.")
def RerandomizeLightConeSolution(self,newSeed,recycle=True):
"""
@@ -383,32 +390,46 @@
output['ProjectionAxis'],output['ProjectionCenter'][0],output['ProjectionCenter'][1],output['ProjectionCenter'][2]))
f.close()
- def SaveLightConeStack(self,filename=None):
- "Save the light cone projection stack as a 3d array in and hdf5 file."
- if (filename is None):
- filename = "%s/%s_data" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
- if not(filename.endswith('.h5')):
- filename += ".h5"
+ 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'])
- if (len(self.projectionStack) == 0):
- if self.verbose: mylog.error("SaveLightConeStack: no projection data loaded.")
- return
+ d_Tolerance = 1e-4
+ max_Iterations = 100
- if self.verbose: mylog.info("Writing light cone data to %s." % filename)
+ targetDistance = self.enzoParameters['CosmologyComovingBoxSize']
- output = h5.openFile(filename, "a")
+ for output in self.allOutputs:
+ z = output['redshift']
- self.projectionStack = na.array(self.projectionStack)
- output.createArray("/", "data",self.projectionStack)
+ # 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
- if (len(self.projectionWeightFieldStack) > 0):
- self.projectionWeightFieldStack = na.array(self.projectionWeightFieldStack)
- output.createArray("/","weight",self.projectionWeightFieldStack)
+ # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
+ distance2 = co.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
- output.close()
+ 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)
- def _CalculateDeltaZMax(self):
- "Calculate delta z that corresponds to full box length going from z to (z - delta 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'])
@@ -416,30 +437,32 @@
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 the box at a given 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.1 # just an initial guess
+ 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-self.enzoParameters['CosmologyComovingBoxSize'])/distance2) > d_Tolerance):
+ while ((na.fabs(distance2 - targetDistance) / distance2) > d_Tolerance):
m = (distance2 - distance1) / (z2 - z1)
z1 = z2
distance1 = distance2
- z2 = ((self.enzoParameters['CosmologyComovingBoxSize'] - distance2) / m) + z2
+ 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['deltaz'] = na.fabs(z2-z)
+ output['deltazMin'] = na.fabs(z2-z)
del co
@@ -481,6 +504,8 @@
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."
@@ -540,6 +565,53 @@
else:
self.lightConeParameters[param] = t
+ def _SaveLightConeStack(self,field=None,weight_field=None,filename=None):
+ "Save the light cone projection stack as a 3d array in and hdf5 file."
+
+ field_node = "%s_%s" % (field,weight_field)
+ weight_field_node = "weight_field_%s" % weight_field
+
+ import tables
+ if (filename is None):
+ filename = "%s/%s_data" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
+ if not(filename.endswith('.h5')):
+ filename += ".h5"
+
+ if (len(self.projectionStack) == 0):
+ if self.verbose: mylog.error("SaveLightConeStack: no projection data loaded.")
+ return
+
+ if self.verbose: mylog.info("Writing light cone data to %s." % filename)
+
+ output = tables.openFile(filename, "a")
+
+ try:
+ node_exists = output.isVisibleNode("/%s" % field_node)
+ except tables.exceptions.NoSuchNodeError:
+ node_exists = False
+
+ if node_exists:
+ mylog.error("Dataset, %s, already exists in %s, not saving." % (field_node,filename))
+ else:
+ mylog.info("Saving %s to %s." % (field_node, filename))
+ self.projectionStack = na.array(self.projectionStack)
+ output.createArray("/",field_node,self.projectionStack)
+
+ if (len(self.projectionWeightFieldStack) > 0):
+ try:
+ node_exists = output.isVisibleNode("/%s" % weight_field_node)
+ except tables.exceptions.NoSuchNodeError:
+ node_exists = False
+
+ if node_exists:
+ mylog.error("Dataset, %s, already exists in %s, not saving." % (weight_field_node,filename))
+ else:
+ mylog.info("Saving %s to %s." % (weight_field_node, filename))
+ self.projectionWeightFieldStack = na.array(self.projectionWeightFieldStack)
+ output.createArray("/",weight_field_node,self.projectionWeightFieldStack)
+
+ output.close()
+
def _SetParameterDefaults(self):
"Set some default parameters to avoid problems if they are not in the parameter file."
self.enzoParameters['GlobalDir'] = ""
@@ -558,6 +630,7 @@
"CosmologyHubbleConstantNow": float,
"CosmologyInitialRedshift": float,
"CosmologyFinalRedshift": float,
+ "TopGridDimensions": float,
"dtDataDump": float,
"RedshiftDumpName": str,
"RedshiftDumpDir": str,
@@ -571,6 +644,5 @@
"ImageResolutionInArcSeconds": float,
"RandomSeed": int,
"UseMinimumNumberOfProjections": int,
- "SaveLightConeSlices": int,
"OutputDir": str,
"OutputPrefix": str}
Modified: trunk/yt/extensions/lightcone/LightConeProjection.py
==============================================================================
--- trunk/yt/extensions/lightcone/LightConeProjection.py (original)
+++ trunk/yt/extensions/lightcone/LightConeProjection.py Tue Dec 16 16:41:22 2008
@@ -25,6 +25,7 @@
from yt.extensions.lightcone import *
from yt.logger import lagosLogger as mylog
+from yt.config import ytcfg
import yt.lagos as lagos
import yt.raven as raven
import numpy as na
@@ -47,7 +48,7 @@
# ProjectionCenter[ProjectionAxis]
# DepthBoxFraction
- node_name = "LightCone_%s_%s_%d_%f_%f" % (field,weight_field,lightConeSlice['ProjectionAxis'],
+ node_name = "LightCone_%d_%f_%f" % (lightConeSlice['ProjectionAxis'],
lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']],
lightConeSlice['DepthBoxFraction'])
@@ -59,7 +60,7 @@
# Make plot collection.
region_center = [0.5 * (dataset_object.parameters['DomainRightEdge'][q] +
dataset_object.parameters['DomainLeftEdge'][q]) \
- for q in range(len(dataset_object.parameters['DomainLeftEdge']))]
+ for q in range(dataset_object.parameters['TopGridRank'])]
pc = raven.PlotCollection(dataset_object,center=region_center)
# 1. The Depth Problem
@@ -85,9 +86,6 @@
pc.add_projection(field,lightConeSlice['ProjectionAxis'],weight_field=weight_field,field_cuts=field_cuts,use_colorbar=True,
node_name=node_name,**kwargs)
- # Serialize projection data.
- pc.plots[0].data._serialize(node_name,force=True)
-
# 2. The Tile Problem
# Tile projection to specified width.
@@ -203,12 +201,16 @@
pc.set_ylim(0,lightConeSlice['WidthBoxFraction'])
# Save an image if requested.
- if (save_image):
+ if (save_image and (ytcfg.getint("yt","__parallel_rank") == 0)):
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(pc.plots[0].data,(0,lightConeSlice['WidthBoxFraction'],0,lightConeSlice['WidthBoxFraction']),
- (pixels,pixels),antialias=False)
-
- return frb
+ if ytcfg.getint("yt","__parallel_rank") == 0:
+ # 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(pc.plots[0].data,(0,lightConeSlice['WidthBoxFraction'],0,lightConeSlice['WidthBoxFraction']),
+ (pixels,pixels),antialias=False)
+
+ return frb
+ else:
+ # If running in parallel and this is not the root process, return None.
+ return None
More information about the yt-svn
mailing list