[Yt-svn] yt-commit r1413 - trunk/yt/lagos
skillman at wrangler.dreamhost.com
skillman at wrangler.dreamhost.com
Thu Aug 27 09:19:36 PDT 2009
Author: skillman
Date: Thu Aug 27 09:19:35 2009
New Revision: 1413
URL: http://yt.spacepope.org/changeset/1413
Log:
This is a fairly large amount of code changes that mainly change bin
positions in the 1,2, and 3d profiles, and how the data is written
out. All your previous scripts should work the same other than
left_collect in BinnedProfiles has now been changed to end_collect
(see below).
Removed:
left_collect parameter - was used to accumulate all data less than the
minimum bound in the first bin. Replaced with end_collect, which does
the same thing but on both ends of the profile.
Added:
end_collect - see above.
bin_style - In write_out and write_out_h5, we now have the option to
choose which bin values are written. Options are:
'both' - stores n+1 bin edges
'left' - n bin left edges
'right' - n bin right edges
'center' - n bin centers (both in log or linear)
write_out and write_out_h5 default to 'left'. Note that if 'both' is
used in write_out, the field will always have a value of 0 in the last
row since there is one extra bin value. In write_out_h5, this is not
true since the bins and data can have different sizes.
Changed:
2D PlotCollection Profiles: If a 2d profile(phase diagram) is made in
PlotCollection, the resulting image now have bins that fill the x,y
bounds. Previously the bins were shifted by 1 row&column due to the
way digitize works. However, since pcolormesh(the generator of the
phase plot images) accepts an (n+1)x(m+1) mesh of points to be filled
with an (n)x(m) dataset, we now are able to fill the entire range
without shifting. Most people won't notice the change.
1D Profile Plots: Plots are explicitly set to draw on left edges.
This could be changed if someone has a good reason for it.
Updated docs for each function call.
Please let me know if you run into any problems. I have tested a
large number of cases, but there may be corner cases that I haven't
thought of.
Cheers,
Sam
Modified:
trunk/yt/lagos/Profiles.py
Modified: trunk/yt/lagos/Profiles.py
==============================================================================
--- trunk/yt/lagos/Profiles.py (original)
+++ trunk/yt/lagos/Profiles.py Thu Aug 27 09:19:35 2009
@@ -3,6 +3,8 @@
Author: Matthew Turk <matthewturk at gmail.com>
Affiliation: KIPAC/SLAC/Stanford
+Author: Samuel Skillman <samskillman at gmail.com>
+Affiliation: CASA, University of Colorado at Boulder
Homepage: http://yt.enzotools.org/
License:
Copyright (C) 2007-2009 Matthew Turk. All Rights Reserved.
@@ -173,32 +175,38 @@
def __init__(self, data_source, n_bins, bin_field,
lower_bound, upper_bound,
log_space = True, lazy_reader=False,
- left_collect = False):
+ end_collect=False):
"""
A 'Profile' produces either a weighted (or unweighted) average or a
straight sum of a field in a bin defined by another field. In the case
of a weighted average, we have: p_i = sum( w_i * v_i ) / sum(w_i)
- We accept a *data_source*, which will be binned into *n_bins* by the
- field *bin_field* between the *lower_bound* and the *upper_bound*.
- These bins may or may not be equally divided in *log_space*, and the
- *lazy_reader* flag controls whether we use a memory conservative
- approach.
+ We accept a *data_source*, which will be binned into *n_bins*
+ by the field *bin_field* between the *lower_bound* and the
+ *upper_bound*. These bins may or may not be equally divided
+ in *log_space*, and the *lazy_reader* flag controls whether we
+ use a memory conservative approach. If *end_collect* is True,
+ take all values outside the given bounds and store them in the
+ 0 and *n_bins*-1 values.
"""
BinnedProfile.__init__(self, data_source, lazy_reader)
self.bin_field = bin_field
self._x_log = log_space
- self.left_collect = left_collect
+ self.end_collect = end_collect
+ self.n_bins = n_bins
+
# Get our bins
if log_space:
func = na.logspace
lower_bound, upper_bound = na.log10(lower_bound), na.log10(upper_bound)
else:
func = na.linspace
+
# These are the bin *edges*
self._bins = func(lower_bound, upper_bound, n_bins + 1)
- # These are the bin *left edges*
+ # These are the bin *left edges*. These are the x-axis values
+ # we plot in the PlotCollection
self[bin_field] = self._bins[:-1]
# If we are not being memory-conservative, grab all the bins
@@ -243,8 +251,8 @@
if source_data.size == 0: # Nothing for us here.
raise EmptyProfileData()
# Truncate at boundaries.
- if self.left_collect:
- mi = na.where(source_data < self._bins.max())
+ if self.end_collect:
+ mi = na.arange(source_data.size)
else:
mi = na.where( (source_data > self._bins.min())
& (source_data < self._bins.max()))
@@ -252,8 +260,12 @@
if sd.size == 0:
raise EmptyProfileData()
# Stick the bins into our fixed bins, set at initialization
- bin_indices = na.digitize(sd, self._bins) - 1
- if self.left_collect: bin_indices = na.maximum(0, bin_indices)
+ bin_indices = na.digitize(sd, self._bins)
+ if self.end_collect: #limit the range of values to 0 and n_bins-1
+ bin_indices = na.minimum(na.maximum(1, bin_indices), self.n_bins) - 1
+ else: #throw away outside values
+ bin_indices -= 1
+
# Now we set up our inverse bin indices
inv_bin_indices = {}
for bin in range(self[self.bin_field].size):
@@ -261,25 +273,55 @@
inv_bin_indices[bin] = na.where(bin_indices == bin)
return (mi, inv_bin_indices)
- def write_out(self, filename, format="%0.16e"):
+ def choose_bins(self, bin_style):
+ # Depending on the bin_style, choose from bin edges 0...N either:
+ # both: 0...N, left: 0...N-1, right: 1...N
+ # center: N bins that are the average (both in linear or log
+ # space) of each pair of left/right edges
+ x = self._data[self.bin_field]
+ if bin_style is 'both': pass
+ elif bin_style is 'left': x = x[:-1]
+ elif bin_style is 'right': x = x[1:]
+ elif bin_style is 'center':
+ if self._x_log: x=na.log10(x)
+ x = 0.5*(x[:-1] + x[1:])
+ if self._x_log: x=10**x
+ else:
+ mylog.error('Did not recognize bin_style')
+ raise ValueError
+ return x
+
+ def write_out(self, filename, format="%0.16e", bin_style='left'):
+ '''
+ Write out data in ascii file, using *format* and
+ *bin_style* (left, right, center, both).
+ '''
fid = open(filename,"w")
fields = [field for field in sorted(self._data.keys()) if field != "UsedBins"]
fid.write("\t".join(["#"] + fields + ["\n"]))
- field_data = na.array([self._data[field] for field in fields])
+ fields.remove(self.bin_field)
+
+ field_data = na.array(self.choose_bins(bin_style))
+ if bin_style is 'both':
+ field_data = na.append([field_data], na.array([self._data[field] for field in fields]), axis=0)
+ else:
+ field_data = na.append([field_data], na.array([self._data[field][:-1] for field in fields]), axis=0)
+
for line in range(field_data.shape[1]):
field_data[:,line].tofile(fid, sep="\t", format=format)
fid.write("\n")
fid.close()
- def write_out_h5(self, filename, group_prefix=None):
+ def write_out_h5(self, filename, group_prefix=None, bin_style='left'):
"""
- Write out data in an hdf5 file. Each profile is put into a
- group, named by the axis fields. Optionally a group_prefix
- can be prepended to the group name. If the group already
- exists, it will delete and replace. However, due to hdf5
- functionality, in only unlinks the data, so an h5repack may be
- necessary to conserve space. Axes values are saved in group
- attributes.
+ Write out data in an hdf5 file *filename*. Each profile is
+ put into a group, named by the axis fields. Optionally a
+ *group_prefix* can be prepended to the group name. If the
+ group already exists, it will delete and replace. However,
+ due to hdf5 functionality, in only unlinks the data, so an
+ h5repack may be necessary to conserve space. Axes values are
+ saved in group attributes. Bins will be saved based on
+ *bin_style* (left, right, center, both).
"""
fid = h5py.File(filename)
fields = [field for field in sorted(self._data.keys()) if (field != "UsedBins" and field != self.bin_field)]
@@ -292,9 +334,9 @@
mylog.info("Profile file is getting larger since you are attempting to overwrite a profile. You may want to repack")
del fid[name]
group = fid.create_group(name)
- group.attrs["x-axis-%s" % self.bin_field] = self._data[self.bin_field]
+ group.attrs["x-axis-%s" % self.bin_field] = self.choose_bins(bin_style)
for field in fields:
- dset = group.create_dataset("%s" % field, data=self._data[field])
+ dset = group.create_dataset("%s" % field, data=self._data[field][:-1])
fid.close()
def _get_bin_fields(self):
@@ -304,37 +346,42 @@
def __init__(self, data_source,
x_n_bins, x_bin_field, x_lower_bound, x_upper_bound, x_log,
y_n_bins, y_bin_field, y_lower_bound, y_upper_bound, y_log,
- lazy_reader=False, left_collect=False):
+ lazy_reader=False, end_collect=False):
"""
- A 'Profile' produces either a weighted (or unweighted) average or a
- straight sum of a field in a bin defined by two other fields. In the case
- of a weighted average, we have: p_i = sum( w_i * v_i ) / sum(w_i)
-
- We accept a *data_source*, which will be binned into *x_n_bins* by the
- field *x_bin_field* between the *x_lower_bound* and the *x_upper_bound*
- and then again binned into *y_n_bins* by the field *y_bin_field*
- between the *y_lower_bound* and the *y_upper_bound*. These bins may or
- may not be equally divided in log-space as specified by *x_log* and *y_log*,
- and the *lazy_reader* flag controls whether we use a memory
- conservative approach.
+ A 'Profile' produces either a weighted (or unweighted) average
+ or a straight sum of a field in a bin defined by two other
+ fields. In the case of a weighted average, we have: p_i =
+ sum( w_i * v_i ) / sum(w_i)
+
+ We accept a *data_source*, which will be binned into
+ *x_n_bins* by the field *x_bin_field* between the
+ *x_lower_bound* and the *x_upper_bound* and then again binned
+ into *y_n_bins* by the field *y_bin_field* between the
+ *y_lower_bound* and the *y_upper_bound*. These bins may or
+ may not be equally divided in log-space as specified by
+ *x_log* and *y_log*, and the *lazy_reader* flag controls
+ whether we use a memory conservative approach. If
+ *end_collect* is True, take all values outside the given
+ bounds and store them in the 0 and *n_bins*-1 values.
"""
BinnedProfile.__init__(self, data_source, lazy_reader)
self.x_bin_field = x_bin_field
self.y_bin_field = y_bin_field
self._x_log = x_log
self._y_log = y_log
- self.left_collect = left_collect
+ self.end_collect = end_collect
+ self.x_n_bins = x_n_bins
+ self.y_n_bins = y_n_bins
func = {True:na.logspace, False:na.linspace}[x_log]
- bounds = fix_bounds(x_upper_bound, x_lower_bound, x_log)
+ bounds = fix_bounds(x_lower_bound, x_upper_bound, x_log)
self._x_bins = func(bounds[0], bounds[1], x_n_bins + 1)
- self[x_bin_field] = self._x_bins[:-1]
+ self[x_bin_field] = self._x_bins
func = {True:na.logspace, False:na.linspace}[y_log]
- bounds = fix_bounds(y_upper_bound, y_lower_bound, y_log)
+ bounds = fix_bounds(y_lower_bound, y_upper_bound, y_log)
self._y_bins = func(bounds[0], bounds[1], y_n_bins + 1)
- self[y_bin_field] = self._y_bins[:-1]
-
+ self[y_bin_field] = self._y_bins
if na.any(na.isnan(self[x_bin_field])) \
or na.any(na.isnan(self[y_bin_field])):
@@ -382,47 +429,85 @@
source_data_y = self._get_field(source, self.y_bin_field, check_cut)
if source_data_x.size == 0:
raise EmptyProfileData()
- if self.left_collect:
- mi = na.where( (source_data_x < self._x_bins.max())
- & (source_data_y < self._y_bins.max()))
+
+ if self.end_collect:
+ mi = na.arange(source_data_x.size)
else:
mi = na.where( (source_data_x > self._x_bins.min())
- & (source_data_x < self._x_bins.max())
- & (source_data_y > self._y_bins.min())
- & (source_data_y < self._y_bins.max()))
+ & (source_data_x < self._x_bins.max())
+ & (source_data_y > self._y_bins.min())
+ & (source_data_y < self._y_bins.max()))
sd_x = source_data_x[mi]
sd_y = source_data_y[mi]
if sd_x.size == 0 or sd_y.size == 0:
raise EmptyProfileData()
+
bin_indices_x = na.digitize(sd_x, self._x_bins) - 1
bin_indices_y = na.digitize(sd_y, self._y_bins) - 1
- if self.left_collect:
- bin_indices_x = na.maximum(bin_indices_x, 0)
- bin_indices_y = na.maximum(bin_indices_y, 0)
+ if self.end_collect:
+ bin_indices_x = na.minimum(na.maximum(1, bin_indices_x), self.x_n_bins) - 1
+ bin_indices_y = na.minimum(na.maximum(1, bin_indices_y), self.y_n_bins) - 1
+
# Now we set up our inverse bin indices
return (mi, bin_indices_x, bin_indices_y)
- def write_out(self, filename, format="%0.16e"):
+ def choose_bins(self, bin_style):
+ # Depending on the bin_style, choose from bin edges 0...N either:
+ # both: 0...N, left: 0...N-1, right: 1...N
+ # center: N bins that are the average (both in linear or log
+ # space) of each pair of left/right edges
+
+ x = self._data[self.x_bin_field]
+ y = self._data[self.y_bin_field]
+ if bin_style is 'both':
+ pass
+ elif bin_style is 'left':
+ x = x[:-1]
+ y = y[:-1]
+ elif bin_style is 'right':
+ x = x[1:]
+ y = y[1:]
+ elif bin_style is 'center':
+ if self._x_log: x=na.log10(x)
+ if self._y_log: y=na.log10(y)
+ x = 0.5*(x[:-1] + x[1:])
+ y = 0.5*(y[:-1] + y[1:])
+ if self._x_log: x=10**x
+ if self._y_log: y=10**y
+ else:
+ mylog.error('Did not recognize bin_style')
+ raise ValueError
+
+ return x,y
+
+ def write_out(self, filename, format="%0.16e", bin_style='left'):
"""
- Write out the values of x,y,v in ascii to *filename* for every field in
- the profile. Optionally a *format* can be specified.
+ Write out the values of x,y,v in ascii to *filename* for every
+ field in the profile. Optionally a *format* can be specified.
+ Bins will be saved based on *bin_style* (left, right, center,
+ both).
"""
fid = open(filename,"w")
fields = [field for field in sorted(self._data.keys()) if field != "UsedBins"]
fid.write("\t".join(["#"] + [self.x_bin_field, self.y_bin_field]
+ fields + ["\n"]))
- x,y = na.meshgrid(self._data[self.y_bin_field],
- self._data[self.y_bin_field])
+ x,y = self.choose_bins(bin_style)
+ x,y = na.meshgrid(x,y)
field_data = [x.ravel(), y.ravel()]
- field_data += [self._data[field].ravel() for field in fields
- if field not in [self.x_bin_field, self.y_bin_field]]
+ if bin_style is not 'both':
+ field_data += [self._data[field][:-1,:-1].ravel() for field in fields
+ if field not in [self.x_bin_field, self.y_bin_field]]
+ else:
+ field_data += [self._data[field].ravel() for field in fields
+ if field not in [self.x_bin_field, self.y_bin_field]]
+
field_data = na.array(field_data)
for line in range(field_data.shape[1]):
field_data[:,line].tofile(fid, sep="\t", format=format)
fid.write("\n")
fid.close()
- def write_out_h5(self, filename, group_prefix=None):
+ def write_out_h5(self, filename, group_prefix=None, bin_style='left'):
"""
Write out data in an hdf5 file. Each profile is put into a
group, named by the axis fields. Optionally a group_prefix
@@ -430,7 +515,8 @@
exists, it will delete and replace. However, due to hdf5
functionality, in only unlinks the data, so an h5repack may be
necessary to conserve space. Axes values are saved in group
- attributes.
+ attributes. Bins will be saved based on *bin_style* (left,
+ right, center, both).
"""
fid = h5py.File(filename)
fields = [field for field in sorted(self._data.keys()) if (field != "UsedBins" and field != self.x_bin_field and field != self.y_bin_field)]
@@ -442,10 +528,12 @@
mylog.info("Profile file is getting larger since you are attempting to overwrite a profile. You may want to repack")
del fid[name]
group = fid.create_group(name)
- group.attrs["x-axis-%s" % self.x_bin_field] = self._data[self.x_bin_field]
- group.attrs["y-axis-%s" % self.y_bin_field] = self._data[self.y_bin_field]
+
+ xbins, ybins = self.choose_bins(bin_style)
+ group.attrs["x-axis-%s" % self.x_bin_field] = xbins
+ group.attrs["y-axis-%s" % self.y_bin_field] = ybins
for field in fields:
- dset = group.create_dataset("%s" % field, data=self._data[field])
+ dset = group.create_dataset("%s" % field, data=self._data[field][:-1,:-1])
fid.close()
def _get_bin_fields(self):
@@ -456,11 +544,26 @@
return upper, lower
class BinnedProfile3D(BinnedProfile):
+ """
+ A 'Profile' produces either a weighted (or unweighted) average
+ or a straight sum of a field in a bin defined by two other
+ fields. In the case of a weighted average, we have: p_i =
+ sum( w_i * v_i ) / sum(w_i)
+
+ We accept a *data_source*, which will be binned into
+ *(x,y,z)_n_bins* by the field *(x,y,z)_bin_field* between the
+ *(x,y,z)_lower_bound* and the *(x,y,z)_upper_bound*. These bins may or
+ may not be equally divided in log-space as specified by
+ *(x,y,z)_log*, and the *lazy_reader* flag controls
+ whether we use a memory conservative approach. If
+ *end_collect* is True, take all values outside the given
+ bounds and store them in the 0 and *n_bins*-1 values.
+ """
def __init__(self, data_source,
x_n_bins, x_bin_field, x_lower_bound, x_upper_bound, x_log,
y_n_bins, y_bin_field, y_lower_bound, y_upper_bound, y_log,
z_n_bins, z_bin_field, z_lower_bound, z_upper_bound, z_log,
- lazy_reader=False):
+ lazy_reader=False, end_collect=False):
BinnedProfile.__init__(self, data_source, lazy_reader)
self.x_bin_field = x_bin_field
self.y_bin_field = y_bin_field
@@ -468,21 +571,25 @@
self._x_log = x_log
self._y_log = y_log
self._z_log = z_log
+ self.end_collect = end_collect
+ self.x_n_bins = x_n_bins
+ self.y_n_bins = y_n_bins
+ self.z_n_bins = z_n_bins
func = {True:na.logspace, False:na.linspace}[x_log]
- bounds = fix_bounds(x_upper_bound, x_lower_bound, x_log)
+ bounds = fix_bounds(x_lower_bound, x_upper_bound, x_log)
self._x_bins = func(bounds[0], bounds[1], x_n_bins + 1)
- self[x_bin_field] = self._x_bins[:-1]
+ self[x_bin_field] = self._x_bins
func = {True:na.logspace, False:na.linspace}[y_log]
- bounds = fix_bounds(y_upper_bound, y_lower_bound, y_log)
+ bounds = fix_bounds(y_lower_bound, y_upper_bound, y_log)
self._y_bins = func(bounds[0], bounds[1], y_n_bins + 1)
- self[y_bin_field] = self._y_bins[:-1]
+ self[y_bin_field] = self._y_bins
func = {True:na.logspace, False:na.linspace}[z_log]
- bounds = fix_bounds(z_upper_bound, z_lower_bound, z_log)
+ bounds = fix_bounds(z_lower_bound, z_upper_bound, z_log)
self._z_bins = func(bounds[0], bounds[1], z_n_bins + 1)
- self[z_bin_field] = self._z_bins[:-1]
+ self[z_bin_field] = self._z_bins
if na.any(na.isnan(self[x_bin_field])) \
or na.any(na.isnan(self[y_bin_field])) \
@@ -537,27 +644,71 @@
source_data_z = self._get_field(source, self.z_bin_field, check_cut)
if source_data_x.size == 0:
raise EmptyProfileData()
- mi = ( (source_data_x > self._x_bins.min())
- & (source_data_x < self._x_bins.max())
- & (source_data_y > self._y_bins.min())
- & (source_data_y < self._y_bins.max())
- & (source_data_z > self._z_bins.min())
- & (source_data_z < self._z_bins.max()))
+ if self.end_collect:
+ mi = na.arange(source_data_x.size)
+ else:
+ mi = ( (source_data_x > self._x_bins.min())
+ & (source_data_x < self._x_bins.max())
+ & (source_data_y > self._y_bins.min())
+ & (source_data_y < self._y_bins.max())
+ & (source_data_z > self._z_bins.min())
+ & (source_data_z < self._z_bins.max()))
sd_x = source_data_x[mi]
sd_y = source_data_y[mi]
sd_z = source_data_z[mi]
if sd_x.size == 0 or sd_y.size == 0 or sd_z.size == 0:
raise EmptyProfileData()
- bin_indices_x = na.digitize(sd_x, self._x_bins)
- bin_indices_y = na.digitize(sd_y, self._y_bins)
- bin_indices_z = na.digitize(sd_z, self._z_bins)
+
+ bin_indices_x = na.digitize(sd_x, self._x_bins) - 1
+ bin_indices_y = na.digitize(sd_y, self._y_bins) - 1
+ bin_indices_z = na.digitize(sd_z, self._z_bins) - 1
+ if self.end_collect:
+ bin_indices_x = na.minimum(na.maximum(1, bin_indices_x), self.x_n_bins) - 1
+ bin_indices_y = na.minimum(na.maximum(1, bin_indices_y), self.y_n_bins) - 1
+ bin_indices_z = na.minimum(na.maximum(1, bin_indices_z), self.z_n_bins) - 1
+
# Now we set up our inverse bin indices
return (mi, bin_indices_x, bin_indices_y, bin_indices_z)
+ def choose_bins(self, bin_style):
+ # Depending on the bin_style, choose from bin edges 0...N either:
+ # both: 0...N, left: 0...N-1, right: 1...N
+ # center: N bins that are the average (both in linear or log
+ # space) of each pair of left/right edges
+
+ x = self._data[self.x_bin_field]
+ y = self._data[self.y_bin_field]
+ z = self._data[self.z_bin_field]
+ if bin_style is 'both':
+ pass
+ elif bin_style is 'left':
+ x = x[:-1]
+ y = y[:-1]
+ z = z[:-1]
+ elif bin_style is 'right':
+ x = x[1:]
+ y = y[1:]
+ z = z[1:]
+ elif bin_style is 'center':
+ if self._x_log: x=na.log10(x)
+ if self._y_log: y=na.log10(y)
+ if self._z_log: z=na.log10(z)
+ x = 0.5*(x[:-1] + x[1:])
+ y = 0.5*(y[:-1] + y[1:])
+ z = 0.5*(z[:-1] + z[1:])
+ if self._x_log: x=10**x
+ if self._y_log: y=10**y
+ if self._z_log: y=10**z
+ else:
+ mylog.error('Did not recognize bin_style')
+ raise ValueError
+
+ return x,y,z
+
def write_out(self, filename, format="%0.16e"):
pass # Will eventually dump HDF5
- def write_out_h5(self, filename, group_prefix=None):
+ def write_out_h5(self, filename, group_prefix=None, bin_style='left'):
"""
Write out data in an hdf5 file. Each profile is put into a
group, named by the axis fields. Optionally a group_prefix
@@ -579,12 +730,14 @@
mylog.info("Profile file is getting larger since you are attempting to overwrite a profile. You may want to repack")
del fid[name]
group = fid.create_group(name)
- group.attrs["x-axis-%s" % self.x_bin_field] = self._data[self.x_bin_field]
- group.attrs["y-axis-%s" % self.y_bin_field] = self._data[self.y_bin_field]
- group.attrs["z-axis-%s" % self.z_bin_field] = self._data[self.z_bin_field]
+
+ xbins, ybins, zbins= self.choose_bins(bin_style)
+ group.attrs["x-axis-%s" % self.x_bin_field] = xbins
+ group.attrs["y-axis-%s" % self.y_bin_field] = ybins
+ group.attrs["z-axis-%s" % self.z_bin_field] = zbins
for field in fields:
- dset = group.create_dataset("%s" % field, data=self._data[field])
+ dset = group.create_dataset("%s" % field, data=self._data[field][:-1,:-1,:-1])
fid.close()
More information about the yt-svn
mailing list