[Yt-svn] yt-commit r938 - in branches/yt-non-3d: doc examples yt yt/extensions/lightcone yt/lagos yt/lagos/hop yt/raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Sat Nov 15 15:40:08 PST 2008
Author: mturk
Date: Sat Nov 15 15:40:05 2008
New Revision: 938
URL: http://yt.spacepope.org/changeset/938
Log:
Merged back changes from trunk, to continue with the 2D and 1D data. (912:HEAD)
Added:
branches/yt-non-3d/yt/extensions/lightcone/Common_nVolume.py
- copied unchanged from r937, /trunk/yt/extensions/lightcone/Common_nVolume.py
Removed:
branches/yt-non-3d/yt/extensions/lightcone/sample_lightcone.par
Modified:
branches/yt-non-3d/doc/install_script.sh
branches/yt-non-3d/examples/makeLightCone.py
branches/yt-non-3d/yt/config.py
branches/yt-non-3d/yt/extensions/lightcone/LightCone.py
branches/yt-non-3d/yt/extensions/lightcone/LightConeProjection.py
branches/yt-non-3d/yt/lagos/BaseDataTypes.py
branches/yt-non-3d/yt/lagos/BaseGridType.py
branches/yt-non-3d/yt/lagos/ContourFinder.py
branches/yt-non-3d/yt/lagos/DataReadingFuncs.py
branches/yt-non-3d/yt/lagos/DerivedFields.py
branches/yt-non-3d/yt/lagos/HelperFunctions.py
branches/yt-non-3d/yt/lagos/HierarchyType.py
branches/yt-non-3d/yt/lagos/ParallelTools.py
branches/yt-non-3d/yt/lagos/PointCombine.c
branches/yt-non-3d/yt/lagos/hop/NewHopOutput.py
branches/yt-non-3d/yt/raven/Callbacks.py
branches/yt-non-3d/yt/raven/PlotTypes.py
branches/yt-non-3d/yt/recipes.py
Modified: branches/yt-non-3d/doc/install_script.sh
==============================================================================
--- branches/yt-non-3d/doc/install_script.sh (original)
+++ branches/yt-non-3d/doc/install_script.sh Sat Nov 15 15:40:05 2008
@@ -53,7 +53,7 @@
# Individual processes
if [ -z "$HDF5_DIR" ]
then
- [ ! -e hdf5-1.6.7.tar.gz ] && wget ftp://ftp.hdfgroup.org/HDF5/current16/src/hdf5-1.6.7.tar.gz
+ [ ! -e hdf5-1.6.8.tar.gz ] && wget ftp://ftp.hdfgroup.org/HDF5/current16/src/hdf5-1.6.8.tar.gz
fi
[ $INST_ZLIB -eq 1 ] && [ ! -e zlib-1.2.3.tar.bz2 ] && wget http://www.zlib.net/zlib-1.2.3.tar.bz2
@@ -98,11 +98,11 @@
if [ -z "$HDF5_DIR" ]
then
- if [ ! -e hdf5-1.6.7/done ]
+ if [ ! -e hdf5-1.6.8/done ]
then
- [ ! -e hdf5-1.6.7 ] && tar xvfz hdf5-1.6.7.tar.gz
+ [ ! -e hdf5-1.6.8 ] && tar xvfz hdf5-1.6.8.tar.gz
echo "Doing HDF5"
- cd hdf5-1.6.7
+ cd hdf5-1.6.8
./configure --prefix=${DEST_DIR}/ || exit 1
make install || exit 1
touch done
Modified: branches/yt-non-3d/examples/makeLightCone.py
==============================================================================
--- branches/yt-non-3d/examples/makeLightCone.py (original)
+++ branches/yt-non-3d/examples/makeLightCone.py Sat Nov 15 15:40:05 2008
@@ -8,6 +8,10 @@
q.CalculateLightConeSolution()
+# If random seed was not provided in the parameter file, it can be given
+# straight to the routine.
+q.CalculateLightConeSolution(seed=123456789)
+
# Save a text file detailing the light cone solution.
q.SaveLightConeSolution()
@@ -30,7 +34,7 @@
# This will allow the projection objects that have already been made
# to be re-used.
# Just don't use the same random seed as the original.
-q.RecycleLightConeSolution(987654321)
+q.RerandomizeLightConeSolution(987654321,recycle=True)
# Save the recycled solution.
q.SaveLightConeSolution(file='light_cone_recycled.out')
@@ -40,3 +44,6 @@
# Save the stack.
q.SaveLightConeStack(file='light_cone_recycled.h5')
+
+# Rerandomize the light cone solution with an entirely new solution.
+q.RerandomizeLightConeSolution(8675309,recycle=False)
Modified: branches/yt-non-3d/yt/config.py
==============================================================================
--- branches/yt-non-3d/yt/config.py (original)
+++ branches/yt-non-3d/yt/config.py Sat Nov 15 15:40:05 2008
@@ -27,6 +27,7 @@
import ConfigParser, os, os.path, types
+
ytcfgDefaults = {
"fido":{
'RunDir': os.path.join(os.getenv("HOME"),'.yt/EnzoRuns/'),
@@ -107,8 +108,12 @@
raise KeyError
self.set(item[0], item[1], val)
-ytcfg = YTConfigParser(['yt.cfg', os.path.expanduser('~/.yt/config')],
- ytcfgDefaults)
+if os.path.exists(os.path.expanduser("~/.yt/config")):
+ ytcfg = YTConfigParser(['yt.cfg', os.path.expanduser('~/.yt/config')],
+ ytcfgDefaults)
+else:
+ ytcfg = YTConfigParser(['yt.cfg'],
+ ytcfgDefaults)
# Now we have parsed the config file. Overrides come from the command line.
Modified: branches/yt-non-3d/yt/extensions/lightcone/LightCone.py
==============================================================================
--- branches/yt-non-3d/yt/extensions/lightcone/LightCone.py (original)
+++ branches/yt-non-3d/yt/extensions/lightcone/LightCone.py Sat Nov 15 15:40:05 2008
@@ -28,6 +28,7 @@
import numpy as na
import random as rand
import tables as h5
+from Common_nVolume import *
class LightCone(object):
def __init__(self,EnzoParameterFile,LightConeParameterFile):
@@ -63,12 +64,15 @@
# Calculate maximum delta z for each data dump.
self._CalculateDeltaZMax()
- def CalculateLightConeSolution(self):
+ def CalculateLightConeSolution(self,seed=None):
"Create list of projections to be added together to make the light cone."
# Make sure recycling flag is off.
self.recycleSolution = False
+ if seed is not None:
+ self.lightConeParameters['RandomSeed'] = seed
+
# Make light cone using minimum number of projections.
if (self.lightConeParameters['UseMinimumNumberOfProjections']):
@@ -215,7 +219,7 @@
# Return the plot collection so the user can remake the plot if they want.
return pc
- def RecycleLightConeSolution(self,newSeed):
+ def RerandomizeLightConeSolution(self,newSeed,recycle=True):
"""
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.
@@ -224,32 +228,69 @@
This routine will take in a new random seed and rerandomize the parts of the light cone that do not contribute
to creating a unique projection object. Additionally, this routine is built such that if the same random
seed is given for the rerandomizing, the solution will be identical to the original.
- """
- mylog.info("Recycling solution made with %s with new seed %s." % (str(self.lightConeParameters['RandomSeed']),
- str(newSeed)))
+ 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.
+ """
- self.recycleSolution = True
- self.recycleRandomSeed = newSeed
+ if recycle:
+ mylog.info("Recycling solution made with %s with new seed %s." % (self.lightConeParameters['RandomSeed'],
+ newSeed))
+ self.recycleRandomSeed = newSeed
+ else:
+ mylog.info("Creating new solution with random seed %s." % newSeed)
+ self.lightConeParameters['RandomSeed'] = newSeed
+ self.recycleRandomSeed = 0
+
+ self.recycleSolution = recycle
+
+ # Keep track of fraction of volume in common between the original and recycled solution.
+ commonVolume = 0.0
+ totalVolume = 0.0
# Seed random number generator with new seed.
- rand.seed(self.recycleRandomSeed)
+ rand.seed(newSeed)
for q,output in enumerate(self.lightConeSolution):
# 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 random projection axis and center.
- # Axis will get thrown away since it is used in creating a unique projection object.
+ # If recyclsing, axis will get thrown away since it is used in creating a unique projection object.
newAxis = rand.randint(0,2)
+ if not recycle:
+ output['ProjectionAxis'] = newAxis
newCenter = [rand.random(),rand.random(),rand.random()]
+ # Make list of rectangle corners to calculate common volume.
+ newCube = na.zeros(shape=(len(newCenter),2))
+ oldCube = na.zeros(shape=(len(newCenter),2))
+ for w in range(len(newCenter)):
+ if (w == output['ProjectionAxis']):
+ oldCube[w] = [output['ProjectionCenter'][w] - 0.5 * output['DepthBoxFraction'],
+ output['ProjectionCenter'][w] + 0.5 * output['DepthBoxFraction']]
+ # If recycling old solution, then keep center for axis in line of sight.
+ if recycle:
+ newCube[w] = oldCube[w]
+ else:
+ newCube[w] = [newCenter[w] - 0.5 * output['DepthBoxFraction'],
+ newCenter[w] + 0.5 * output['DepthBoxFraction']]
+ else:
+ oldCube[w] = [output['ProjectionCenter'][w] - 0.5 * output['WidthBoxFraction'],
+ output['ProjectionCenter'][w] + 0.5 * output['WidthBoxFraction']]
+ newCube[w] = [newCenter[w] - 0.5 * output['WidthBoxFraction'],
+ newCenter[w] + 0.5 * output['WidthBoxFraction']]
+
+ commonVolume += commonNVolume(oldCube,newCube,periodic=na.array([[0,1],[0,1],[0,1]]))
+ totalVolume += output['DepthBoxFraction'] * output['WidthBoxFraction']**2
+
# Replaces centers for every axis except the line of sight axis.
for w in range(len(newCenter)):
- if (w != self.lightConeSolution[q]['ProjectionAxis']):
+ if not(recycle and (w == self.lightConeSolution[q]['ProjectionAxis'])):
self.lightConeSolution[q]['ProjectionCenter'][w] = newCenter[w]
+ mylog.info("Fraction of total volume in common with old solution is %.2e." % (commonVolume/totalVolume))
def SaveLightConeSolution(self,file="light_cone.dat"):
"Write out a text file with information on light cone solution."
Modified: branches/yt-non-3d/yt/extensions/lightcone/LightConeProjection.py
==============================================================================
--- branches/yt-non-3d/yt/extensions/lightcone/LightConeProjection.py (original)
+++ branches/yt-non-3d/yt/extensions/lightcone/LightConeProjection.py Sat Nov 15 15:40:05 2008
@@ -23,7 +23,7 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
-from yt.lagos.lightcone import *
+from yt.extensions.lightcone import *
def LightConeProjection(lightConeSlice,field,pixels,weight_field=None,save_image=False,name="",**kwargs):
"Create a single projection to be added into the light cone stack."
Modified: branches/yt-non-3d/yt/lagos/BaseDataTypes.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/BaseDataTypes.py (original)
+++ branches/yt-non-3d/yt/lagos/BaseDataTypes.py Sat Nov 15 15:40:05 2008
@@ -72,14 +72,17 @@
class FakeGridForParticles(object):
"""
Mock up a grid to insert particle positions and radii
- into for purposes of confinement in an :class:`Enzo3DData`.
+ into for purposes of confinement in an :class:`AMR3DData`.
"""
def __init__(self, grid):
self._corners = grid._corners
self.field_parameters = {}
self.data = {'x':grid['particle_position_x'],
'y':grid['particle_position_y'],
- 'z':grid['particle_position_z']}
+ 'z':grid['particle_position_z'],
+ 'dx':grid['dx'],
+ 'dy':grid['dy'],
+ 'dz':grid['dz']}
self.real_grid = grid
self.child_mask = 1
def __getitem__(self, field):
@@ -94,9 +97,9 @@
else: tr = self.data[field]
return tr
-class EnzoData:
+class AMRData:
"""
- Generic EnzoData container. By itself, will attempt to
+ Generic AMRData container. By itself, will attempt to
generate field, read fields (method defined by derived classes)
and deal with passing back and forth field parameters.
"""
@@ -106,7 +109,7 @@
def __init__(self, pf, fields, **kwargs):
"""
Typically this is never called directly, but only due to inheritance.
- It associates a :class:`~yt.lagos.EnzoStaticOutput` with the class,
+ It associates a :class:`~yt.lagos.StaticOutput` with the class,
sets its initial set of fields, and the remainder of the arguments
are passed as field_parameters.
"""
@@ -159,7 +162,7 @@
def clear_data(self):
"""
- Clears out all data from the EnzoData instance, freeing memory.
+ Clears out all data from the AMRData instance, freeing memory.
"""
for key in self.data.keys():
del self.data[key]
@@ -304,10 +307,10 @@
__del_gridLevels)
-class Enzo1DData(EnzoData, GridPropertiesMixin):
+class AMR1DData(AMRData, GridPropertiesMixin):
_spatial = False
def __init__(self, pf, fields, **kwargs):
- EnzoData.__init__(self, pf, fields, **kwargs)
+ AMRData.__init__(self, pf, fields, **kwargs)
self._grids = None
def _generate_field_in_grids(self, field, num_ghost_zones=0):
@@ -351,14 +354,14 @@
[self._get_data_from_grid(grid, field)
for grid in self._grids])
-class EnzoOrthoRayBase(Enzo1DData):
+class AMROrthoRayBase(AMR1DData):
_key_fields = ['x','y','z','dx','dy','dz']
def __init__(self, axis, coords, fields=None, pf=None, **kwargs):
"""
Dimensionality is reduced to one, and an ordered list of points at an
(x,y) tuple along *axis* are available.
"""
- Enzo1DData.__init__(self, pf, fields, **kwargs)
+ AMR1DData.__init__(self, pf, fields, **kwargs)
self.axis = axis
self.px_ax = x_dict[self.axis]
self.py_ax = y_dict[self.axis]
@@ -396,14 +399,14 @@
gf = grid[field][sl]
return gf[na.where(grid.child_mask[sl])]
-class EnzoRayBase(Enzo1DData):
+class AMRRayBase(AMR1DData):
def __init__(self, start_point, end_point, fields=None, pf=None, **kwargs):
"""
We accept a start point and an end point and then get all the data
between those two.
"""
mylog.warning("This code is poorly tested. It may give bad data!")
- Enzo1DData.__init__(self, pf, fields, **kwargs)
+ AMR1DData.__init__(self, pf, fields, **kwargs)
self.start_point = na.array(start_point)
self.end_point = na.array(end_point)
self.vec = self.end_point - self.start_point
@@ -452,20 +455,20 @@
self.center, self.vec)
return mask
-class Enzo2DData(EnzoData, GridPropertiesMixin):
+class AMR2DData(AMRData, GridPropertiesMixin):
_key_fields = ['px','py','pdx','pdy']
"""
- Class to represent a set of :class:`EnzoData` that's 2-D in nature, and
+ Class to represent a set of :class:`AMRData` that's 2-D in nature, and
thus does not have as many actions as the 3-D data types.
"""
_spatial = False
def __init__(self, axis, fields, pf=None, **kwargs):
"""
- Prepares the Enzo2DData, normal to *axis*. If *axis* is 4, we are not
+ Prepares the AMR2DData, normal to *axis*. If *axis* is 4, we are not
aligned with any axis.
"""
self.axis = axis
- EnzoData.__init__(self, pf, fields, **kwargs)
+ AMRData.__init__(self, pf, fields, **kwargs)
self.set_field_parameter("axis",axis)
@@ -566,9 +569,9 @@
self[f] = array[i,:]
return True
-class EnzoSliceBase(Enzo2DData):
+class AMRSliceBase(AMR2DData):
"""
- EnzoSlice is an orthogonal slice through the data, taking all the points
+ AMRSlice is an orthogonal slice through the data, taking all the points
at the finest resolution available and then indexing them. It is more
appropriately thought of as a slice 'operator' than an object,
however, as its field and coordinate can both change.
@@ -582,7 +585,7 @@
Slice along *axis*:ref:`axis-specification`, at the coordinate *coord*.
Optionally supply fields.
"""
- Enzo2DData.__init__(self, axis, fields, pf, **kwargs)
+ AMR2DData.__init__(self, axis, fields, pf, **kwargs)
self.center = center
self.coord = coord
if node_name is False:
@@ -693,16 +696,16 @@
def _gen_node_name(self):
return "%s_%s_%s" % (self.fields[0], self.axis, self.coord)
-class EnzoCuttingPlaneBase(Enzo2DData):
+class AMRCuttingPlaneBase(AMR2DData):
"""
- EnzoCuttingPlane is an oblique plane through the data,
+ AMRCuttingPlane is an oblique plane through the data,
defined by a normal vector and a coordinate. It attempts to guess
an 'up' vector, which cannot be overridden, and then it pixelizes
the appropriate data onto the plane without interpolation.
"""
_plane = None
_top_node = "/CuttingPlanes"
- _key_fields = Enzo2DData._key_fields + ['pz','pdz']
+ _key_fields = AMR2DData._key_fields + ['pz','pdz']
def __init__(self, normal, center, fields = None, node_name = None,
**kwargs):
"""
@@ -710,7 +713,7 @@
the *normal* vector and the *center* to define the viewing plane.
The 'up' direction is guessed at automatically.
"""
- Enzo2DData.__init__(self, 4, fields, **kwargs)
+ AMR2DData.__init__(self, 4, fields, **kwargs)
self.center = center
self.set_field_parameter('center',center)
# Let's set up our plane equation
@@ -821,19 +824,19 @@
L_name = ("%s" % self._norm_vec).replace(" ","_")[1:-1]
return "%s_c%s_L%s" % (self.fields[0], cen_name, L_name)
-class EnzoProjBase(Enzo2DData, ParallelAnalysisInterface):
+class AMRProjBase(AMR2DData, ParallelAnalysisInterface):
_top_node = "/Projections"
- _key_fields = Enzo2DData._key_fields + ['weight_field']
+ _key_fields = AMR2DData._key_fields + ['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):
"""
- EnzoProj is a projection of a *field* along an *axis*. The field
+ AMRProj is a projection of a *field* along an *axis*. The field
can have an associated *weight_field*, in which case the values are
multiplied by a weight before being summed, and then divided by the sum
of that weight.
"""
- Enzo2DData.__init__(self, axis, field, pf, node_name = None, **kwargs)
+ AMR2DData.__init__(self, axis, field, pf, node_name = None, **kwargs)
if field_cuts is not None:
field_cuts = ['grid%s' % cut for cut in ensure_list(field_cuts)]
self._field_cuts = field_cuts
@@ -989,8 +992,9 @@
args += self.__retval_coords[grid2.id] + [self.__retval_fields[grid2.id]]
args += self.__retval_coords[grid1.id] + [self.__retval_fields[grid1.id]]
args.append(1) # Refinement factor
+ args.append(na.ones(args[0].shape, dtype='int64'))
kk = PointCombine.CombineGrids(*args)
- goodI = na.where(self.__retval_coords[grid2.id][0] > -1)
+ goodI = args[-1].astype('bool')
self.__retval_coords[grid2.id] = \
[coords[goodI] for coords in self.__retval_coords[grid2.id]]
self.__retval_fields[grid2.id] = \
@@ -1004,18 +1008,20 @@
% (level, self._max_level), len(grids))
for pi, grid1 in enumerate(grids):
pbar.update(pi)
- for grid2 in self.source._grids[grids_up][self.__overlap_masks[grid1.Parent.id]]:
- if self.__retval_coords[grid2.id][0].shape[0] == 0: continue
- args = []
- args += self.__retval_coords[grid2.id] + [self.__retval_fields[grid2.id]]
- args += self.__retval_coords[grid1.id] + [self.__retval_fields[grid1.id]]
- args.append(int(grid2.dx / grid1.dx))
- kk = PointCombine.CombineGrids(*args)
- goodI = (self.__retval_coords[grid2.id][0] > -1)
- self.__retval_coords[grid2.id] = \
- [coords[goodI] for coords in self.__retval_coords[grid2.id]]
- self.__retval_fields[grid2.id] = \
- [fields[goodI] for fields in self.__retval_fields[grid2.id]]
+ for parent in ensure_list(grid1.Parent):
+ for grid2 in self.source._grids[grids_up][self.__overlap_masks[parent.id]]:
+ if self.__retval_coords[grid2.id][0].shape[0] == 0: continue
+ args = []
+ args += self.__retval_coords[grid2.id] + [self.__retval_fields[grid2.id]]
+ args += self.__retval_coords[grid1.id] + [self.__retval_fields[grid1.id]]
+ args.append(int(grid2.dx / grid1.dx))
+ args.append(na.ones(args[0].shape, dtype='int64'))
+ kk = PointCombine.CombineGrids(*args)
+ goodI = args[-1].astype('bool')
+ self.__retval_coords[grid2.id] = \
+ [coords[goodI] for coords in self.__retval_coords[grid2.id]]
+ self.__retval_fields[grid2.id] = \
+ [fields[goodI] for fields in self.__retval_fields[grid2.id]]
for grid1 in self.source.select_grids(level-1):
if not self._check_region and self.__retval_coords[grid1.id][0].size != 0:
mylog.error("Something messed up, and %s still has %s points of data",
@@ -1131,7 +1137,7 @@
def _gen_node_name(self):
return "%s_%s_%s" % (self.fields[0], self._weight, self.axis)
-class Enzo3DData(EnzoData, GridPropertiesMixin):
+class AMR3DData(AMRData, GridPropertiesMixin):
_key_fields = ['x','y','z','dx','dy','dz']
"""
Class describing a cluster of data points, not necessarily sharing any
@@ -1141,11 +1147,11 @@
_num_ghost_zones = 0
def __init__(self, center, fields, pf = None, **kwargs):
"""
- Returns an instance of Enzo3DData, or prepares one. Usually only
+ Returns an instance of AMR3DData, or prepares one. Usually only
used as a base class. Note that *center* is supplied, but only used
for fields and quantities that require it.
"""
- EnzoData.__init__(self, pf, fields, **kwargs)
+ AMRData.__init__(self, pf, fields, **kwargs)
self.center = center
self.set_field_parameter("center",center)
self.coords = None
@@ -1332,14 +1338,14 @@
grid[field][self._get_point_indices()] = value
-class ExtractedRegionBase(Enzo3DData):
+class ExtractedRegionBase(AMR3DData):
"""
ExtractedRegions are arbitrarily defined containers of data, useful
for things like selection along a baryon field.
"""
def __init__(self, base_region, indices, **kwargs):
cen = base_region.get_field_parameter("center")
- Enzo3DData.__init__(self, center=cen,
+ AMR3DData.__init__(self, center=cen,
fields=None, pf=base_region.pf, **kwargs)
self._base_region = base_region # We don't weakly reference because
# It is not cyclic
@@ -1386,7 +1392,7 @@
# Yeah, if it's not true, we don't care.
return self._indices[grid.id-1]
-class EnzoCylinderBase(Enzo3DData):
+class AMRCylinderBase(AMR3DData):
"""
We can define a cylinder (or disk) to act as a data object.
"""
@@ -1397,7 +1403,7 @@
can define a cylinder of any proportion. Only cells whose centers are
within the cylinder will be selected.
"""
- Enzo3DData.__init__(self, na.array(center), fields, pf, **kwargs)
+ AMR3DData.__init__(self, na.array(center), fields, pf, **kwargs)
self._norm_vec = na.array(normal)/na.sqrt(na.dot(normal,normal))
self.set_field_parameter("height_vector", self._norm_vec)
self._height = height
@@ -1447,9 +1453,9 @@
& (r < self._radius))
return cm
-class EnzoRegionBase(Enzo3DData):
+class AMRRegionBase(AMR3DData):
"""
- EnzoRegions are rectangular prisms of data.
+ AMRRegions are rectangular prisms of data.
"""
def __init__(self, center, left_edge, right_edge, fields = None,
pf = None, **kwargs):
@@ -1458,7 +1464,7 @@
three *right_edge* coordinates, and a *center* that need not be the
center.
"""
- Enzo3DData.__init__(self, center, fields, pf, **kwargs)
+ AMR3DData.__init__(self, center, fields, pf, **kwargs)
self.left_edge = left_edge
self.right_edge = right_edge
self._refresh_data()
@@ -1484,9 +1490,9 @@
& (grid['z'] + 0.5*grid['dz'] >= self.left_edge[2]) )
return cm
-class EnzoPeriodicRegionBase(Enzo3DData):
+class AMRPeriodicRegionBase(AMR3DData):
"""
- EnzoRegions are rectangular prisms of data.
+ AMRRegions are rectangular prisms of data.
"""
def __init__(self, center, left_edge, right_edge, fields = None,
pf = None, **kwargs):
@@ -1495,7 +1501,7 @@
three *right_edge* coordinates, and a *center* that need not be the
center.
"""
- Enzo3DData.__init__(self, center, fields, pf, **kwargs)
+ AMR3DData.__init__(self, center, fields, pf, **kwargs)
self.left_edge = na.array(left_edge)
self.right_edge = na.array(right_edge)
self._refresh_data()
@@ -1536,7 +1542,7 @@
& (grid['z'] + grid['dz'] + off_z >= self.left_edge[2]) )
return cm
-class EnzoGridCollection(Enzo3DData):
+class AMRGridCollection(AMR3DData):
"""
An arbitrary selection of grids, within which we accept all points.
"""
@@ -1546,7 +1552,7 @@
By selecting an arbitrary *grid_list*, we can act on those grids.
Child cells are not returned.
"""
- Enzo3DData.__init__(self, center, fields, pf, **kwargs)
+ AMR3DData.__init__(self, center, fields, pf, **kwargs)
self._grids = na.array(grid_list)
self.fields = fields
self.connection_pool = True
@@ -1568,7 +1574,7 @@
pointI = na.where(k == True)
return pointI
-class EnzoSphereBase(Enzo3DData):
+class AMRSphereBase(AMR3DData):
"""
A sphere of points
"""
@@ -1577,7 +1583,7 @@
The most famous of all the data objects, we define it via a
*center* and a *radius*.
"""
- Enzo3DData.__init__(self, center, fields, pf, **kwargs)
+ AMR3DData.__init__(self, center, fields, pf, **kwargs)
if radius < self.hierarchy.get_smallest_dx():
raise SyntaxError("Your radius is smaller than your finest cell!")
self.set_field_parameter('radius',radius)
@@ -1609,7 +1615,7 @@
self._cut_masks[grid.id] = cm
return cm
-class EnzoCoveringGrid(Enzo3DData):
+class AMRCoveringGrid(AMR3DData):
"""
Covering grids represent fixed-resolution data over a given region.
In order to achieve this goal -- for instance in order to obtain ghost
@@ -1625,7 +1631,7 @@
generating fixed resolution data between *left_edge* and *right_edge*
that is *dims* (3-values) on a side.
"""
- Enzo3DData.__init__(self, center=None, fields=fields, pf=pf, **kwargs)
+ AMR3DData.__init__(self, center=None, fields=fields, pf=pf, **kwargs)
self.left_edge = na.array(left_edge)
self.right_edge = na.array(right_edge)
self.level = level
@@ -1644,6 +1650,7 @@
na.any(self.right_edge > self.pf["DomainRightEdge"]):
grids,ind = self.pf.hierarchy.get_periodic_box_grids(
self.left_edge, self.right_edge)
+ ind = slice(None)
else:
grids,ind = self.pf.hierarchy.get_box_grids(
self.left_edge, self.right_edge)
@@ -1655,7 +1662,7 @@
mylog.error("Sorry, dude, do it yourself, it's already in 3-D.")
def _refresh_data(self):
- Enzo3DData._refresh_data(self)
+ AMR3DData._refresh_data(self)
self['dx'] = self.dx * na.ones(self.ActiveDimensions, dtype='float64')
self['dy'] = self.dy * na.ones(self.ActiveDimensions, dtype='float64')
self['dz'] = self.dz * na.ones(self.ActiveDimensions, dtype='float64')
@@ -1731,24 +1738,30 @@
self.left_edge, self.right_edge, c_dx, c_fields,
ll, self.pf["DomainLeftEdge"], self.pf["DomainRightEdge"])
-class EnzoSmoothedCoveringGrid(EnzoCoveringGrid):
+class AMRSmoothedCoveringGrid(AMRCoveringGrid):
def __init__(self, *args, **kwargs):
dlog2 = na.log10(kwargs['dims'])/na.log10(2)
if not na.all(na.floor(dlog2) == na.ceil(dlog2)):
mylog.warning("Must be power of two dimensions")
#raise ValueError
kwargs['num_ghost_zones'] = 0
- EnzoCoveringGrid.__init__(self, *args, **kwargs)
- if na.any(self.left_edge == 0):
+ AMRCoveringGrid.__init__(self, *args, **kwargs)
+ if na.any(self.left_edge == self.pf["DomainLeftEdge"]):
self.left_edge += self.dx
self.ActiveDimensions -= 1
- if na.any(self.right_edge == 0):
+ if na.any(self.right_edge == self.pf["DomainRightEdge"]):
self.right_edge -= self.dx
self.ActiveDimensions -= 1
def _get_list_of_grids(self):
- grids, ind = self.pf.hierarchy.get_box_grids(self.left_edge-self.dx,
- self.right_edge+self.dx)
+ if na.any(self.left_edge - self.dx < self.pf["DomainLeftEdge"]) or \
+ na.any(self.right_edge + self.dx > self.pf["DomainRightEdge"]):
+ grids,ind = self.pf.hierarchy.get_periodic_box_grids(
+ self.left_edge - self.dx, self.right_edge + self.dx)
+ ind = slice(None)
+ else:
+ grids,ind = self.pf.hierarchy.get_box_grids(
+ self.left_edge - self.dx, self.right_edge + self.dx)
level_ind = na.where(self.pf.hierarchy.gridLevels.ravel()[ind] <= self.level)
sort_ind = na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind])
self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)]
@@ -1831,3 +1844,18 @@
def flush_data(self, *args, **kwargs):
raise KeyError("Can't do this")
+
+
+class EnzoOrthoRayBase(AMROrthoRayBase): pass
+class EnzoRayBase(AMRRayBase): pass
+class EnzoSliceBase(AMRSliceBase): pass
+class EnzoCuttingPlaneBase(AMRCuttingPlaneBase): pass
+class EnzoProjBase(AMRProjBase): pass
+class EnzoCylinderBase(AMRCylinderBase): pass
+class EnzoRegionBase(AMRRegionBase): pass
+class EnzoPeriodicRegionBase(AMRPeriodicRegionBase): pass
+class EnzoGridCollection(AMRGridCollection): pass
+class EnzoSphereBase(AMRSphereBase): pass
+class EnzoCoveringGrid(AMRCoveringGrid): pass
+class EnzoSmoothedCoveringGrid(AMRSmoothedCoveringGrid): pass
+
Modified: branches/yt-non-3d/yt/lagos/BaseGridType.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/BaseGridType.py (original)
+++ branches/yt-non-3d/yt/lagos/BaseGridType.py Sat Nov 15 15:40:05 2008
@@ -27,35 +27,25 @@
#import yt.enki, gc
from yt.funcs import *
-class EnzoGridBase(EnzoData):
+class AMRGridPatch(AMRData):
_spatial = True
_num_ghost_zones = 0
- """
- Class representing a single Enzo Grid instance.
- """
_grids = None
+ _id_offset = 1
def __init__(self, id, filename=None, hierarchy = None):
- """
- Returns an instance of EnzoGrid with *id*, associated with
- *filename* and *hierarchy*.
- """
- #EnzoData's init function is a time-burglar.
- #All of the field parameters will be passed to us as needed.
- #EnzoData.__init__(self, None, [])
self.data = {}
self.field_parameters = {}
self.fields = []
self.start_index = None
self.id = id
+ if (id % 1e4) == 0: mylog.debug("Prepared grid %s", id)
if hierarchy: self.hierarchy = weakref.proxy(hierarchy)
if filename: self.set_filename(filename)
self.overlap_masks = [None, None, None]
self._overlap_grids = [None, None, None]
self._file_access_pooling = False
-
- def __len__(self):
- return na.prod(self.ActiveDimensions)
+ self.pf = self.hierarchy.parameter_file # weakref already
def _generate_field(self, field):
if self.pf.field_info.has_key(field):
@@ -74,7 +64,7 @@
self[field] = self.pf.field_info[field](self)
else: # Can't find the field, try as it might
raise exceptions.KeyError, field
-
+
def get_data(self, field):
"""
Returns a field or set of fields for a key or set of keys
@@ -101,111 +91,17 @@
self._generate_field(field)
return self.data[field]
- def clear_all_grid_references(self):
- """
- This clears out all references this grid has to any others, as
- well as the hierarchy. It's like extra-cleaning after clear_data.
- """
- self.clear_all_derived_quantities()
- if hasattr(self, 'hierarchy'):
- del self.hierarchy
- if hasattr(self, 'Parent'):
- if self.Parent != None:
- self.Parent.clear_all_grid_references()
- del self.Parent
- if hasattr(self, 'Children'):
- for i in self.Children:
- if i != None:
- del i
- del self.Children
-
- def _prepare_grid(self):
- """
- Copies all the appropriate attributes from the hierarchy
- """
- # This is definitely the slowest part of generating the hierarchy
- # Now we give it pointers to all of its attributes
- # Note that to keep in line with Enzo, we have broken PEP-8
- h = self.hierarchy # cache it
- self.Dimensions = h.gridDimensions[self.id-1]
- self.StartIndices = h.gridStartIndices[self.id-1]
- self.EndIndices = h.gridEndIndices[self.id-1]
- self.LeftEdge = h.gridLeftEdge[self.id-1]
- self.RightEdge = h.gridRightEdge[self.id-1]
- self.Level = h.gridLevels[self.id-1,0]
- self.Time = h.gridTimes[self.id-1,0]
- self.NumberOfParticles = h.gridNumberOfParticles[self.id-1,0]
- self.ActiveDimensions = (self.EndIndices - self.StartIndices + 1)
- self.Children = h.gridTree[self.id-1]
- pID = h.gridReverseTree[self.id-1]
- if pID != None and pID != -1:
- self.Parent = weakref.proxy(h.grids[pID - 1])
- else:
- self.Parent = None
-
def _setup_dx(self):
# So first we figure out what the index is. We don't assume
# that dx=dy=dz , at least here. We probably do elsewhere.
- self.dx = self.hierarchy.gridDxs[self.id-1,0]
- self.dy = self.hierarchy.gridDys[self.id-1,0]
- self.dz = self.hierarchy.gridDzs[self.id-1,0]
+ id = self.id - self._id_offset
+ self.dx = self.hierarchy.gridDxs[id,0]
+ self.dy = self.hierarchy.gridDys[id,0]
+ self.dz = self.hierarchy.gridDzs[id,0]
self.data['dx'] = self.dx
self.data['dy'] = self.dy
self.data['dz'] = self.dz
- self._corners = self.hierarchy.gridCorners[:,:,self.id-1]
-
- def _guess_properties_from_parent(self):
- """
- We know that our grid boundary occurs on the cell boundary of our
- parent. This can be a very expensive process, but it is necessary
- in some hierarchys, where yt is unable to generate a completely
- space-filling tiling of grids, possibly due to the finite accuracy in a
- standard Enzo hierarchy file.
- """
- le = self.LeftEdge
- self.dx = self.Parent.dx/2.0
- self.dy = self.Parent.dy/2.0
- self.dz = self.Parent.dz/2.0
- ParentLeftIndex = na.rint((self.LeftEdge-self.Parent.LeftEdge)/self.Parent.dx)
- self.start_index = 2*(ParentLeftIndex + self.Parent.get_global_startindex()).astype('int64')
- self.LeftEdge = self.Parent.LeftEdge + self.Parent.dx * ParentLeftIndex
- self.RightEdge = self.LeftEdge + \
- self.ActiveDimensions*na.array([self.dx,self.dy,self.dz])
- self.hierarchy.gridDxs[self.id-1,0] = self.dx
- self.hierarchy.gridDys[self.id-1,0] = self.dy
- self.hierarchy.gridDzs[self.id-1,0] = self.dz
- self.hierarchy.gridLeftEdge[self.id-1,:] = self.LeftEdge
- self.hierarchy.gridRightEdge[self.id-1,:] = self.RightEdge
- self.hierarchy.gridCorners[:,:,self.id-1] = na.array([ # Unroll!
- [self.LeftEdge[0], self.LeftEdge[1], self.LeftEdge[2]],
- [self.RightEdge[0], self.LeftEdge[1], self.LeftEdge[2]],
- [self.RightEdge[0], self.RightEdge[1], self.LeftEdge[2]],
- [self.RightEdge[0], self.RightEdge[1], self.RightEdge[2]],
- [self.LeftEdge[0], self.RightEdge[1], self.RightEdge[2]],
- [self.LeftEdge[0], self.LeftEdge[1], self.RightEdge[2]],
- [self.RightEdge[0], self.LeftEdge[1], self.RightEdge[2]],
- [self.LeftEdge[0], self.RightEdge[1], self.LeftEdge[2]],
- ], dtype='float64')
- self.__child_mask = None
- self.__child_index_mask = None
- self.__child_indices = None
- self._setup_dx()
-
- def get_global_startindex(self):
- """
- Return the integer starting index for each dimension at the current
- level.
- """
- if self.start_index != None:
- return self.start_index
- if self.Parent == None:
- start_index = self.LeftEdge / na.array([self.dx, self.dy, self.dz])
- return na.rint(start_index).astype('int64').ravel()
- pdx = na.array([self.Parent.dx, self.Parent.dy, self.Parent.dz]).ravel()
- start_index = (self.Parent.get_global_startindex()) + \
- na.rint((self.LeftEdge - self.Parent.LeftEdge)/pdx)
- self.start_index = (start_index*2).astype('int64').ravel()
- return self.start_index
+ self._corners = self.hierarchy.gridCorners[:,:,id]
def _generate_overlap_masks(self, axis, LE, RE):
"""
@@ -227,6 +123,24 @@
def __int__(self):
return self.id
+ def clear_all_grid_references(self):
+ """
+ This clears out all references this grid has to any others, as
+ well as the hierarchy. It's like extra-cleaning after clear_data.
+ """
+ self.clear_all_derived_quantities()
+ if hasattr(self, 'hierarchy'):
+ del self.hierarchy
+ if hasattr(self, 'Parent'):
+ if self.Parent != None:
+ self.Parent.clear_all_grid_references()
+ del self.Parent
+ if hasattr(self, 'Children'):
+ for i in self.Children:
+ if i != None:
+ del i
+ del self.Children
+
def clear_data(self):
"""
Clear out the following things: child_mask, child_indices,
@@ -238,21 +152,35 @@
del self.coarseData
if hasattr(self, 'retVal'):
del self.retVal
- EnzoData.clear_data(self)
+ AMRData.clear_data(self)
self._setup_dx()
- def set_filename(self, filename):
+ def _prepare_grid(self):
"""
- Intelligently set the filename.
+ Copies all the appropriate attributes from the hierarchy
"""
- if self.hierarchy._strip_path:
- self.filename = os.path.join(self.hierarchy.directory,
- os.path.basename(filename))
- elif filename[0] == os.path.sep:
- self.filename = filename
+ # This is definitely the slowest part of generating the hierarchy
+ # Now we give it pointers to all of its attributes
+ # Note that to keep in line with Enzo, we have broken PEP-8
+ h = self.hierarchy # cache it
+ self.Dimensions = h.gridDimensions[self.id-1]
+ self.StartIndices = h.gridStartIndices[self.id-1]
+ self.EndIndices = h.gridEndIndices[self.id-1]
+ self.LeftEdge = h.gridLeftEdge[self.id-1]
+ self.RightEdge = h.gridRightEdge[self.id-1]
+ self.Level = h.gridLevels[self.id-1,0]
+ self.Time = h.gridTimes[self.id-1,0]
+ self.NumberOfParticles = h.gridNumberOfParticles[self.id-1,0]
+ self.ActiveDimensions = (self.EndIndices - self.StartIndices + 1)
+ self.Children = h.gridTree[self.id-1]
+ pID = h.gridReverseTree[self.id-1]
+ if pID != None and pID != -1:
+ self.Parent = weakref.proxy(h.grids[pID - 1])
else:
- self.filename = os.path.join(self.hierarchy.directory, filename)
- return
+ self.Parent = None
+
+ def __len__(self):
+ return na.prod(self.ActiveDimensions)
def find_max(self, field):
"""
@@ -300,28 +228,6 @@
del self.child_mask
del self.child_ind
- def __get_enzo_grid(self):
- """
- **DO NOT USE**
-
- This attempts to get an instance of this particular grid from the SWIG
- interface. Note that it first checks to see if the ParameterFile has
- been instantiated.
- """
- if self.hierarchy.eiTopGrid == None:
- self.hierarchy.initializeEnzoInterface()
- p=re.compile("Grid = %s\n" % (self.id))
- h=open(self.hierarchyFilename,"r").read()
- m=re.search(p,h)
- h=open(self.hierarchyFilename,"r")
- retVal = yt.enki.EnzoInterface.fseek(h, long(m.end()), 0)
- self.eiGrid=yt.enki.EnzoInterface.grid()
- cwd = os.getcwd() # Hate doing this, need to for relative pathnames
- os.chdir(self.hierarchy.directory)
- self.eiGrid.ReadGrid(h, 1)
- os.chdir(cwd)
- mylog.debug("Grid read with SWIG")
-
def _set_child_mask(self, newCM):
if self.__child_mask != None:
mylog.warning("Overriding child_mask attribute! This is probably unwise!")
@@ -369,6 +275,16 @@
self.__child_index_mask = None
#@time_execution
+ def __fill_child_mask(self, child, mask, tofill):
+ startIndex = na.maximum(0, na.rint((child.LeftEdge - self.LeftEdge)/self.dx))
+ endIndex = na.minimum(na.rint((child.RightEdge - self.LeftEdge)/self.dx),
+ self.ActiveDimensions)
+ #startIndex + self.ActiveDimensions)
+ startIndex = na.maximum(0, startIndex)
+ mask[startIndex[0]:endIndex[0],
+ startIndex[1]:endIndex[1],
+ startIndex[2]:endIndex[2]] = tofill
+
def __generate_child_mask(self):
"""
Generates self.child_mask, which is zero where child grids exist (and
@@ -376,12 +292,7 @@
"""
self.__child_mask = na.ones(self.ActiveDimensions, 'int32')
for child in self.Children:
- # Now let's get our overlap
- startIndex = na.rint((child.LeftEdge - self.LeftEdge)/self.dx)
- endIndex = na.rint((child.RightEdge - self.LeftEdge)/self.dx)
- self.__child_mask[startIndex[0]:endIndex[0],
- startIndex[1]:endIndex[1],
- startIndex[2]:endIndex[2]] = 0
+ self.__fill_child_mask(child, self.__child_mask, 0)
self.__child_indices = (self.__child_mask==0) # bool, possibly redundant
def __generate_child_index_mask(self):
@@ -391,12 +302,7 @@
"""
self.__child_index_mask = na.zeros(self.ActiveDimensions, 'int32') - 1
for child in self.Children:
- # Now let's get our overlap
- startIndex = na.rint((child.LeftEdge - self.LeftEdge)/self.dx)
- endIndex = na.rint((child.RightEdge - self.LeftEdge)/self.dx)
- self.__child_index_mask[startIndex[0]:endIndex[0],
- startIndex[1]:endIndex[1],
- startIndex[2]:endIndex[2]] = child.id
+ self.__fill_child_mask(child, self.__child_mask, child.id)
def _get_coords(self):
if self.__coords == None: self._generate_coords()
@@ -468,3 +374,83 @@
for key in data.keys():
if key not in self.__current_data_keys:
del self.data[key]
+
+class EnzoGridBase(AMRGridPatch):
+ """
+ Class representing a single Enzo Grid instance.
+ """
+ def __init__(self, id, filename=None, hierarchy = None):
+ """
+ Returns an instance of EnzoGrid with *id*, associated with
+ *filename* and *hierarchy*.
+ """
+ #All of the field parameters will be passed to us as needed.
+ AMRGridPatch.__init__(self, id, filename, hierarchy)
+ self._file_access_pooling = False
+
+ def _guess_properties_from_parent(self):
+ """
+ We know that our grid boundary occurs on the cell boundary of our
+ parent. This can be a very expensive process, but it is necessary
+ in some hierarchys, where yt is unable to generate a completely
+ space-filling tiling of grids, possibly due to the finite accuracy in a
+ standard Enzo hierarchy file.
+ """
+ le = self.LeftEdge
+ self.dx = self.Parent.dx/2.0
+ self.dy = self.Parent.dy/2.0
+ self.dz = self.Parent.dz/2.0
+ ParentLeftIndex = na.rint((self.LeftEdge-self.Parent.LeftEdge)/self.Parent.dx)
+ self.start_index = 2*(ParentLeftIndex + self.Parent.get_global_startindex()).astype('int64')
+ self.LeftEdge = self.Parent.LeftEdge + self.Parent.dx * ParentLeftIndex
+ self.RightEdge = self.LeftEdge + \
+ self.ActiveDimensions*na.array([self.dx,self.dy,self.dz])
+ self.hierarchy.gridDxs[self.id-1,0] = self.dx
+ self.hierarchy.gridDys[self.id-1,0] = self.dy
+ self.hierarchy.gridDzs[self.id-1,0] = self.dz
+ self.hierarchy.gridLeftEdge[self.id-1,:] = self.LeftEdge
+ self.hierarchy.gridRightEdge[self.id-1,:] = self.RightEdge
+ self.hierarchy.gridCorners[:,:,self.id-1] = na.array([ # Unroll!
+ [self.LeftEdge[0], self.LeftEdge[1], self.LeftEdge[2]],
+ [self.RightEdge[0], self.LeftEdge[1], self.LeftEdge[2]],
+ [self.RightEdge[0], self.RightEdge[1], self.LeftEdge[2]],
+ [self.RightEdge[0], self.RightEdge[1], self.RightEdge[2]],
+ [self.LeftEdge[0], self.RightEdge[1], self.RightEdge[2]],
+ [self.LeftEdge[0], self.LeftEdge[1], self.RightEdge[2]],
+ [self.RightEdge[0], self.LeftEdge[1], self.RightEdge[2]],
+ [self.LeftEdge[0], self.RightEdge[1], self.LeftEdge[2]],
+ ], dtype='float64')
+ self.__child_mask = None
+ self.__child_index_mask = None
+ self.__child_indices = None
+ self._setup_dx()
+
+ def get_global_startindex(self):
+ """
+ Return the integer starting index for each dimension at the current
+ level.
+ """
+ if self.start_index != None:
+ return self.start_index
+ if self.Parent == None:
+ start_index = self.LeftEdge / na.array([self.dx, self.dy, self.dz])
+ return na.rint(start_index).astype('int64').ravel()
+ pdx = na.array([self.Parent.dx, self.Parent.dy, self.Parent.dz]).ravel()
+ start_index = (self.Parent.get_global_startindex()) + \
+ na.rint((self.LeftEdge - self.Parent.LeftEdge)/pdx)
+ self.start_index = (start_index*2).astype('int64').ravel()
+ return self.start_index
+
+ def set_filename(self, filename):
+ """
+ Intelligently set the filename.
+ """
+ if self.hierarchy._strip_path:
+ self.filename = os.path.join(self.hierarchy.directory,
+ os.path.basename(filename))
+ elif filename[0] == os.path.sep:
+ self.filename = filename
+ else:
+ self.filename = os.path.join(self.hierarchy.directory, filename)
+ return
+
Modified: branches/yt-non-3d/yt/lagos/ContourFinder.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/ContourFinder.py (original)
+++ branches/yt-non-3d/yt/lagos/ContourFinder.py Sat Nov 15 15:40:05 2008
@@ -103,7 +103,8 @@
my_queue.add(data_source._grids)
for i,grid in enumerate(my_queue):
max_before = grid["tempContours"].max()
- if na.all(grid.LeftEdge == 0.0) and na.all(grid.RightEdge == 1.0):
+ if na.all(grid.LeftEdge == grid.pf["DomainLeftEdge"]) and \
+ na.all(grid.RightEdge == grid.pf["DomainRightEdge"]):
cg = grid.retrieve_ghost_zones(0,["tempContours","GridIndices"])
else:
cg = grid.retrieve_ghost_zones(1,["tempContours","GridIndices"])
Modified: branches/yt-non-3d/yt/lagos/DataReadingFuncs.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/DataReadingFuncs.py (original)
+++ branches/yt-non-3d/yt/lagos/DataReadingFuncs.py Sat Nov 15 15:40:05 2008
@@ -233,7 +233,7 @@
sets = list(sets)
for g in grids: files_keys[g.filename].append(g)
for file in files_keys:
- mylog.debug("Starting read %s", file)
+ mylog.debug("Starting read %s (%s)", file, sets)
nodes = [g.id for g in files_keys[file]]
nodes.sort()
data = HDF5LightReader.ReadMultipleGrids(file, nodes, sets)
Modified: branches/yt-non-3d/yt/lagos/DerivedFields.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/DerivedFields.py (original)
+++ branches/yt-non-3d/yt/lagos/DerivedFields.py Sat Nov 15 15:40:05 2008
@@ -746,10 +746,11 @@
units=r"\rm{km}\rm{Mpc}/\rm{s}", validators=[ValidateParameter('center')])
def _AngularMomentum(field, data):
return data["CellMass"] * data["SpecificAngularMomentum"]
-add_field("AngularMomentum", units=r"\rm{g}\/\rm{cm}^2/\rm{s}")
+add_field("AngularMomentum", units=r"\rm{g}\/\rm{cm}^2/\rm{s}",
+ vector_field=True)
def _AngularMomentumMSUNKMSMPC(field, data):
return data["CellMassMsun"] * data["SpecificAngularMomentumKMSMPC"]
-add_field("AngularMomentumMSUNKMSMPC",
+add_field("AngularMomentumMSUNKMSMPC", vector_field=True,
units=r"M_{\odot}\rm{km}\rm{Mpc}/\rm{s}")
def _ParticleSpecificAngularMomentum(field, data):
@@ -782,10 +783,11 @@
units=r"\rm{km}\rm{Mpc}/\rm{s}", validators=[ValidateParameter('center')])
def _ParticleAngularMomentum(field, data):
return data["ParticleMass"] * data["ParticleSpecificAngularMomentum"]
-add_field("ParticleAngularMomentum", units=r"\rm{g}\/\rm{cm}^2/\rm{s}")
+add_field("ParticleAngularMomentum", units=r"\rm{g}\/\rm{cm}^2/\rm{s}",
+ vector_field=True)
def _ParticleAngularMomentumMSUNKMSMPC(field, data):
return data["ParticleMass"] * data["ParticleSpecificAngularMomentumKMSMPC"]
-add_field("ParticleAngularMomentumMSUNKMSMPC",
+add_field("ParticleAngularMomentumMSUNKMSMPC", vector_field=True,
units=r"M_{\odot}\rm{km}\rm{Mpc}/\rm{s}")
Modified: branches/yt-non-3d/yt/lagos/HelperFunctions.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/HelperFunctions.py (original)
+++ branches/yt-non-3d/yt/lagos/HelperFunctions.py Sat Nov 15 15:40:05 2008
@@ -38,7 +38,7 @@
orig_shape = data_object[self.x_name].shape
x_vals = data_object[self.x_name].ravel().astype('float64')
- x_i = na.digitize(x_vals, self.x_bins) - 1
+ x_i = (na.digitize(x_vals, self.x_bins) - 1).astype('int32')
if na.any((x_i == -1) | (x_i == len(self.x_bins)-1)):
mylog.error("Sorry, but your values are outside" + \
" the table! Dunno what to do, so dying.")
@@ -63,8 +63,8 @@
x_vals = data_object[self.x_name].ravel().astype('float64')
y_vals = data_object[self.y_name].ravel().astype('float64')
- x_i = na.digitize(x_vals, self.x_bins) - 1
- y_i = na.digitize(y_vals, self.y_bins) - 1
+ x_i = (na.digitize(x_vals, self.x_bins) - 1).astype('int32')
+ y_i = (na.digitize(y_vals, self.y_bins) - 1).astype('int32')
if na.any((x_i == -1) | (x_i == len(self.x_bins)-1)) \
or na.any((y_i == -1) | (y_i == len(self.y_bins)-1)):
if not self.truncate:
Modified: branches/yt-non-3d/yt/lagos/HierarchyType.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/HierarchyType.py (original)
+++ branches/yt-non-3d/yt/lagos/HierarchyType.py Sat Nov 15 15:40:05 2008
@@ -165,18 +165,18 @@
self._data_file = None
def _setup_classes(self, dd):
- self.proj = classobj("EnzoProj",(EnzoProjBase,), dd)
- self.slice = classobj("EnzoSlice",(EnzoSliceBase,), dd)
- self.region = classobj("EnzoRegion",(EnzoRegionBase,), dd)
- self.periodic_region = classobj("EnzoPeriodicRegion",(EnzoPeriodicRegionBase,), dd)
- self.covering_grid = classobj("EnzoCoveringGrid",(EnzoCoveringGrid,), dd)
- self.smoothed_covering_grid = classobj("EnzoSmoothedCoveringGrid",(EnzoSmoothedCoveringGrid,), dd)
- self.sphere = classobj("EnzoSphere",(EnzoSphereBase,), dd)
- self.cutting = classobj("EnzoCuttingPlane",(EnzoCuttingPlaneBase,), dd)
- self.ray = classobj("EnzoRay",(EnzoRayBase,), dd)
- self.ortho_ray = classobj("EnzoOrthoRay",(EnzoOrthoRayBase,), dd)
- self.disk = classobj("EnzoCylinder",(EnzoCylinderBase,), dd)
- self.grid_collection = classobj("EnzoGridCollection",(EnzoGridCollection,), dd)
+ self.proj = classobj("AMRProj",(AMRProjBase,), dd)
+ self.slice = classobj("AMRSlice",(AMRSliceBase,), dd)
+ self.region = classobj("AMRRegion",(AMRRegionBase,), dd)
+ self.periodic_region = classobj("AMRPeriodicRegion",(AMRPeriodicRegionBase,), dd)
+ self.covering_grid = classobj("AMRCoveringGrid",(AMRCoveringGrid,), dd)
+ self.smoothed_covering_grid = classobj("AMRSmoothedCoveringGrid",(AMRSmoothedCoveringGrid,), dd)
+ self.sphere = classobj("AMRSphere",(AMRSphereBase,), dd)
+ self.cutting = classobj("AMRCuttingPlane",(AMRCuttingPlaneBase,), dd)
+ self.ray = classobj("AMRRay",(AMRRayBase,), dd)
+ self.ortho_ray = classobj("AMROrthoRay",(AMROrthoRayBase,), dd)
+ self.disk = classobj("AMRCylinder",(AMRCylinderBase,), dd)
+ self.grid_collection = classobj("AMRGridCollection",(AMRGridCollection,), dd)
def _deserialize_hierarchy(self, harray):
mylog.debug("Cached entry found.")
@@ -537,8 +537,7 @@
self.data_style = data_style
self.hierarchy_filename = os.path.abspath(pf.parameter_filename) \
+ ".hierarchy"
- self.__hierarchy_lines = open(self.hierarchy_filename).readlines()
- if len(self.__hierarchy_lines) == 0:
+ if os.path.getsize(self.hierarchy_filename) == 0:
raise IOError(-1,"File empty", self.hierarchy_filename)
self.boundary_filename = os.path.abspath(pf.parameter_filename) \
+ ".boundary"
@@ -546,11 +545,14 @@
# Now we search backwards from the end of the file to find out how many
# grids we have, which allows us to preallocate memory
self.__hierarchy_string = open(self.hierarchy_filename).read()
- for line in reversed(self.__hierarchy_lines):
- if line.startswith("Grid ="):
- self.num_grids = int(line.split("=")[-1])
+ for line in rlines(open(self.hierarchy_filename, "rb")):
+ if line.startswith("BaryonFileName") or \
+ line.startswith("FileName "):
+ testGrid = line.split("=")[-1].strip().rstrip()
+ if line.startswith("Grid "):
+ self.num_grids = testGridID = int(line.split("=")[-1])
break
- self.__guess_data_style(pf["TopGridRank"])
+ self.__guess_data_style(pf["TopGridRank"], testGrid, testGridID)
# For some reason, r8 seems to want Float64
if pf.has_key("CompilerPrecision") \
and pf["CompilerPrecision"] == "r4":
@@ -560,22 +562,15 @@
AMRHierarchy.__init__(self, pf)
- del self.__hierarchy_string, self.__hierarchy_lines
+ del self.__hierarchy_string
def _setup_classes(self):
dd = self._get_data_reader_dict()
self.grid = classobj("EnzoGrid",(EnzoGridBase,), dd)
AMRHierarchy._setup_classes(self, dd)
- def __guess_data_style(self, rank):
+ def __guess_data_style(self, rank, testGrid, testGridID):
if self.data_style: return
- for line in reversed(self.__hierarchy_lines):
- if line.startswith("BaryonFileName") or \
- line.startswith("FileName "):
- testGrid = line.split("=")[-1].strip().rstrip()
- if line.startswith("Grid "):
- testGridID = int(line.split("=")[-1])
- break
if testGrid[0] != os.path.sep:
testGrid = os.path.join(self.directory, testGrid)
if not os.path.exists(testGrid):
@@ -658,14 +653,13 @@
for v in vals.split():
toAdd[curGrid-1,j] = func(v)
j+=1
- for line_index, line in enumerate(self.__hierarchy_lines):
+ for line_index, line in enumerate(open(self.hierarchy_filename)):
# We can do this the slow, 'reliable' way by stripping
# or we can manually pad all our strings, which speeds it up by a
# factor of about ten
#param, vals = map(strip,line.split("="))
if (line_index % 1e5) == 0:
- mylog.debug("Parsing line % 9i / % 9i",
- line_index, len(self.__hierarchy_lines))
+ mylog.debug("Parsing line % 9i", line_index)
if len(line) < 2:
continue
param, vals = line.split("=")
@@ -985,3 +979,47 @@
rs += "(%s)\s*" % (scanf_regex[t])
rs +="$"
return re.compile(rs,re.M)
+
+# These next two functions are taken from
+# http://www.reddit.com/r/Python/comments/6hj75/reverse_file_iterator/c03vms4
+# Credit goes to "Brian" on Reddit
+
+def rblocks(f, blocksize=4096*256):
+ """Read file as series of blocks from end of file to start.
+
+ The data itself is in normal order, only the order of the blocks is reversed.
+ ie. "hello world" -> ["ld","wor", "lo ", "hel"]
+ Note that the file must be opened in binary mode.
+ """
+ if 'b' not in f.mode.lower():
+ raise Exception("File must be opened using binary mode.")
+ size = os.stat(f.name).st_size
+ fullblocks, lastblock = divmod(size, blocksize)
+
+ # The first(end of file) block will be short, since this leaves
+ # the rest aligned on a blocksize boundary. This may be more
+ # efficient than having the last (first in file) block be short
+ f.seek(-lastblock,2)
+ yield f.read(lastblock)
+
+ for i in range(fullblocks-1,-1, -1):
+ f.seek(i * blocksize)
+ yield f.read(blocksize)
+
+def rlines(f, keepends=False):
+ """Iterate through the lines of a file in reverse order.
+
+ If keepends is true, line endings are kept as part of the line.
+ """
+ buf = ''
+ for block in rblocks(f):
+ buf = block + buf
+ lines = buf.splitlines(keepends)
+ # Return all lines except the first (since may be partial)
+ if lines:
+ lines.reverse()
+ buf = lines.pop() # Last line becomes end of new first line.
+ for line in lines:
+ yield line
+ yield buf # First line.
+
Modified: branches/yt-non-3d/yt/lagos/ParallelTools.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/ParallelTools.py (original)
+++ branches/yt-non-3d/yt/lagos/ParallelTools.py Sat Nov 15 15:40:05 2008
@@ -28,7 +28,8 @@
import yt.logger
import itertools, sys
-if os.path.basename(sys.argv[0]) == "mpi4py":
+if os.path.basename(sys.executable) in ["mpi4py"] \
+ or "--parallel" in sys.argv:
from mpi4py import MPI
parallel_capable = (MPI.COMM_WORLD.size > 1)
if parallel_capable:
@@ -108,6 +109,17 @@
if not self.just_list: self.pobj._finalize_parallel()
raise StopIteration
+def parallel_simple_proxy(func):
+ if not parallel_capable: return func
+ @wraps(func)
+ def single_proc_results(self, *args, **kwargs):
+ retval = None
+ if self._owned:
+ retval = func(self, *args, **kwargs)
+ MPI.COMM_WORLD.Bcast(retval)
+ return retval
+ return single_proc_results
+
def parallel_passthrough(func):
@wraps(func)
def passage(self, data):
@@ -152,7 +164,8 @@
LE[y_dict[axis]] = y[0]
RE[y_dict[axis]] = y[1]
- return True, self.hierarchy.region(self.center, LE, RE)
+ reg = self.hierarchy.region(self.center, LE, RE)
+ return True, reg
def _partition_hierarchy_3d(self, padding=0.0):
if not parallel_capable:
@@ -179,12 +192,21 @@
def _mpi_catdict(self, data):
mylog.debug("Opening MPI Barrier on %s", MPI.COMM_WORLD.rank)
MPI.COMM_WORLD.Barrier()
- if MPI.COMM_WORLD.rank == 0:
- data = self.__mpi_recvdict(data)
- else:
- MPI.COMM_WORLD.Send(data, dest=0, tag=0)
- mylog.debug("Opening MPI Broadcast on %s", MPI.COMM_WORLD.rank)
- data = MPI.COMM_WORLD.Bcast(data, root=0)
+ field_keys = data.keys()
+ field_keys.sort()
+ np = MPI.COMM_WORLD.size
+ for key in field_keys:
+ mylog.debug("Joining %s (%s) on %s", key, type(data[key]),
+ MPI.COMM_WORLD.rank)
+ if MPI.COMM_WORLD.rank == 0:
+ data[key] = na.concatenate([data[key]] +
+ [MPI.COMM_WORLD.Recv(source=i, tag=0) for i in range(1, np)],
+ axis=-1)
+ else:
+ MPI.COMM_WORLD.Send(data[key], dest=0, tag=0)
+ MPI.COMM_WORLD.Barrier()
+ data[key] = MPI.COMM_WORLD.Bcast(data[key], root=0)
+ mylog.debug("Done joining dictionary on %s", MPI.COMM_WORLD.rank)
MPI.COMM_WORLD.Barrier()
return data
Modified: branches/yt-non-3d/yt/lagos/PointCombine.c
==============================================================================
--- branches/yt-non-3d/yt/lagos/PointCombine.c (original)
+++ branches/yt-non-3d/yt/lagos/PointCombine.c Sat Nov 15 15:40:05 2008
@@ -43,23 +43,29 @@
Py_CombineGrids(PyObject *obj, PyObject *args)
{
PyObject *ogrid_src_x, *ogrid_src_y, *ogrid_src_vals,
- *ogrid_src_mask, *ogrid_src_wgt;
+ *ogrid_src_mask, *ogrid_src_wgt, *ogrid_used_mask;
PyObject *ogrid_dst_x, *ogrid_dst_y, *ogrid_dst_vals,
*ogrid_dst_mask, *ogrid_dst_wgt;
PyArrayObject *grid_src_x, *grid_src_y, **grid_src_vals,
- *grid_src_mask, *grid_src_wgt;
+ *grid_src_mask, *grid_src_wgt, *grid_used_mask;
PyArrayObject *grid_dst_x, *grid_dst_y, **grid_dst_vals,
*grid_dst_mask, *grid_dst_wgt;
+ grid_src_x = grid_src_y = //grid_src_vals =
+ grid_src_mask = grid_src_wgt = grid_used_mask =
+ grid_dst_x = grid_dst_y = //grid_dst_vals =
+ grid_dst_mask = grid_dst_wgt = NULL;
+
int NumArrays, src_len, dst_len, refinement_factor;
+ NumArrays = 0;
- if (!PyArg_ParseTuple(args, "OOOOOOOOOOi",
+ if (!PyArg_ParseTuple(args, "OOOOOOOOOOiO",
&ogrid_src_x, &ogrid_src_y,
&ogrid_src_mask, &ogrid_src_wgt, &ogrid_src_vals,
&ogrid_dst_x, &ogrid_dst_y,
&ogrid_dst_mask, &ogrid_dst_wgt, &ogrid_dst_vals,
- &refinement_factor))
+ &refinement_factor, &ogrid_used_mask))
return PyErr_Format(_combineGridsError,
"CombineGrids: Invalid parameters.");
@@ -97,6 +103,15 @@
goto _fail;
}
+ grid_used_mask = (PyArrayObject *) PyArray_FromAny(ogrid_used_mask,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((grid_used_mask == NULL) || (PyArray_SIZE(grid_used_mask) != src_len)) {
+ PyErr_Format(_combineGridsError,
+ "CombineGrids: src_x and used_mask must be the same shape.");
+ goto _fail;
+ }
+
grid_dst_x = (PyArrayObject *) PyArray_FromAny(ogrid_dst_x,
PyArray_DescrFromType(NPY_INT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
@@ -172,6 +187,7 @@
npy_int64 *src_y = (npy_int64 *) PyArray_GETPTR1(grid_src_y,0);
npy_float64 *src_wgt = (npy_float64 *) PyArray_GETPTR1(grid_src_wgt,0);
npy_int64 *src_mask = (npy_int64 *) PyArray_GETPTR1(grid_src_mask,0);
+ npy_int64 *src_used_mask = (npy_int64 *) PyArray_GETPTR1(grid_used_mask,0);
npy_int64 *dst_x = (npy_int64 *) PyArray_GETPTR1(grid_dst_x,0);
npy_int64 *dst_y = (npy_int64 *) PyArray_GETPTR1(grid_dst_y,0);
@@ -183,7 +199,7 @@
int num_found = 0;
for (si = 0; si < src_len; si++) {
- if (src_x[si] < 0) continue;
+ if (src_used_mask[si] == 0) continue;
init_x = refinement_factor * src_x[si];
init_y = refinement_factor * src_y[si];
for (x_off = 0; x_off < refinement_factor; x_off++) {
@@ -191,7 +207,6 @@
fine_x = init_x + x_off;
fine_y = init_y + y_off;
for (di = 0; di < dst_len; di++) {
- if (dst_x[di] < 0) continue;
if ((fine_x == dst_x[di]) &&
(fine_y == dst_y[di])) {
num_found++;
@@ -200,7 +215,7 @@
((refinement_factor != 1) && (dst_mask[di])));
// So if they are on the same level, then take the logical and
// otherwise, set it to the destination mask
- src_x[si] = -2;
+ src_used_mask[si] = 0;
for (i = 0; i < NumArrays; i++) {
dst_vals[i][di] += src_vals[i][si];
}
@@ -215,6 +230,7 @@
Py_DECREF(grid_src_y);
Py_DECREF(grid_src_mask);
Py_DECREF(grid_src_wgt);
+ Py_DECREF(grid_used_mask);
Py_DECREF(grid_dst_x);
Py_DECREF(grid_dst_y);
@@ -241,6 +257,7 @@
Py_XDECREF(grid_src_y);
Py_XDECREF(grid_src_wgt);
Py_XDECREF(grid_src_mask);
+ Py_XDECREF(grid_used_mask);
Py_XDECREF(grid_dst_x);
Py_XDECREF(grid_dst_y);
@@ -703,22 +720,14 @@
for(i=0;i<3;i++) {
itc = 1;
if (ac_le[i] < adl_edge[i]) {
- ac_le_p[i][itc ] = adr_edge[i] + ac_le[i];
- ac_re_p[i][itc++] = adr_edge[i] + ac_re[i];
- //fprintf(stderr, "le axis: %d (%d) le: %0.5f re: %0.5f\n",
- // i, itc-1, ac_le_p[i][itc-1], ac_re_p[i][itc-1]);
+ ac_le_p[i][itc ] = adr_edge[i] - (adl_edge[i] - ac_le[i]);
+ ac_re_p[i][itc++] = adr_edge[i] + (ac_re[i] - adl_edge[i]);
}
if (ac_re[i] > adr_edge[i]) {
- ac_le_p[i][itc ] = adl_edge[i] - (adr_edge[i]-adl_edge[i])
- + ac_le[i];
- ac_re_p[i][itc++] = adl_edge[i] - (adr_edge[i]-adl_edge[i])
- + ac_re[i];
- //fprintf(stderr, "re axis: %d (%d) le: %0.5f re: %0.5f\n",
- // i, itc-1, ac_le_p[i][itc-1], ac_re_p[i][itc-1]);
+ ac_le_p[i][itc ] = ac_le[i] - (adr_edge[i] - adl_edge[i]);
+ ac_re_p[i][itc++] = adl_edge[i] + (ac_re[i] - adr_edge[i]);
}
p_niter[i] = itc;
- //fprintf(stderr, "Iterating on %d %d (%d) times\n",
- // i, itc, p_niter[i]);
}
for (xg = 0; xg < g_data[0]->dimensions[0]; xg++) {
Modified: branches/yt-non-3d/yt/lagos/hop/NewHopOutput.py
==============================================================================
--- branches/yt-non-3d/yt/lagos/hop/NewHopOutput.py (original)
+++ branches/yt-non-3d/yt/lagos/hop/NewHopOutput.py Sat Nov 15 15:40:05 2008
@@ -3,9 +3,13 @@
Author: Matthew Turk <matthewturk at gmail.com>
Affiliation: KIPAC/SLAC/Stanford
+
+Author: Stephen Skory <sskory at physics.ucsd.edu>
+Affiliation: UC San Diego
+
Homepage: http://yt.enzotools.org/
License:
- Copyright (C) 2008 Matthew Turk. All Rights Reserved.
+ Copyright (C) 2008 Matthew Turk, Stephen Skory. All Rights Reserved.
This file is part of yt.
@@ -70,6 +74,10 @@
self._base_indices = na.arange(tot_part)[ii]
def __enlarge_data(self):
+ """
+ Take the data fields and reflect them around the sides of the
+ box, using periodicity.
+ """
sizetemp = self.particle_fields["particle_position_x"].size
self.tempx = [i for i in range(3*sizetemp)]
self.tempy = [i for i in range(3*sizetemp)]
@@ -98,6 +106,11 @@
self.tempm[3*sizetemp -1 - part] = i
def __slice_data(self):
+ """
+ Cut up the particle data into subboxes, depending on *cuts* and the over-sampling
+ parameter *padding*. This also enlarges the box such that the particles run
+ from 0 to 1, which may not be neccesary.
+ """
cuts = 3
padding = 0.2
self.temp2x = {}
@@ -143,6 +156,9 @@
self.tracking2[((i*cuts) + j)*cuts + k] = t_t
def __reduce_data(self):
+ """
+ Take the boxes after HOP has handled them and make them small, to their true size.
+ """
cuts = 3
padding = 0.2
# loop over the sub-boxes
@@ -162,6 +178,9 @@
self.temp2z[((i*cuts) + j)*cuts + k] = self.temp2z[((i*cuts) + j)*cuts + k]*(zright-zleft)+zleft
def __run_hops(self):
+ """
+ Run HOP on the different boxes.
+ """
cuts = 3
padding = 0.2
nParts = self.particle_fields["particle_position_x"].size
Modified: branches/yt-non-3d/yt/raven/Callbacks.py
==============================================================================
--- branches/yt-non-3d/yt/raven/Callbacks.py (original)
+++ branches/yt-non-3d/yt/raven/Callbacks.py Sat Nov 15 15:40:05 2008
@@ -3,9 +3,13 @@
Author: Matthew Turk <matthewturk at gmail.com>
Affiliation: KIPAC/SLAC/Stanford
+Author: J. S. Oishi <jsoishi at astro.berkeley.edu>
+Affiliation: UC Berkeley
+Author: Stephen Skory <sskory at physics.ucsd.edu>
+Affiliation: UC San Diego
Homepage: http://yt.enzotools.org/
License:
- Copyright (C) 2008 Matthew Turk. All Rights Reserved.
+ Copyright (C) 2008 Matthew Turk, JS Oishi, Stephen Skory. All Rights Reserved.
This file is part of yt.
@@ -539,3 +543,69 @@
plot._axes.text(center_x, center_y, "%s" % i,
fontsize=self.font_size)
+class CoordAxesCallback(PlotCallback):
+ """Creates x and y axes for a VMPlot. In the future, it will
+ attempt to guess the proper units to use.
+
+ """
+ def __init__(self,unit=None,coords=False):
+ PlotCallback.__init__(self)
+ self.unit = unit
+ self.coords = coords
+
+ def __call__(self,plot):
+ # 1. find out what the domain is
+ # 2. pick a unit for it
+ # 3. run self._axes.set_xlabel & self._axes.set_ylabel to actually lay shit down.
+ # 4. adjust extent information to make sure labels are visable.
+
+ # put the plot into data coordinates
+ nx,ny = plot.image._A.shape
+ dx = (plot.xlim[1] - plot.xlim[0])/nx
+ dy = (plot.ylim[1] - plot.ylim[0])/ny
+
+ unit_conversion = plot.data.hierarchy[plot.im["Unit"]]
+ aspect = (plot.xlim[1]-plot.xlim[0])/(plot.ylim[1]-plot.ylim[0])
+
+ print "aspect ratio = ", aspect
+
+ # if coords is False, label axes relative to the center of the
+ # display. if coords is True, label axes with the absolute
+ # coordinates of the region.
+ xcenter = 0.
+ ycenter = 0.
+ if not self.coords:
+ center = plot.data.center
+ if plot.data.axis == 0:
+ xcenter = center[1]
+ ycenter = center[2]
+ elif plot.data.axis == 1:
+ xcenter = center[0]
+ ycenter = center[2]
+ else:
+ xcenter = center[0]
+ ycenter = center[1]
+
+
+ xformat_function = lambda a,b: '%7.1e' %((a*dx + plot.xlim[0] - xcenter)*unit_conversion)
+ yformat_function = lambda a,b: '%7.1e' %((a*dy + plot.ylim[0] - ycenter)*unit_conversion)
+ else:
+ xformat_function = lambda a,b: '%7.1e' %((a*dx + plot.xlim[0])*unit_conversion)
+ yformat_function = lambda a,b: '%7.1e' %((a*dy + plot.ylim[0])*unit_conversion)
+
+ xticker = matplotlib.ticker.FuncFormatter(xformat_function)
+ yticker = matplotlib.ticker.FuncFormatter(yformat_function)
+ plot._axes.xaxis.set_major_formatter(xticker)
+ plot._axes.yaxis.set_major_formatter(yticker)
+
+ xlabel = '%s (%s)' % (lagos.axis_labels[plot.data.axis][0],plot.im["Unit"])
+ ylabel = '%s (%s)' % (lagos.axis_labels[plot.data.axis][1],plot.im["Unit"])
+ xticksize = nx/4.
+ yticksize = ny/4.
+ plot._axes.xaxis.set_major_locator(matplotlib.ticker.FixedLocator([i*xticksize for i in range(0,5)]))
+ plot._axes.yaxis.set_major_locator(matplotlib.ticker.FixedLocator([i*yticksize for i in range(0,5)]))
+
+ plot._axes.set_xlabel(xlabel,visible=True)
+ plot._axes.set_ylabel(ylabel,visible=True)
+ plot._figure.subplots_adjust(left=0.1,right=0.8)
+
Modified: branches/yt-non-3d/yt/raven/PlotTypes.py
==============================================================================
--- branches/yt-non-3d/yt/raven/PlotTypes.py (original)
+++ branches/yt-non-3d/yt/raven/PlotTypes.py Sat Nov 15 15:40:05 2008
@@ -229,7 +229,7 @@
shrink=0.95)
else:
self.colorbar = None
- self.set_width(1,'1')
+ self.set_width(1,'unitary')
def _get_buff(self, width=None):
x0, x1 = self.xlim
Modified: branches/yt-non-3d/yt/recipes.py
==============================================================================
--- branches/yt-non-3d/yt/recipes.py (original)
+++ branches/yt-non-3d/yt/recipes.py Sat Nov 15 15:40:05 2008
@@ -63,6 +63,15 @@
# yt-generalization : needs to be changed to 'unitary'
return 0.1 / pf["1"]
+def _fix_width(pf, width):
+ if width is not None:
+ if iterable(width):
+ return width[0] / pf[width[1]]
+ return width
+ mylog.info("Setting width to be the full box")
+ # yt-generalization : needs to be changed to 'unitary'
+ return 1.0 / pf["1"]
+
def _fix_axis(pf, axis):
if axis is None:
raise ValueError("You need to specify an Axis!")
@@ -78,9 +87,15 @@
_arg_fixer = {
'center' : (True, _fix_center),
'radius' : (True, _fix_radius),
+ 'width' : (True, _fix_width),
'axis' : (True, _fix_axis),
}
+#
+# * Need to change filename to be an argument in the descriptor
+# * Need to add cmap to the descriptor
+#
+
def fix_plot_args(plot_func):
@wraps(plot_func)
def call_func(self, pf, **kwargs):
@@ -126,4 +141,24 @@
p.switch_x("CellMassMsun")
return pc
+ @fix_plot_args
+ def slice_axis(self, pf, field="Density", center=None, width=None, filename=None, axis=None, coord=None):
+ pc = which_pc(pf, center=center)
+ p = pc.add_slice(field, axis, coord=coord)
+ pc.set_width(width, '1')
+ return pc
+
+ @fix_plot_args
+ def slice(self, pf, field="Density", center=None, width=None, coord=None, filename=None):
+ pc = which_pc(pf, center=center)
+ for axis in range(3): pc.add_slice(field, axis, coord=coord)
+ pc.set_width(width, '1')
+ return pc
+
+ def dispatch(self):
+ pass
+
RecipeBook = _RecipeBook()
+
+if __name__ == "__main__":
+ RecipeBook.dispatch()
More information about the yt-svn
mailing list