[Yt-svn] yt-commit r1376 - branches/yt-1.5/yt branches/yt-1.5/yt/extensions branches/yt-1.5/yt/extensions/lightcone branches/yt-1.5/yt/lagos branches/yt-1.5/yt/raven trunk/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri Jul 17 07:30:17 PDT 2009
Author: mturk
Date: Fri Jul 17 07:30:13 2009
New Revision: 1376
URL: http://yt.spacepope.org/changeset/1376
Log:
Backporting a bunch of changes to yt-1.5, including the halo profiler command
line (closes #95 when it's documented), some new mods to the plotting command,
fixes to the light cone generator for parallelization, a bunch of
parallelization fixes, inline (not a target for this release but doesn't hurt)
fixes, h5py version fixes and upgrade path for pytables fixes, fixes to the
non-ortho ray object (but not the .c file, which will be in the next commit),
some field fixes, the plot collection mods from Adam Ginsburg and Sam Skillman,
and that's about it!
Modified:
branches/yt-1.5/yt/commands.py
branches/yt-1.5/yt/extensions/HaloProfiler.py
branches/yt-1.5/yt/extensions/lightcone/LightCone.py
branches/yt-1.5/yt/lagos/BaseDataTypes.py
branches/yt-1.5/yt/lagos/DataReadingFuncs.py
branches/yt-1.5/yt/lagos/EnzoFields.py
branches/yt-1.5/yt/lagos/HierarchyType.py
branches/yt-1.5/yt/lagos/OutputTypes.py
branches/yt-1.5/yt/lagos/ParallelTools.py
branches/yt-1.5/yt/lagos/RTIntegrator.pyx
branches/yt-1.5/yt/lagos/UniversalFields.py
branches/yt-1.5/yt/raven/FixedResolution.py
branches/yt-1.5/yt/raven/PlotCollection.py
branches/yt-1.5/yt/raven/PlotTypes.py
trunk/yt/lagos/BaseDataTypes.py
Modified: branches/yt-1.5/yt/commands.py
==============================================================================
--- branches/yt-1.5/yt/commands.py (original)
+++ branches/yt-1.5/yt/commands.py Fri Jul 17 07:30:13 2009
@@ -27,7 +27,7 @@
from yt.funcs import *
from yt.recipes import _fix_pf
import yt.cmdln as cmdln
-import optparse, os, os.path, math
+import optparse, os, os.path, math, sys
_common_options = dict(
axis = dict(short="-a", long="--axis",
@@ -38,6 +38,10 @@
action="store_true",
dest="takelog", default=True,
help="Take the log of the field?"),
+ text = dict(short="-t", long="--text",
+ action="store", type="string",
+ dest="text", default=None,
+ help="Textual annotation"),
field = dict(short="-f", long="--field",
action="store", type="string",
dest="field", default="Density",
@@ -55,6 +59,11 @@
dest="zlim", default=None,
nargs=2,
help="Color limits (min, max)"),
+ dex = dict(short="", long="--dex",
+ action="store", type="float",
+ dest="dex", default=None,
+ nargs=1,
+ help="Number of dex above min to display"),
width = dict(short="-w", long="--width",
action="store", type="float",
dest="width", default=1.0,
@@ -133,6 +142,32 @@
action="store_true",
dest="grids", default=False,
help="Show the grid boundaries"),
+ halos = dict(short="", long="--halos",
+ action="store", type="string",
+ dest="halos",default="multiple",
+ help="Run halo profiler on a 'single' halo or 'multiple' halos."),
+ halo_radius = dict(short="", long="--halo_radius",
+ action="store", type="float",
+ dest="halo_radius",default=0.1,
+ help="Constant radius for profiling halos if using hop output files with no radius entry. Default: 0.1."),
+ halo_radius_units = dict(short="", long="--halo_radius_units",
+ action="store", type="string",
+ dest="halo_radius_units",default="1",
+ help="Units for radius used with --halo_radius flag. Default: '1' (code units)."),
+ halo_hop_style = dict(short="", long="--halo_hop_style",
+ action="store", type="string",
+ dest="halo_hop_style",default="new",
+ help="Style of hop output file. 'new' for yt_hop files and 'old' for enzo_hop files."),
+ halo_parameter_file = dict(short="", long="--halo_parameter_file",
+ action="store", type="string",
+ dest="halo_parameter_file",default=None,
+ help="HaloProfiler parameter file."),
+ make_profiles = dict(short="", long="--make_profiles",
+ action="store_true", default=False,
+ help="Make profiles with halo profiler."),
+ make_projections = dict(short="", long="--make_projections",
+ action="store_true", default=False,
+ help="Make projections with halo profiler.")
)
def _add_options(parser, *options):
@@ -201,8 +236,29 @@
else: fn = opts.output
hop_list.write_out(fn)
+ @add_cmd_options(['make_profiles','make_projections','halo_parameter_file',
+ 'halos','halo_hop_style','halo_radius','halo_radius_units'])
+ def do_halos(self, subcmd, opts, arg):
+ """
+ Run HaloProfiler on one dataset.
+
+ ${cmd_option_list}
+ """
+ import yt.extensions.HaloProfiler as HP
+ kwargs = {'halos': opts.halos,
+ 'hop_style': opts.halo_hop_style,
+ 'radius': opts.halo_radius,
+ 'radius_units': opts.halo_radius_units}
+
+ hp = HP.HaloProfiler(arg,opts.halo_parameter_file,**kwargs)
+ if opts.make_profiles:
+ hp.makeProfiles()
+ if opts.make_projections:
+ hp.makeProjections()
+
@add_cmd_options(["maxw", "minw", "proj", "axis", "field", "weight",
- "zlim", "nframes", "output", "cmap", "uboxes"])
+ "zlim", "nframes", "output", "cmap", "uboxes", "dex",
+ "text"])
def do_zoomin(self, subcmd, opts, arg):
"""
Create a set of zoomin frames
@@ -217,11 +273,14 @@
axes = [opts.axis]
pc = PlotCollection(pf)
for ax in axes:
- if opts.projection: pc.add_projection(opts.field, ax,
+ if opts.projection: p = pc.add_projection(opts.field, ax,
weight_field=opts.weight)
- else: pc.add_slice(opts.field, ax)
- if opts.unit_boxes: pc.plots[-1].add_callback(
- UnitBoundaryCallback(factor=8))
+ else: p = pc.add_slice(opts.field, ax)
+ if opts.unit_boxes: p.modify["units"](factor=8)
+ if opts.text is not None:
+ p.modify["text"](
+ (0.02, 0.05), opts.text.replace(r"\n", "\n"),
+ text_args=dict(size="medium", color="w"))
pc.set_width(opts.max_width,'1')
# Check the output directory
if not os.path.isdir(opts.output):
@@ -238,10 +297,13 @@
mylog.info("Setting width to %0.3e", w)
mylog.info("Saving frame %06i",i)
pc.set_width(w,"1")
- if opts.zlim: pc.set_zlim(*opts.zlim)
+ if opts.zlim:
+ pc.set_zlim(*opts.zlim)
+ if opts.dex:
+ pc.set_zlim('min', None, opts.dex)
pc.set_cmap(opts.cmap)
pc.save(os.path.join(opts.output,"%s_frame%06i" % (pf,i)))
- w *= factor
+ w = factor**i
@add_cmd_options(["width", "unit", "bn", "proj", "center",
"zlim", "axis", "field", "weight", "skip",
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 Fri Jul 17 07:30:13 2009
@@ -195,7 +195,6 @@
def makeProjections(self,save_images=True,save_cube=True,**kwargs):
"Make projections of all halos using specified fields."
-
# Get virial quantities.
self._LoadVirialData()
@@ -252,8 +251,9 @@
y_axis = coords[1]
for field in self.projectionFields.keys():
- pc.add_projection(field,w,weight_field=self.projectionFields[field],source=region,**kwargs)
-
+ pc.add_projection(field,w,weight_field=self.projectionFields[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']
@@ -281,16 +281,16 @@
frb = raven.FixedResolutionBuffer(pc.plots[e].data,(proj_left[0],proj_right[0],proj_left[1],proj_right[1]),
(projectionResolution,projectionResolution),
antialias=False)
- output.create_dataset("/%s" % field,data=frb[field])
+ dataset_name = "%s_%s" % (field,self.projectionFields[field])
+ if dataset_name in output.listnames(): del output[dataset_name]
+ output.create_dataset(dataset_name,data=frb[field])
output.close()
if save_images:
- pc.save("%s/Halo_%04d" % (outputDir,halo['id']))
+ pc.save("%s/Halo_%04d" % (outputDir,halo['id']),force_save=True)
pc.clear_plots()
-
del region
-
del pc
@lagos.parallel_root_only
@@ -314,7 +314,7 @@
rho_crit_now = 1.8788e-29 * self.pf['CosmologyHubbleConstantNow']**2.0 # g cm^-3
Msun2g = 1.989e33
- rho_crit = rho_crit_now * ((1 + self.pf['CosmologyCurrentRedshift'])**3.0)
+ rho_crit = rho_crit_now * ((1.0 + self.pf['CosmologyCurrentRedshift'])**3.0)
profile['ActualOverdensity'] = (Msun2g * profile['TotalMassMsun']) / \
profile['CellVolume'] / rho_crit
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 Fri Jul 17 07:30:13 2009
@@ -33,6 +33,7 @@
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):
@@ -79,8 +80,10 @@
# Calculate redshifts for dt data dumps.
if (self.enzoParameters.has_key('dtDataDump')):
- self._CalculateTimeDumpRedshifts()
-
+ if self.lightConeParameters['dtSplit'] is None:
+ self._CalculateTimeDumpRedshifts()
+ else:
+ self._CalculateSplitTimeDumpRedshifts(self.lightConeParameters['dtSplit'],self.lightConeParameters['dt_new'])
# Combine all data dumps.
self._CombineDataOutputs()
@@ -250,7 +253,7 @@
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,**kwargs):
+ save_stack=True,save_slice_images=False,flatten_stack=False,photon_field=False,**kwargs):
"Create projections for light cone, then add them together."
# Clear projection stack.
@@ -273,6 +276,19 @@
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:
+ if photon_field:
+ #need to 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
+ boxSizeProper = self.enzoParameters['CosmologyComovingBoxSize'] / (self.enzoParameters['CosmologyHubbleConstantNow'] *
+ (1.0 + output['redshift']))
+ pixelarea = (boxSizeProper/self.pixels)**2 #in proper cm^2
+ factor = pixelarea/(4.0*na.pi*dL**2)
+ 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):
@@ -334,7 +350,7 @@
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.add_fixed_resolution_plot(frb,field,**kwargs)
pc.save(filename)
# Return the plot collection so the user can remake the plot if they want.
@@ -582,6 +598,35 @@
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
@@ -686,6 +731,7 @@
if node_exists:
if over_write:
mylog.info("Dataset, %s, already exists, overwriting." % field_node)
+ write_data = True
del output[field_node]
else:
mylog.info("Dataset, %s, already exists in %s, not saving." % (field_node,filename))
@@ -726,7 +772,7 @@
def _SetParameterDefaults(self):
"Set some default parameters to avoid problems if they are not in the parameter file."
- self.enzoParameters['GlobalDir'] = ""
+ self.enzoParameters['GlobalDir'] = "./"
self.enzoParameters['RedshiftDumpName'] = "RD"
self.enzoParameters['RedshiftDumpDir'] = "RD"
self.enzoParameters['DataDumpName'] = "DD"
@@ -737,6 +783,8 @@
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,
@@ -762,4 +810,6 @@
"MinimumCoherentBoxFraction": float,
"MinimumSliceDeltaz": float,
"OutputDir": str,
- "OutputPrefix": str}
+ "OutputPrefix": str,
+ "dtSplit": float,
+ "dt_new": float}
Modified: branches/yt-1.5/yt/lagos/BaseDataTypes.py
==============================================================================
--- branches/yt-1.5/yt/lagos/BaseDataTypes.py (original)
+++ branches/yt-1.5/yt/lagos/BaseDataTypes.py Fri Jul 17 07:30:13 2009
@@ -376,7 +376,7 @@
continue
mylog.info("Getting field %s from %s", field, len(self._grids))
if field not in self.hierarchy.field_list and not in_grids:
- if self._generate_field(field):
+ if field != "dts" and self._generate_field(field):
continue # True means we already assigned it
self[field] = na.concatenate(
[self._get_data_from_grid(grid, field)
@@ -439,11 +439,13 @@
"""
mylog.warning("This code is poorly tested. It may give bad data!")
AMR1DData.__init__(self, pf, fields, **kwargs)
- self.start_point = na.array(start_point)
- self.end_point = na.array(end_point)
+ self.start_point = na.array(start_point, dtype='float64')
+ self.end_point = na.array(end_point, dtype='float64')
self.vec = self.end_point - self.start_point
+ #self.vec /= na.sqrt(na.dot(self.vec, self.vec))
self.center = self.start_point
self.set_field_parameter('center', self.start_point)
+ self._dts = {}
#self._refresh_data()
def _get_list_of_grids(self):
@@ -456,16 +458,16 @@
i1 = (i+1) % 3
i2 = (i+2) % 3
vs = self._get_line_at_coord(LE[:,i], i)
- p = p | ( ( (LE[:,i1] < vs[:,i1]) & (RE[:,i1] > vs[:,i1]) ) \
- & ( (LE[:,i2] < vs[:,i2]) & (RE[:,i2] > vs[:,i2]) ) )
+ p = p | ( ( (LE[:,i1] <= vs[:,i1]) & (RE[:,i1] >= vs[:,i1]) ) \
+ & ( (LE[:,i2] <= vs[:,i2]) & (RE[:,i2] >= vs[:,i2]) ) )
vs = self._get_line_at_coord(RE[:,i], i)
- p = p | ( ( (LE[:,i1] < vs[:,i1]) & (RE[:,i1] > vs[:,i1]) ) \
- & ( (LE[:,i2] < vs[:,i2]) & (RE[:,i2] > vs[:,i2]) ) )
- p = p | ( na.all( LE < self.start_point, axis=1 )
- & na.all( RE > self.start_point, axis=1 ) )
- p = p | ( na.all( LE < self.end_point, axis=1 )
- & na.all( RE > self.end_point, axis=1 ) )
- self._grids = self.hierarchy.grids.copy()#[p]
+ p = p | ( ( (LE[:,i1] <= vs[:,i1]) & (RE[:,i1] >= vs[:,i1]) ) \
+ & ( (LE[:,i2] <= vs[:,i2]) & (RE[:,i2] >= vs[:,i2]) ) )
+ p = p | ( na.all( LE <= self.start_point, axis=1 )
+ & na.all( RE >= self.start_point, axis=1 ) )
+ p = p | ( na.all( LE <= self.end_point, axis=1 )
+ & na.all( RE >= self.end_point, axis=1 ) )
+ self._grids = self.hierarchy.grids[p]
def _get_line_at_coord(self, v, index):
# t*self.vec + self.start_point = self.end_point
@@ -476,14 +478,17 @@
def _get_data_from_grid(self, grid, field):
mask = na.logical_and(self._get_cut_mask(grid),
grid.child_mask)
+ if field == 'dts': return self._dts[grid.id][mask]
return grid[field][mask]
@cache_mask
def _get_cut_mask(self, grid):
mask = na.zeros(grid.ActiveDimensions, dtype='int')
+ dts = na.zeros(grid.ActiveDimensions, dtype='float64')
import RTIntegrator as RT
- RT.VoxelTraversal(mask, grid.LeftEdge, grid.RightEdge,
+ RT.VoxelTraversal(mask, dts, grid.LeftEdge, grid.RightEdge,
grid.dds, self.center, self.vec)
+ self._dts[grid.id] = na.abs(dts)
return mask
class AMR2DData(AMRData, GridPropertiesMixin, ParallelAnalysisInterface):
@@ -528,9 +533,11 @@
continue # A "True" return means we did it
# To ensure that we use data from this object as much as possible,
# we're going to have to set the same thing several times
- temp_data[field] = na.concatenate(
- [self._get_data_from_grid(grid, field)
- for grid in self._get_grids()])
+ data = [self._get_data_from_grid(grid, field)
+ for grid in self._get_grids()]
+ if len(data) == 0: data = None
+ else: data = na.concatenate(data)
+ temp_data[field] = data
# Now the next field can use this field
self[field] = temp_data[field]
# We finalize
@@ -675,7 +682,9 @@
points = []
for grid in self._get_grids():
points.append(self._generate_grid_coords(grid))
- t = self._mpi_catarray(na.concatenate(points))
+ if len(points) == 0: points = None
+ else: points = na.concatenate(points)
+ t = self._mpi_catarray(points)
self['px'] = t[:,0]
self['py'] = t[:,1]
self['pz'] = t[:,2]
@@ -755,6 +764,13 @@
return "%s/%s_%s" % \
(self._top_node, self.axis, self.coord)
+ def __get_quantities(self):
+ if self.__quantities is None:
+ self.__quantities = DerivedQuantityCollection(self)
+ return self.__quantities
+ __quantities = None
+ quantities = property(__get_quantities)
+
class AMRCuttingPlaneBase(AMR2DData):
"""
AMRCuttingPlane is an oblique plane through the data,
@@ -847,7 +863,9 @@
points = []
for grid in self._get_grids():
points.append(self._generate_grid_coords(grid))
- t = self._mpi_catarray(na.concatenate(points))
+ if len(points) == 0: points = None
+ else: points = na.concatenate(points)
+ t = self._mpi_catarray(points)
pos = (t[:,0:3] - self.center)
self['px'] = na.dot(pos, self._x_vec)
self['py'] = na.dot(pos, self._y_vec)
@@ -883,7 +901,7 @@
return na.where(k)
def _gen_node_name(self):
- cen_name = ("%s" % self.center).replace(" ","_")[1:-1]
+ cen_name = ("%s" % (self.center,)).replace(" ","_")[1:-1]
L_name = ("%s" % self._norm_vec).replace(" ","_")[1:-1]
return "%s/c%s_L%s" % \
(self._top_node, cen_name, L_name)
@@ -895,7 +913,8 @@
_con_args = ('axis', 'field', 'weight_field')
def __init__(self, axis, field, weight_field = None,
max_level = None, center = None, pf = None,
- source=None, node_name = None, field_cuts = None, **kwargs):
+ source=None, node_name = None, field_cuts = None,
+ serialize=True,**kwargs):
"""
AMRProj is a projection of a *field* along an *axis*. The field
can have an associated *weight_field*, in which case the values are
@@ -904,6 +923,7 @@
"""
AMR2DData.__init__(self, axis, field, pf, node_name = None, **kwargs)
self._field_cuts = field_cuts
+ self.serialize = serialize
self.center = center
if center is not None: self.set_field_parameter('center',center)
self._node_name = node_name
@@ -923,7 +943,7 @@
self._temp = {}
self._deserialize(node_name)
self._refresh_data()
- if self._okay_to_serialize: self._serialize(node_name=self._node_name)
+ if self._okay_to_serialize and self.serialize: self._serialize(node_name=self._node_name)
def _convert_field_name(self, field):
if field == "weight_field": return "weight_field_%s" % self._weight
@@ -1137,7 +1157,7 @@
field_data = na.vsplit(data.pop('fields'), len(fields))
for fi, field in enumerate(fields):
self[field] = field_data[fi].ravel()
- self._store_fields(field, self._node_name)
+ if self.serialize: self._store_fields(field, self._node_name)
for i in data.keys(): self[i] = data.pop(i)
def add_fields(self, fields, weight = "CellMassMsun"):
Modified: branches/yt-1.5/yt/lagos/DataReadingFuncs.py
==============================================================================
--- branches/yt-1.5/yt/lagos/DataReadingFuncs.py (original)
+++ branches/yt-1.5/yt/lagos/DataReadingFuncs.py Fri Jul 17 07:30:13 2009
@@ -156,12 +156,10 @@
def readDataSliceInMemory(self, grid, field, axis, coord):
import enzo
- sl = [slice(None), slice(None), slice(None)]
- sl[axis] = slice(coord, coord + 1)
- #sl = tuple(reversed(sl))
- sl = tuple(sl)
- bsl = (slice(3,-3), slice(3,-3), slice(3,-3))
- return enzo.grid_data[grid.id][field][bsl][sl]#.swapaxes(0,2)
+ sl = [slice(3,-3), slice(3,-3), slice(3,-3)]
+ sl[axis] = slice(coord + 3, coord + 4)
+ sl = tuple(reversed(sl))
+ return enzo.grid_data[grid.id][field][sl].swapaxes(0,2)
def getExceptionInMemory():
return KeyError
@@ -288,11 +286,11 @@
def _read_set(self, grid, field):
import enzo
if grid.id not in self.grids_in_memory: raise KeyError
- t1 = enzo.yt_parameter_file["InitialTime"]
- t2 = enzo.hierarchy_information["GridOldTimes"][grid.id - 1]
+ return self.grids_in_memory[grid.id][field].swapaxes(0,2)[self.my_slice]
coef1 = max((grid.Time - t1)/(grid.Time - t2), 0.0)
coef2 = 1.0 - coef1
- return self.grids_in_memory[grid.id][field][self.my_slice]
+ t1 = enzo.yt_parameter_file["InitialTime"]
+ t2 = enzo.hierarchy_information["GridOldTimes"][grid.id]
return (coef1*self.grids_in_memory[grid.id][field] + \
coef2*self.old_grids_in_memory[grid.id][field])\
[self.my_slice]
Modified: branches/yt-1.5/yt/lagos/EnzoFields.py
==============================================================================
--- branches/yt-1.5/yt/lagos/EnzoFields.py (original)
+++ branches/yt-1.5/yt/lagos/EnzoFields.py Fri Jul 17 07:30:13 2009
@@ -38,7 +38,7 @@
_speciesList = ["HI","HII","Electron",
"HeI","HeII","HeIII",
"H2I","H2II","HM",
- "DI","DII","HDI","Metal"]
+ "DI","DII","HDI","Metal","PreShock"]
def _SpeciesFraction(field, data):
sp = field.name.split("_")[0] + "_Density"
return data[sp]/data["Density"]
@@ -59,7 +59,7 @@
return data["Total_Energy"]
else:
if data.pf["DualEnergyFormalism"]:
- return data["Gas_Energy"]
+ return data["GasEnergy"]
else:
return data["Total_Energy"] - 0.5*(
data["x-velocity"]**2.0
Modified: branches/yt-1.5/yt/lagos/HierarchyType.py
==============================================================================
--- branches/yt-1.5/yt/lagos/HierarchyType.py (original)
+++ branches/yt-1.5/yt/lagos/HierarchyType.py Fri Jul 17 07:30:13 2009
@@ -157,14 +157,18 @@
"""
if self._data_mode != 'a': return
+ if "ArgsError" in dir(h5py.h5):
+ exception = h5py.h5.ArgsError
+ else:
+ exception = h5py.h5.H5Error
try:
node_loc = self._data_file[node]
if name in node_loc.listnames() and force:
mylog.info("Overwriting node %s/%s", node, name)
del self._data_file[node][name]
elif name in node_loc.listnames() and passthrough:
- return
- except h5py.h5.ArgsError:
+ return
+ except exception:
pass
myGroup = self._data_file['/']
for q in node.split('/'):
@@ -175,10 +179,17 @@
self._data_file.flush()
def _reload_data_file(self, *args, **kwargs):
+ if self._data_file is None: return
self._data_file.close()
del self._data_file
self._data_file = h5py.File(self.__data_filename, self._data_mode)
+ def _reset_save_data(self,round_robin=False):
+ if round_robin:
+ self.save_data = self._save_data
+ else:
+ self.save_data = parallel_splitter(self._save_data, self._reload_data_file)
+
save_data = parallel_splitter(_save_data, _reload_data_file)
def save_object(self, obj, name):
@@ -203,10 +214,18 @@
if self._data_file == None:
return None
if node[0] != "/": node = "/%s" % node
- full_name = "%s/%s" % (node, name)
- if full_name not in self._data_file:
+
+ myGroup = self._data_file['/']
+ for group in node.split('/'):
+ if group:
+ if group not in myGroup.listnames():
+ return None
+ myGroup = myGroup[group]
+ if name not in myGroup.listnames():
return None
- return self._data_file["%s/%s" % (node, name)]
+
+ full_name = "%s/%s" % (node, name)
+ return self._data_file[full_name]
def _close_data_file(self):
if self._data_file:
@@ -840,9 +859,9 @@
self.gridLevels[secondGrid] = self.gridLevels[firstGrid]
pTree = [ [ grid.id - 1 for grid in self.gridTree[i] ] for i in range(self.num_grids) ]
self.gridReverseTree[0] = -1
- self.save_data(cPickle.dumps(pTree, protocol=-1), "/", "Tree")
- self.save_data(na.array(self.gridReverseTree), "/", "ReverseTree")
- self.save_data(self.gridLevels, "/", "Levels")
+ self.save_data(cPickle.dumps(pTree, protocol=-1), "/", "Tree", force=True)
+ self.save_data(na.array(self.gridReverseTree), "/", "ReverseTree", force=True)
+ self.save_data(self.gridLevels, "/", "Levels", force=True)
@parallel_blocking_call
def _populate_hierarchy(self):
@@ -866,12 +885,15 @@
self.__setup_grid_tree()
else:
mylog.debug("Grabbing serialized tree data")
- pTree = cPickle.loads(treeArray.value)
- self.gridReverseTree = list(self.get_data("/","ReverseTree"))
- self.gridTree = [ [ weakref.proxy(self.grids[i]) for i in pTree[j] ]
- for j in range(self.num_grids) ]
- self.gridLevels = self.get_data("/","Levels")[:]
- mylog.debug("Grabbed")
+ try:
+ pTree = cPickle.loads(treeArray.value)
+ self.gridReverseTree = list(self.get_data("/","ReverseTree"))
+ self.gridTree = [ [ weakref.proxy(self.grids[i]) for i in pTree[j] ]
+ for j in range(self.num_grids) ]
+ self.gridLevels = self.get_data("/","Levels")[:]
+ mylog.debug("Grabbed")
+ except EOFError:
+ self.__setup_grid_tree()
for i,v in enumerate(self.gridReverseTree):
# For multiple grids on the root level
if v == -1: self.gridReverseTree[i] = None
@@ -962,8 +984,14 @@
return self.grids[(random_sample,)]
class EnzoHierarchyInMemory(EnzoHierarchy):
- def __init__(self, pf, data_style = 8):
- import enzo
+ _data_style = 8
+ def _obtain_enzo(self):
+ import enzo; return enzo
+
+ def __init__(self, pf, data_style = None):
+ self.parameter_file = weakref.proxy(pf) # for _obtain_enzo
+ if data_style is None: data_style = self._data_style
+ enzo = self._obtain_enzo()
self.float_type = 'float64'
self.data_style = data_style # Mandated
self.directory = os.getcwd()
@@ -980,16 +1008,16 @@
data = list(field_list)
if MPI.COMM_WORLD.rank == 0:
for i in range(1, MPI.COMM_WORLD.size):
- data += MPI.COMM_WORLD.Recv(source=i, tag=0)
+ data += MPI.COMM_WORLD.recv(source=i, tag=0)
data = list(set(data))
else:
- MPI.COMM_WORLD.Send(data, dest=0, tag=0)
+ MPI.COMM_WORLD.send(data, dest=0, tag=0)
MPI.COMM_WORLD.Barrier()
- return MPI.COMM_WORLD.Bcast(data, root=0)
+ return MPI.COMM_WORLD.bcast(data, root=0)
def _populate_hierarchy(self):
self._copy_hierarchy_structure()
- import enzo
+ enzo = self._obtain_enzo()
mylog.debug("Copying reverse tree")
self.gridReverseTree = enzo.hierarchy_information["GridParentIDs"].ravel().tolist()
# Initial setup:
@@ -1024,7 +1052,7 @@
mylog.debug("Hierarchy fully populated.")
def _copy_hierarchy_structure(self):
- import enzo
+ enzo = self._obtain_enzo()
self.gridDimensions[:] = enzo.hierarchy_information["GridDimensions"][:]
self.gridStartIndices[:] = enzo.hierarchy_information["GridStartIndices"][:]
self.gridEndIndices[:] = enzo.hierarchy_information["GridEndIndices"][:]
@@ -1046,6 +1074,9 @@
random_sample = na.mgrid[0:max(len(gg)-1,1)].astype("int32")
return gg[(random_sample,)]
+ def save_data(self, *args, **kwargs):
+ pass
+
class EnzoHierarchy1D(EnzoHierarchy):
def __init__(self, *args, **kwargs):
EnzoHierarchy.__init__(self, *args, **kwargs)
Modified: branches/yt-1.5/yt/lagos/OutputTypes.py
==============================================================================
--- branches/yt-1.5/yt/lagos/OutputTypes.py (original)
+++ branches/yt-1.5/yt/lagos/OutputTypes.py Fri Jul 17 07:30:13 2009
@@ -43,7 +43,7 @@
output_type_registry[name]=cls
mylog.debug("Registering: %s as %s", name, cls)
- def __new__(cls, filename, *args, **kwargs):
+ def __new__(cls, filename=None, *args, **kwargs):
apath = os.path.abspath(filename)
if not os.path.exists(apath): raise IOError(filename)
if apath not in _cached_pfs:
@@ -84,10 +84,13 @@
return self.basename
def _hash(self):
- import hashlib
s = "%s;%s;%s" % (self.basename,
self["InitialTime"], self["CurrentTimeIdentifier"])
- return hashlib.md5(s).hexdigest()
+ try:
+ import hashlib
+ return hashlib.md5(s).hexdigest()
+ except ImportError:
+ return s.replace(";", "*")
@classmethod
def _is_valid(cls, *args, **kwargs):
@@ -393,16 +396,27 @@
class EnzoStaticOutputInMemory(EnzoStaticOutput):
_hierarchy_class = EnzoHierarchyInMemory
+ _data_style = 8
+
+ def __new__(cls, *args, **kwargs):
+ obj = object.__new__(cls)
+ obj.__init__(*args, **kwargs)
+ return obj
+
def __init__(self, parameter_override=None, conversion_override=None):
if parameter_override is None: parameter_override = {}
self.__parameter_override = parameter_override
if conversion_override is None: conversion_override = {}
self.__conversion_override = conversion_override
- StaticOutput.__init__(self, "InMemoryParameterFile", 8)
+ StaticOutput.__init__(self, "InMemoryParameterFile", self._data_style)
+
+ self.field_info = self._fieldinfo_class()
def _parse_parameter_file(self):
- import enzo
+ enzo = self._obtain_enzo()
+ self.basename = "cycle%08i" % (
+ enzo.yt_parameter_file["NumberOfPythonCalls"])
self.parameters['CurrentTimeIdentifier'] = time.time()
self.parameters.update(enzo.yt_parameter_file)
self.conversion_factors.update(enzo.conversion_factors)
@@ -417,6 +431,9 @@
for p, v in self.__conversion_override.items():
self.conversion_factors[p] = v
+ def _obtain_enzo(self):
+ import enzo; return enzo
+
@classmethod
def _is_valid(cls, *args, **kwargs):
return False
Modified: branches/yt-1.5/yt/lagos/ParallelTools.py
==============================================================================
--- branches/yt-1.5/yt/lagos/ParallelTools.py (original)
+++ branches/yt-1.5/yt/lagos/ParallelTools.py Fri Jul 17 07:30:13 2009
@@ -276,9 +276,11 @@
MPI.COMM_WORLD.rank)
if MPI.COMM_WORLD.rank == 0:
temp_data = []
+ if data[key] is not None: temp_data.append(data[key])
for i in range(1,np):
- temp_data.append(_recv_array(source=i, tag=0))
- data[key] = na.concatenate([data[key]] + temp_data, axis=-1)
+ buf = _recv_array(source=i, tag=0)
+ if buf is not None: temp_data.append(buf)
+ data[key] = na.concatenate(temp_data, axis=-1)
else:
_send_array(data[key], dest=0, tag=0)
self._barrier()
@@ -324,7 +326,7 @@
# First we receive, then we make a new list.
for i in range(1,MPI.COMM_WORLD.size):
buf = _recv_array(source=i, tag=0)
- data = na.concatenate([data, buf])
+ if buf is not None: data = na.concatenate([data, buf])
return data
@parallel_passthrough
Modified: branches/yt-1.5/yt/lagos/RTIntegrator.pyx
==============================================================================
--- branches/yt-1.5/yt/lagos/RTIntegrator.pyx (original)
+++ branches/yt-1.5/yt/lagos/RTIntegrator.pyx Fri Jul 17 07:30:13 2009
@@ -70,8 +70,10 @@
i_s = o_s[i]
return i_s
+ at cython.wraparound(False)
@cython.boundscheck(False)
def VoxelTraversal(np.ndarray[np.int_t, ndim=3] grid_mask,
+ np.ndarray[np.float64_t, ndim=3] grid_t,
np.ndarray[np.float64_t, ndim=1] left_edge,
np.ndarray[np.float64_t, ndim=1] right_edge,
np.ndarray[np.float64_t, ndim=1] dx,
@@ -81,60 +83,84 @@
# Find the first place the ray hits the grid on its path
# Do left edge then right edge in each dim
cdef int i, x, y
- cdef double tl, tr, intersect_t, intersect
- cdef np.ndarray step = np.ones(3, dtype=np.float64) # maybe just ints?
- cdef np.ndarray cur_ind = np.zeros(3, dtype=np.int64) # maybe just ints?
- cdef np.ndarray tdelta = np.zeros(3, dtype=np.float64)
- cdef np.ndarray tmax = np.zeros(3, dtype=np.float64)
- intersect_t = 1e30
+ cdef np.float64_t tl, tr, intersect_t, enter_t, exit_t, dt_tolerance
+ cdef np.ndarray[np.int64_t, ndim=1] step = np.empty((3,), dtype=np.int64)
+ cdef np.ndarray[np.int64_t, ndim=1] cur_ind = np.empty((3,), dtype=np.int64)
+ cdef np.ndarray[np.float64_t, ndim=1] tdelta = np.empty((3,), dtype=np.float64)
+ cdef np.ndarray[np.float64_t, ndim=1] tmax = np.empty((3,), dtype=np.float64)
+ cdef np.ndarray[np.float64_t, ndim=1] intersect = np.empty((3,), dtype=np.float64)
+ intersect_t = 1
+ dt_tolerance = 1e-6
+ # recall p = v * t + u
+ # where p is position, v is our vector, u is the start point
for i in range(3):
# As long as we're iterating, set some other stuff, too
- if (v[i] < 0): step[i] = -1
+ if(v[i] < 0): step[i] = -1
+ else: step[i] = 1
x = (i+1)%3
y = (i+2)%3
tl = (left_edge[i] - u[i])/v[i]
tr = (right_edge[i] - u[i])/v[i]
- if (left_edge[x] <= (u[x] + tl*v[x]) < right_edge[x]) and \
- (left_edge[y] <= (u[y] + tl*v[y]) < right_edge[y]) and \
- (0 <= tl < intersect_t):
+ if (left_edge[x] <= (u[x] + tl*v[x]) <= right_edge[x]) and \
+ (left_edge[y] <= (u[y] + tl*v[y]) <= right_edge[y]) and \
+ (0.0 <= tl < intersect_t):
intersect_t = tl
- if (left_edge[x] <= (u[x] + tr*v[x]) < right_edge[x]) and \
- (left_edge[y] <= (u[y] + tr*v[y]) < right_edge[y]) and \
- (0 <= tr < intersect_t):
+ if (left_edge[x] <= (u[x] + tr*v[x]) <= right_edge[x]) and \
+ (left_edge[y] <= (u[y] + tr*v[y]) <= right_edge[y]) and \
+ (0.0 <= tr < intersect_t):
intersect_t = tr
+ # if fully enclosed
if (left_edge[0] <= u[0] <= right_edge[0]) and \
(left_edge[1] <= u[1] <= right_edge[1]) and \
(left_edge[2] <= u[2] <= right_edge[2]):
- intersect_t = 0
- if intersect_t > 1e29: return
+ intersect_t = 0.0
+ if not (0 <= intersect_t <= 1): return
# Now get the indices of the intersection
+ intersect = u + intersect_t * v
+ cdef int ncells = 0
for i in range(3):
- intersect = u[i] + intersect_t * v[i]
- cur_ind[i] = np.floor((intersect - left_edge[i])/dx[i])
- if step[i] > 0: tdelta[i] = intersect - ((cur_ind[i]+1) * dx[i])
- if step[i] < 0: tdelta[i] = intersect - (cur_ind[i] * dx[i])
- tmax[i] = dx[i]/v[i]
- # We now know where we have pierced the grid initially.
- cdef int in_cells = 1
+ cur_ind[i] = np.floor((intersect[i] + 1e-8*dx[i] - left_edge[i])/dx[i])
+ tmax[i] = (((cur_ind[i]+step[i])*dx[i])+left_edge[i]-u[i])/v[i]
+ if cur_ind[i] == grid_mask.shape[i] and step[i] < 0:
+ cur_ind[i] = grid_mask.shape[i] - 1
+ if step[i] > 0: tmax[i] = (((cur_ind[i]+1)*dx[i])+left_edge[i]-u[i])/v[i]
+ if step[i] < 0: tmax[i] = (((cur_ind[i]+0)*dx[i])+left_edge[i]-u[i])/v[i]
+ tdelta[i] = np.abs((dx[i]/v[i]))
+ # The variable intersect contains the point we first pierce the grid
+ enter_t = intersect_t
while 1:
- if cur_ind[0] >= grid_mask.shape[0] or \
- cur_ind[1] >= grid_mask.shape[1] or \
- cur_ind[2] >= grid_mask.shape[2]:
+ if (not (0 <= cur_ind[0] < grid_mask.shape[0])) or \
+ (not (0 <= cur_ind[1] < grid_mask.shape[1])) or \
+ (not (0 <= cur_ind[2] < grid_mask.shape[2])):
break
- else:
- grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1
+ # Note that we are calculating t on the fly, but we get *negative* t
+ # values from what they should be.
+ # If we've reached t = 1, we are done.
+ grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1
+ if (tmax[0] > 1.0) and (tmax[1] > 1.0) and (tmax[2] > 1.0):
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = 1.0 - enter_t
+ break
+ ncells += 1
if tmax[0] < tmax[1]:
if tmax[0] < tmax[2]:
- cur_ind[0] += step[0]
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[0] - enter_t
+ enter_t = tmax[0]
tmax[0] += tdelta[0]
+ cur_ind[0] += step[0]
else:
- cur_ind[2] += step[2]
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[2] - enter_t
+ enter_t = tmax[2]
tmax[2] += tdelta[2]
+ cur_ind[2] += step[2]
else:
if tmax[1] < tmax[2]:
- cur_ind[1] += step[1]
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[1] - enter_t
+ enter_t = tmax[1]
tmax[1] += tdelta[1]
+ cur_ind[1] += step[1]
else:
- cur_ind[2] += step[2]
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[2] - enter_t
+ enter_t = tmax[2]
tmax[2] += tdelta[2]
+ cur_ind[2] += step[2]
return
Modified: branches/yt-1.5/yt/lagos/UniversalFields.py
==============================================================================
--- branches/yt-1.5/yt/lagos/UniversalFields.py (original)
+++ branches/yt-1.5/yt/lagos/UniversalFields.py Fri Jul 17 07:30:13 2009
@@ -356,6 +356,10 @@
function=_TotalMass,
convert_function=_convertCellMassMsun)
+def _Matter_Density(field,data):
+ return (data['Density'] + data['Dark_Matter_Density'])
+add_field("Matter_Density",function=_Matter_Density,units=r"\rm{g}/\rm{cm^3}")
+
def _CellVolume(field, data):
if data['dx'].size == 1:
try:
Modified: branches/yt-1.5/yt/raven/FixedResolution.py
==============================================================================
--- branches/yt-1.5/yt/raven/FixedResolution.py (original)
+++ branches/yt-1.5/yt/raven/FixedResolution.py Fri Jul 17 07:30:13 2009
@@ -26,7 +26,6 @@
from yt.raven import *
import _MPL
-
class FixedResolutionBuffer(object):
def __init__(self, data_source, bounds, buff_size, antialias = True):
"""
@@ -92,10 +91,17 @@
provided by STSci to interface with FITS-format files.
"""
import pyfits
+ extra_fields = ['x','y','z','px','py','pz','pdx','pdy','pdz','weight_field']
if filename_prefix.endswith('.fits'): filename_prefix=filename_prefix[:-5]
- if fields is None: fields = self.data.keys()
+ if fields is None:
+ fields = [field for field in self.data_source.fields
+ if field not in extra_fields]
for field in fields:
hdu = pyfits.PrimaryHDU(self[field])
+ if self.data_source.has_key('weight_field'):
+ weightname = self.data_source._weight
+ if weightname is None: weightname = 'None'
+ field = field +'_'+weightname
hdu.writeto("%s_%s.fits" % (filename_prefix, field))
def open_in_ds9(self, field, take_log=True):
Modified: branches/yt-1.5/yt/raven/PlotCollection.py
==============================================================================
--- branches/yt-1.5/yt/raven/PlotCollection.py (original)
+++ branches/yt-1.5/yt/raven/PlotCollection.py Fri Jul 17 07:30:13 2009
@@ -67,7 +67,7 @@
for p in self.plots:
yield p
- def save(self, basename, format="png", override=False):
+ def save(self, basename, format="png", override=False, force_save=False):
"""
Same plots with automatically generated names, prefixed with *basename*
(including directory path) unless *override* is specified, and in
@@ -77,7 +77,7 @@
for plot in self.plots:
fn.append(plot.save_image(basename, \
format=format, submit=self._run_id,
- override=override))
+ override=override, force_save=force_save))
if self.submit:
im = plot.im.copy()
im["Filename"] = self._http_prefix + "/" \
@@ -102,13 +102,23 @@
for plot in self.plots:
plot.set_ylim(ymin, ymax)
- def set_zlim(self, zmin, zmax):
+ def set_zlim(self, zmin, zmax, **kwargs):
"""
- Set the limits of the colorbar.
+ Set the limits of the colorbar. 'min' or 'max' are possible inputs
+ when combined with dex=value, where value gives the maximum number of
+ dex to go above/below the min/max. If value is larger than the true
+ range of values, min/max are limited to true range.
+
+ Only ONE of the following options can be specified. If all 3 are
+ specified, they will be used in the following precedence order:
+ ticks - a list of floating point numbers at which to put ticks
+ minmaxtick - display DEFAULT ticks with min & max also displayed
+ nticks - if ticks not specified, can automatically determine a
+ number of ticks to be evenly spaced in log space
"""
for plot in self.plots:
plot.set_autoscale(False)
- plot.set_zlim(zmin, zmax)
+ plot.set_zlim(zmin, zmax, **kwargs)
def set_lim(self, lim):
"""
Modified: branches/yt-1.5/yt/raven/PlotTypes.py
==============================================================================
--- branches/yt-1.5/yt/raven/PlotTypes.py (original)
+++ branches/yt-1.5/yt/raven/PlotTypes.py Fri Jul 17 07:30:13 2009
@@ -101,7 +101,7 @@
def __getitem__(self, item):
return self.data[item] # Should be returned in CGS
- def save_image(self, prefix, format="png", submit=None, override=True):
+ def save_image(self, prefix, format="png", submit=None, override=True, force_save=False):
"""
Save this plot image. Will generate a filename based on the *prefix*,
*format*. *submit* will govern the submission to the Deliverator and
@@ -115,7 +115,10 @@
my_prefix = prefix
fn = ".".join([my_prefix, format])
canvas = engineVals["canvas"](self._figure)
- only_on_root(canvas.print_figure, fn)
+ if force_save:
+ canvas.print_figure(fn)
+ else:
+ only_on_root(canvas.print_figure, fn)
self["Type"] = self._type_name
self["GeneratedAt"] = self.data.pf["CurrentTimeIdentifier"]
return fn
@@ -138,12 +141,47 @@
"""
self._axes.set_ylim(ymin, ymax)
- def set_zlim(self, zmin, zmax):
+ def set_zlim(self, zmin, zmax, dex=None, nticks=None, ticks=None, minmaxtick=False):
"""
Set the z boundaries of this plot.
+
+ Only ONE of the following options can be specified. If all 3 are
+ specified, they will be used in the following precedence order:
+ ticks - a list of floating point numbers at which to put ticks
+ minmaxtick - display DEFAULT ticks with min & max also displayed
+ nticks - if ticks not specified, can automatically determine a
+ number of ticks to be evenly spaced in log space
"""
# This next call fixes some things, but is slower...
#self._redraw_image()
+ if (zmin in (None,'min')) or (zmax in (None,'max')):
+ imbuff = self._axes.images[-1]._A
+ if zmin == 'min':
+ zmin = na.nanmin(imbuff[na.nonzero(imbuff)])
+ if dex is not None:
+ zmax = min(zmin*10**(dex),na.nanmax(imbuff))
+ if zmax == 'max':
+ zmax = na.nanmax(imbuff)
+ if dex is not None:
+ zmin = max(zmax/(10**(dex)),na.nanmin(imbuff))
+ if ticks is not None:
+ ticks = na.sort(ticks)
+ self.colorbar.locator = matplotlib.ticker.FixedLocator(ticks)
+ self.colorbar.formatter = matplotlib.ticker.FixedFormatter(["%0.2e" % (x) for x in ticks])
+ elif minmaxtick:
+ ticks = na.array(self.colorbar._ticker()[1],dtype='float')
+ ticks = [zmin] + ticks.tolist() + [zmax]
+ self.colorbar.locator = matplotlib.ticker.FixedLocator(ticks)
+ self.colorbar.formatter = matplotlib.ticker.FixedFormatter(["%0.2e" % (x) for x in ticks])
+ elif nticks is not None:
+ lin = na.linspace(na.log10(zmin),na.log10(zmax),nticks)
+ self.colorbar.locator = matplotlib.ticker.FixedLocator(10**lin)
+ self.colorbar.formatter = matplotlib.ticker.FixedFormatter(["%0.2e" % (10**x) for x in lin])
+ else:
+ if hasattr(self,'_old_locator'):
+ self.colorbar.locator = self._old_locator
+ if hasattr(self,'_old_formatter'):
+ self.colorbar.formatter = self._old_formatter
self.norm.autoscale(na.array([zmin,zmax]))
self.image.changed()
if self.colorbar is not None:
@@ -264,6 +302,8 @@
self.colorbar = self._figure.colorbar(self._axes.images[-1], \
extend='neither', \
shrink=0.95)
+ self._old_locator = self.colorbar.locator
+ self._old_formatter = self.colorbar.formatter
else:
self.colorbar = None
self.set_width(1,'unitary')
@@ -730,10 +770,21 @@
def set_ylim(self, ymin, ymax):
self._ylim = (ymin,ymax)
- def set_zlim(self, zmin, zmax):
+ def set_zlim(self, zmin, zmax, dex=None):
"""
Set the z boundaries of this plot.
"""
+ # This next call fixes some things, but is slower...
+ #self._redraw_image()
+ if (zmin is None) or (zmax is None):
+ if zmin == 'min':
+ zmin = na.nanmin(self._axes.images[-1]._A)
+ if dex is not None:
+ zmax = min(zmin*10**(dex),na.nanmax(self._axes.images[-1]._A))
+ if zmax == 'max':
+ zmax = na.nanmax(self._axes.images[-1]._A)
+ if dex is not None:
+ zmin = max(zmax/(10**(dex)),na.nanmin(self._axes.images[-1]._A))
self._zlim = (zmin, zmax)
def set_log_field(self, val):
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Fri Jul 17 07:30:13 2009
@@ -467,8 +467,7 @@
& na.all( RE >= self.start_point, axis=1 ) )
p = p | ( na.all( LE <= self.end_point, axis=1 )
& na.all( RE >= self.end_point, axis=1 ) )
- #self._grids = self.hierarchy.grids[p]
- self._grids = self.hierarchy.grids.copy()
+ self._grids = self.hierarchy.grids[p]
def _get_line_at_coord(self, v, index):
# t*self.vec + self.start_point = self.end_point
More information about the yt-svn
mailing list