[Yt-svn] yt-commit r1533 - in trunk: tests yt yt/extensions yt/extensions/volume_rendering yt/lagos yt/lagos/hop yt/raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri Nov 20 20:27:34 PST 2009
Author: mturk
Date: Fri Nov 20 20:27:29 2009
New Revision: 1533
URL: http://yt.enzotools.org/changeset/1533
Log:
The big merge back from mercurial. This includes a lot of work from Stephen
Skory and Britton Smith, as well as stuff I've done for the refactoring of the
code. See the message to the mailing list for more information on all the
changes.
Added:
trunk/yt/extensions/volume_rendering/CUDARayCast.py
trunk/yt/extensions/volume_rendering/_cuda_caster.cu
trunk/yt/lagos/ObjectFindingMixin.py
trunk/yt/lagos/ParticleIO.py
Modified:
trunk/tests/test_lagos.py
trunk/yt/config.py
trunk/yt/extensions/HierarchySubset.py
trunk/yt/extensions/volume_rendering/grid_partitioner.py
trunk/yt/extensions/volume_rendering/software_sampler.py
trunk/yt/funcs.py
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/BaseGridType.py
trunk/yt/lagos/DataReadingFuncs.py
trunk/yt/lagos/DerivedQuantities.py
trunk/yt/lagos/EnzoFields.py
trunk/yt/lagos/FieldInfoContainer.py
trunk/yt/lagos/HDF5LightReader.c
trunk/yt/lagos/HaloFinding.py
trunk/yt/lagos/HierarchyType.py
trunk/yt/lagos/OutputTypes.py
trunk/yt/lagos/ParallelTools.py
trunk/yt/lagos/Profiles.py
trunk/yt/lagos/UniversalFields.py
trunk/yt/lagos/__init__.py
trunk/yt/lagos/hop/EnzoHop.c
trunk/yt/lagos/hop/hop_hop.c
trunk/yt/lagos/hop/hop_kd.c
trunk/yt/performance_counters.py
trunk/yt/raven/Callbacks.py
trunk/yt/raven/PlotTypes.py
Modified: trunk/tests/test_lagos.py
==============================================================================
--- trunk/tests/test_lagos.py (original)
+++ trunk/tests/test_lagos.py Fri Nov 20 20:27:29 2009
@@ -99,7 +99,7 @@
self.assert_(child.Parent.id == self.hierarchy.grids[0].id)
def testGetSelectLevels(self):
- for level in range(self.hierarchy.maxLevel+1):
+ for level in range(self.hierarchy.max_level+1):
for grid in self.hierarchy.select_grids(level):
self.assert_(grid.Level == level)
@@ -181,8 +181,8 @@
def _returnProfile2DFunction(field, weight, accumulation, lazy):
def add_field_function(self):
self.data.set_field_parameter("center",[.5,.5,.5])
- cv_min = self.hierarchy.gridDxs.min()**3.0
- cv_max = self.hierarchy.gridDxs.max()**3.0
+ cv_min = self.hierarchy.get_smallest_dx()**3.0
+ cv_max = 1.0 / max(self.OutputFile["TopGridDimensions"])
profile = yt.lagos.BinnedProfile2D(self.data,
8, "RadiusCode", 1e-3, 1.0, True,
8, "CellVolumeCode", cv_min, cv_max, True, lazy)
Modified: trunk/yt/config.py
==============================================================================
--- trunk/yt/config.py (original)
+++ trunk/yt/config.py Fri Nov 20 20:27:29 2009
@@ -58,6 +58,7 @@
'LogLevel': '20',
'timefunctions':'False',
'inGui':'False',
+ 'inline':'False',
'__parallel':'False',
'__parallel_rank':'0',
'__parallel_size':'1',
Modified: trunk/yt/extensions/HierarchySubset.py
==============================================================================
--- trunk/yt/extensions/HierarchySubset.py (original)
+++ trunk/yt/extensions/HierarchySubset.py Fri Nov 20 20:27:29 2009
@@ -25,6 +25,7 @@
from yt.mods import *
+from yt.lagos import AMRGridPatch, StaticOutput, AMRHierarchy
import h5py, os.path
import yt.commands as commands
@@ -32,66 +33,63 @@
class DummyHierarchy(object):
pass
-class ConstructedRootGrid(object):
- id = -1
- def __init__(self, pf, level, left_edge, right_edge):
+class ConstructedRootGrid(AMRGridPatch):
+ __slots__ = ['base_grid', 'id', 'base_pf']
+ _id_offset = 1
+ def __init__(self, base_pf, pf, hierarchy, level, left_edge, right_edge):
"""
This is a fake root grid, constructed by creating a
:class:`yt.lagos.CoveringGridBase` at a given *level* between
*left_edge* and *right_edge*.
"""
self.pf = pf
- self.hierarchy = DummyHierarchy()
- self.hierarchy.data_style = -1
+ self.base_pf = base_pf
+ self.field_parameters = {}
+ self.fields = []
+ self.NumberOfParticles = 0
+ self.id = 1
+ self.hierarchy = hierarchy
+ self._child_mask = self._child_indices = self._child_index_mask = None
self.Level = level
self.LeftEdge = left_edge
self.RightEdge = right_edge
- self.index = na.min([grid.get_global_startindex() for grid in
- pf.h.select_grids(level)], axis=0).astype('int64')
- self.dds = pf.h.select_grids(level)[0].dds.copy()
+ self.start_index = na.min([grid.get_global_startindex() for grid in
+ base_pf.h.select_grids(level)], axis=0).astype('int64')
+ self.dds = base_pf.h.select_grids(level)[0].dds.copy()
dims = (self.RightEdge-self.LeftEdge)/self.dds
self.ActiveDimensions = dims
print "Constructing base grid of size %s" % (self.ActiveDimensions)
- self.cg = pf.h.smoothed_covering_grid(level, self.LeftEdge,
+ self.base_grid = base_pf.h.smoothed_covering_grid(level, self.LeftEdge,
self.RightEdge, dims=dims)
- self._calculate_child_masks()
+ self.base_grid.Level = self.base_grid.level
+ self.data = {}
+ #self._calculate_child_masks()
+ self.Parent = None
+ self.Children = []
- def _calculate_child_masks(self):
- # This might be slow
- grids, grid_ind = self.pf.hierarchy.get_box_grids(
- self.LeftEdge, self.RightEdge)
- self.Children = [g for g in grids if g.Level == self.Level + 1]
- self.child_mask = na.ones(self.ActiveDimensions, dtype='int32')
- for c in self.Children:
- si = na.maximum(0, na.rint((c.LeftEdge - self.LeftEdge)/self.dds))
- ei = na.minimum(self.ActiveDimensions,
- na.rint((c.RightEdge - self.LeftEdge)/self.dds))
- self.child_mask[si[0]:ei[0], si[1]:ei[1], si[2]:ei[2]] = 0
-
- def __getitem__(self, field):
- return self.cg[field]
+ def get_vertex_centered_data(self, field, smoothed=True):
+ vc = self.base_pf.h.smoothed_covering_grid(self.base_grid.Level,
+ self.base_grid.LeftEdge - self.base_grid.dds*0.5,
+ self.base_grid.RightEdge + self.base_grid.dds*0.5,
+ dims = self.ActiveDimensions + 1)
+ return vc[field]
+
+class AMRExtractedGridProxy(AMRGridPatch):
+ __slots__ = ['base_grid']
+ _id_offset = 1
+ def __init__(self, grid_id, base_grid, hierarchy):
+ # We make a little birdhouse in our soul for the base_grid
+ # (they're the only bee in our bonnet!)
+ self.base_grid = base_grid
+ AMRGridPatch.__init__(self, grid_id, filename = None, hierarchy=hierarchy)
+ self.Parent = None
+ self.Children = []
+ self.Level = -1
- def get_global_startindex(self):
- return self.index
+ def get_vertex_centered_data(self, *args, **kwargs):
+ return self.base_grid.get_vertex_centered_data(*args, **kwargs)
- def get_vertex_centered_data(self, field, smoothed=True):
- # We discard the 'smoothed' keyword argument here
- cg = self.pf.h.smoothed_covering_grid(self.Level,
- self.LeftEdge - self.dds,
- self.RightEdge + self.dds,
- dims = self.ActiveDimensions + 2,
- num_ghost_zones = 1, fields=[field])
- bds = na.array(zip(cg.left_edge+cg.dds/2.0, cg.right_edge-cg.dds/2.0)).ravel()
- interp = lagos.TrilinearFieldInterpolator(cg[field], bds, ['x','y','z'])
- ad = self.ActiveDimensions + 1
- x,y,z = na.mgrid[self.LeftEdge[0]:self.RightEdge[0]:ad[0]*1j,
- self.LeftEdge[1]:self.RightEdge[1]:ad[1]*1j,
- self.LeftEdge[2]:self.RightEdge[2]:ad[2]*1j]
- dd = {'x':x,'y':y,'z':z}
- scalars = interp(dd)
- return scalars
-
-class ExtractedHierarchy(object):
+class OldExtractedHierarchy(object):
def __init__(self, pf, min_level, max_level = -1, offset = None,
always_copy=False):
@@ -196,3 +194,160 @@
def _convert_coords(self, val):
return (val - self.left_edge_offset)*self.mult_factor
+
+class ExtractedHierarchy(AMRHierarchy):
+
+ grid = AMRExtractedGridProxy
+
+ def __init__(self, pf, data_style):
+ # First we set up our translation between original and extracted
+ self.data_style = data_style
+ self.min_level = pf.min_level
+ self.int_offset = na.min([grid.get_global_startindex() for grid in
+ pf.base_pf.h.select_grids(pf.min_level)], axis=0).astype('float64')
+ min_left = na.min([grid.LeftEdge for grid in
+ pf.base_pf.h.select_grids(pf.min_level)], axis=0).astype('float64')
+ max_right = na.max([grid.RightEdge for grid in
+ pf.base_pf.h.select_grids(pf.min_level)], axis=0).astype('float64')
+ level_dx = pf.base_pf.h.select_grids(pf.min_level)[0].dds[0]
+ dims = ((max_right-min_left)/level_dx)
+ max_right += (dims.max() - dims) * level_dx
+ offset = pf.offset
+ if offset is None: offset = min_left
+ self.left_edge_offset = offset
+ pf.offset = offset
+ self.mult_factor = 2**pf.min_level
+ self.min_left_edge = self._convert_coords(min_left)
+ self.max_right_edge = self._convert_coords(max_right)
+ self.min_left, self.max_right = min_left, max_right
+ if pf.max_level == -1: max_level = pf.base_pf.h.max_level
+ self.max_level = min(max_level, pf.base_pf.h.max_level)
+ self.final_level = self.max_level - self.min_level
+
+ # Now we utilize the existing machinery for generating the appropriate
+ # arrays of grids, etc etc.
+ self.base_pf = pf.base_pf
+ AMRHierarchy.__init__(self, pf, data_style)
+
+ # Now a few cleanups
+ self.pf.override["DomainRightEdge"] = self.max_right_edge
+ self.pf.override["DomainLeftEdge"] = self.min_left_edge
+ for u,v in self.base_pf.units.items():
+ self.pf.override[u] = v / self.mult_factor
+ self.pf.override['unitary'] = 1.0 / (self.pf["DomainRightEdge"] -
+ self.pf["DomainLeftEdge"]).max()
+
+ def _count_grids(self):
+ self.num_grids = 1 + sum( ( # 1 is the base grid
+ len(self.base_pf.h.select_grids(level))
+ for level in range(self.min_level+1, self.max_level)) )
+
+ def _parse_hierarchy(self):
+ # Here we need to set up the grid info, which for the Enzo hierarchy
+ # is done like:
+ # self.grid_dimensions.flat[:] = ei
+ # self.grid_dimensions -= na.array(si, self.float_type)
+ # self.grid_dimensions += 1
+ # self.grid_left_edge.flat[:] = LE
+ # self.grid_right_edge.flat[:] = RE
+ # self.grid_particle_count.flat[:] = np
+ # self.grids = na.array(self.grids, dtype='object')
+ #
+ # For now, we make the presupposition that all of our grids are
+ # strictly nested and we are not doing any cuts. However, we do
+ # construct a root grid!
+ root_level_grids = self.base_pf.h.select_grids(self.min_level)
+ base_grid = ConstructedRootGrid(self.base_pf, self.pf, self,
+ self.min_level, self.min_left, self.max_right)
+ self._fill_grid_arrays(base_grid, 0)
+ grids = [base_grid]
+ # We need to ensure we have the correct parentage relationships
+ # However, we want the parent/child to be to the new proxy grids
+ # so we need to map between the old ids and the new ids
+ self.id_map = {}
+ grid_id = 2 # id 0 is the base grid
+ for level in range(self.min_level+1, self.max_level):
+ for grid in self.base_pf.h.select_grids(level):
+ # This next little bit will have to be changed if we ever move to
+ # not-strictly-nested AMR hierarchies
+ parent = self.id_map.get(grid.Parent.id, base_grid)
+ grids.append(self.grid(grid_id, grid, self))
+ parent.Children.append(grids[-1])
+ grids[-1].Parent = parent
+ self.id_map[grid.id] = grids[-1]
+ # Now we fill in our arrays of values -- note that we
+ # are filling in values from the base grids, not the newly
+ # extracted grids. We will perform bulk changes after we
+ # finish.
+ self._fill_grid_arrays(grid, grid_id-1)
+ grid_id += 1
+
+ self.grid_left_edge = self._convert_coords(self.grid_left_edge)
+ self.grid_right_edge = self._convert_coords(self.grid_right_edge)
+ self.grids = na.array(grids, dtype='object')
+
+ def _fill_grid_arrays(self, grid, i):
+ # This just fills in the grid arrays for a single grid --
+ # note that we presuppose here that we are being handed a grid object
+ # that has these defined; this means we are being handed the *base*
+ # grid, not the newly extracted one
+ self.grid_dimensions[i,:] = grid.ActiveDimensions
+ self.grid_left_edge[i,:] = grid.LeftEdge
+ self.grid_right_edge[i,:] = grid.RightEdge
+ self.grid_particle_count[i] = grid.NumberOfParticles
+
+ def _populate_grid_objects(self):
+ for grid in self.grids:
+ grid.Level = grid.base_grid.Level - self.pf.min_level
+ grid._prepare_grid()
+ grid._setup_dx()
+ grid.start_index = None
+ self.max_level -= self.pf.min_level
+ print "New max level:", self.max_level
+
+ def _convert_coords(self, val):
+ return (val - self.left_edge_offset)*self.mult_factor
+
+ def _detect_fields(self):
+ self.field_list = self.base_pf.h.field_list[:]
+
+ def _setup_unknown_fields(self):
+ pass # Done in the base_h
+
+ def _setup_derived_fields(self):
+ self.derived_field_list = self.base_pf.h.derived_field_list[:]
+
+ def _initialize_data_storage(self):
+ self._data_file = None
+
+ def _setup_classes(self):
+ dd = self._get_data_reader_dict()
+ AMRHierarchy._setup_classes(self, dd)
+ self.object_types.sort()
+
+class ExtractedParameterFile(StaticOutput):
+ _hierarchy_class = ExtractedHierarchy
+ data_style = "extracted"
+
+ def __init__(self, base_pf, min_level, max_level = -1, offset = None):
+ self.base_pf = base_pf
+ self.min_level = min_level
+ self.max_level = max_level
+ self.offset = offset
+ self.override = {}
+
+ def __repr__(self):
+ return "extracted_%s" % self.base_pf
+
+ def __getattr__(self, name):
+ # This won't get called if 'name' is found already
+ # and we'd like it to raise AttributeError if it's not anywhere
+ if name in ['h', 'hierarchy']:
+ return StaticOutput._get_hierarchy(self)
+ return getattr(self.base_pf, name)
+
+ def __getitem__(self, key):
+ if key not in self.override:
+ return self.base_pf[key]
+ return self.override[key]
+
Added: trunk/yt/extensions/volume_rendering/CUDARayCast.py
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/volume_rendering/CUDARayCast.py Fri Nov 20 20:27:29 2009
@@ -0,0 +1,168 @@
+"""
+An attempt at putting the ray-casting operation into CUDA
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2009 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import sys;sys.path.insert(0,'.')
+
+from yt.mods import *
+import yt.extensions.HierarchySubset as hs
+import numpy as na
+import h5py, time
+
+import matplotlib;matplotlib.use("Agg");import pylab
+
+from yt.extensions.volume_rendering.TransferFunction import ColorTransferFunction
+
+if __name__ == "__main__":
+
+ # This is boilerplate code for setting up pycuda
+ import pycuda.driver as cuda
+ import pycuda.compiler as compiler
+ import pycuda.autoinit
+ import pycuda.gpuarray as gpuarray
+ cuda.init()
+ assert (cuda.Device.count() >= 1)
+
+ print "Extracting hierarchy."
+ opf = load("/u/ki/mturk/ki05/MSM96-SIM3-restart-J64/DataDump0081.dir/DataDump0081")
+ pf = hs.ExtractedParameterFile(opf, 20)
+
+ cpu = {}
+ gpu = {}
+
+ print "Reading data."
+ #fn = "DataDump0081_partitioned.h5"
+ fn = "RedshiftOutput0005_partitioned.h5"
+ f = h5py.File("/u/ki/mturk/ki05/%s" % fn)
+ cpu['grid_data'] = f["/PGrids/Data"][:].astype("float32")
+ cpu['dims'] = f["/PGrids/Dims"][:].astype("int32") - 1
+ cpu['left_edge'] = f["/PGrids/LeftEdges"][:].astype("float32")
+ cpu['right_edge'] = f["/PGrids/RightEdges"][:].astype("float32")
+
+ print "Constructing transfer function."
+ if "Data" in fn:
+ mh = na.log10(1.67e-24)
+ tf = ColorTransferFunction((7.5+mh, 14.0+mh))
+ tf.add_gaussian( 8.25+mh, 0.002, [0.2, 0.2, 0.4, 0.1])
+ tf.add_gaussian( 9.75+mh, 0.002, [0.0, 0.0, 0.3, 0.1])
+ tf.add_gaussian(10.25+mh, 0.004, [0.0, 0.3, 0.0, 0.1])
+ tf.add_gaussian(11.50+mh, 0.005, [1.0, 0.0, 0.0, 0.7])
+ tf.add_gaussian(12.75+mh, 0.010, [1.0, 1.0, 1.0, 1.0])
+ elif "Red" in fn:
+ tf = ColorTransferFunction((-31, -27))
+ tf.add_gaussian(-30.0, 0.05, [1.0, 0.0, 0.0, 0.1])
+ tf.add_gaussian(-29.5, 0.03, [0.0, 1.0, 0.0, 0.3])
+ tf.add_gaussian(-29.0, 0.05, [0.0, 0.0, 1.0, 0.5])
+ tf.add_gaussian(-28.5, 0.05, [1.0, 1.0, 1.0, 1.0])
+ else: raise RuntimeError
+
+ cpu['ngrids'] = na.array([cpu['dims'].shape[0]], dtype='int32')
+ cpu['tf_r'] = tf.red.y.astype("float32")
+ cpu['tf_g'] = tf.green.y.astype("float32")
+ cpu['tf_b'] = tf.blue.y.astype("float32")
+ cpu['tf_a'] = tf.alpha.y.astype("float32")
+
+ cpu['tf_bounds'] = na.array(tf.x_bounds, dtype='float32')
+
+ cpu['v_dir'] = na.array([0.3, 0.5, 0.6], dtype='float32')
+
+ c = na.array([0.47284317, 0.48062515, 0.58282089], dtype='float32')
+
+ print "Getting cutting plane."
+ cp = pf.h.cutting(cpu['v_dir'], c)
+
+ W = 2000.0/pf['au']
+ W = 0.25
+ Nvec = 128
+ back_c = c - cp._norm_vec * W
+ front_c = c + cp._norm_vec * W
+
+ px, py = na.mgrid[-W:W:Nvec*1j,-W:W:Nvec*1j]
+ xv = cp._inv_mat[0,0]*px + cp._inv_mat[0,1]*py + cp.center[0]
+ yv = cp._inv_mat[1,0]*px + cp._inv_mat[1,1]*py + cp.center[1]
+ zv = cp._inv_mat[2,0]*px + cp._inv_mat[2,1]*py + cp.center[2]
+ cpu['v_pos'] = na.array([xv, yv, zv], dtype='float32').transpose()
+
+ cpu['image_r'] = na.zeros((Nvec, Nvec), dtype='float32').ravel()
+ cpu['image_g'] = na.zeros((Nvec, Nvec), dtype='float32').ravel()
+ cpu['image_b'] = na.zeros((Nvec, Nvec), dtype='float32').ravel()
+ cpu['image_a'] = na.zeros((Nvec, Nvec), dtype='float32').ravel()
+
+ print "Generating module"
+ source = open("yt/extensions/volume_rendering/_cuda_caster.cu").read()
+ mod = compiler.SourceModule(source)
+ func = mod.get_function("ray_cast")
+
+ for n, a in cpu.items():
+ ss = a.size * a.dtype.itemsize
+ print "Allocating %0.3e megabytes for %s" % (ss/(1024*1024.), n)
+ gpu[n] = cuda.to_device(a.ravel('F'))
+ #pycuda.autoinit.context.synchronize()
+
+ BLOCK_SIZE = 8
+ grid_size = Nvec / BLOCK_SIZE
+
+ print "Running ray_cast function."
+ t1 = time.time()
+ ret = func(gpu['ngrids'],
+ gpu['grid_data'],
+ gpu['dims'],
+ gpu['left_edge'],
+ gpu['right_edge'],
+ gpu['tf_r'],
+ gpu['tf_g'],
+ gpu['tf_b'],
+ gpu['tf_a'],
+ gpu['tf_bounds'],
+ gpu['v_dir'],
+ gpu['v_pos'],
+ gpu['image_r'],
+ gpu['image_g'],
+ gpu['image_b'],
+ gpu['image_a'],
+ block=(BLOCK_SIZE,BLOCK_SIZE,1),
+ grid=(grid_size, grid_size), time_kernel=True)
+ t2 = time.time()
+ print "BACK: %0.3e" % (t2-t1)
+
+ mi, ma = 1e300, -1e300
+ image = []
+ for im in 'rgb':
+ ii = 'image_%s' % im
+ sh, dtype = cpu[ii].shape, cpu[ii].dtype
+ del cpu[ii]
+ cpu[ii] = cuda.from_device(gpu[ii], sh, dtype).reshape((Nvec,Nvec))
+ mi, ma = min(cpu[ii].min(),mi), max(cpu[ii].max(), ma)
+ image.append(cpu[ii])
+ print "Min/max of %s %0.3e %0.3e" % (
+ im, image[-1].min(), image[-1].max())
+ pylab.clf()
+ pylab.imshow(image[-1], interpolation='nearest')
+ pylab.savefig("/u/ki/mturk/public_html/vr6/%s.png" % (ii))
+
+ image = na.array(image).transpose()
+ image = (image - mi) / (ma - mi)
+ pylab.clf()
+ pylab.imshow(image, interpolation='nearest')
+ pylab.savefig("/u/ki/mturk/public_html/vr6/image_rgb.png")
Added: trunk/yt/extensions/volume_rendering/_cuda_caster.cu
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/volume_rendering/_cuda_caster.cu Fri Nov 20 20:27:29 2009
@@ -0,0 +1,314 @@
+/***********************************************************************
+An attempt at putting the ray-casting operation into CUDA
+An attempt at putting the ray-casting operation into CUDA
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2009 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+
+//extern __shared__ float array[];
+
+
+#define NUM_SAMPLES 25
+#define VINDEX(A,B,C) tg.data[((((A)+ci[0])*(tg.dims[1]+1)+((B)+ci[1]))*(tg.dims[2]+1)+ci[2]+(C))]
+
+#define fmin(A, B) ( A * (A < B) + B * (B < A) )
+#define fmax(A, B) ( A * (A > B) + B * (B > A) )
+#define fclip(A, B, C) ( fmax( fmin(A, C), B) )
+
+struct transfer_function
+{
+ float vs[4][256];
+ float dbin;
+ float bounds[2];
+};
+
+__shared__ struct transfer_function tf;
+
+struct grid
+{
+ float left_edge[3];
+ float right_edge[3];
+ float dds[3];
+ int dims[3];
+ float *data;
+};
+
+__shared__ struct grid tg;
+
+__device__ float interpolate(float *cache, int *ds, int *ci, float *dp)
+{
+ int i;
+ float dv, dm[3];
+ for(i=0;i<3;i++)dm[i] = (1.0 - dp[i]);
+ dv = 0.0;
+ dv += cache[0] * (dm[0]*dm[1]*dm[2]);
+ dv += cache[1] * (dm[0]*dm[1]*dp[2]);
+ dv += cache[2] * (dm[0]*dp[1]*dm[2]);
+ dv += cache[3] * (dm[0]*dp[1]*dp[2]);
+ dv += cache[4] * (dp[0]*dm[1]*dm[2]);
+ dv += cache[5] * (dp[0]*dm[1]*dp[2]);
+ dv += cache[6] * (dp[0]*dp[1]*dm[2]);
+ dv += cache[7] * (dp[0]*dp[1]*dp[2]);
+ return dv;
+}
+
+__device__ void eval_transfer(float dt, float dv, float *rgba,
+ transfer_function &tf)
+{
+ int i, bin_id;
+ float temp, bv, dy, dd, ta;
+
+ bin_id = (int) ((dv - tf.bounds[0]) / tf.dbin);
+ bv = tf.vs[3][bin_id ];
+ dy = tf.vs[3][bin_id+1] - bv;
+ dd = dv - (tf.bounds[0] + bin_id*tf.dbin);
+ temp = bv+dd*(dy/tf.dbin);
+ ta = temp;
+ for (i = 0; i < 3; i++)
+ {
+ bv = tf.vs[i][bin_id ];
+ dy = tf.vs[i][bin_id+1];
+ dd = dv - (tf.bounds[0] + bin_id*tf.dbin);
+ temp = bv+dd*(dy/tf.dbin);
+ rgba[i] += (1.0 - rgba[3])*ta*temp*dt;
+ }
+ //rgba[3] += (1.0 - rgba[3])*ta*dt;
+}
+
+__device__ void sample_values(float v_pos[3], float v_dir[3],
+ float enter_t, float exit_t, int ci[3], float rgba[4],
+ transfer_function &tf, grid &tg)
+{
+ float cp[3], dp[3], dt, t, dv;
+ int dti, i;
+ float cache[8];
+
+ cache[0] = VINDEX(0,0,0);
+ cache[1] = VINDEX(0,0,1);
+ cache[2] = VINDEX(0,1,0);
+ cache[3] = VINDEX(0,1,1);
+ cache[4] = VINDEX(1,0,0);
+ cache[5] = VINDEX(1,0,1);
+ cache[6] = VINDEX(1,1,0);
+ cache[7] = VINDEX(1,1,1);
+
+ dt = (exit_t - enter_t) / (NUM_SAMPLES-1);
+ for (dti = 0; dti < NUM_SAMPLES - 1; dti++)
+ {
+ t = enter_t + dt*dti;
+ for (i = 0; i < 3; i++)
+ {
+ cp[i] = v_pos[i] + t * v_dir[i];
+ dp[i] = fclip(fmod(cp[i], tg.dds[i])/tg.dds[i], 0.0, 1.0);
+ }
+ dv = interpolate(cache, tg.dims, ci, dp);
+ eval_transfer(dt, dv, rgba, tf);
+ }
+}
+
+
+/* We need to know several things if we want to ray cast through a grid.
+ We need the grid spatial information, as well as its values. We also need
+ the transfer function, which defines what our image will look like. */
+
+__global__ void ray_cast(int ngrids,
+ float *grid_data,
+ int *dims,
+ float *left_edge,
+ float *right_edge,
+ float *tf_r,
+ float *tf_g,
+ float *tf_b,
+ float *tf_a,
+ float *tf_bounds,
+ float *v_dir,
+ float *av_pos,
+ float *image_r,
+ float *image_g,
+ float *image_b,
+ float *image_a)
+{
+
+ int cur_ind[3], step[3], x, y, i, direction;
+ float intersect_t = 1.0, intersect_ts[3];
+ float tmax[3];
+ float tl, tr, temp_xl, temp_yl, temp_xr, temp_yr;
+
+ int offset;
+
+ //transfer_function tf;
+ for (i = 0; i < 4; i++)
+ {
+ x = 4 * (8 * threadIdx.x + threadIdx.y) + i;
+ tf.vs[0][x] = tf_r[x];
+ tf.vs[1][x] = tf_g[x];
+ tf.vs[2][x] = tf_b[x];
+ tf.vs[3][x] = tf_a[x];
+ }
+
+ tf.bounds[0] = tf_bounds[0]; tf.bounds[1] = tf_bounds[1];
+ tf.dbin = (tf.bounds[1] - tf.bounds[0])/255.0;
+
+ /* Set up the grid, just for convenience */
+ //grid tg;
+ int grid_i;
+
+ int tidx = (blockDim.x * gridDim.x) * (
+ blockDim.y * blockIdx.y + threadIdx.y)
+ + (blockDim.x * blockIdx.x + threadIdx.x);
+
+ float rgba[4];
+ //rgba[0] = image_r[tidx];
+ //rgba[1] = image_g[tidx];
+ //rgba[2] = image_b[tidx];
+ //rgba[3] = image_a[tidx];
+
+ float v_pos[3];
+ v_pos[0] = av_pos[tidx + 0];
+ v_pos[1] = av_pos[tidx + 1];
+ v_pos[2] = av_pos[tidx + 2];
+
+ tg.data = grid_data;
+ int skip;
+ for (i = 0; i < 3; i++)
+ {
+ step[i] = 0;
+ step[i] += (v_dir[i] > 0);
+ step[i] += -1 * (v_dir[i] < 0);
+ }
+
+ for(grid_i = 0; grid_i < ngrids; grid_i++) {
+ skip = 0;
+
+ if (threadIdx.x == 0)
+ {
+ if (threadIdx.y == 0) tg.dims[0] = dims[3*grid_i + 0];
+ if (threadIdx.y == 1) tg.dims[1] = dims[3*grid_i + 1];
+ if (threadIdx.y == 2) tg.dims[2] = dims[3*grid_i + 2];
+ }
+
+ if (threadIdx.x == 1)
+ {
+ if (threadIdx.y == 0) tg.left_edge[0] = left_edge[3*grid_i + 0];
+ if (threadIdx.y == 1) tg.left_edge[1] = left_edge[3*grid_i + 1];
+ if (threadIdx.y == 2) tg.left_edge[2] = left_edge[3*grid_i + 2];
+ }
+
+ if (threadIdx.x == 2) {
+ if (threadIdx.y == 0) tg.right_edge[0] = right_edge[3*grid_i + 0];
+ if (threadIdx.y == 1) tg.right_edge[1] = right_edge[3*grid_i + 1];
+ if (threadIdx.y == 2) tg.right_edge[2] = right_edge[3*grid_i + 2];
+ }
+
+ if (threadIdx.x == 3) {
+ if (threadIdx.y == 0) tg.dds[0] = (tg.right_edge[0] - tg.left_edge[0])/tg.dims[0];
+ if (threadIdx.y == 1) tg.dds[1] = (tg.right_edge[1] - tg.left_edge[1])/tg.dims[1];
+ if (threadIdx.y == 2) tg.dds[2] = (tg.right_edge[2] - tg.left_edge[2])/tg.dims[2];
+ }
+
+ /* We integrate our ray */
+
+ for (i = 0; i < 3; i++)
+ {
+ x = (i + 1) % 3;
+ y = (i + 2) % 3;
+
+ tl = (tg.left_edge[i] - v_pos[i])/v_dir[i];
+ temp_xl = (v_pos[i] + tl*v_dir[x]);
+ temp_yr = (v_pos[i] + tl*v_dir[y]);
+
+ tr = (tg.right_edge[i] - v_pos[i])/v_dir[i];
+ temp_xr = (v_pos[x] + tr*v_dir[x]);
+ temp_yr = (v_pos[y] + tr*v_dir[y]);
+
+ intersect_ts[i] = 1.0;
+
+ intersect_ts[i] +=
+ ( (tg.left_edge[x] <= temp_xl) &&
+ (temp_xl <= tg.right_edge[x]) &&
+ (tg.left_edge[y] <= temp_yl) &&
+ (temp_yl <= tg.right_edge[y]) &&
+ (0.0 <= tl) && (tl < intersect_ts[i]) && (tl < tr) ) * tl;
+
+ intersect_ts[i] +=
+ ( (tg.left_edge[x] <= temp_xr) &&
+ (temp_xr <= tg.right_edge[x]) &&
+ (tg.left_edge[y] <= temp_yr) &&
+ (temp_yr <= tg.right_edge[y]) &&
+ (0.0 <= tr) && (tr < intersect_ts[i]) && (tr < tl) ) * tr;
+
+ intersect_t = ( intersect_ts[i] < intersect_t) * intersect_ts[i];
+
+ }
+
+ intersect_t *= (!( (tg.left_edge[0] <= v_pos[0]) &&
+ (v_pos[0] <= tg.right_edge[0]) &&
+ (tg.left_edge[0] <= v_pos[0]) &&
+ (v_pos[0] <= tg.right_edge[0]) &&
+ (tg.left_edge[0] <= v_pos[0]) &&
+ (v_pos[0] <= tg.right_edge[0])));
+
+ skip = ((intersect_t < 0) || (intersect_t > 1.0));
+
+ for (i = 0; i < 3; i++)
+ {
+ cur_ind[i] = (int) floor(((v_pos[i] + intersect_t * v_dir[i]) +
+ step[i]*1e-7*tg.dds[i] -
+ tg.left_edge[i])/tg.dds[i]);
+ tmax[i] = (((cur_ind[i]+step[i])*tg.dds[i])+
+ tg.left_edge[i]-v_pos[i])/v_dir[i];
+ cur_ind[i] -= ((cur_ind[i] == tg.dims[i]) && (step[i] < 0));
+ skip = ((cur_ind[i] < 0) || (cur_ind[i] >= tg.dims[i]));
+ offset = (step[i] > 0);
+ tmax[i] = (((cur_ind[i]+offset)*tg.dds[i])+tg.left_edge[i]-v_pos[i])/v_dir[i];
+ }
+
+ /* This is the primary grid walking loop */
+ while(!( (skip)
+ ||((cur_ind[0] < 0) || (cur_ind[0] >= tg.dims[0])
+ || (cur_ind[1] < 0) || (cur_ind[1] >= tg.dims[1])
+ || (cur_ind[2] < 0) || (cur_ind[2] >= tg.dims[2]))))
+ {
+ direction = 0;
+ direction += 2 * (tmax[0] < tmax[1]) * (tmax[0] >= tmax[2]);
+ direction += 1 * (tmax[0] >= tmax[1]) * (tmax[1] < tmax[2]);
+ direction += 2 * (tmax[0] >= tmax[1]) * (tmax[1] >= tmax[2]);
+ sample_values(v_pos, v_dir, intersect_t, tmax[direction],
+ cur_ind, rgba, tf, tg);
+ cur_ind[direction] += step[direction];
+ intersect_t = tmax[direction];
+ tmax[direction] += abs(tg.dds[direction]/v_dir[direction]);
+ }
+
+ tg.data += (tg.dims[0]+1) * (tg.dims[1]+1) * (tg.dims[2]+1);
+
+ }
+
+ int iy = threadIdx.y + blockDim.y * blockIdx.y;
+ int ix = threadIdx.x + blockDim.x * blockIdx.x;
+ __syncthreads();
+
+ image_r[tidx] = rgba[0];
+ image_g[tidx] = rgba[1];
+ image_b[tidx] = rgba[2];
+ image_a[tidx] = rgba[3];
+}
Modified: trunk/yt/extensions/volume_rendering/grid_partitioner.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/grid_partitioner.py (original)
+++ trunk/yt/extensions/volume_rendering/grid_partitioner.py Fri Nov 20 20:27:29 2009
@@ -27,16 +27,15 @@
from yt.funcs import *
import h5py
-from VolumeIntegrator import PartitionedGrid
+from yt.utils import PartitionedGrid
def partition_grid(start_grid, field, log_field = True, threshold = None):
if threshold is not None:
if start_grid[field].max() < threshold[0] or \
start_grid[field].min() > threshold[1]: return None
- to_cut_up = start_grid.get_vertex_centered_data(field).astype('float64')
+ to_cut_up = start_grid.get_vertex_centered_data(field, smoothed=True).astype('float64')
if log_field: to_cut_up = na.log10(to_cut_up)
- assert(na.any(na.isnan(to_cut_up)) == False)
if len(start_grid.Children) == 0:
pg = PartitionedGrid(
@@ -50,27 +49,27 @@
y_vert = [0, start_grid.ActiveDimensions[1]]
z_vert = [0, start_grid.ActiveDimensions[2]]
+ gi = start_grid.get_global_startindex()
for grid in start_grid.Children:
- gi = start_grid.get_global_startindex()
si = grid.get_global_startindex()/2 - gi
ei = si + grid.ActiveDimensions/2
x_vert += [si[0], ei[0]]
y_vert += [si[1], ei[1]]
z_vert += [si[2], ei[2]]
- cim = start_grid.child_index_mask
-
# Now we sort by our vertices, in axis order
x_vert.sort()
y_vert.sort()
z_vert.sort()
- grids = []
+ return [g for g in _partition(start_grid, to_cut_up, x_vert, y_vert, z_vert)]
- covered = na.zeros(start_grid.ActiveDimensions)
+def _partition(grid, grid_data, x_vert, y_vert, z_vert):
+ grids = []
+ cim = grid.child_index_mask
for xs, xe in zip(x_vert[:-1], x_vert[1:]):
- for ys, ye in zip(y_vert[:-1], y_vert[1:]):
+ for ys, ye in zip(y_vert[:-1:-1], y_vert[1::-1]):
for zs, ze in zip(z_vert[:-1], z_vert[1:]):
sl = (slice(xs, xe), slice(ys, ye), slice(zs, ze))
dd = cim[sl]
@@ -78,30 +77,28 @@
uniq = na.unique(dd)
if uniq.size > 1: continue
if uniq[0] > -1: continue
- data = to_cut_up[xs:xe+1,ys:ye+1,zs:ze+1].copy()
+ data = grid_data[xs:xe+1,ys:ye+1,zs:ze+1].copy()
dims = na.array(dd.shape, dtype='int64')
start_index = na.array([xs,ys,zs], dtype='int64')
- left_edge = start_grid.LeftEdge + start_index * start_grid.dds
- right_edge = left_edge + dims * start_grid.dds
- grids.append(PartitionedGrid(
- data, left_edge, right_edge, dims))
- covered[xs:xe,ys:ye,zs:ze] += 1
- assert(na.all(covered == start_grid.child_mask))
- assert(covered.max() <= 1)
-
- return grids
+ left_edge = grid.LeftEdge + start_index * grid.dds
+ right_edge = left_edge + dims * grid.dds
+ yield PartitionedGrid(
+ data, left_edge, right_edge, dims)
def partition_all_grids(grid_list, field = "Density",
threshold = (-1e300, 1e300), eval_func = None):
new_grids = []
pbar = get_pbar("Partitioning ", len(grid_list))
if eval_func is None: eval_func = lambda a: True
+ dx = 1e300
for i, g in enumerate(grid_list):
if not eval_func(g): continue
pbar.update(i)
+ if g.dds[0] < dx: dx = g.dds[0]
to_add = partition_grid(g, field, True, threshold)
if to_add is not None: new_grids += to_add
pbar.finish()
+ for g in new_grids: g.min_dds = dx
return na.array(new_grids, dtype='object')
def export_partitioned_grids(grid_list, fn):
@@ -110,6 +107,7 @@
nelem = sum((grid.my_data.size for grid in grid_list))
ngrids = len(grid_list)
group = f.create_group("/PGrids")
+ group.attrs["min_dds"] = grid_list[0].min_dds
left_edge = na.concatenate([[grid.LeftEdge,] for grid in grid_list])
f.create_dataset("/PGrids/LeftEdges", data=left_edge); del left_edge
right_edge = na.concatenate([[grid.RightEdge,] for grid in grid_list])
@@ -131,14 +129,117 @@
data = f["/PGrids/Data"][:]
pbar = get_pbar("Reading Grids", dims.shape[0])
curpos = 0
+ dx = f["/PGrids"].attrs["min_dds"]
for i in xrange(dims.shape[0]):
gd = dims[i,:]
gle, gre = left_edges[i,:], right_edges[i,:]
gdata = data[curpos:curpos+gd.prod()].reshape(gd)
# Vertex -> Grid, so we -1 from dims in this
grid_list.append(PartitionedGrid(gdata, gle, gre, gd - 1))
+ grid_list[-1].min_dds = dx
curpos += gd.prod()
pbar.update(i)
pbar.finish()
f.close()
return na.array(grid_list, dtype='object')
+
+class PartitionRegion(object):
+ _count = 0
+ def __init__(self, dims, source_offset, source_vertices, cim_base):
+ self.source_offset = source_offset
+ self.dims = dims
+ cv = []
+ self._cim = cim_base
+ for vertex in source_vertices:
+ if na.any(vertex - source_offset >= 0) or \
+ na.any(vertex - source_offset < dims):
+ cv.append(vertex)
+ self.child_vertices = na.array(cv)
+
+ @property
+ def cim(self):
+ sls = self.source_offset
+ sle = self.source_offset + self.dims
+ return self._cim[sls[0]:sle[0],
+ sls[1]:sle[1],
+ sls[2]:sle[2]]
+
+ def split(self, axis, coord):
+ dims_left = self.dims.copy()
+ dims_left[axis] = coord - self.source_offset[axis]
+ off_left = self.source_offset.copy()
+ left_region = PartitionRegion(dims_left, off_left,
+ self.child_vertices, self._cim)
+ dims_right = self.dims.copy()
+ dims_right[axis] = self.dims[axis] - coord + self.source_offset[axis]
+ off_right = self.source_offset.copy()
+ off_right[axis] = coord
+ right_region = PartitionRegion(dims_right, off_right,
+ self.child_vertices, self._cim)
+ return left_region, right_region
+
+ def find_hyperplane(self, axis):
+ # Our axis is the normal to the hyperplane
+ # Region boundaries is [2][3]
+ # child_vertices is flat 3D array
+ min_balance = 1e30
+ considered = set([self.source_offset[axis]])
+ considered.add(self.source_offset[axis] + self.dims[axis])
+ best_coord = self.source_offset[axis] + self.dims[axis]
+ for v in self.child_vertices:
+ coord = v[axis]
+ sc = coord - self.source_offset[axis]
+ if coord in considered: continue
+ if sc >= self.dims[axis]: continue
+ if sc < 0: continue
+ eff = self.evaluate_hyperplane(axis, coord)
+ if eff < min_balance:
+ min_balance = eff
+ best_coord = coord
+ considered.add(coord)
+ return best_coord
+
+ def evaluate_hyperplane(self, axis, coord):
+ # We check that we're roughly evenly balanced on either side of the grid
+ # Calculate how well balanced it is...
+ vert = self.child_vertices[:,axis]
+ n_left = (vert <= coord).sum()
+ n_right = (vert > coord).sum()
+ eff = abs(0.5 - (n_left / float(vert.shape[0])))
+ return eff
+
+def partition_region(region, axis=0):
+ # region_boundaries is in ints
+ split_coord = region.find_hyperplane(axis)
+ sc = split_coord - region.source_offset[axis]
+ if sc == 0 or sc == region.dims[axis]:
+ rc = na.unique(region.cim)
+ if rc.size > 1 and rc[0] == -1:
+ region._count += 1
+ if region._count > 3:
+ import pdb;pdb.set_trace()
+ return partition_region(region, (axis+1)%3)
+ elif rc.size > 1 and rc[0] > -1:
+ return []
+ left_region, right_region = region.split(axis, split_coord)
+ lrc = na.unique(left_region.cim)
+ rrc = na.unique(right_region.cim)
+ if lrc.size > 1 and lrc[0] == -1:
+ #print axis, split_coord, "Splitting left region", lrc
+ left_region = partition_region(left_region, (axis + 1) % 3)
+ elif lrc.size > 1 and lrc[0] > -1:
+ left_region = []
+ #print axis, split_coord, "Not splitting left region", lrc
+ else:
+ left_region = [left_region]
+
+ if rrc.size > 1 and rrc[0] == -1:
+ #print axis, split_coord, "Splitting right region", rrc
+ right_region = partition_region(right_region, (axis + 1) % 3)
+ elif rrc.size > 1 and rrc[0] > -1:
+ right_region = []
+ #print axis, split_coord, "Not splitting right region", rrc
+ else:
+ right_region = [right_region]
+
+ return left_region + right_region
Modified: trunk/yt/extensions/volume_rendering/software_sampler.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/software_sampler.py (original)
+++ trunk/yt/extensions/volume_rendering/software_sampler.py Fri Nov 20 20:27:29 2009
@@ -82,4 +82,4 @@
pos = g.cast_plane(tfp, vp)
pbar.finish()
- return partitioned_grids, image, vectors, norm_vec, pos
+ return partitioned_grids[ind], image, vectors, norm_vec, pos
Modified: trunk/yt/funcs.py
==============================================================================
--- trunk/yt/funcs.py (original)
+++ trunk/yt/funcs.py Fri Nov 20 20:27:29 2009
@@ -118,7 +118,7 @@
pid = os.getpid()
status_file = "/proc/%s/statm" % (pid)
if not os.path.isfile(status_file):
- return None
+ return 0.0
line = open(status_file).read()
size, resident, share, text, library, data, dt = [int(i) for i in line.split()]
return resident * pagesize / (1024 * 1024) # return in megs
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Fri Nov 20 20:27:29 2009
@@ -203,18 +203,15 @@
This will attempt to convert a given unit to cgs from code units.
It either returns the multiplicative factor or throws a KeyError.
"""
- return self.hierarchy[datatype]
+ return self.pf[datatype]
def clear_data(self):
"""
Clears out all data from the AMRData instance, freeing memory.
"""
- for key in self.data.keys():
- del self.data[key]
- del self.data
+ self.data.clear()
if self._grids is not None:
for grid in self._grids: grid.clear_data()
- self.data = {}
def clear_cache(self):
"""
@@ -238,6 +235,9 @@
self.clear_data()
self.get_data()
+ def keys(self):
+ return self.data.keys()
+
def __getitem__(self, key):
"""
Returns a single field. Will add if necessary.
@@ -320,75 +320,56 @@
grids = [g for g in self._grids if g.Level == level]
return grids
- def __get_levelIndices(self):
- if self.__levelIndices: return self.__levelIndices
- # Otherwise, generate
- # We only have to do this once, so it's not terribly expensive:
- ll = {}
- for level in range(MAXLEVEL):
- t = [i for i in range(len(self._grids)) if self._grids[i].Level == level]
- ll[level] = na.array(t)
- self.__levelIndices = ll
- return self.__levelIndices
-
- def __set_levelIndices(self, val):
- self.__levelIndices = val
-
- def __del_levelIndices(self):
- del self.__levelIndices
- self.__levelIndices = None
-
- __levelIndices = None
- levelIndices = property(__get_levelIndices, __set_levelIndices,
- __del_levelIndices)
-
- def __get_gridLeftEdge(self):
- if self.__gridLeftEdge == None:
- self.__gridLeftEdge = na.array([g.LeftEdge for g in self._grids])
- return self.__gridLeftEdge
-
- def __del_gridLeftEdge(self):
- del self.__gridLeftEdge
- self.__gridLeftEdge = None
-
- def __set_gridLeftEdge(self, val):
- self.__gridLeftEdge = val
-
- __gridLeftEdge = None
- gridLeftEdge = property(__get_gridLeftEdge, __set_gridLeftEdge,
- __del_gridLeftEdge)
-
- def __get_gridRightEdge(self):
- if self.__gridRightEdge == None:
- self.__gridRightEdge = na.array([g.RightEdge for g in self._grids])
- return self.__gridRightEdge
-
- def __del_gridRightEdge(self):
- del self.__gridRightEdge
- self.__gridRightEdge = None
-
- def __set_gridRightEdge(self, val):
- self.__gridRightEdge = val
-
- __gridRightEdge = None
- gridRightEdge = property(__get_gridRightEdge, __set_gridRightEdge,
- __del_gridRightEdge)
-
- def __get_gridLevels(self):
- if self.__gridLevels == None:
- self.__gridLevels = na.array([g.Level for g in self._grids])
- return self.__gridLevels
-
- def __del_gridLevels(self):
- del self.__gridLevels
- self.__gridLevels = None
-
- def __set_gridLevels(self, val):
- self.__gridLevels = val
-
- __gridLevels = None
- gridLevels = property(__get_gridLevels, __set_gridLevels,
- __del_gridLevels)
+ def select_grid_indices(self, level):
+ return na.where(self.grid_levels == level)
+
+ def __get_grid_left_edge(self):
+ if self.__grid_left_edge == None:
+ self.__grid_left_edge = na.array([g.LeftEdge for g in self._grids])
+ return self.__grid_left_edge
+
+ def __del_grid_left_edge(self):
+ del self.__grid_left_edge
+ self.__grid_left_edge = None
+
+ def __set_grid_left_edge(self, val):
+ self.__grid_left_edge = val
+
+ __grid_left_edge = None
+ grid_left_edge = property(__get_grid_left_edge, __set_grid_left_edge,
+ __del_grid_left_edge)
+
+ def __get_grid_right_edge(self):
+ if self.__grid_right_edge == None:
+ self.__grid_right_edge = na.array([g.RightEdge for g in self._grids])
+ return self.__grid_right_edge
+
+ def __del_grid_right_edge(self):
+ del self.__grid_right_edge
+ self.__grid_right_edge = None
+
+ def __set_grid_right_edge(self, val):
+ self.__grid_right_edge = val
+
+ __grid_right_edge = None
+ grid_right_edge = property(__get_grid_right_edge, __set_grid_right_edge,
+ __del_grid_right_edge)
+
+ def __get_grid_levels(self):
+ if self.__grid_levels == None:
+ self.__grid_levels = na.array([g.Level for g in self._grids])
+ return self.__grid_levels
+
+ def __del_grid_levels(self):
+ del self.__grid_levels
+ self.__grid_levels = None
+
+ def __set_grid_levels(self, val):
+ self.__grid_levels = val
+
+ __grid_levels = None
+ grid_levels = property(__get_grid_levels, __set_grid_levels,
+ __del_grid_levels)
class AMR1DData(AMRData, GridPropertiesMixin):
@@ -468,10 +449,10 @@
def _get_list_of_grids(self):
# This bugs me, but we will give the tie to the LeftEdge
- y = na.where( (self.px >= self.pf.hierarchy.gridLeftEdge[:,self.px_ax])
- & (self.px < self.pf.hierarchy.gridRightEdge[:,self.px_ax])
- & (self.py >= self.pf.hierarchy.gridLeftEdge[:,self.py_ax])
- & (self.py < self.pf.hierarchy.gridRightEdge[:,self.py_ax]))
+ y = na.where( (self.px >= self.pf.hierarchy.grid_left_edge[:,self.px_ax])
+ & (self.px < self.pf.hierarchy.grid_right_edge[:,self.px_ax])
+ & (self.py >= self.pf.hierarchy.grid_left_edge[:,self.py_ax])
+ & (self.py < self.pf.hierarchy.grid_right_edge[:,self.py_ax]))
self._grids = self.hierarchy.grids[y]
def _get_data_from_grid(self, grid, field):
@@ -516,8 +497,8 @@
def _get_list_of_grids(self):
# Get the value of the line at each LeftEdge and RightEdge
- LE = self.pf.h.gridLeftEdge
- RE = self.pf.h.gridRightEdge
+ LE = self.pf.h.grid_left_edge
+ RE = self.pf.h.grid_right_edge
p = na.zeros(self.pf.h.num_grids, dtype='bool')
# Check left faces first
for i in range(3):
@@ -740,7 +721,7 @@
# Here we assume that the grid is the max level
level = self.hierarchy.max_level
self.coord
- dx = self.hierarchy.gridDxs[self.hierarchy.levelIndices[level][0]]
+ dx = self.hierarchy.select_grids(level)[0].dds[self.axis]
self.coord += dx * val
else:
raise ValueError(val)
@@ -768,8 +749,8 @@
self.ActiveDimensions = (t.shape[0], 1, 1)
def _get_list_of_grids(self):
- goodI = ((self.pf.h.gridRightEdge[:,self.axis] > self.coord)
- & (self.pf.h.gridLeftEdge[:,self.axis] <= self.coord ))
+ goodI = ((self.pf.h.grid_right_edge[:,self.axis] > self.coord)
+ & (self.pf.h.grid_left_edge[:,self.axis] <= self.coord ))
self._grids = self.pf.h.grids[goodI] # Using sources not hierarchy
def __cut_mask_child_mask(self, grid):
@@ -819,7 +800,7 @@
conv_factor = 1.0
if self.pf.field_info.has_key(field):
conv_factor = self.pf.field_info[field]._convert_function(self)
- dv = self._read_data_slice(grid, field, self.axis, wantedIndex) * conv_factor
+ dv = self.hierarchy.io._read_data_slice(grid, field, self.axis, wantedIndex) * conv_factor
else:
dv = grid[field]
if dv.size == 1: dv = na.ones(grid.ActiveDimensions)*dv
@@ -893,8 +874,8 @@
# onto the normal vector of a plane is:
# D = (a x_0 + b y_0 + c z_0 + d)/sqrt(a^2+b^2+c^2)
# @todo: Convert to using corners
- LE = self.pf.h.gridLeftEdge
- RE = self.pf.h.gridRightEdge
+ LE = self.pf.h.grid_left_edge
+ RE = self.pf.h.grid_right_edge
vertices = na.array([[LE[:,0],LE[:,1],LE[:,2]],
[RE[:,0],RE[:,1],RE[:,2]],
[LE[:,0],LE[:,1],RE[:,2]],
@@ -1168,8 +1149,8 @@
def _get_point_indices(self, grid):
if self._pixelmask.max() == 0: return []
k = amr_utils.PointsInVolume(self._coord, self._pixelmask,
- grid.LeftEdge, grid.RightEdge,
- grid.child_mask, just_one(grid['dx']))
+ grid.LeftEdge, grid.RightEdge,
+ grid.child_mask, just_one(grid['dx']))
return k
def _gen_node_name(self):
@@ -1187,7 +1168,7 @@
def __init__(self, axis, field, weight_field = None,
max_level = None, center = None, pf = None,
source=None, node_name = None, field_cuts = None,
- serialize=True,**kwargs):
+ preload_style='level', 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
@@ -1203,11 +1184,12 @@
self._initialize_source(source)
self._grids = self.source._grids
if max_level == None:
- max_level = self.hierarchy.maxLevel
+ max_level = self.hierarchy.max_level
if self.source is not None:
- max_level = min(max_level, self.source.gridLevels.max())
+ max_level = min(max_level, self.source.grid_levels.max())
self._max_level = max_level
self._weight = weight_field
+ self.preload_style = preload_style
self.func = na.sum # for the future
self.__retval_coords = {}
self.__retval_fields = {}
@@ -1241,31 +1223,15 @@
self._okay_to_serialize = True
#@time_execution
- def __cache_data(self):
- rdf = self.hierarchy.grid.readDataFast
- self.hierarchy.grid.readDataFast = readDataPackedHandle
- for fn,g_list in self.hierarchy.cpu_map.items():
- to_read = na.intersect1d(g_list, self.source._grids)
- if len(to_read) == 0: continue
- fh = h5py.File(to_read[0].filename,'r')
- for g in to_read:
- g.handle = fh
- for field in ensure_list(self.fields):
- g[field]
- del g.handle
- fh.close()
- self.hierarchy.grid.readDataFast = readDataPackedHandle
-
- #@time_execution
def __calculate_overlap(self, level):
s = self.source
mylog.info("Generating overlap masks for level %s", level)
i = 0
pbar = get_pbar("Reading and masking grids ", len(s._grids))
mylog.debug("Examining level %s", level)
- grids = s.levelIndices[level]
- RE = s.gridRightEdge[grids]
- LE = s.gridLeftEdge[grids]
+ grids = s.select_grid_indices(level)
+ RE = s.grid_right_edge[grids]
+ LE = s.grid_left_edge[grids]
for grid in s._grids[grids]:
pbar.update(i)
self.__overlap_masks[grid.id] = \
@@ -1325,7 +1291,7 @@
def __combine_grids_on_level(self, level):
grids = self.source.select_grids(level)
- grids_i = self.source.levelIndices[level]
+ grids_i = self.source.select_grid_indices(level)
pbar = get_pbar('Combining level % 2i / % 2i ' \
% (level, self._max_level), len(grids))
# We have an N^2 check, so we try to be as quick as possible
@@ -1352,7 +1318,7 @@
def __refine_to_level(self, level):
grids = self.source.select_grids(level)
- grids_up = self.source.levelIndices[level-1]
+ grids_up = self.source.select_grid_indices(level - 1)
pbar = get_pbar('Refining to level % 2i / % 2i ' \
% (level, self._max_level), len(grids))
for pi, grid1 in enumerate(grids):
@@ -1365,7 +1331,7 @@
args += self.__retval_coords[grid2.id] + [self.__retval_fields[grid2.id]]
args += self.__retval_coords[grid1.id] + [self.__retval_fields[grid1.id]]
# Refinement factor, which is same in all directions
- args.append(int(grid2['dx'] / grid1['dx']))
+ args.append(int(grid2.dds[0] / grid1.dds[0]))
args.append(na.ones(args[0].shape, dtype='int64'))
kk = PointCombine.CombineGrids(*args)
goodI = args[-1].astype('bool')
@@ -1394,11 +1360,15 @@
# We do this here, but I am not convinced it should be done here
# It is probably faster, as it consolidates IO, but if we did it in
# _project_level, then it would be more memory conservative
+ if self.preload_style == 'all':
+ print "Preloading %s grids and getting %s" % (
+ len(self.source._grids), self._get_dependencies(fields))
+ self._preload(self.source._grids,
+ self._get_dependencies(fields), self.hierarchy.io)
for level in range(0, self._max_level+1):
- level_ind = self.source.levelIndices[level]
- if len(level_ind) == 0: continue
- self._preload(self.source._grids[level_ind],
- self._get_dependencies(fields), self.hierarchy.queue)
+ if self.preload_style == 'level':
+ self._preload(self.source.select_grids(level),
+ self._get_dependencies(fields), self.hierarchy.io)
self.__calculate_overlap(level)
my_coords, my_dx, my_fields = self.__project_level(level, fields)
coord_data.append(my_coords)
@@ -1409,10 +1379,11 @@
if len(check) > 0: all_data.append(check)
# Now, we should clean up after ourselves...
for grid in self.source.select_grids(level - 1):
- grid.clear_data()
del self.__retval_coords[grid.id]
del self.__retval_fields[grid.id]
del self.__overlap_masks[grid.id]
+ mylog.debug("End of projecting level level %s, memory usage %0.3e",
+ level, get_memory_usage()/1024.)
coord_data = na.concatenate(coord_data, axis=1)
field_data = na.concatenate(field_data, axis=1)
dxs = na.concatenate(dxs, axis=1)
@@ -1433,6 +1404,7 @@
self[field] = field_data[fi].ravel()
if self.serialize: self._store_fields(field, self._node_name)
for i in data.keys(): self[i] = data.pop(i)
+ mylog.info("Projection completed")
def add_fields(self, fields, weight = "CellMassMsun"):
pass
@@ -1523,8 +1495,8 @@
else:
grids,ind = self.pf.hierarchy.get_box_grids(
self.left_edge, self.right_edge)
- level_ind = (self.pf.hierarchy.gridLevels.ravel()[ind] <= self.level)
- sort_ind = na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind])
+ level_ind = (self.pf.hierarchy.grid_levels.ravel()[ind] <= self.level)
+ sort_ind = na.argsort(self.pf.h.grid_levels.ravel()[ind][level_ind])
self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)][::-1]
def _generate_coords(self):
@@ -1802,6 +1774,14 @@
grid[field] = na.ones(grid.ActiveDimensions)*default_value
grid[field][self._get_point_indices(grid)] = value
+ _particle_handler = None
+
+ @property
+ def particles(self):
+ if self._particle_handler is None:
+ self._particle_handler = \
+ particle_handler_registry[self._type_name](self.pf, self)
+ return self._particle_handler
class ExtractedRegionBase(AMR3DData):
"""
@@ -1965,9 +1945,9 @@
self._refresh_data()
def _get_list_of_grids(self):
- H = na.sum(self._norm_vec.reshape((1,3,1)) * self.pf.h.gridCorners,
+ H = na.sum(self._norm_vec.reshape((1,3,1)) * self.pf.h.grid_corners,
axis=1) + self._d
- D = na.sqrt(na.sum((self.pf.h.gridCorners -
+ D = na.sqrt(na.sum((self.pf.h.grid_corners -
self.center.reshape((1,3,1)))**2.0,axis=1))
R = na.sqrt(D**2.0-H**2.0)
self._grids = self.hierarchy.grids[
@@ -2237,8 +2217,8 @@
else:
grids,ind = self.pf.hierarchy.get_box_grids(
self.left_edge, self.right_edge)
- level_ind = na.where(self.pf.hierarchy.gridLevels.ravel()[ind] <= self.level)
- sort_ind = na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind])
+ level_ind = na.where(self.pf.hierarchy.grid_levels.ravel()[ind] <= self.level)
+ sort_ind = na.argsort(self.pf.h.grid_levels.ravel()[ind][level_ind])
self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)][::-1]
def extract_region(self, indices):
@@ -2365,8 +2345,8 @@
grids,ind = self.pf.hierarchy.get_box_grids(
self.left_edge - self.dds,
self.right_edge + self.dds)
- level_ind = na.where(self.pf.hierarchy.gridLevels.ravel()[ind] <= self.level)
- sort_ind = na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind])
+ level_ind = na.where(self.pf.hierarchy.grid_levels.ravel()[ind] <= self.level)
+ sort_ind = na.argsort(self.pf.h.grid_levels.ravel()[ind][level_ind])
self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)]
def _get_level_array(self, level, fields):
@@ -2419,7 +2399,7 @@
# How do we find out the root grid base dx?
idims = na.array([3,3,3])
dx = na.minimum((self.right_edge-self.left_edge)/(idims-2),
- self.pf.h.grids[0]['dx'])
+ self.pf.h.grids[0].dds[0])
idims = na.floor((self.right_edge-self.left_edge)/dx) + 2
for ax in 'xyz': self['cd%s'%ax] = dx[0]
self[field] = na.zeros(idims,dtype='float64')-999
@@ -2479,8 +2459,8 @@
grids,ind = self.pf.hierarchy.get_box_grids(
self.left_edge - buffer,
self.right_edge + buffer)
- level_ind = (self.pf.hierarchy.gridLevels.ravel()[ind] <= self.level)
- sort_ind = na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind])
+ level_ind = (self.pf.hierarchy.grid_levels.ravel()[ind] <= self.level)
+ sort_ind = na.argsort(self.pf.h.grid_levels.ravel()[ind][level_ind])
self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)][::-1]
def _refresh_data(self):
Modified: trunk/yt/lagos/BaseGridType.py
==============================================================================
--- trunk/yt/lagos/BaseGridType.py (original)
+++ trunk/yt/lagos/BaseGridType.py Fri Nov 20 20:27:29 2009
@@ -27,7 +27,7 @@
#import yt.enki, gc
from yt.funcs import *
-class AMRGridPatch(AMRData):
+class AMRGridPatch(object):
_spatial = True
_num_ghost_zones = 0
_grids = None
@@ -37,18 +37,74 @@
_skip_add = True
_con_args = ('id', 'filename')
- filename = None
-
- def __init__(self, id, filename=None, hierarchy = None):
+ __slots__ = ['data', 'field_parameters', 'id', 'hierarchy', 'pf',
+ 'ActiveDimensions', 'LeftEdge', 'RightEdge', 'Level',
+ 'NumberOfParticles', 'Children', 'Parent',
+ 'start_index', 'filename', '__weakref__', 'dds',
+ '_child_mask', '_child_indices', '_child_index_mask',
+ '_parent_id', '_children_ids']
+ def __init__(self, id, filename = None, hierarchy = 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.pf = self.hierarchy.parameter_file # weakref already
+ self._child_mask = self._child_indices = self._child_index_mask = None
+ self.start_index = None
+
+ 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 / self.dds
+ return na.rint(start_index).astype('int64').ravel()
+ pdx = self.Parent.dds
+ start_index = (self.Parent.get_global_startindex()) + \
+ na.rint((self.LeftEdge - self.Parent.LeftEdge)/pdx)
+ self.start_index = (start_index*self.pf["RefineBy"]).astype('int64').ravel()
+ return self.start_index
+
+
+ def get_field_parameter(self, name, default=None):
+ """
+ This is typically only used by derived field functions, but
+ it returns parameters used to generate fields.
+ """
+ if self.field_parameters.has_key(name):
+ return self.field_parameters[name]
+ else:
+ return default
+
+ def set_field_parameter(self, name, val):
+ """
+ Here we set up dictionaries that get passed up and down and ultimately
+ to derived fields.
+ """
+ self.field_parameters[name] = val
+
+ def has_field_parameter(self, name):
+ """
+ Checks if a field parameter is set.
+ """
+ return self.field_parameters.has_key(name)
+
+ def convert(self, datatype):
+ """
+ This will attempt to convert a given unit to cgs from code units.
+ It either returns the multiplicative factor or throws a KeyError.
+ """
+ return self.pf[datatype]
+
+ def __repr__(self):
+ # We'll do this the slow way to be clear what's going on
+ s = "%s (%s): " % (self.__class__.__name__, self.pf)
+ s += ", ".join(["%s=%s" % (i, getattr(self,i))
+ for i in self._con_args])
+ return s
def _generate_field(self, field):
if self.pf.field_info.has_key(field):
@@ -59,7 +115,7 @@
# This is only going to be raised if n_gz > 0
n_gz = ngt_exception.ghost_zones
f_gz = ngt_exception.fields
- gz_grid = self.retrieve_ghost_zones(n_gz, f_gz, smoothed=True)
+ gz_grid = self.retrieve_ghost_zones(n_gz, f_gz, smoothed=False)
temp_array = self.pf.field_info[field](gz_grid)
sl = [slice(n_gz,-n_gz)] * 3
self[field] = temp_array[sl]
@@ -67,6 +123,32 @@
self[field] = self.pf.field_info[field](self)
else: # Can't find the field, try as it might
raise exceptions.KeyError, field
+
+ def has_key(self, key):
+ return (key in self.data)
+
+ def __getitem__(self, key):
+ """
+ Returns a single field. Will add if necessary.
+ """
+ if not self.data.has_key(key):
+ self.get_data(key)
+ return self.data[key]
+
+ def __setitem__(self, key, val):
+ """
+ Sets a field to be some other value.
+ """
+ self.data[key] = val
+
+ def __delitem__(self, key):
+ """
+ Deletes a field
+ """
+ del self.data[key]
+
+ def keys(self):
+ return self.data.keys()
def get_data(self, field):
"""
@@ -83,12 +165,9 @@
self[field] = na.array([],dtype='int64')
return self.data[field]
try:
- if hasattr(self.hierarchy, 'queue'):
- temp = self.hierarchy.queue.pop(self, field)
- else:
- temp = self.readDataFast(field)
- self[field] = temp * conv_factor
- except self._read_exception, exc:
+ temp = self.hierarchy.io.pop(self, field)
+ self[field] = na.multiply(temp, conv_factor, temp)
+ except self.hierarchy.io._read_exception, exc:
if field in self.pf.field_info:
if self.pf.field_info[field].not_in_all:
self[field] = na.zeros(self.ActiveDimensions, dtype='float64')
@@ -103,14 +182,26 @@
# So first we figure out what the index is. We don't assume
# that dx=dy=dz , at least here. We probably do elsewhere.
id = self.id - self._id_offset
- self.dds = na.array([self.hierarchy.gridDxs[id,0],
- self.hierarchy.gridDys[id,0],
- self.hierarchy.gridDzs[id,0]])
+ if self.Parent is not None:
+ self.dds = self.Parent.dds / self.pf["RefineBy"]
+ else:
+ LE, RE = self.hierarchy.grid_left_edge[id,:], \
+ self.hierarchy.grid_right_edge[id,:]
+ self.dds = na.array((RE-LE)/self.ActiveDimensions)
self.data['dx'], self.data['dy'], self.data['dz'] = self.dds
@property
def _corners(self):
- return self.hierarchy.gridCorners[:,:,self.id - self._id_offset]
+ return 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')
def _generate_overlap_masks(self, axis, LE, RE):
"""
@@ -132,24 +223,6 @@
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,
@@ -157,13 +230,12 @@
"""
self._del_child_mask()
self._del_child_indices()
- if hasattr(self, 'coarseData'):
- del self.coarseData
- if hasattr(self, 'retVal'):
- del self.retVal
- AMRData.clear_data(self)
+ self.data.clear()
self._setup_dx()
+ def check_child_masks(self):
+ return self._child_mask, self._child_indices
+
def _prepare_grid(self):
"""
Copies all the appropriate attributes from the hierarchy
@@ -173,20 +245,13 @@
# Note that to keep in line with Enzo, we have broken PEP-8
h = self.hierarchy # cache it
my_ind = self.id - self._id_offset
- self.ActiveDimensions = h.gridEndIndices[my_ind] \
- - h.gridStartIndices[my_ind] + 1
- self.LeftEdge = h.gridLeftEdge[my_ind]
- self.RightEdge = h.gridRightEdge[my_ind]
- self.Level = h.gridLevels[my_ind,0]
+ self.ActiveDimensions = h.grid_dimensions[my_ind]
+ self.LeftEdge = h.grid_left_edge[my_ind]
+ self.RightEdge = h.grid_right_edge[my_ind]
+ h.grid_levels[my_ind, 0] = self.Level
# This might be needed for streaming formats
#self.Time = h.gridTimes[my_ind,0]
- self.NumberOfParticles = h.gridNumberOfParticles[my_ind,0]
- self.Children = h.gridTree[my_ind]
- pID = h.gridReverseTree[my_ind]
- if pID != None and pID != -1:
- self.Parent = weakref.proxy(h.grids[pID - self._id_offset])
- else:
- self.Parent = None
+ self.NumberOfParticles = h.grid_particle_count[my_ind,0]
def __len__(self):
return na.prod(self.ActiveDimensions)
@@ -238,59 +303,58 @@
del self.child_ind
def _set_child_mask(self, newCM):
- if self.__child_mask != None:
+ if self._child_mask != None:
mylog.warning("Overriding child_mask attribute! This is probably unwise!")
- self.__child_mask = newCM
+ self._child_mask = newCM
def _set_child_indices(self, newCI):
- if self.__child_indices != None:
+ if self._child_indices != None:
mylog.warning("Overriding child_indices attribute! This is probably unwise!")
- self.__child_indices = newCI
+ self._child_indices = newCI
def _get_child_mask(self):
- if self.__child_mask == None:
+ if self._child_mask == None:
self.__generate_child_mask()
- return self.__child_mask
+ return self._child_mask
def _get_child_indices(self):
- if self.__child_indices == None:
+ if self._child_indices == None:
self.__generate_child_mask()
- return self.__child_indices
+ return self._child_indices
def _del_child_indices(self):
try:
- del self.__child_indices
+ del self._child_indices
except AttributeError:
pass
- self.__child_indices = None
+ self._child_indices = None
def _del_child_mask(self):
try:
- del self.__child_mask
+ del self._child_mask
except AttributeError:
pass
- self.__child_mask = None
+ self._child_mask = None
def _get_child_index_mask(self):
- if self.__child_index_mask is None:
+ if self._child_index_mask is None:
self.__generate_child_index_mask()
- return self.__child_index_mask
+ return self._child_index_mask
def _del_child_index_mask(self):
try:
- del self.__child_index_mask
+ del self._child_index_mask
except AttributeError:
pass
- self.__child_index_mask = None
+ 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.dds))
- endIndex = na.minimum(na.rint(
- (child.RightEdge - self.LeftEdge)/self.dds),
+ rf = self.pf["RefineBy"]
+ gi, cgi = self.get_global_startindex(), child.get_global_startindex()
+ startIndex = na.maximum(0, cgi/rf - gi)
+ endIndex = na.minimum( (cgi+child.ActiveDimensions)/rf - gi,
self.ActiveDimensions)
- startIndex = na.maximum(0, startIndex)
mask[startIndex[0]:endIndex[0],
startIndex[1]:endIndex[1],
startIndex[2]:endIndex[2]] = tofill
@@ -300,19 +364,19 @@
Generates self.child_mask, which is zero where child grids exist (and
thus, where higher resolution data is available.)
"""
- self.__child_mask = na.ones(self.ActiveDimensions, 'int32')
+ self._child_mask = na.ones(self.ActiveDimensions, 'int32')
for child in self.Children:
- self.__fill_child_mask(child, self.__child_mask, 0)
- self.__child_indices = (self.__child_mask==0) # bool, possibly redundant
+ 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):
"""
Generates self.child_index_mask, which is -1 where there is no child,
and otherwise has the ID of the grid that resides there.
"""
- self.__child_index_mask = na.zeros(self.ActiveDimensions, 'int32') - 1
+ self._child_index_mask = na.zeros(self.ActiveDimensions, 'int32') - 1
for child in self.Children:
- self.__fill_child_mask(child, self.__child_index_mask,
+ self.__fill_child_mask(child, self._child_index_mask,
child.id)
def _get_coords(self):
@@ -337,10 +401,6 @@
LE = na.reshape(self.LeftEdge,(3,1,1,1))
self['x'], self['y'], self['z'] = (ind+0.5)*self.dds+LE
- __child_mask = None
- __child_indices = None
- __child_index_mask = None
-
child_mask = property(fget=_get_child_mask, fdel=_del_child_mask)
child_index_mask = property(fget=_get_child_index_mask, fdel=_del_child_index_mask)
child_indices = property(fget=_get_child_indices, fdel = _del_child_indices)
@@ -382,38 +442,22 @@
scalars = 10**interp(dict(x=x,y=y,z=z))
return scalars
- def _save_data_state(self):
- self.__current_data_keys = self.data.keys()
- if self.__child_mask != None:
- self.__current_child_mask == True
- else:
- self.__current_child_mask = False
-
- if self.__child_indices != None:
- self.__current_child_indices == True
- else:
- self.__current_child_indices = False
-
- def _restore_data_state(self):
- if not self.__current_child_mask:
- self._del_child_mask()
- if not self.__current_child_indices:
- self._del_child_indices()
- for key in data.keys():
- if key not in self.__current_data_keys:
- del self.data[key]
-
-class EnzoGridBase(AMRGridPatch):
+class EnzoGrid(AMRGridPatch):
"""
Class representing a single Enzo Grid instance.
"""
- def __init__(self, id, filename=None, hierarchy = None):
+
+ __slots__ = []
+ def __init__(self, id, hierarchy):
"""
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)
+ AMRGridPatch.__init__(self, id, filename = None, hierarchy = hierarchy)
+ self._children_ids = []
+ self._parent_id = -1
+ self.Level = -1
def _guess_properties_from_parent(self):
"""
@@ -433,42 +477,13 @@
self.start_index = rf*(ParentLeftIndex + self.Parent.get_global_startindex()).astype('int64')
self.LeftEdge = self.Parent.LeftEdge + self.Parent.dds * ParentLeftIndex
self.RightEdge = self.LeftEdge + self.ActiveDimensions*self.dds
- self.hierarchy.gridDxs[my_ind,0] = self['dx']
- self.hierarchy.gridDys[my_ind,0] = self['dy']
- self.hierarchy.gridDzs[my_ind,0] = self['dz']
- self.hierarchy.gridLeftEdge[my_ind,:] = self.LeftEdge
- self.hierarchy.gridRightEdge[my_ind,:] = self.RightEdge
- self.hierarchy.gridCorners[:,:,my_ind] = 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.hierarchy.grid_left_edge[my_ind,:] = self.LeftEdge
+ self.hierarchy.grid_right_edge[my_ind,:] = self.RightEdge
+ 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 / self.dds
- return na.rint(start_index).astype('int64').ravel()
- pdx = self.Parent.dds
- start_index = (self.Parent.get_global_startindex()) + \
- na.rint((self.LeftEdge - self.Parent.LeftEdge)/pdx)
- self.start_index = (start_index*self.pf["RefineBy"]).astype('int64').ravel()
- return self.start_index
-
def set_filename(self, filename):
"""
Intelligently set the filename.
@@ -485,7 +500,22 @@
def __repr__(self):
return "EnzoGrid_%04i" % (self.id)
-class OrionGridBase(AMRGridPatch):
+ @property
+ def Parent(self):
+ if self._parent_id == -1: return None
+ return self.hierarchy.grids[self._parent_id - self._id_offset]
+
+ @property
+ def Children(self):
+ return [self.hierarchy.grids[cid - self._id_offset]
+ for cid in self._children_ids]
+
+class EnzoGridInMemory(EnzoGrid):
+ __slots__ = ['proc_num']
+ def set_filename(self, filename):
+ pass
+
+class OrionGrid(AMRGridPatch):
_id_offset = 0
def __init__(self, LeftEdge, RightEdge, index, level, filename, offset, dimensions,start,stop,paranoia=False):
AMRGridPatch.__init__(self, index)
@@ -516,8 +546,8 @@
self.StartIndices = h.gridStartIndices[self.id]
self.EndIndices = h.gridEndIndices[self.id]
h.gridLevels[self.id,0] = self.Level
- h.gridLeftEdge[self.id,:] = self.LeftEdge[:]
- h.gridRightEdge[self.id,:] = self.RightEdge[:]
+ h.grid_left_edge[self.id,:] = self.LeftEdge[:]
+ h.grid_right_edge[self.id,:] = self.RightEdge[:]
self.Time = h.gridTimes[self.id,0]
self.NumberOfParticles = h.gridNumberOfParticles[self.id,0]
self.Children = h.gridTree[self.id]
Modified: trunk/yt/lagos/DataReadingFuncs.py
==============================================================================
--- trunk/yt/lagos/DataReadingFuncs.py (original)
+++ trunk/yt/lagos/DataReadingFuncs.py Fri Nov 20 20:27:29 2009
@@ -25,186 +25,20 @@
from yt.lagos import *
import exceptions
-def getFieldsHDF4(self):
- """
- Returns a list of fields associated with the filename
- Should *only* be called as EnzoGridInstance.getFields, never as getFields(object)
- """
- return SD.SD(self.filename).datasets().keys()
-
-def getFieldsHDF5(self):
- """
- Returns a list of fields associated with the filename
- Should *only* be called as EnzoGridInstance.getFields, never as getFields(object)
- """
- return HDF5LightReader.ReadListOfDatasets(self.filename, "/")
-
-def readDataHDF4(self, field):
- """
- Returns after having obtained or generated a field. Should throw an
- exception. Should only be called as EnzoGridInstance.readData()
-
- @param field: field to read
- @type field: string
- """
- t = SD.SD(self.filename).select(field).get()
- return t.swapaxes(0,2)
-
-def readAllDataHDF4(self):
- """
- Reads all fields inside an HDF4 file. Should only be called as
- EnzoGridInstance.readAllData() .
- """
- sets = SD.SD(self.filename).datasets()
- for set in sets:
- self[set] = self.readDataFast(set)
-
-def readDataHDF5(self, field):
- return HDF5LightReader.ReadData(self.filename, "/%s" % field).swapaxes(0,2)
-
-def readAllDataHDF5(self):
- """
- Not implemented. Fix me!
- """
- pass
-
-def readAllDataPacked(self):
- """
- Not implemented. Fix me!
- """
- pass
-
-def readDataSliceHDF5(self, grid, field, axis, coord):
- """
- Reads a slice through the HDF5 data
-
- @param grid: Grid to slice
- @type grid: L{EnzoGrid<EnzoGrid>}
- @param field: field to get
- @type field: string
- @param axis: axis to slice along
- @param coord: coord to slice at
- """
- axis = {0:2,1:1,2:0}[axis]
- t = HDF5LightReader.ReadDataSlice(grid.filename, "/%s" %
- (field), axis, coord).transpose()
- return t
-
-def readDataSliceHDF4(self, grid, field, axis, coord):
- """
- Reads a slice through the HDF4 data
-
- @param grid: Grid to slice
- @type grid: L{EnzoGrid<EnzoGrid>}
- @param field: field to get
- @type field: string
- @param sl: region to get
- @type sl: SliceType
- """
- sl = [slice(None), slice(None), slice(None)]
- sl[axis] = slice(coord, coord + 1)
- sl = tuple(reversed(sl))
- return SD.SD(grid.filename).select(field)[sl].swapaxes(0,2)
-
-def readDataPackedHandle(self, field):
- t = self.handle.getNode("/Grid%08i" % (self.id), field).read().astype('float64')
- t = t.swapaxes(0,2)
- return t
-
-def readDataPacked(self, field, dspace=None):
- if dspace is not None:
- dspace="/Grid%08i/%s" % (self.id, dspace)
- return HDF5LightReader.ReadData(self.filename,
- "/Grid%08i/%s" % (self.id, field), dspace).swapaxes(0,2)
- return HDF5LightReader.ReadData(self.filename,
- "/Grid%08i/%s" % (self.id, field)).swapaxes(0,2)
-
-def readDataSlicePacked(self, grid, field, axis, coord):
- """
- Reads a slice through the HDF5 data
-
- @param grid: Grid to slice
- @type grid: L{EnzoGrid<EnzoGrid>}
- @param field: field to get
- @type field: string
- @param sl: region to get
- @type sl: SliceType
- """
- axis = {0:2,1:1,2:0}[axis]
- t = HDF5LightReader.ReadDataSlice(grid.filename, "/Grid%08i/%s" %
- (grid.id, field), axis, coord).transpose()
- return t
-
-def getFieldsPacked(self):
- """
- Returns a list of fields associated with the filename
- Should *only* be called as EnzoGridInstance.getFields, never as getFields(object)
- """
- return HDF5LightReader.ReadListOfDatasets(self.filename, "/Grid%08i" % self.id)
-
-def getExceptionHDF4():
- return SD.HDF4Error
-
-def getExceptionHDF5():
- return (exceptions.KeyError, HDF5LightReader.ReadingError)
-
-def readDataInMemory(self, field):
- import enzo
- return enzo.grid_data[self.id][field].swapaxes(0,2)
-
-def readAllDataInMemory(self):
- pass
-
-def getFieldsInMemory(self):
- import enzo
- return enzo.grid_data[self.id].keys()
-
-def readDataSliceInMemory(self, grid, field, axis, coord):
- import enzo
- 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
-
-def readDataSliceHDF4_2D(self, grid, field, axis, coord):
- t = SD.SD(grid.filename).select(field).get()
- return t.transpose()
-
-
-def readDataSlicePacked2D(self, grid, field, axis, coord):
- """
- Reads a slice through the HDF5 data
+_axis_ids = {0:2,1:1,2:0}
- @param grid: Grid to slice
- @type grid: L{EnzoGrid<EnzoGrid>}
- @param field: field to get
- @type field: string
- @param sl: region to get
- @type sl: SliceType
- """
- t = HDF5LightReader.ReadData(grid.filename, "/Grid%08i/%s" %
- (grid.id, field)).transpose()
- return t
+io_registry = {}
-def readDataSlicePacked1D(self, grid, field, axis, coord):
- """
- Reads a slice through the HDF5 data
+class BaseIOHandler(object):
- @param grid: Grid to slice
- @type grid: L{EnzoGrid<EnzoGrid>}
- @param field: field to get
- @type field: string
- @param sl: region to get
- @type sl: SliceType
- """
- t = HDF5LightReader.ReadData(grid.filename, "/Grid%08i/%s" %
- (grid.id, field))
- return t
+ _data_style = None
+ _particle_reader = True
-class BaseDataQueue(object):
+ class __metaclass__(type):
+ def __init__(cls, name, b, d):
+ type.__init__(cls, name, b, d)
+ if hasattr(cls, "_data_style"):
+ io_registry[cls._data_style] = cls
def __init__(self):
self.queue = defaultdict(dict)
@@ -220,7 +54,14 @@
return self.modify(self.queue[grid.id].pop(field))
else:
# We only read the one set and do not store it if it isn't pre-loaded
- return self._read_set(grid, field)
+ return self._read_data_set(grid, field)
+
+ def _read_particles(self, fields, rtype, args, grid_list, enclosed,
+ conv_factors):
+ filenames = [g.filename for g in grid_list]
+ ids = [g.id for g in grid_list]
+ return HDF5LightReader.ReadParticles(
+ rtype, fields, filenames, ids, conv_factors, args, 0)
def peek(self, grid, field):
return self.queue[grid.id].get(field, None)
@@ -230,31 +71,139 @@
raise ValueError
self.queue[grid][field] = data
-class DataQueueHDF4(BaseDataQueue):
- def _read_set(self, grid, field):
- return readDataHDF4(grid, field)
+ # Now we define our interface
+ def _read_data_set(self, grid, field):
+ pass
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ pass
+
+ def _read_field_names(self, grid):
+ pass
+
+ @property
+ def _read_exception(self):
+ return None
+
+class IOHandlerHDF4(BaseIOHandler):
+
+ _data_style = "enzo_hdf4"
def modify(self, field):
return field.swapaxes(0,2)
-class DataQueueHDF4_2D(BaseDataQueue):
- def _read_set(self, grid, field):
+ def _read_field_names(self, grid):
+ """
+ Returns a list of fields associated with the filename
+ Should *only* be called as EnzoGridInstance.getFields, never as getFields(object)
+ """
+ return SD.SD(grid.filename).datasets().keys()
+
+ def _read_data_set(self, grid, field):
+ """
+ Returns after having obtained or generated a field. Should throw an
+ exception. Should only be called as EnzoGridInstance.readData()
+
+ @param field: field to read
+ @type field: string
+ """
+ return SD.SD(grid.filename).select(field).get()
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ """
+ Reads a slice through the HDF4 data
+
+ @param grid: Grid to slice
+ @type grid: L{EnzoGrid<EnzoGrid>}
+ @param field: field to get
+ @type field: string
+ @param sl: region to get
+ @type sl: SliceType
+ """
+ sl = [slice(None), slice(None), slice(None)]
+ sl[axis] = slice(coord, coord + 1)
+ sl = tuple(reversed(sl))
+ return SD.SD(grid.filename).select(field)[sl].swapaxes(0,2)
+
+ @property
+ def _read_exception(self):
+ return SD.HDF4Error
+
+class IOHandlerHDF4_2D(IOHandlerHDF4):
+
+ _data_style = "enzo_hdf4_2d"
+
+ def _read_data_set(self, grid, field):
t = SD.SD(grid.filename).select(field).get()[:,:,None]
return t.swapaxes(0,1)
+ def _read_data_slice(self, grid, field, axis, coord):
+ t = SD.SD(grid.filename).select(field).get()
+ return t.transpose()
+
def modify(self, field):
- pass
+ return field
-class DataQueueHDF5(BaseDataQueue):
- def _read_set(self, grid, field):
- return readDataHDF5(grid, field)
+class IOHandlerHDF5(BaseIOHandler):
+
+ _data_style = "enzo_hdf5"
+
+ def _read_field_names(self, grid):
+ """
+ Returns a list of fields associated with the filename
+ Should *only* be called as EnzoGridInstance.getFields, never as getFields(object)
+ """
+ return HDF5LightReader.ReadListOfDatasets(grid.filename, "/")
+
+ def _read_data_set(self, grid, field):
+ return HDF5LightReader.ReadData(grid.filename, "/%s" % field).swapaxes(0,2)
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ """
+ Reads a slice through the HDF5 data
+
+ @param grid: Grid to slice
+ @type grid: L{EnzoGrid<EnzoGrid>}
+ @param field: field to get
+ @type field: string
+ @param axis: axis to slice along
+ @param coord: coord to slice at
+ """
+ axis = {0:2,1:1,2:0}[axis]
+ t = HDF5LightReader.ReadDataSlice(grid.filename, "/%s" %
+ (field), axis, coord).transpose()
+ return t
def modify(self, field):
return field.swapaxes(0,2)
-class DataQueuePackedHDF5(BaseDataQueue):
- def _read_set(self, grid, field):
- return readDataPacked(grid, field)
+ @property
+ def _read_exception(self):
+ return (exceptions.KeyError, HDF5LightReader.ReadingError)
+
+class IOHandlerExtracted(BaseIOHandler):
+
+ _data_style = 'extracted'
+
+ def _read_data_set(self, grid, field):
+ return (grid.base_grid[field] / grid.base_grid.convert(field))
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ sl = [slice(None), slice(None), slice(None)]
+ sl[axis] = slice(coord, coord + 1)
+ return grid.base_grid[field][tuple(sl)] / grid.base_grid.convert(field)
+
+class IOHandlerPackedHDF5(BaseIOHandler):
+
+ _data_style = "enzo_packed_3d"
+ _particle_reader = True
+
+ def _read_particles(self, fields, rtype, args, grid_list, enclosed,
+ conv_factors):
+ filenames = [g.filename for g in grid_list]
+ ids = [g.id for g in grid_list]
+ return HDF5LightReader.ReadParticles(
+ rtype, fields, filenames, ids, conv_factors, args, 1)
def modify(self, field):
return field.swapaxes(0,2)
@@ -265,7 +214,7 @@
pf_field_list = grids[0].pf.h.field_list
sets = [dset for dset in list(sets) if dset in pf_field_list]
for g in grids: files_keys[g.filename].append(g)
- exc = getExceptionHDF5()
+ exc = self._read_exception
for file in files_keys:
mylog.debug("Starting read %s (%s)", file, sets)
nodes = [g.id for g in files_keys[file]]
@@ -278,20 +227,42 @@
for gid in data: self.queue[gid].update(data[gid])
mylog.debug("Finished read of %s", sets)
-class DataQueueInMemory(BaseDataQueue):
+ def _read_data_set(self, grid, field):
+ return HDF5LightReader.ReadData(grid.filename,
+ "/Grid%08i/%s" % (grid.id, field)).swapaxes(0,2)
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ axis = _axis_ids[axis]
+ return HDF5LightReader.ReadDataSlice(grid.filename, "/Grid%08i/%s" %
+ (grid.id, field), axis, coord).transpose()
+
+ def _read_field_names(self, grid):
+ return HDF5LightReader.ReadListOfDatasets(
+ grid.filename, "/Grid%08i" % grid.id)
+
+
+class IOHandlerInMemory(BaseIOHandler):
+
+ _data_style = "enzo_inline"
+
def __init__(self, ghost_zones=3):
import enzo
+ self.enzo = enzo
self.grids_in_memory = enzo.grid_data
self.old_grids_in_memory = enzo.old_grid_data
self.my_slice = (slice(ghost_zones,-ghost_zones),
slice(ghost_zones,-ghost_zones),
slice(ghost_zones,-ghost_zones))
- BaseDataQueue.__init__(self)
+ BaseIOHandler.__init__(self)
- def _read_set(self, grid, field):
- import enzo
+ def _read_data_set(self, grid, field):
if grid.id not in self.grids_in_memory: raise KeyError
- return self.grids_in_memory[grid.id][field].swapaxes(0,2)[self.my_slice]
+ tr = self.grids_in_memory[grid.id][field]
+ # If it's particles, we copy.
+ if len(tr.shape) == 1: return tr.copy()
+ # New in-place unit conversion breaks if we don't copy first
+ return tr.swapaxes(0,2)[self.my_slice].copy()
+ # We don't do this, because we currently do not interpolate
coef1 = max((grid.Time - t1)/(grid.Time - t2), 0.0)
coef2 = 1.0 - coef1
t1 = enzo.yt_parameter_file["InitialTime"]
@@ -303,32 +274,64 @@
def modify(self, field):
return field.swapaxes(0,2)
- def preload(self, grids, sets):
- pass
+ def _read_field_names(self, grid):
+ return self.grids_in_memory[grid.id].keys()
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ sl = [slice(3,-3), slice(3,-3), slice(3,-3)]
+ sl[axis] = slice(coord + 3, coord + 4)
+ sl = tuple(reversed(sl))
+ tr = self.grids_in_memory[grid.id][field][sl].swapaxes(0,2)
+ # In-place unit conversion requires we return a copy
+ return tr.copy()
+
+ @property
+ def _read_exception(self):
+ return KeyError
+
+class IOHandlerNative(BaseIOHandler):
-class DataQueueNative(BaseDataQueue):
- def _read_set(self, grid, field):
+ _data_style = "orion_native"
+
+ def _read_data_set(self, grid, field):
return readDataNative(grid, field)
def modify(self, field):
return field.swapaxes(0,2)
-class DataQueuePacked2D(BaseDataQueue):
- def _read_set(self, grid, field):
+class IOHandlerPacked2D(IOHandlerPackedHDF5):
+
+ _data_style = "enzo_packed_2d"
+
+ def _read_data_set(self, grid, field):
return HDF5LightReader.ReadData(grid.filename,
"/Grid%08i/%s" % (grid.id, field)).transpose()[:,:,None]
def modify(self, field):
pass
-class DataQueuePacked1D(BaseDataQueue):
- def _read_set(self, grid, field):
+ def _read_data_slice(self, grid, field, axis, coord):
+ t = HDF5LightReader.ReadData(grid.filename, "/Grid%08i/%s" %
+ (grid.id, field)).transpose()
+ return t
+
+
+class IOHandlerPacked1D(IOHandlerPackedHDF5):
+
+ _data_style = "enzo_packed_1d"
+
+ def _read_data_set(self, grid, field):
return HDF5LightReader.ReadData(grid.filename,
"/Grid%08i/%s" % (grid.id, field)).transpose()[:,None,None]
def modify(self, field):
pass
+ def _read_data_slice(self, grid, field, axis, coord):
+ t = HDF5LightReader.ReadData(grid.filename, "/Grid%08i/%s" %
+ (grid.id, field))
+ return t
+
#
# BoxLib/Orion data readers follow
#
Modified: trunk/yt/lagos/DerivedQuantities.py
==============================================================================
--- trunk/yt/lagos/DerivedQuantities.py (original)
+++ trunk/yt/lagos/DerivedQuantities.py Fri Nov 20 20:27:29 2009
@@ -67,7 +67,7 @@
self.func(e, *args, **kwargs)
mylog.debug("Preloading %s", e.requested)
self._preload([g for g in self._get_grid_objs()], e.requested,
- self._data_source.pf.h.queue)
+ self._data_source.pf.h.io)
if lazy_reader and not self.force_unlazy:
return self._call_func_lazy(args, kwargs)
else:
@@ -316,6 +316,9 @@
int idx2 = blockIdx.y * blockDim.x;
int offset = threadIdx.x;
+ /* Here we're just setting up convenience pointers to our
+ shared array */
+
float* x_data1 = (float*) array;
float* y_data1 = (float*) &x_data1[blockDim.x];
float* z_data1 = (float*) &y_data1[blockDim.x];
Modified: trunk/yt/lagos/EnzoFields.py
==============================================================================
--- trunk/yt/lagos/EnzoFields.py (original)
+++ trunk/yt/lagos/EnzoFields.py Fri Nov 20 20:27:29 2009
@@ -24,6 +24,7 @@
"""
from UniversalFields import *
+from yt.amr_utils import CICDeposit_3
rho_crit_now = 1.8788e-29 # times h^2
@@ -228,6 +229,21 @@
add_field("particle_density", function=_pdensity,
validators=[ValidateSpatial(0)], convert_function=_convertDensity)
+def _pdensity_pyx(field, data):
+ blank = na.zeros(data.ActiveDimensions, dtype='float32')
+ if data.NumberOfParticles == 0: return blank
+ CICDeposit_3(data["particle_position_x"].astype(na.float64),
+ data["particle_position_y"].astype(na.float64),
+ data["particle_position_z"].astype(na.float64),
+ data["particle_mass"].astype(na.float32),
+ na.int64(data.NumberOfParticles),
+ blank, na.array(data.LeftEdge).astype(na.float64),
+ na.array(data.ActiveDimensions).astype(na.int32),
+ na.float64(data['dx']))
+ return blank
+add_field("particle_density_pyx", function=_pdensity_pyx,
+ validators=[ValidateSpatial(0)], convert_function=_convertDensity)
+
def _spdensity(field, data):
blank = na.zeros(data.ActiveDimensions, dtype='float32', order="FORTRAN")
if data.NumberOfParticles == 0: return blank
@@ -242,6 +258,23 @@
add_field("star_density", function=_spdensity,
validators=[ValidateSpatial(0)], convert_function=_convertDensity)
+def _spdensity_pyx(field, data):
+ blank = na.zeros(data.ActiveDimensions, dtype='float32')
+ if data.NumberOfParticles == 0: return blank
+ filter = data['creation_time'] > 0.0
+ if not filter.any(): return blank
+ CICDeposit_3(data["particle_position_x"][filter].astype(na.float64),
+ data["particle_position_y"][filter].astype(na.float64),
+ data["particle_position_z"][filter].astype(na.float64),
+ data["particle_mass"][filter].astype(na.float32),
+ na.int64(na.where(filter)[0].size),
+ blank, na.array(data.LeftEdge).astype(na.float64),
+ na.array(data.ActiveDimensions).astype(na.int32),
+ na.float64(data['dx']))
+ return blank
+add_field("star_density_pyx", function=_spdensity_pyx,
+ validators=[ValidateSpatial(0)], convert_function=_convertDensity)
+
EnzoFieldInfo["Temperature"].units = r"K"
#
Modified: trunk/yt/lagos/FieldInfoContainer.py
==============================================================================
--- trunk/yt/lagos/FieldInfoContainer.py (original)
+++ trunk/yt/lagos/FieldInfoContainer.py Fri Nov 20 20:27:29 2009
@@ -148,7 +148,7 @@
class FieldDetector(defaultdict):
Level = 1
- NumberOfParticles = 0
+ NumberOfParticles = 1
_read_exception = None
def __init__(self, nd = 16, pf = None):
self.nd = nd
@@ -156,10 +156,16 @@
self.LeftEdge = [0.0,0.0,0.0]
self.RightEdge = [1.0,1.0,1.0]
self['dx'] = self['dy'] = self['dz'] = na.array([1.0])
- self.fields = []
if pf is None:
pf = defaultdict(lambda: 1)
self.pf = pf
+ class fake_hierarchy(object):
+ class fake_io(object):
+ def _read_data_set(io_self, data, field):
+ return self._read_data(field)
+ _read_exception = RuntimeError
+ io = fake_io()
+ self.hierarchy = fake_hierarchy()
self.requested = []
self.requested_parameters = []
defaultdict.__init__(self, lambda: na.ones((nd,nd,nd)))
@@ -184,6 +190,10 @@
def _read_data(self, field_name):
self.requested.append(field_name)
+ if FieldInfo.has_key(field_name) and \
+ FieldInfo[field_name].particle_type:
+ self.requested.append(field_name)
+ return na.ones(self.NumberOfParticles)
return defaultdict.__missing__(self, field_name)
def get_field_parameter(self, param):
@@ -201,6 +211,7 @@
class DerivedField(object):
def __init__(self, name, function,
convert_function = None,
+ particle_convert_function = None,
units = "", projected_units = "",
take_log = True, validators = None,
particle_type = False, vector_field=False,
@@ -237,6 +248,7 @@
if not convert_function:
convert_function = lambda a: 1.0
self._convert_function = convert_function
+ self._particle_convert_function = particle_convert_function
self.particle_type = particle_type
self.vector_field = vector_field
self.projection_conversion = projection_conversion
@@ -282,10 +294,10 @@
Return the value of the field in a given *data* object.
"""
ii = self.check_available(data)
- original_fields = data.fields[:] # Copy
+ original_fields = data.keys() # Copy
dd = self._function(self, data)
dd *= self._convert_function(data)
- for field_name in data.fields:
+ for field_name in data.keys():
if field_name not in original_fields:
del data[field_name]
return dd
@@ -309,6 +321,11 @@
data_label += r"$"
return data_label
+ def particle_convert(self, data):
+ if self._particle_convert_function is not None:
+ return self._particle_convert_function(data)
+ return None
+
class FieldValidator(object):
pass
Modified: trunk/yt/lagos/HDF5LightReader.c
==============================================================================
--- trunk/yt/lagos/HDF5LightReader.c (original)
+++ trunk/yt/lagos/HDF5LightReader.c Fri Nov 20 20:27:29 2009
@@ -38,6 +38,46 @@
static PyObject *_hdf5ReadError;
herr_t iterate_dataset(hid_t loc_id, const char *name, void *nodelist);
+/* Structures for particle reading */
+
+typedef struct particle_validation_ {
+ int total_valid_particles;
+ int particles_to_check;
+ int nread;
+ int stride_size;
+ int *mask;
+ int update_count;
+ int nfields;
+ char **field_names;
+ PyArrayObject *conv_factors;
+ PyArrayObject **return_values;
+ int (*count_func)(struct particle_validation_ *data);
+ int (*count_func_float)(struct particle_validation_ *data);
+ int (*count_func_double)(struct particle_validation_ *data);
+ int (*count_func_longdouble)(struct particle_validation_ *data);
+ void *validation_reqs;
+ void *particle_position[3];
+} particle_validation;
+
+typedef struct region_validation_ {
+ /* These cannot contain any pointers */
+ npy_float64 left_edge[3];
+ npy_float64 right_edge[3];
+ npy_float64 period[3];
+ int periodic;
+} region_validation;
+
+/* Forward declarations */
+int setup_validator_region(particle_validation *data, PyObject *InputData);
+int run_validators(particle_validation *pv, char *filename,
+ int grid_id, const int read, const int packed,
+ int grid_index);
+
+/* Different data type validators */
+
+int count_particles_region_FLOAT(particle_validation *data);
+int count_particles_region_DOUBLE(particle_validation *data);
+
int get_my_desc_type(hid_t native_type_id){
/*
@@ -531,6 +571,12 @@
long id;
char grid_node_name[13]; // Grid + 8 + \0
+ /* Similar to the way Enzo does it, we're going to set the file access
+ property to store bigger bits in RAM. */
+
+ int memory_increment = 1024*1024; /* in bytes */
+ int dump_flag = 0;
+
file_id = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0) {
@@ -658,11 +704,612 @@
return NULL;
}
+static PyObject *
+Py_ReadParticles(PyObject *obj, PyObject *args)
+{
+ // Process:
+ // - Interpret arguments:
+ // - List of filenames
+ // - List of grid ids (this is Enzo HDF5-specific)
+ // - List of conversion factors for particles
+ // - Validation arguments
+ // - List of true/false for validation
+ // - Set up validation scheme
+ // - Iterate over grids, counting
+ // - Allocate array
+ // - Iterate over grids, turning on and off dataspace
+
+ int source_type, ig, id, i, ifield, packed;
+ Py_ssize_t ngrids, nfields;
+
+ PyObject *field_list, *filename_list, *grid_ids, *oconv_factors, *vargs;
+ PyArrayObject *conv_factors = NULL;
+ int stride_size = 10000000;
+ particle_validation pv;
+
+ /* We explicitly initialize pointers to NULL to aid with failure handling */
+
+ vargs = field_list = filename_list = grid_ids = NULL;
+ pv.mask = NULL;
+ pv.field_names = NULL;
+ pv.validation_reqs = NULL;
+ pv.particle_position[0] = pv.particle_position[1] = pv.particle_position[2] = NULL;
+ pv.return_values = NULL;
+
+ /* Set initial values for pv */
+ pv.stride_size = stride_size;
+ pv.total_valid_particles = pv.particles_to_check = pv.nread = pv.nfields = 0;
+
+ if (!PyArg_ParseTuple(args, "iOOOOOi",
+ &source_type, /* This is a non-public API, so we just take ints */
+ &field_list, &filename_list, &grid_ids, &oconv_factors, &vargs,
+ &packed))
+ return PyErr_Format(_hdf5ReadError,
+ "ReadParticles: Invalid parameters.");
+
+ if (!PyList_Check(field_list)) {
+ PyErr_Format(_hdf5ReadError,
+ "ReadParticles: field_list is not a list!\n");
+ goto _fail;
+ }
+ nfields = PyList_Size(field_list);
+ pv.nfields = nfields;
+
+ if (!PyList_Check(filename_list)) {
+ PyErr_Format(_hdf5ReadError,
+ "ReadParticles: filename_list is not a list!\n");
+ goto _fail;
+ }
+ ngrids = PyList_Size(filename_list);
+
+ if (!PyList_Check(grid_ids)
+ || (PyList_Size(grid_ids) != ngrids)) {
+ PyErr_Format(_hdf5ReadError,
+ "ReadParticles: grid_ids is not a list of correct length!\n");
+ goto _fail;
+ }
+
+ conv_factors = (PyArrayObject *) PyArray_FromAny(oconv_factors,
+ PyArray_DescrFromType(NPY_FLOAT64), 2, 2,
+ 0, NULL);
+ if( (conv_factors == NULL) ||
+ !(PyArray_DIM(conv_factors, 0) == ngrids) ||
+ !(PyArray_DIM(conv_factors, 1) == nfields) ) {
+ PyErr_Format(_hdf5ReadError,
+ "ReadParticles: conv_factors is not an array of (ngrids, nfields)");
+ goto _fail;
+ }
+
+ if (!PyTuple_Check(vargs)) {
+ PyErr_Format(_hdf5ReadError,
+ "ReadParticles: vargs is not a tuple!\n");
+ goto _fail;
+ }
+
+ /* We've now parsed all our arguments and it is time to set up the
+ validator */
+
+ /* First initialize our particle_validation structure */
+
+ pv.mask = (int*) malloc(sizeof(int) * stride_size);
+
+ switch(source_type) {
+ case 0:
+ /* Region type */
+ setup_validator_region(&pv, vargs);
+ break;
+ default:
+ PyErr_Format(_hdf5ReadError,
+ "Unrecognized data source.\n");
+ goto _fail;
+ break;
+ }
+
+ /* Okay, now we open our files and stride over each grid in the files. */
+
+ PyObject *temp = NULL;
+ char *filename = NULL;
+ pv.update_count = 1;
+ pv.conv_factors = conv_factors;
+
+ for (ig = 0; ig < ngrids ; ig++) {
+ temp = PyList_GetItem(filename_list, ig);
+ filename = PyString_AsString(temp);
+ temp = PyList_GetItem(grid_ids, ig);
+ id = PyInt_AsLong(temp);
+ //fprintf(stderr, "Counting from grid %d\n", id);
+ if(run_validators(&pv, filename, id, 0, packed, ig) < 0) {
+ goto _fail;
+ }
+ }
+ /* Now we know how big to make our array, hooray. */
+ pv.update_count = 0;
+
+ pv.return_values = (PyArrayObject**) malloc(
+ sizeof(PyArrayObject*) * nfields);
+ pv.field_names = (char **) malloc(
+ sizeof(char *) * nfields);
+ for (ifield = 0; ifield < nfields; ifield++) {
+ pv.return_values[ifield] = NULL;
+ pv.field_names[ifield] = PyString_AsString(PyList_GetItem(field_list, ifield));
+ }
+
+ /* Now we know how many particles we want. */
+
+ for (ig = 0; ig < ngrids ; ig++) {
+ /* We should have some logic here to read the entire thing
+ in a stride, without checking particle positions,
+ if it's fully-enclosed. */
+ temp = PyList_GetItem(filename_list, ig);
+ filename = PyString_AsString(temp);
+ temp = PyList_GetItem(grid_ids, ig);
+ id = PyInt_AsLong(temp);
+ //fprintf(stderr, "Reading from grid %d\n", id);
+ if(run_validators(&pv, filename, id, 1, packed, ig) < 0) {
+ goto _fail;
+ }
+ }
+
+ /* Let's pack up our return values */
+ PyObject *my_list = PyList_New(0);
+ for (i = 0; i < pv.nfields ; i++){
+ PyList_Append(my_list, (PyObject *) pv.return_values[i]);
+ Py_DECREF(pv.return_values[i]);
+ }
+ PyObject *return_value = Py_BuildValue("N", my_list);
+
+ /* Now we do some finalization */
+ free(pv.mask);
+ free(pv.field_names);
+ free(pv.return_values); /* Has to happen after packing our return value */
+ for (i = 0; i<3; i++) {
+ free(pv.particle_position[i]);
+ }
+ Py_DECREF(conv_factors);
+ free(pv.validation_reqs);
+ /* We don't need to free pv */
+
+ return return_value;
+
+ _fail:
+
+ if(pv.mask != NULL) free(pv.mask);
+ if(pv.field_names != NULL) {
+ for (i = 0; i<pv.nfields; i++) {
+ free(pv.field_names[i]);
+ }
+ free(pv.field_names);
+ }
+ if(conv_factors != NULL) { Py_DECREF(conv_factors); }
+ if(pv.return_values != NULL){
+ for (i = 0; i < pv.nfields; i++) {
+ if(pv.return_values[i] != NULL) { Py_DECREF(pv.return_values[i]); }
+ }
+ free(pv.return_values);
+ }
+ for(i = 0; i < 3; i++) {
+ if(pv.particle_position[i] != NULL) free(pv.particle_position[i]);
+ }
+ if(pv.validation_reqs != NULL) free(pv.validation_reqs);
+
+ return NULL;
+}
+
+int setup_validator_region(particle_validation *data, PyObject *InputData)
+{
+ int i;
+ /* These are borrowed references */
+ PyArrayObject *left_edge = (PyArrayObject *) PyTuple_GetItem(InputData, 0);
+ PyArrayObject *right_edge = (PyArrayObject *) PyTuple_GetItem(InputData, 1);
+ PyObject *operiodic = PyTuple_GetItem(InputData, 2);
+ npy_float64 DW;
+
+ /* This will get freed in the finalization of particle validation */
+ region_validation *rv = (region_validation *)
+ malloc(sizeof(region_validation));
+ data->validation_reqs = (void *) rv;
+
+ for (i = 0; i < 3; i++){
+ rv->left_edge[i] = *(npy_float64*) PyArray_GETPTR1(left_edge, i);
+ rv->right_edge[i] = *(npy_float64*) PyArray_GETPTR1(right_edge, i);
+ }
+
+ rv->periodic = PyInt_AsLong(operiodic);
+ if(rv->periodic == 1) {
+ PyArrayObject *domain_left_edge = (PyArrayObject *) PyTuple_GetItem(InputData, 3);
+ PyArrayObject *domain_right_edge = (PyArrayObject *) PyTuple_GetItem(InputData, 4);
+ for (i = 0; i < 3; i++){
+ DW = (*(npy_float64*) PyArray_GETPTR1(domain_right_edge, i))
+ - (*(npy_float64*) PyArray_GETPTR1(domain_left_edge, i));
+ rv->period[i] = DW;
+ fprintf(stderr, "Setting period equal to %lf\n", rv->period[i]);
+ }
+ }
+
+ data->count_func = NULL;
+ data->count_func_float = count_particles_region_FLOAT;
+ data->count_func_double = count_particles_region_DOUBLE;
+
+ /* We need to insert more periodic logic here */
+
+ return 1;
+}
+
+int run_validators(particle_validation *pv, char *filename,
+ int grid_id, const int read, const int packed,
+ int grid_index)
+{
+ int i, ifield, p_ind;
+ hid_t file_id;
+ hid_t dataset_x, dataset_y, dataset_z;
+ hid_t dataspace, memspace;
+ hid_t datatype_id, native_type_id;
+ hid_t rdatatype_id, rnative_type_id;
+ char name_x[255], name_y[255], name_z[255];
+ hsize_t num_part_this_grid;
+ hsize_t current_pos = 0;
+ hsize_t num_particles_to_read = 0;
+ hid_t *dataset_read;
+ dataset_read = NULL;
+ npy_float64 *cfactors = NULL;
+
+ cfactors = (npy_float64*) malloc(pv->nfields * sizeof(npy_float64));
+
+ for (ifield = 0; ifield < pv->nfields; ifield++){
+ cfactors[ifield] = *(npy_float64 *) PyArray_GETPTR2(
+ pv->conv_factors, grid_index, ifield);
+ }
+
+ /* We set these to -1 to identify which haven't been used */
+ file_id = dataset_x = dataset_y = dataset_z = -1;
+ dataspace = memspace = datatype_id = native_type_id = -1;
+ rdatatype_id = rnative_type_id = -1;
+
+ int npy_type;
+
+ if (packed == 1) {
+ snprintf(name_x, 254, "/Grid%08d/particle_position_x", grid_id);
+ snprintf(name_y, 254, "/Grid%08d/particle_position_y", grid_id);
+ snprintf(name_z, 254, "/Grid%08d/particle_position_z", grid_id);
+ } else {
+ snprintf(name_x, 254, "/particle_position_x");
+ snprintf(name_y, 254, "/particle_position_y");
+ snprintf(name_z, 254, "/particle_position_z");
+ }
+
+ /* First we open the file */
+
+ file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+ if (file_id < 0) {
+ PyErr_Format(_hdf5ReadError,
+ "ReadHDF5DataSet: Unable to open %s", filename);
+ goto _fail;
+ }
+
+ /* Figure out the size of our array */
+
+ /* Note that we are using one dimension, not multiple. If you're using
+ Enzo data with multiple particle dimensions, well, I'm afraid I can't help you.
+ Just now, at least. */
+
+ dataset_x = H5Dopen(file_id, name_x);
+ dataset_y = H5Dopen(file_id, name_y);
+ dataset_z = H5Dopen(file_id, name_z);
+
+ dataspace = H5Dget_space(dataset_x);
+ num_part_this_grid = H5Sget_simple_extent_npoints(dataspace);
+
+ /* Let's get the information about the datatype now */
+ datatype_id = H5Dget_type(dataset_x);
+ native_type_id = H5Tget_native_type(datatype_id, H5T_DIR_ASCEND);
+
+ /* If the count function is not set, set it and allocate arrays */
+ if(pv->count_func == NULL){
+ int typesize = 0;
+ if(H5Tequal(native_type_id, H5T_NATIVE_FLOAT) > 0) {
+ pv->count_func = pv->count_func_float;
+ typesize = sizeof(float);
+ } else if (H5Tequal(native_type_id, H5T_NATIVE_DOUBLE) > 0 ) {
+ pv->count_func = pv->count_func_double;
+ typesize = sizeof(double);
+ } else if (H5Tequal(native_type_id, H5T_NATIVE_LDOUBLE) > 0 ) {
+ pv->count_func = pv->count_func_longdouble;
+ typesize = sizeof(long double);
+ } else {
+ H5Tclose(datatype_id); H5Tclose(native_type_id);
+ PyErr_Format(_hdf5ReadError,
+ "ReadHDF5DataSet: Unrecognized particle position array type");
+ goto _fail;
+ }
+ /* We allocate arrays here */
+ for(i = 0; i < 3; i++)
+ pv->particle_position[i] = malloc(pv->stride_size * typesize);
+ }
+
+ if(read == 1) {
+ /* We allocate space for the to-read datasets */
+ dataset_read = (hid_t*) malloc(pv->nfields * sizeof(hid_t));
+ for (i = 0; i < pv->nfields; i++) {
+ char toread[255];
+ if (packed == 1) {
+ snprintf(toread, 255, "/Grid%08d/%s", grid_id, pv->field_names[i]);
+ } else {
+ snprintf(toread, 255, "/%s", pv->field_names[i]);
+ }
+
+ dataset_read[i] = H5Dopen(file_id, toread);
+ /* We know how many particles we will want, so we allocate our
+ output arrays */
+ if(pv->return_values[i] == NULL){
+ /* Get the data type */
+ npy_intp dims = pv->total_valid_particles;
+ rdatatype_id = H5Dget_type(dataset_read[i]);
+ rnative_type_id = H5Tget_native_type(rdatatype_id, H5T_DIR_ASCEND);
+ npy_type = get_my_desc_type(rnative_type_id);
+ //fprintf(stderr, "Allocating array of size %d\n", (int) dims);
+ pv->return_values[i] = (PyArrayObject *)
+ PyArray_SimpleNewFromDescr(
+ 1, &dims, PyArray_DescrFromType(npy_type));
+ H5Tclose(rnative_type_id);
+ H5Tclose(rdatatype_id);
+ }
+ }
+ }
+
+ /* Now we create an in-memory dataspace */
+
+ /* We begin our iteration over the strides here */
+ while(current_pos < num_part_this_grid) {
+ num_particles_to_read = (
+ (current_pos + pv->stride_size >= num_part_this_grid) ?
+ num_part_this_grid - current_pos : pv->stride_size );
+ pv->particles_to_check = num_particles_to_read;
+
+ memspace = H5Screate_simple(1, &num_particles_to_read, NULL);
+
+ H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, ¤t_pos, NULL,
+ &num_particles_to_read, NULL);
+ H5Dread(dataset_x, native_type_id, memspace, dataspace, H5P_DEFAULT,
+ pv->particle_position[0]);
+ H5Dread(dataset_y, native_type_id, memspace, dataspace, H5P_DEFAULT,
+ pv->particle_position[1]);
+ H5Dread(dataset_z, native_type_id, memspace, dataspace, H5P_DEFAULT,
+ pv->particle_position[2]);
+
+ hsize_t read_here = pv->count_func(pv);
+
+ if((read == 1) && (read_here > 0)) {
+ /* First we select the dataspace */
+ H5Sselect_none(dataspace);
+ int num_copied = 0;
+ hsize_t *coords = malloc(sizeof(hsize_t) * read_here);
+ for(i = 0 ; i < num_particles_to_read ; i++) {
+ /* This might be quite slow */
+ if(pv->mask[i] == 1){
+ coords[num_copied++] = current_pos + i;
+ }
+ }
+ H5Sselect_elements(dataspace, H5S_SELECT_SET, read_here, coords);
+ H5Sset_extent_simple(memspace, 1, &read_here, &read_here);
+ free(coords);
+ for (ifield = 0; ifield < pv->nfields; ifield++){
+ rdatatype_id = H5Dget_type(dataset_read[ifield]);
+ rnative_type_id = H5Tget_native_type(rdatatype_id, H5T_DIR_ASCEND);
+ H5Dread(dataset_read[ifield], rnative_type_id,
+ memspace, dataspace, H5P_DEFAULT,
+ (void*) PyArray_GETPTR1(pv->return_values[ifield], pv->nread));
+ H5Tclose(rnative_type_id);
+ H5Tclose(rdatatype_id);
+ /* Now we multiply our fields by the appropriate conversion factor */
+ if (cfactors[ifield] != 1.0) {
+ for(p_ind = 0; p_ind < read_here; p_ind++)
+ if (npy_type == 11) { // floats
+ *(npy_float32 *) PyArray_GETPTR1(
+ pv->return_values[ifield], p_ind + pv->nread)
+ *= cfactors[ifield];
+ } else { // doubles
+ *(npy_float64 *) PyArray_GETPTR1(
+ pv->return_values[ifield], p_ind + pv->nread)
+ *= cfactors[ifield];
+ }
+ }
+ }
+ pv->nread += read_here;
+ }
+
+ current_pos += num_particles_to_read;
+ }
+
+ H5Dclose(dataset_x);
+ H5Dclose(dataset_y);
+ H5Dclose(dataset_z);
+ H5Sclose(dataspace);
+ H5Sclose(memspace);
+ H5Tclose(datatype_id);
+ H5Tclose(native_type_id);
+ free(cfactors);
+
+ if (read == 1) {
+ for (i = 0; i < pv->nfields; i++) {
+ H5Dclose(dataset_read[i]);
+ }
+ free(dataset_read);
+ }
+
+ H5Fclose(file_id);
+
+ return 1;
+
+ _fail:
+ /* Now the nasty business of closing our HDF5 references */
+ if(cfactors!=NULL)free(cfactors);
+ if(!(dataset_x <= 0)&&(H5Iget_ref(dataset_x))) H5Dclose(dataset_x);
+ if(!(dataset_y <= 0)&&(H5Iget_ref(dataset_y))) H5Dclose(dataset_y);
+ if(!(dataset_z <= 0)&&(H5Iget_ref(dataset_z))) H5Dclose(dataset_z);
+ if(!(dataspace <= 0)&&(H5Iget_ref(dataspace))) H5Sclose(dataspace);
+ if(!(memspace <= 0)&&(H5Iget_ref(memspace))) H5Sclose(memspace);
+ if(!(native_type_id <= 0)&&(H5Iget_ref(native_type_id))) H5Tclose(native_type_id);
+ if(!(datatype_id <= 0)&&(H5Iget_ref(datatype_id))) H5Tclose(datatype_id);
+ if ((read == 1) && (dataset_read != NULL)) {
+ for (i = 0; i < pv->nfields; i++) {
+ if(!(dataset_read[i] <= 0)&&(H5Iget_ref(dataset_read[i])))
+ H5Dclose(dataset_read[i]);
+ }
+ free(dataset_read);
+ }
+ if(!(file_id <= 0)&&(H5Iget_ref(file_id))) H5Fclose(file_id);
+
+ return 0;
+}
+
+/* Hate to do copy-pasta here, but I think it is necessary */
+
+int count_particles_region_FLOAT(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n;
+ region_validation *vdata;
+
+ vdata = (region_validation*) data->validation_reqs;
+
+ float **particle_data = (float **) data->particle_position;
+
+ float *particle_position_x = particle_data[0];
+ float *particle_position_y = particle_data[1];
+ float *particle_position_z = particle_data[2];
+ float tempx, tempy, tempz;
+
+ if (vdata->periodic == 0) {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ if ( (particle_position_x[ind] >= vdata->left_edge[0])
+ && (particle_position_x[ind] <= vdata->right_edge[0])
+ && (particle_position_y[ind] >= vdata->left_edge[1])
+ && (particle_position_y[ind] <= vdata->right_edge[1])
+ && (particle_position_z[ind] >= vdata->left_edge[2])
+ && (particle_position_z[ind] <= vdata->right_edge[2])) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ } else {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ tempx = particle_position_x[ind];
+ tempy = particle_position_y[ind];
+ tempz = particle_position_z[ind];
+ if ( (tempx < vdata->left_edge[0]) && (tempx < vdata->right_edge[0]) ) {
+ tempx += vdata->period[0];
+ } else if ( (tempx > vdata->left_edge[0]) && (tempx > vdata->right_edge[0]) ) {
+ tempx -= vdata->period[0];
+ }
+ if ( (tempy < vdata->left_edge[1]) && (tempx < vdata->right_edge[1]) ) {
+ tempy += vdata->period[1];
+ } else if ( (tempy > vdata->left_edge[1]) && (tempx > vdata->right_edge[1]) ) {
+ tempy -= vdata->period[1];
+ }
+ if ( (tempz < vdata->left_edge[2]) && (tempx < vdata->right_edge[2]) ) {
+ tempz += vdata->period[2];
+ } else if ( (tempz > vdata->left_edge[2]) && (tempx > vdata->right_edge[2]) ) {
+ tempz -= vdata->period[2];
+ }
+ if ( (tempx >= vdata->left_edge[0])
+ && (tempx <= vdata->right_edge[0])
+ && (tempy >= vdata->left_edge[1])
+ && (tempy <= vdata->right_edge[1])
+ && (tempz >= vdata->left_edge[2])
+ && (tempz <= vdata->right_edge[2])) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ }
+ return n;
+}
+
+int count_particles_region_DOUBLE(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n=0;
+ region_validation *vdata;
+
+ vdata = (region_validation*) data->validation_reqs;
+
+ double **particle_data = (double **) data->particle_position;
+
+ double *particle_position_x = particle_data[0];
+ double *particle_position_y = particle_data[1];
+ double *particle_position_z = particle_data[2];
+ double tempx, tempy, tempz;
+
+ if (vdata->periodic == 0) {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ if ( (particle_position_x[ind] >= vdata->left_edge[0])
+ && (particle_position_x[ind] <= vdata->right_edge[0])
+ && (particle_position_y[ind] >= vdata->left_edge[1])
+ && (particle_position_y[ind] <= vdata->right_edge[1])
+ && (particle_position_z[ind] >= vdata->left_edge[2])
+ && (particle_position_z[ind] <= vdata->right_edge[2])) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ } else {
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ tempx = particle_position_x[ind];
+ tempy = particle_position_y[ind];
+ tempz = particle_position_z[ind];
+ if ( (tempx < vdata->left_edge[0]) && (tempx < vdata->right_edge[0]) ) {
+ tempx += vdata->period[0];
+ } else if ( (tempx > vdata->left_edge[0]) && (tempx > vdata->right_edge[0]) ) {
+ tempx -= vdata->period[0];
+ }
+ if ( (tempy < vdata->left_edge[1]) && (tempx < vdata->right_edge[1]) ) {
+ tempy += vdata->period[1];
+ } else if ( (tempy > vdata->left_edge[1]) && (tempx > vdata->right_edge[1]) ) {
+ tempy -= vdata->period[1];
+ }
+ if ( (tempz < vdata->left_edge[2]) && (tempx < vdata->right_edge[2]) ) {
+ tempz += vdata->period[2];
+ } else if ( (tempz > vdata->left_edge[2]) && (tempx > vdata->right_edge[2]) ) {
+ tempz -= vdata->period[2];
+ }
+ if ( (tempx >= vdata->left_edge[0])
+ && (tempx <= vdata->right_edge[0])
+ && (tempy >= vdata->left_edge[1])
+ && (tempy <= vdata->right_edge[1])
+ && (tempz >= vdata->left_edge[2])
+ && (tempz <= vdata->right_edge[2])) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ }
+ return n;
+}
+
static PyMethodDef _hdf5LightReaderMethods[] = {
{"ReadData", Py_ReadHDF5DataSet, METH_VARARGS},
{"ReadDataSlice", Py_ReadHDF5DataSetSlice, METH_VARARGS},
{"ReadListOfDatasets", Py_ReadListOfDatasets, METH_VARARGS},
{"ReadMultipleGrids", Py_ReadMultipleGrids, METH_VARARGS},
+ {"ReadParticles", Py_ReadParticles, METH_VARARGS},
{NULL, NULL}
};
Modified: trunk/yt/lagos/HaloFinding.py
==============================================================================
--- trunk/yt/lagos/HaloFinding.py (original)
+++ trunk/yt/lagos/HaloFinding.py Fri Nov 20 20:27:29 2009
@@ -28,12 +28,19 @@
from yt.lagos import *
from yt.lagos.hop.EnzoHop import RunHOP
try:
+ from yt.lagos.parallelHOP.parallelHOP import *
+except ImportError:
+ mylog.debug("ParallelHOP not imported.")
+
+try:
from yt.lagos.fof.EnzoFOF import RunFOF
except ImportError:
pass
+from yt.performance_counters import yt_counters, time_function
from kd import *
-import math
+import math, sys
+from collections import defaultdict
class Halo(object):
"""
@@ -48,13 +55,22 @@
dont_wrap = ["get_sphere", "write_particle_list"]
extra_wrap = ["__getitem__"]
- def __init__(self, halo_list, id, indices = None):
+ def __init__(self, halo_list, id, indices = None, size=None, CoM=None,
+ max_dens_point=None, group_total_mass=None, max_radius=None, bulk_vel=None,
+ tasks=None):
self.halo_list = halo_list
self.id = id
self.data = halo_list._data_source
if indices is not None: self.indices = halo_list._base_indices[indices]
# We assume that if indices = None, the instantiator has OTHER plans
# for us -- i.e., setting it somehow else
+ self.size = size
+ self.CoM = CoM
+ self.max_dens_point = max_dens_point
+ self.group_total_mass = group_total_mass
+ self.max_radius = max_radius
+ self.bulk_vel = bulk_vel
+ self.tasks = tasks
def center_of_mass(self):
"""
@@ -116,7 +132,10 @@
return r.max()
def __getitem__(self, key):
- return self.data[key][self.indices]
+ if ytcfg.getboolean("yt","inline") == False:
+ return self.data.particles[key][self.indices]
+ else:
+ return self.data[key][self.indices]
def get_sphere(self, center_of_mass=True):
"""
@@ -140,7 +159,7 @@
handle.create_group("/%s" % gn)
for field in ["particle_position_%s" % ax for ax in 'xyz'] \
+ ["particle_velocity_%s" % ax for ax in 'xyz'] \
- + ["particle_index"]:
+ + ["particle_index"] + ["ParticleMassMsun"]:
handle.create_dataset("/%s/%s" % (gn, field), data=self[field])
n = handle["/%s" % gn]
# set attributes on n
@@ -149,6 +168,135 @@
class HOPHalo(Halo):
pass
+class parallelHOPHalo(Halo,ParallelAnalysisInterface):
+ dont_wrap = ["maximum_density","maximum_density_location",
+ "center_of_mass","total_mass","bulk_velocity","maximum_radius",
+ "get_size","get_sphere", "write_particle_list","__getitem__"]
+
+ def maximum_density(self):
+ """
+ Return the HOP-identified maximum density.
+ """
+ if self.max_dens_point is not None:
+ return self.halo_list._max_dens[self.id][0]
+ max = self._mpi_allmax(self.halo_list._max_dens[self.id][0])
+ return max
+
+ def maximum_density_location(self):
+ """
+ Return the location HOP identified as maximally dense.
+ """
+ if self.max_dens_point is not None:
+ return na.array([self.halo_list._max_dens[self.id][1], self.halo_list._max_dens[self.id][2],
+ self.halo_list._max_dens[self.id][3]])
+ # If I own the maximum density, my location is globally correct.
+ max_dens = self.maximum_density()
+ if self.halo_list._max_dens[self.id][0] == max_dens:
+ value = na.array([
+ self.halo_list._max_dens[self.id][1],
+ self.halo_list._max_dens[self.id][2],
+ self.halo_list._max_dens[self.id][3]])
+ else:
+ value = na.array([0,0,0])
+ # This works, and isn't appropriate but for now will be fine...
+ value = self._mpi_allsum(value)
+ return value
+
+ def center_of_mass(self):
+ """
+ Calculate and return the center of mass.
+ """
+ # If it's precomputed, we save time!
+ if self.CoM is not None:
+ return self.CoM
+ # This need to be called by all tasks, but not all will end up using
+ # it.
+ c_vec = self.maximum_density_location() - na.array([0.5,0.5,0.5])
+ if self.indices is not None:
+ pm = self["ParticleMassMsun"]
+ cx = (self["particle_position_x"] - c_vec[0])
+ cy = (self["particle_position_y"] - c_vec[1])
+ cz = (self["particle_position_z"] - c_vec[2])
+ com = na.array([v-na.floor(v) for v in [cx,cy,cz]])
+ my_mass = pm.sum()
+ my_com = ((com*pm).sum(axis=1)/my_mass + c_vec) * my_mass
+ else:
+ my_mass = 0.
+ my_com = na.array([0.,0.,0.])
+ global_mass = self._mpi_allsum(my_mass)
+ global_com = self._mpi_allsum(my_com)
+ return global_com / global_mass
+
+ def total_mass(self):
+ """
+ Returns the total mass in solar masses of the halo.
+ """
+ if self.group_total_mass is not None:
+ return self.group_total_mass
+ if self.indices is not None:
+ my_mass = self["ParticleMassMsun"].sum()
+ else:
+ my_mass = 0.
+ global_mass = self._mpi_allsum(float(my_mass))
+ return global_mass
+
+ def bulk_velocity(self):
+ """
+ Returns the mass-weighted average velocity.
+ """
+ if self.bulk_vel is not None:
+ return self.bulk_vel
+ # Unf. this cannot be reasonably computed inside of parallelHOP because
+ # we don't pass velocities in.
+ if self.indices is not None:
+ pm = self["ParticleMassMsun"]
+ vx = (self["particle_velocity_x"] * pm).sum()
+ vy = (self["particle_velocity_y"] * pm).sum()
+ vz = (self["particle_velocity_z"] * pm).sum()
+ pm = pm.sum()
+ else:
+ pm = 0.
+ vx = 0.
+ vy = 0.
+ vz = 0.
+ bv = na.array([vx,vy,vz,pm])
+ global_bv = self._mpi_allsum(bv)
+ return global_bv[:3]/global_bv[3]
+
+ def maximum_radius(self, center_of_mass=True):
+ """
+ Returns the maximum radius in the halo for all particles,
+ either from the point of maximum density or from the (default)
+ *center_of_mass*.
+ """
+ if self.max_radius is not None:
+ return self.max_radius
+ if center_of_mass: center = self.center_of_mass()
+ else: center = self.maximum_density_location()
+ if self.indices is not None:
+ rx = na.abs(self["particle_position_x"]-center[0])
+ ry = na.abs(self["particle_position_y"]-center[1])
+ rz = na.abs(self["particle_position_z"]-center[2])
+ r = na.sqrt(na.minimum(rx, 1.0-rx)**2.0
+ + na.minimum(ry, 1.0-ry)**2.0
+ + na.minimum(rz, 1.0-rz)**2.0)
+ my_max = r.max()
+
+ else:
+ my_max = 0.
+ return self._mpi_allmax(my_max)
+
+ def get_size(self):
+ if self.size is not None:
+ return self.size
+ if self.indices is not None:
+ my_size = self.indices.size
+ else:
+ my_size = 0
+ global_size = self._mpi_allsum(my_size)
+ return global_size
+
+
class FOFHalo(Halo):
def center_of_mass(self):
@@ -198,8 +346,18 @@
else: ii = slice(None)
self.particle_fields = {}
for field in self._fields:
- tot_part = self._data_source[field].size
- self.particle_fields[field] = self._data_source[field][ii]
+ if ytcfg.getboolean("yt","inline") == False:
+ tot_part = self._data_source.particles[field].size
+ if field == "particle_index":
+ self.particle_fields[field] = self._data_source.particles[field][ii].astype('int64')
+ else:
+ self.particle_fields[field] = self._data_source.particles[field][ii].astype('float64')
+ else:
+ tot_part = self._data_source[field].size
+ if field == "particle_index":
+ self.particle_fields[field] = self._data_source[field][ii].astype('int64')
+ else:
+ self.particle_fields[field] = self._data_source[field][ii].astype('float64')
self._base_indices = na.arange(tot_part)[ii]
def __get_dm_indices(self):
@@ -212,6 +370,7 @@
else:
mylog.warning("No particle_type, no creation_time, so not distinguishing.")
return slice(None)
+
def _parse_output(self):
unique_ids = na.unique(self.tags)
@@ -341,6 +500,28 @@
f.flush()
f.close()
+ def write_particle_lists_txt(self, prefix, fp=None):
+ """
+ Write out the location of halo data in hdf5 files to *prefix*.
+ """
+ if hasattr(fp, 'write'):
+ f = fp
+ else:
+ f = open("%s.txt" % prefix,"w")
+ for group in self:
+ if group.tasks is not None:
+ fn = ""
+ for task in group.tasks:
+ fn += "%s.h5 " % self._get_filename(prefix, rank=task)
+ elif self._distributed:
+ fn = "%s.h5" % self._get_filename(prefix, rank=group._owner)
+ else:
+ fn = "%s.h5" % self._get_filename(prefix)
+ gn = "Halo%08i" % (group.id)
+ f.write("%s %s\n" % (gn, fn))
+ f.flush()
+ f.close()
+
class HOPHaloList(HaloList):
_name = "HOP"
@@ -393,8 +574,168 @@
def write_out(self, filename="FOFAnalysis.out"):
HaloList.write_out(self, filename)
-class GenericHaloFinder(ParallelAnalysisInterface):
- def __init__(self, pf, dm_only=True, padding=0.02):
+class parallelHOPHaloList(HaloList,ParallelAnalysisInterface):
+ _name = "parallelHOP"
+ _halo_class = parallelHOPHalo
+ _fields = ["particle_position_%s" % ax for ax in 'xyz'] + \
+ ["ParticleMassMsun", "particle_index"]
+
+ def __init__(self, data_source, padding, num_neighbors, bounds, total_mass,
+ period, threshold=160.0, dm_only=True, rearrange=True):
+ """
+ Run hop on *data_source* with a given density *threshold*. If
+ *dm_only* is set, only run it on the dark matter particles, otherwise
+ on all particles. Returns an iterable collection of *HopGroup* items.
+ """
+ self.threshold = threshold
+ self.num_neighbors = num_neighbors
+ self.bounds = bounds
+ self.total_mass = total_mass
+ self.rearrange = rearrange
+ self.period = period
+ self._data_source = data_source
+ mylog.info("Initializing HOP")
+ HaloList.__init__(self, data_source, dm_only)
+
+ def _run_finder(self):
+ yt_counters("Reading Data")
+ obj = RunParallelHOP(self.period, self.padding,
+ self.num_neighbors, self.bounds,
+ self.particle_fields["particle_position_x"],
+ self.particle_fields["particle_position_y"],
+ self.particle_fields["particle_position_z"],
+ self.particle_fields["particle_index"],
+ self.particle_fields["ParticleMassMsun"]/self.total_mass,
+ self.threshold, rearrange=self.rearrange)
+ self.densities, self.tags = obj.density, obj.chainID
+ # I'm going to go ahead and delete self.densities because it's not
+ # actually being used. I'm not going to remove it altogether because
+ # it may be useful to someone someday.
+ del self.densities
+ self.group_count = obj.group_count
+ self.group_sizes = obj.group_sizes
+ self.CoM = obj.CoM
+ self.Tot_M = obj.Tot_M * self.total_mass
+ self.max_dens_point = obj.max_dens_point
+ self.max_radius = obj.max_radius
+ # Precompute the bulk velocity in parallel.
+ yt_counters("Precomp bulk vel.")
+ self.bulk_vel = na.zeros((self.group_count, 3), dtype='float64')
+ yt_counters("bulk vel. reading data")
+ pm = self.particle_fields["ParticleMassMsun"]
+ if ytcfg.getboolean("yt","inline") == False:
+ xv = self._data_source.particles["particle_velocity_x"][self._base_indices] * \
+ self._data_source.convert("x-velocity")
+ yv = self._data_source.particles["particle_velocity_y"][self._base_indices] * \
+ self._data_source.convert("y-velocity")
+ zv = self._data_source.particles["particle_velocity_z"][self._base_indices] * \
+ self._data_source.convert("z-velocity")
+ else:
+ xv = self._data_source["particle_velocity_x"][self._base_indices]
+ yv = self._data_source["particle_velocity_y"][self._base_indices]
+ zv = self._data_source["particle_velocity_z"][self._base_indices]
+ yt_counters("bulk vel. reading data")
+ yt_counters("bulk vel. computing")
+ select = (self.tags >= 0)
+ calc = len(na.where(select == True)[0])
+ if calc:
+ vel = na.empty((calc, 3), dtype='float64')
+ ms = pm[select]
+ vel[:,0] = xv[select] * ms
+ vel[:,1] = yv[select] * ms
+ vel[:,2] = zv[select] * ms
+ subchain = self.tags[select]
+ sort = subchain.argsort()
+ vel = vel[sort]
+ sort_subchain = subchain[sort]
+ uniq_subchain = na.unique(sort_subchain)
+ diff_subchain = na.ediff1d(sort_subchain)
+ marks = (diff_subchain > 0)
+ marks = na.arange(calc)[marks] + 1
+ marks = na.concatenate(([0], marks, [calc]))
+ for i, u in enumerate(uniq_subchain):
+ self.bulk_vel[u] = na.sum(vel[marks[i]:marks[i+1]], axis=0)
+ del select, vel, ms, subchain, sort, sort_subchain, uniq_subchain
+ del diff_subchain, marks
+ # Bring it together, and divide by the previously computed total mass
+ # of each halo.
+ self.bulk_vel = self._mpi_Allsum_double(self.bulk_vel)
+ for groupID in xrange(self.group_count):
+ self.bulk_vel[groupID] = self.bulk_vel[groupID] / self.Tot_M[groupID]
+ yt_counters("bulk vel. computing")
+ self.taskID = obj.mine
+ self.halo_taskmap = obj.halo_taskmap # A defaultdict.
+ del obj
+ yt_counters("Precomp bulk vel.")
+
+ def _parse_output(self):
+ yt_counters("Final Grouping")
+ """
+ Each task will make an entry for all groups, but it may be empty.
+ """
+ unique_ids = na.unique(self.tags)
+ counts = na.bincount((self.tags+1).tolist())
+ sort_indices = na.argsort(self.tags)
+ grab_indices = na.indices(self.tags.shape).ravel()[sort_indices]
+ del sort_indices
+ cp = 0
+ index = 0
+ # We want arrays for parallel HOP
+ self._groups = na.empty(self.group_count, dtype='object')
+ self._max_dens = na.empty((self.group_count, 4), dtype='float64')
+ for i in unique_ids:
+ if i == -1:
+ cp += counts[i+1]
+ continue
+ # If there is a gap in the unique_ids, make empty groups to
+ # fill it in.
+ while index < i:
+ self._groups[index] = self._halo_class(self, index, \
+ size=self.group_sizes[index], CoM=self.CoM[index], \
+ max_dens_point=self.max_dens_point[index], \
+ group_total_mass=self.Tot_M[index], max_radius=self.max_radius[index],
+ bulk_vel=self.bulk_vel[index], tasks=self.halo_taskmap[index])
+ # I don't own this halo
+ self._do_not_claim_object(self._groups[index])
+ self._max_dens[index] = [self.max_dens_point[index][0], self.max_dens_point[index][1], \
+ self.max_dens_point[index][2], self.max_dens_point[index][3]]
+ index += 1
+ cp_c = cp + counts[i+1]
+ group_indices = grab_indices[cp:cp_c]
+ self._groups[index] = self._halo_class(self, i, group_indices, \
+ size=self.group_sizes[i], CoM=self.CoM[i], \
+ max_dens_point=self.max_dens_point[i], \
+ group_total_mass=self.Tot_M[i], max_radius=self.max_radius[i],
+ bulk_vel=self.bulk_vel[i], tasks=self.halo_taskmap[index])
+ # This halo may be owned by many, including this task
+ self._claim_object(self._groups[index])
+ self._max_dens[index] = [self.max_dens_point[i][0], self.max_dens_point[i][1], \
+ self.max_dens_point[i][2], self.max_dens_point[i][3]]
+ cp += counts[i+1]
+ index += 1
+ # If there are missing groups at the end, add them.
+ while index < self.group_count:
+ self._groups[index] = self._halo_class(self, index, \
+ size=self.group_sizes[index], CoM=self.CoM[index], \
+ max_dens_point=self.max_dens_point[i], \
+ group_total_mass=self.Tot_M[index], max_radius=self.max_radius[index],
+ bulk_vel=self.bulk_vel[index], tasks=self.halo_taskmap[index])
+ self._do_not_claim_object(self._groups[index])
+ self._max_dens[index] = [self.max_dens_point[index][0], self.max_dens_point[index][1], \
+ self.max_dens_point[index][2], self.max_dens_point[index][3]]
+ index += 1
+ # Clean up
+ del self.max_dens_point, self.Tot_M, self.max_radius, self.bulk_vel
+ del self.halo_taskmap, self.tags
+
+ def __len__(self):
+ return self.group_count
+
+ def write_out(self, filename="parallelHopAnalysis.out"):
+ HaloList.write_out(self, filename)
+
+class GenericHaloFinder(HaloList, ParallelAnalysisInterface):
+ def __init__(self, pf, dm_only=True, padding=0.0):
self.pf = pf
self.hierarchy = pf.h
self.center = (pf["DomainRightEdge"] + pf["DomainLeftEdge"])/2.0
@@ -475,10 +816,14 @@
arr[arr > RE[i]+self.padding] -= dw[i]
def write_out(self, filename):
- self._data_source.get_data(["particle_velocity_%s" % ax for ax in 'xyz'])
+ #self._data_source.get_data(["particle_velocity_%s" % ax for ax in 'xyz'])
f = self._write_on_root(filename)
HaloList.write_out(self, f)
+ def write_particle_lists_txt(self, prefix):
+ f = self._write_on_root("%s.txt" % prefix)
+ HaloList.write_particle_lists_txt(self, prefix, fp=f)
+
@parallel_blocking_call
def write_particle_lists(self, prefix):
fn = "%s.h5" % self._get_filename(prefix)
@@ -487,15 +832,161 @@
if not self._is_mine(halo): continue
halo.write_particle_list(f)
+class parallelHF(GenericHaloFinder, parallelHOPHaloList):
+ def __init__(self, pf, threshold=160, dm_only=True, resize=True, rearrange=True,\
+ fancy_padding=True, safety=1.5):
+ GenericHaloFinder.__init__(self, pf, dm_only, padding=0.0)
+ self.padding = 0.0
+ self.num_neighbors = 65
+ self.safety = safety
+ period = pf["DomainRightEdge"] - pf["DomainLeftEdge"]
+ # Cut up the volume evenly initially, with no padding.
+ padded, LE, RE, self._data_source = self._partition_hierarchy_3d(padding=self.padding)
+ # also get the total mass of particles
+ yt_counters("Reading Data")
+ # Adaptive subregions by bisection.
+ ds_names = ["particle_position_x","particle_position_y","particle_position_z"]
+ if ytcfg.getboolean("yt","inline") == False and \
+ resize and self._mpi_get_size() is not None:
+ cut_list = self._partition_hierarchy_3d_bisection_list()
+ for i,cut in enumerate(cut_list):
+ dim = cut[0]
+ if i == 0:
+ width = pf["DomainRightEdge"][dim] - pf["DomainLeftEdge"][dim]
+ new_LE = pf["DomainLeftEdge"]
+ else:
+ new_LE, new_RE = new_top_bounds
+ width = new_RE[dim] - new_LE[dim]
+ if ytcfg.getboolean("yt","inline") == False:
+ data = self._data_source.particles[ds_names[dim]]
+ else:
+ data = self._data_source[ds_names[dim]]
+ if i == 0:
+ local_parts = data.size
+ n_parts = self._mpi_allsum(local_parts)
+ min = self._mpi_allmin(local_parts)
+ max = self._mpi_allmax(local_parts)
+ mylog.info("Initial distribution: (min,max,tot) (%d,%d,%d)" \
+ % (min, max, n_parts))
+ num_bins = 1000
+ bin_width = float(width)/float(num_bins)
+ bins = na.arange(num_bins+1, dtype='float64') * bin_width + new_LE[dim]
+ counts, bins = na.histogram(data, bins, new=True)
+ if i == 0:
+ new_group, new_comm, LE, RE, new_top_bounds, new_cc, self._data_source= \
+ self._partition_hierarchy_3d_bisection(dim, bins, counts, top_bounds = None,\
+ old_group = None, old_comm = None, cut=cut)
+ else:
+ new_group, new_comm, LE, RE, new_top_bounds, new_cc, self._data_source = \
+ self._partition_hierarchy_3d_bisection(dim, bins, counts, top_bounds = new_top_bounds,\
+ old_group = new_group, old_comm = new_comm, cut=cut, old_cc=new_cc)
+ del bins, counts
+ # If this isn't parallel, define the region as an AMRRegionStrict so
+ # particle IO works.
+ if self._mpi_get_size() == None or self._mpi_get_size() == 1:
+ self._data_source = self.hierarchy.periodic_region_strict([0.5]*3, LE, RE)
+ # get the average spacing between particles for this region
+ # The except is for the serial case, where the full box is what we want.
+ if ytcfg.getboolean("yt","inline") == False:
+ data = self._data_source.particles["particle_position_x"]
+ else:
+ data = self._data_source["particle_position_x"]
+ try:
+ l = self._data_source.right_edge - self._data_source.left_edge
+ except AttributeError:
+ l = pf["DomainRightEdge"] - pf["DomainLeftEdge"]
+ vol = l[0] * l[1] * l[2]
+ full_vol = vol
+ if not fancy_padding:
+ avg_spacing = (float(vol) / data.size)**(1./3.)
+ # padding is a function of inter-particle spacing, this is an
+ # approximation, but it's OK with the safety factor
+ padding = (self.num_neighbors)**(1./3.) * self.safety * avg_spacing
+ self.padding = (na.ones(3,dtype='float64')*padding, na.ones(3,dtype='float64')*padding)
+ mylog.info('padding %s avg_spacing %f vol %f local_parts %d' % \
+ (str(self.padding), avg_spacing, vol, data.size))
+ # Another approach to padding, perhaps more accurate.
+ elif fancy_padding and self._distributed:
+ LE_padding, RE_padding = na.empty(3,dtype='float64'), na.empty(3,dtype='float64')
+ for dim in xrange(3):
+ if ytcfg.getboolean("yt","inline") == False:
+ data = self._data_source.particles[ds_names[dim]]
+ else:
+ data = self._data_source[ds_names[dim]]
+ num_bins = 1000
+ width = self._data_source.right_edge[dim] - self._data_source.left_edge[dim]
+ area = (self._data_source.right_edge[(dim+1)%3] - self._data_source.left_edge[(dim+1)%3]) * \
+ (self._data_source.right_edge[(dim+2)%3] - self._data_source.left_edge[(dim+2)%3])
+ bin_width = float(width)/float(num_bins)
+ bins = na.arange(num_bins+1, dtype='float64') * bin_width + self._data_source.left_edge[dim]
+ counts, bins = na.histogram(data, bins, new=True)
+ # left side.
+ start = 0
+ count = counts[0]
+ while count < self.num_neighbors:
+ start += 1
+ count += counts[start]
+ # Get the avg spacing in just this boundary.
+ vol = area * (bins[start+1] - bins[0])
+ avg_spacing = (float(vol) / count)**(1./3.)
+ LE_padding[dim] = (self.num_neighbors)**(1./3.) * self.safety * avg_spacing
+ # right side.
+ start = -1
+ count = counts[-1]
+ while count < self.num_neighbors:
+ start -= 1
+ count += counts[start]
+ vol = area * (bins[-1] - bins[start-1])
+ avg_spacing = (float(vol) / count)**(1./3.)
+ RE_padding[dim] = (self.num_neighbors)**(1./3.) * self.safety * avg_spacing
+ self.padding = (LE_padding, RE_padding)
+ del bins, counts
+ mylog.info('fancy_padding %s avg_spacing %f full_vol %f local_parts %d %s' % \
+ (str(self.padding), avg_spacing, full_vol, data.size, str(self._data_source)))
+ # Now we get the full box mass after we have the final composition of
+ # subvolumes.
+ if ytcfg.getboolean("yt","inline") == False:
+ total_mass = self._mpi_allsum((self._data_source.particles["ParticleMassMsun"].astype('float64')).sum())
+ else:
+ total_mass = self._mpi_allsum((self._data_source["ParticleMassMsun"].astype('float64')).sum())
+ if not self._distributed:
+ self.padding = (na.zeros(3,dtype='float64'), na.zeros(3,dtype='float64'))
+ self.bounds = (LE, RE)
+ (LE_padding, RE_padding) = self.padding
+ parallelHOPHaloList.__init__(self, self._data_source, self.padding, \
+ self.num_neighbors, self.bounds, total_mass, period, \
+ threshold=threshold, dm_only=dm_only, rearrange=rearrange)
+ self._join_halolists()
+ yt_counters("Final Grouping")
+
+ def _join_halolists(self):
+ gs = -self.group_sizes.copy()
+ Cx = self.CoM[:,0].copy()
+ indexes = na.arange(self.group_count)
+ sorted = na.asarray([indexes[i] for i in na.lexsort([indexes, Cx, gs])])
+ del indexes, Cx, gs
+ self._groups = self._groups[sorted]
+ self._max_dens = self._max_dens[sorted]
+ for i in xrange(self.group_count):
+ self._groups[i].id = i
+ del sorted, self.group_sizes, self.CoM
+
+
class HOPHaloFinder(GenericHaloFinder, HOPHaloList):
def __init__(self, pf, threshold=160, dm_only=True, padding=0.02):
GenericHaloFinder.__init__(self, pf, dm_only, padding)
# do it once with no padding so the total_mass is correct (no duplicated particles)
self.padding = 0.0
- padded, LE, RE, data_source = self._partition_hierarchy_3d(padding=self.padding)
+ padded, LE, RE, self._data_source = self._partition_hierarchy_3d(padding=self.padding)
# For scaling the threshold, note that it's a passthrough
- total_mass = self._mpi_allsum(data_source["ParticleMassMsun"].sum())
+ if dm_only:
+ select = self._data_source["creation_time"] > 0
+ total_mass = self._mpi_allsum((self._data_source["ParticleMassMsun"][select]).sum())
+ sub_mass = (self._data_source["ParticleMassMsun"][select]).sum()
+ else:
+ total_mass = self._mpi_allsum(self._data_source["ParticleMassMsun"].sum())
+ sub_mass = self._data_source["ParticleMassMsun"].sum()
# MJT: Note that instead of this, if we are assuming that the particles
# are all on different processors, we should instead construct an
# object representing the entire domain and sum it "lazily" with
@@ -504,8 +995,8 @@
padded, LE, RE, self._data_source = self._partition_hierarchy_3d(padding=self.padding)
self.bounds = (LE, RE)
# reflect particles around the periodic boundary
- self._reposition_particles((LE, RE))
- sub_mass = self._data_source["ParticleMassMsun"].sum()
+ #self._reposition_particles((LE, RE))
+ #sub_mass = self._data_source["ParticleMassMsun"].sum()
HOPHaloList.__init__(self, self._data_source, threshold*total_mass/sub_mass, dm_only)
self._parse_halolist(total_mass/sub_mass)
self._join_halolists()
@@ -527,7 +1018,7 @@
padded, LE, RE, self._data_source = self._partition_hierarchy_3d(padding=self.padding)
self.bounds = (LE, RE)
# reflect particles around the periodic boundary
- self._reposition_particles((LE, RE))
+ #self._reposition_particles((LE, RE))
# here is where the FOF halo finder is run
FOFHaloList.__init__(self, self._data_source, link * avg_spacing, dm_only)
self._parse_halolist(1.)
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Fri Nov 20 20:27:29 2009
@@ -27,93 +27,114 @@
from yt.funcs import *
import string, re, gc, time
import cPickle
-#import yt.enki
+from itertools import chain, izip
-_data_style_funcs = \
- { 'enzo_hdf4': (readDataHDF4,readAllDataHDF4, getFieldsHDF4, readDataSliceHDF4,
- getExceptionHDF4, DataQueueHDF4),
- 'enzo_hdf4_2d': (readDataHDF4, readAllDataHDF4, getFieldsHDF4, readDataSliceHDF4_2D,
- getExceptionHDF4, DataQueueHDF4_2D),
- 'enzo_hdf5': (readDataHDF5, readAllDataHDF5, getFieldsHDF5, readDataSliceHDF5,
- getExceptionHDF5, DataQueueHDF5),
- 'enzo_packed_3d': (readDataPacked, readAllDataPacked, getFieldsPacked, readDataSlicePacked,
- getExceptionHDF5, DataQueuePackedHDF5),
- 'orion_native': (readDataNative, readAllDataNative, None, readDataSliceNative,
- getExceptionHDF5, DataQueueNative), \
- 'enzo_inline': (readDataInMemory, readAllDataInMemory, getFieldsInMemory, readDataSliceInMemory,
- getExceptionInMemory, DataQueueInMemory),
- 'enzo_packed_2d': (readDataPacked, readAllDataPacked, getFieldsPacked, readDataSlicePacked2D,
- getExceptionHDF5, DataQueuePacked2D),
- 'enzo_packed_1d': (readDataPacked, readAllDataPacked, getFieldsPacked, readDataSlicePacked1D,
- getExceptionHDF5, DataQueuePacked1D),
- }
-
-class AMRHierarchy:
- _data_mode = None # Default
- def __init__(self, pf):
+class AMRHierarchy(ObjectFindingMixin, ParallelAnalysisInterface):
+ float_type = 'float64'
+
+ def __init__(self, pf, data_style):
self.parameter_file = weakref.proxy(pf)
- self._max_locations = {}
- self._data_file = None
+ self.pf = self.parameter_file
+
+ self._initialize_state_variables()
+
+ mylog.debug("Initializing data storage.")
+ self._initialize_data_storage()
+
+ mylog.debug("Counting grids.")
+ self._count_grids()
+
+ # Must be defined in subclass
+ mylog.debug("Setting up classes.")
self._setup_classes()
- self._initialize_grids()
- # For use with derived quantities depending on centers
- # Although really, I think perhaps we should take a closer look
- # at how "center" is used.
- self.center = None
+ mylog.debug("Counting grids.")
+ self._initialize_grid_arrays()
+
+ mylog.debug("Parsing hierarchy.")
+ self._parse_hierarchy()
+
+ mylog.debug("Constructing grid objects.")
+ self._populate_grid_objects()
+ mylog.debug("Initializing data grid data IO")
+ self._setup_data_io()
+
+ mylog.debug("Detecting fields.")
+ self.field_list = []
+ self._detect_fields()
+
+ mylog.debug("Adding unknown detected fields")
+ self._setup_unknown_fields()
+
+ mylog.debug("Setting up derived fields")
+ self._setup_derived_fields()
+
+ mylog.debug("Re-examining hierarchy")
self._initialize_level_stats()
- mylog.debug("Initializing data file")
- self._initialize_data_file()
- mylog.debug("Populating hierarchy")
- self._populate_hierarchy()
- mylog.debug("Done populating hierarchy")
- self._add_detected_fields()
- mylog.debug("Done adding auto-detected fields")
+ def _get_parameters(self):
+ return self.parameter_file.parameters
+ parameters=property(_get_parameters)
- def _initialize_grids(self):
- mylog.debug("Allocating memory for %s grids", self.num_grids)
- self.gridDimensions = na.zeros((self.num_grids,3), 'int32')
- self.gridStartIndices = na.zeros((self.num_grids,3), 'int32')
- self.gridEndIndices = na.zeros((self.num_grids,3), 'int32')
- self.gridLeftEdge = na.zeros((self.num_grids,3), self.float_type)
- self.gridRightEdge = na.zeros((self.num_grids,3), self.float_type)
- self.gridLevels = na.zeros((self.num_grids,1), 'int32')
- self.gridDxs = na.zeros((self.num_grids,1), self.float_type)
- self.gridDys = na.zeros((self.num_grids,1), self.float_type)
- self.gridDzs = na.zeros((self.num_grids,1), self.float_type)
- self.gridTimes = na.zeros((self.num_grids,1), 'float64')
- self.gridNumberOfParticles = na.zeros((self.num_grids,1))
- mylog.debug("Done allocating")
- mylog.debug("Creating grid objects")
- self.grids = na.array([self.grid(i+1) for i in xrange(self.num_grids)])
- self.gridReverseTree = [-1] * self.num_grids
- self.gridTree = [ [] for i in range(self.num_grids)]
- mylog.debug("Done creating grid objects")
+ def select_grids(self, level):
+ """
+ Returns an array of grids at *level*.
+ """
+ return self.grids[self.grid_levels.flat == level]
- def _initialize_level_stats(self):
- # Now some statistics:
- # 0 = number of grids
- # 1 = number of cells
- # 2 = blank
- desc = {'names': ['numgrids','numcells','level'],
- 'formats':['Int32']*3}
- self.level_stats = blankRecordArray(desc, MAXLEVEL)
- self.level_stats['level'] = [i for i in range(MAXLEVEL)]
- self.level_stats['numgrids'] = [0 for i in range(MAXLEVEL)]
- self.level_stats['numcells'] = [0 for i in range(MAXLEVEL)]
+ def _initialize_state_variables(self):
+ self._parallel_locking = False
+ self._data_file = None
+ self._data_mode = None
+ self._max_locations = {}
+ self.num_grids = None
- def __setup_filemap(self, grid):
- if not self.data_style == 'enzo_packed_3d':
- return
- self.cpu_map[grid.filename].append(grid)
+ def _initialize_grid_arrays(self):
+ mylog.debug("Allocating arrays for %s grids", self.num_grids)
+ self.grid_dimensions = na.ones((self.num_grids,3), 'int32')
+ self.grid_left_edge = na.zeros((self.num_grids,3), self.float_type)
+ self.grid_right_edge = na.ones((self.num_grids,3), self.float_type)
+ self.grid_levels = na.zeros((self.num_grids,1), 'int32')
+ self.grid_particle_count = na.zeros((self.num_grids,1), 'int32')
+
+ def _setup_classes(self, dd):
+ # Called by subclass
+ self.object_types = []
+ self.objects = []
+ for name, cls in sorted(data_object_registry.items()):
+ cname = cls.__name__
+ if cname.endswith("Base"): cname = cname[:-4]
+ self._add_object_class(name, cname, cls, dd)
+ self.object_types.sort()
+
+ # Now all the object related stuff
+
+ def all_data(self, find_max=False):
+ pf = self.parameter_file
+ if find_max: c = self.find_max("Density")[1]
+ else: c = (pf["DomainRightEdge"] + pf["DomainLeftEdge"])/2.0
+ return self.region(c,
+ pf["DomainLeftEdge"], pf["DomainRightEdge"])
+
+ def clear_all_data(self):
+ """
+ This routine clears all the data currently being held onto by the grids
+ and the data io handler.
+ """
+ for g in self.grids: g.clear_data()
+ self.io.queue.clear()
- def _initialize_data_file(self):
+ def _get_data_reader_dict(self):
+ dd = { 'pf' : self.parameter_file, # Already weak
+ 'hierarchy': weakref.proxy(self) }
+ return dd
+
+ def _initialize_data_storage(self):
if not ytcfg.getboolean('lagos','serialize'): return
if os.path.isfile(os.path.join(self.directory,
- "%s.yt" % self["CurrentTimeIdentifier"])):
- fn = os.path.join(self.directory,"%s.yt" % self["CurrentTimeIdentifier"])
+ "%s.yt" % self.pf["CurrentTimeIdentifier"])):
+ fn = os.path.join(self.directory,"%s.yt" % self.pf["CurrentTimeIdentifier"])
else:
fn = os.path.join(self.directory,
"%s.yt" % self.parameter_file.basename)
@@ -130,33 +151,13 @@
self.__data_filename = fn
self._data_file = h5py.File(fn, self._data_mode)
- def clear_all_data(self):
- """
- This routine clears all the data currently being held onto by the grids
- and the data queue.
- """
- for g in self.grids: g.clear_data()
- self.queue.queue.clear()
-
@parallel_root_only
def __create_data_file(self, fn):
f = h5py.File(fn, 'a')
f.close()
- def _setup_data_queue(self):
- self.queue = _data_style_funcs[self.data_style][5]()
-
- def _setup_grid_corners(self):
- self.gridCorners = na.array([ # Unroll!
- [self.gridLeftEdge[:,0], self.gridLeftEdge[:,1], self.gridLeftEdge[:,2]],
- [self.gridRightEdge[:,0], self.gridLeftEdge[:,1], self.gridLeftEdge[:,2]],
- [self.gridRightEdge[:,0], self.gridRightEdge[:,1], self.gridLeftEdge[:,2]],
- [self.gridRightEdge[:,0], self.gridRightEdge[:,1], self.gridRightEdge[:,2]],
- [self.gridLeftEdge[:,0], self.gridRightEdge[:,1], self.gridRightEdge[:,2]],
- [self.gridLeftEdge[:,0], self.gridLeftEdge[:,1], self.gridRightEdge[:,2]],
- [self.gridRightEdge[:,0], self.gridLeftEdge[:,1], self.gridRightEdge[:,2]],
- [self.gridLeftEdge[:,0], self.gridRightEdge[:,1], self.gridLeftEdge[:,2]],
- ], dtype='float64')
+ def _setup_data_io(self):
+ self.io = io_registry[self.data_style]()
def _save_data(self, array, node, name, set_attr=None, force=False, passthrough = False):
"""
@@ -242,7 +243,7 @@
return None
full_name = "%s/%s" % (node, name)
- return self._data_file[full_name]
+ return self._data_file[full_name][:]
def _close_data_file(self):
if self._data_file:
@@ -250,28 +251,8 @@
del self._data_file
self._data_file = None
- def _add_object_class(self, name, class_name, base, dd):
- self.object_types.append(name)
- obj = classobj(class_name, (base,), dd)
- setattr(self, name, obj)
-
- def _setup_classes(self, dd):
- self.object_types = []
- self.objects = []
- for name, cls in sorted(data_object_registry.items()):
- cname = cls.__name__
- if cname.endswith("Base"): cname = cname[:-4]
- self._add_object_class(name, cname, cls, dd)
- self.object_types.sort()
-
- def all_data(self, find_max=False):
- pf = self.parameter_file
- if find_max: c = self.find_max("Density")[1]
- else: c = (pf["DomainRightEdge"] + pf["DomainLeftEdge"])/2.0
- return self.region(c,
- pf["DomainLeftEdge"], pf["DomainRightEdge"])
-
def _deserialize_hierarchy(self, harray):
+ # THIS IS BROKEN AND NEEDS TO BE FIXED
mylog.debug("Cached entry found.")
self.gridDimensions[:] = harray[:,0:3]
self.gridStartIndices[:] = harray[:,3:6]
@@ -282,123 +263,57 @@
self.gridTimes[:] = harray[:,16:17]
self.gridNumberOfParticles[:] = harray[:,17:18]
- def _get_data_reader_dict(self):
- dd = { 'readDataFast' : _data_style_funcs[self.data_style][0],
- 'readAllData' : _data_style_funcs[self.data_style][1],
- 'getFields' : _data_style_funcs[self.data_style][2],
- 'readDataSlice' : _data_style_funcs[self.data_style][3],
- '_read_data' : _data_style_funcs[self.data_style][0],
- '_read_all_data' : _data_style_funcs[self.data_style][1],
- '_read_field_names' : _data_style_funcs[self.data_style][2],
- '_read_data_slice' : _data_style_funcs[self.data_style][3],
- '_read_exception' : _data_style_funcs[self.data_style][4](),
- 'pf' : self.parameter_file, # Already weak
- 'hierarchy': weakref.proxy(self) }
- return dd
-
def get_smallest_dx(self):
"""
Returns (in code units) the smallest cell size in the simulation.
"""
- return self.gridDxs.min()
-
- def find_ray_grids(self, coord, axis):
- """
- Returns the (objects, indices) of grids that an (x,y) ray intersects
- along *axis*
- """
- # Let's figure out which grids are on the slice
- mask=na.ones(self.num_grids)
- # So if gRE > coord, we get a mask, if not, we get a zero
- # if gLE > coord, we get a zero, if not, mask
- # Thus, if the coordinate is between the two edges, we win!
- na.choose(na.greater(self.gridRightEdge[:,x_dict[axis]],coord[0]),(0,mask),mask)
- na.choose(na.greater(self.gridLeftEdge[:,x_dict[axis]],coord[0]),(mask,0),mask)
- na.choose(na.greater(self.gridRightEdge[:,y_dict[axis]],coord[1]),(0,mask),mask)
- na.choose(na.greater(self.gridLeftEdge[:,y_dict[axis]],coord[1]),(mask,0),mask)
- ind = na.where(mask == 1)
- return self.grids[ind], ind
-
- @time_execution
- def find_max(self, field, finestLevels = True):
- """
- Returns (value, center) of location of maximum for a given field.
- """
- if (field, finestLevels) in self._max_locations:
- return self._max_locations[(field, finestLevels)]
- mg, mc, mv, pos = self.find_max_cell_location(field, finestLevels)
- self._max_locations[(field, finestLevels)] = (mv, pos)
- return mv, pos
- findMax = deprecate(find_max)
-
- def find_max_cell_location(self, field, finestLevels = True):
- if finestLevels:
- gi = (self.gridLevels >= self.max_level - NUMTOCHECK).ravel()
- source = self.grid_collection([0.0]*3,
- self.grids[gi])
- else:
- source = self.all_data()
- mylog.debug("Searching %s grids for maximum value of %s",
- len(source._grids), field)
- max_val, maxi, mx, my, mz, mg = source.quantities["MaxLocation"](
- field, lazy_reader=True)
- max_grid = self.grids[mg]
- mc = na.unravel_index(maxi, max_grid.ActiveDimensions)
- mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s %s", \
- max_val, mx, my, mz, max_grid, max_grid.Level, mc)
- self.parameters["Max%sValue" % (field)] = max_val
- self.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
- return max_grid, mc, max_val, na.array((mx,my,mz), dtype='float64')
-
- @time_execution
- def find_min(self, field):
- """
- Returns (value, center) of location of minimum for a given field
- """
- gI = na.where(self.gridLevels >= 0) # Slow but pedantic
- minVal = 1e100
- for grid in self.grids[gI[0]]:
- mylog.debug("Checking %s (level %s)", grid.id, grid.Level)
- val, coord = grid.find_min(field)
- if val < minVal:
- minCoord = coord
- minVal = val
- minGrid = grid
- mc = na.array(minCoord)
- pos=minGrid.get_position(mc)
- mylog.info("Min Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s", \
- minVal, pos[0], pos[1], pos[2], minGrid, minGrid.Level)
- self.center = pos
- # This probably won't work for anyone else
- self.binkVelocity = (minGrid["x-velocity"][minCoord], \
- minGrid["y-velocity"][minCoord], \
- minGrid["z-velocity"][minCoord])
- self.parameters["Min%sValue" % (field)] = minVal
- self.parameters["Min%sPos" % (field)] = "%s" % (pos)
- return minVal, pos
+ return self.select_grids(self.grid_levels.max())[0].dds[0]
- def _get_parameters(self):
- return self.parameter_file.parameters
- parameters=property(_get_parameters)
+ def _add_object_class(self, name, class_name, base, dd):
+ self.object_types.append(name)
+ obj = classobj(class_name, (base,), dd)
+ setattr(self, name, obj)
- def __getitem__(self, item):
- return self.parameter_file[item]
- def select_grids(self, level):
- """
- Returns an array of grids at *level*.
- """
- return self.grids[self._select_level(level)]
+ def _initialize_level_stats(self):
+ # Now some statistics:
+ # 0 = number of grids
+ # 1 = number of cells
+ # 2 = blank
+ desc = {'names': ['numgrids','numcells','level'],
+ 'formats':['Int32']*3}
+ self.level_stats = blankRecordArray(desc, MAXLEVEL)
+ self.level_stats['level'] = [i for i in range(MAXLEVEL)]
+ self.level_stats['numgrids'] = [0 for i in range(MAXLEVEL)]
+ self.level_stats['numcells'] = [0 for i in range(MAXLEVEL)]
+ for level in xrange(self.max_level+1):
+ self.level_stats[level]['numgrids'] = na.sum(self.grid_levels == level)
+ li = (self.grid_levels[:,0] == level)
+ self.level_stats[level]['numcells'] = self.grid_dimensions[li,:].prod(axis=1).sum()
+
+ @property
+ def grid_corners(self):
+ return na.array([
+ [self.grid_left_edge[:,0], self.grid_left_edge[:,1], self.grid_left_edge[:,2]],
+ [self.grid_right_edge[:,0], self.grid_left_edge[:,1], self.grid_left_edge[:,2]],
+ [self.grid_right_edge[:,0], self.grid_right_edge[:,1], self.grid_left_edge[:,2]],
+ [self.grid_right_edge[:,0], self.grid_right_edge[:,1], self.grid_right_edge[:,2]],
+ [self.grid_left_edge[:,0], self.grid_right_edge[:,1], self.grid_right_edge[:,2]],
+ [self.grid_left_edge[:,0], self.grid_left_edge[:,1], self.grid_right_edge[:,2]],
+ [self.grid_right_edge[:,0], self.grid_left_edge[:,1], self.grid_right_edge[:,2]],
+ [self.grid_left_edge[:,0], self.grid_right_edge[:,1], self.grid_left_edge[:,2]],
+ ], dtype='float64')
def print_stats(self):
"""
Prints out (stdout) relevant information about the simulation
"""
- for i in xrange(MAXLEVEL):
- if (self.level_stats['numgrids'][i]) == 0:
+ for level in xrange(MAXLEVEL):
+ if (self.level_stats['numgrids'][level]) == 0:
break
print "% 3i\t% 6i\t% 11i" % \
- (i, self.level_stats['numgrids'][i], self.level_stats['numcells'][i])
- dx = self.gridDxs[self.levelIndices[i][0]]
+ (level, self.level_stats['numgrids'][level],
+ self.level_stats['numcells'][level])
+ dx = self.select_grids(level)[0].dds[0]
print "-" * 28
print " \t% 6i\t% 11i" % (self.level_stats['numgrids'].sum(), self.level_stats['numcells'].sum())
print "\n"
@@ -406,9 +321,9 @@
print "z = %0.8f" % (self["CosmologyCurrentRedshift"])
except:
pass
- t_s = self["InitialTime"] * self["Time"]
+ t_s = self.pf["InitialTime"] * self.pf["Time"]
print "t = %0.8e = %0.8e s = %0.8e years" % \
- (self["InitialTime"], \
+ (self.pf["InitialTime"], \
t_s, t_s / (365*24*3600.0) )
print "\nSmallest Cell:"
u=[]
@@ -418,270 +333,25 @@
for unit in u:
print "\tWidth: %0.3e %s" % (dx*unit[0], unit[1])
- def find_point(self, coord):
- """
- Returns the (objects, indices) of grids containing an (x,y,z) point
- """
- mask=na.ones(self.num_grids)
- for i in xrange(len(coord)):
- na.choose(na.greater(self.gridLeftEdge[:,i],coord[i]), (mask,0), mask)
- na.choose(na.greater(self.gridRightEdge[:,i],coord[i]), (0,mask), mask)
- ind = na.where(mask == 1)
- return self.grids[ind], ind
-
- def find_slice_grids(self, coord, axis):
- """
- Returns the (objects, indices) of grids that a slice intersects along
- *axis*
- """
- # Let's figure out which grids are on the slice
- mask=na.ones(self.num_grids)
- # So if gRE > coord, we get a mask, if not, we get a zero
- # if gLE > coord, we get a zero, if not, mask
- # Thus, if the coordinate is between the edges, we win!
- #ind = na.where( na.logical_and(self.gridRightEdge[:,axis] > coord, \
- #self.gridLeftEdge[:,axis] < coord))
- na.choose(na.greater(self.gridRightEdge[:,axis],coord),(0,mask),mask)
- na.choose(na.greater(self.gridLeftEdge[:,axis],coord),(mask,0),mask)
- ind = na.where(mask == 1)
- return self.grids[ind], ind
-
- def find_sphere_grids(self, center, radius):
- """
- Returns objects, indices of grids within a sphere
- """
- centers = (self.gridRightEdge + self.gridLeftEdge)/2.0
- long_axis = na.maximum.reduce(self.gridRightEdge - self.gridLeftEdge, 1)
- t = na.abs(centers - center)
- DW = self.parameter_file["DomainRightEdge"] \
- - self.parameter_file["DomainLeftEdge"]
- na.minimum(t, na.abs(DW-t), t)
- dist = na.sqrt(na.sum((t**2.0), axis=1))
- gridI = na.where(na.logical_and((self.gridDxs<=radius)[:,0],(dist < (radius + long_axis))) == 1)
- return self.grids[gridI], gridI
-
- def get_box_grids(self, left_edge, right_edge):
- """
- Gets back all the grids between a left edge and right edge
- """
- grid_i = na.where((na.all(self.gridRightEdge > left_edge, axis=1)
- & na.all(self.gridLeftEdge < right_edge, axis=1)) == True)
- return self.grids[grid_i], grid_i
-
- def get_periodic_box_grids(self, left_edge, right_edge):
- left_edge = na.array(left_edge)
- right_edge = na.array(right_edge)
- mask = na.zeros(self.grids.shape, dtype='bool')
- dl = self.parameters["DomainLeftEdge"]
- dr = self.parameters["DomainRightEdge"]
- db = right_edge - left_edge
- for off_x in [-1, 0, 1]:
- nle = left_edge.copy()
- nre = left_edge.copy()
- nle[0] = dl[0] + (dr[0]-dl[0])*off_x + left_edge[0]
- for off_y in [-1, 0, 1]:
- nle[1] = dl[1] + (dr[1]-dl[1])*off_y + left_edge[1]
- for off_z in [-1, 0, 1]:
- nle[2] = dl[2] + (dr[2]-dl[2])*off_z + left_edge[2]
- nre = nle + db
- g, gi = self.get_box_grids(nle, nre)
- mask[gi] = True
- return self.grids[mask], na.where(mask)
-
- @time_execution
- def export_particles_pb(self, filename, filter = 1, indexboundary = 0, fields = None, scale=1.0):
- """
- Exports all the star particles, or a subset, to pb-format *filename*
- for viewing in partiview. Filters based on particle_type=*filter*,
- particle_index>=*indexboundary*, and exports *fields*, if supplied.
- Otherwise, index, position(x,y,z). Optionally *scale* by a given
- factor before outputting.
- """
- import struct
- pbf_magic = 0xffffff98
- header_fmt = 'Iii'
- fmt = 'ifff'
- f = open(filename,"w")
- if fields:
- fmt += len(fields)*'f'
- padded_fields = string.join(fields,"\0") + "\0"
- header_fmt += "%ss" % len(padded_fields)
- args = [pbf_magic, struct.calcsize(header_fmt), len(fields), padded_fields]
- fields = ["particle_index","particle_position_x","particle_position_y","particle_position_z"] \
- + fields
- format = 'Int32,Float32,Float32,Float32' + ',Float32'*(len(fields)-4)
- else:
- args = [pbf_magic, struct.calcsize(header_fmt), 0]
- fields = ["particle_index","particle_position_x","particle_position_y","particle_position_z"]
- format = 'Int32,Float32,Float32,Float32'
- f.write(struct.pack(header_fmt, *args))
- tot = 0
- sc = na.array([1.0] + [scale] * 3 + [1.0]*(len(fields)-4))
- gI = na.where(self.gridNumberOfParticles.ravel() > 0)
- for g in self.grids[gI]:
- pI = na.where(na.logical_and((g["particle_type"] == filter),(g["particle_index"] >= indexboundary)) == 1)
- tot += pI[0].shape[0]
- toRec = []
- for field, scale in zip(fields, sc):
- toRec.append(scale*g[field][pI])
- particle_info = rec.array(toRec,formats=format)
- particle_info.tofile(f)
- f.close()
- mylog.info("Wrote %s particles to %s", tot, filename)
-
- @time_execution
- def export_boxes_pv(self, filename):
- """
- Exports the grid structure in partiview text format.
- """
- f=open(filename,"w")
- for l in xrange(self.maxLevel):
- f.write("add object g%s = l%s\n" % (l,l))
- ind = self._select_level(l)
- for i in ind:
- f.write("add box -n %s -l %s %s,%s %s,%s %s,%s\n" % \
- (i+1, self.gridLevels.ravel()[i],
- self.gridLeftEdge[i,0], self.gridRightEdge[i,0],
- self.gridLeftEdge[i,1], self.gridRightEdge[i,1],
- self.gridLeftEdge[i,2], self.gridRightEdge[i,2]))
-
- def _select_level(self, level):
- # We return a numarray of the indices of all the grids on a given level
- indices = na.where(self.gridLevels[:,0] == level)[0]
- return indices
-
- def __initialize_octree_list(self, fields):
- import DepthFirstOctree as dfo
- o_length = r_length = 0
- grids = []
- levels_finest, levels_all = defaultdict(lambda: 0), defaultdict(lambda: 0)
- for g in self.grids:
- ff = na.array([g[f] for f in fields])
- grids.append(dfo.OctreeGrid(
- g.child_index_mask.astype('int32'),
- ff.astype("float64"),
- g.LeftEdge.astype('float64'),
- g.ActiveDimensions.astype('int32'),
- na.ones(1,dtype='float64') * g.dds[0], g.Level))
- levels_all[g.Level] += g.ActiveDimensions.prod()
- levels_finest[g.Level] += g.child_mask.ravel().sum()
- g.clear_data()
- ogl = dfo.OctreeGridList(grids)
- return ogl, levels_finest, levels_all
-
- def _generate_flat_octree(self, fields):
- """
- Generates two arrays, one of the actual values in a depth-first flat
- octree array, and the other of the values describing the refinement.
- This allows for export to a code that understands this. *field* is the
- field used in the data array.
- """
- import DepthFirstOctree as dfo
- fields = ensure_list(fields)
- ogl, levels_finest, levels_all = self.__initialize_octree_list(fields)
- o_length = na.sum(levels_finest.values())
- r_length = na.sum(levels_all.values())
- output = na.zeros((o_length,len(fields)), dtype='float64')
- refined = na.zeros(r_length, dtype='int32')
- position = dfo.position()
- dfo.RecurseOctreeDepthFirst(0, 0, 0,
- ogl[0].dimensions[0],
- ogl[0].dimensions[1],
- ogl[0].dimensions[2],
- position, 1,
- output, refined, ogl)
- dd = {}
- for i,field in enumerate(fields): dd[field] = output[:,i]
- return dd, refined
-
- def _generate_levels_octree(self, fields):
- import DepthFirstOctree as dfo
- fields = ensure_list(fields) + ["Ones", "Ones"]
- ogl, levels_finest, levels_all = self.__initialize_octree_list(fields)
- o_length = na.sum(levels_finest.values())
- r_length = na.sum(levels_all.values())
- output = na.zeros((r_length,len(fields)), dtype='float64')
- genealogy = na.zeros((r_length, 3), dtype='int32') - 1 # init to -1
- corners = na.zeros((r_length, 3), dtype='float64')
- position = na.add.accumulate(
- na.array([0] + [levels_all[v] for v in
- sorted(levels_all)[:-1]], dtype='int32'))
- pp = position.copy()
- dfo.RecurseOctreeByLevels(0, 0, 0,
- ogl[0].dimensions[0],
- ogl[0].dimensions[1],
- ogl[0].dimensions[2],
- position.astype('int32'), 1,
- output, genealogy, corners, ogl)
- return output, genealogy, levels_all, levels_finest, pp, corners
-
- def _add_detected_fields(self):
- """add any extra fields in the
-
-
- """
- for field in self.field_list:
- if field in self.parameter_file.field_info: continue
- mylog.info("Adding %s to list of fields", field)
- cf = None
- if self.parameter_file.has_key(field):
- def external_wrapper(f):
- def _convert_function(data):
- return data.convert(f)
- return _convert_function
- cf = external_wrapper(field)
- add_field(field, lambda a, b: None, convert_function=cf,take_log=False)
- self.derived_field_list = []
- for field in self.parameter_file.field_info:
- try:
- fd = self.parameter_file.field_info[field].get_dependencies(pf = self.parameter_file)
- except:
- continue
- available = na.all([f in self.field_list for f in fd.requested])
- if available: self.derived_field_list.append(field)
- for field in self.field_list:
- if field not in self.derived_field_list:
- self.derived_field_list.append(field)
-
class EnzoHierarchy(AMRHierarchy):
- eiTopGrid = None
- _strip_path = False
- @time_execution
- def __init__(self, pf, data_style=None):
- """
- This is the grid structure as Enzo sees it, with some added bonuses.
- It's primarily used as a class factory, to generate data objects and
- access grids.
- It should never be created directly -- you should always access it via
- calls to an affiliated :class:`~yt.lagos.EnzoStaticOutput`.
+ _strip_path = False
+ grid = EnzoGrid
- On instantiation, it processes the hierarchy and generates the grids.
- """
- # Expect filename to be the name of the parameter file, not the
- # hierarchy
+ def __init__(self, pf, data_style):
+
self.data_style = data_style
- self.hierarchy_filename = os.path.abspath(pf.parameter_filename) \
- + ".hierarchy"
- if os.path.getsize(self.hierarchy_filename) == 0:
+ self.hierarchy_filename = os.path.abspath(
+ "%s.hierarchy" % (pf.parameter_filename))
+ harray_fn = self.hierarchy_filename[:-9] + "harrays"
+ if os.path.exists(harray_fn):
+ harray_fp = h5py.File(harray_fn)
+ self.num_grids = harray_fp["/Level"].len()
+ elif 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"
self.directory = os.path.dirname(self.hierarchy_filename)
- # 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()
- testGrid = testGridID = None
- 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"], testGrid, testGridID)
- self._setup_data_queue()
+
# For some reason, r8 seems to want Float64
if pf.has_key("CompilerPrecision") \
and pf["CompilerPrecision"] == "r4":
@@ -689,34 +359,45 @@
else:
self.float_type = 'float64'
- AMRHierarchy.__init__(self, pf)
-
+ AMRHierarchy.__init__(self, pf, data_style)
# sync it back
self.parameter_file.data_style = self.data_style
- del self.__hierarchy_string
-
def _setup_classes(self):
dd = self._get_data_reader_dict()
AMRHierarchy._setup_classes(self, dd)
- self._add_object_class('grid', "EnzoGrid", EnzoGridBase, dd)
self.object_types.sort()
- def __guess_data_style(self, rank, testGrid, testGridID):
- if self.data_style: return
- if testGrid[0] != os.path.sep:
- testGrid = os.path.join(self.directory, testGrid)
- if not os.path.exists(testGrid):
- testGrid = os.path.join(self.directory,
- os.path.basename(testGrid))
+ def _count_grids(self):
+ if self.num_grids is not None: return
+ test_grid = test_grid_id = None
+ self.num_stars = 0
+ for line in rlines(open(self.hierarchy_filename, "rb")):
+ if line.startswith("BaryonFileName") or \
+ line.startswith("FileName "):
+ test_grid = line.split("=")[-1].strip().rstrip()
+ if line.startswith("NumberOfStarParticles"):
+ self.num_stars = int(line.split("=")[-1])
+ if line.startswith("Grid "):
+ self.num_grids = test_grid_id = int(line.split("=")[-1])
+ break
+ self._guess_data_style(self.pf["TopGridRank"], test_grid, test_grid_id)
+
+ def _guess_data_style(self, rank, test_grid, test_grid_id):
+ if self.data_style is not None: return
+ if test_grid[0] != os.path.sep:
+ test_grid = os.path.join(self.directory, test_grid)
+ if not os.path.exists(test_grid):
+ test_grid = os.path.join(self.directory,
+ os.path.basename(test_grid))
mylog.debug("Your data uses the annoying hardcoded path.")
self._strip_path = True
try:
- a = SD.SD(testGrid)
+ a = SD.SD(test_grid)
self.data_style = 'enzo_hdf4'
mylog.debug("Detected HDF4")
except:
- list_of_sets = HDF5LightReader.ReadListOfDatasets(testGrid, "/")
+ list_of_sets = HDF5LightReader.ReadListOfDatasets(test_grid, "/")
if len(list_of_sets) == 0 and rank == 3:
mylog.debug("Detected packed HDF5")
self.data_style = 'enzo_packed_3d'
@@ -732,257 +413,212 @@
else:
raise TypeError
- def __setup_filemap(self, grid):
- if not self.data_style == 'enzo_packed_3d':
- return
- try:
- self.cpu_map[grid.filename].append(grid)
- except AttributeError:
- pass
+ # Sets are sorted, so that won't work!
+ def _parse_hierarchy(self):
+ def _next_token_line(token, f):
+ line = f.readline()
+ while token not in line:
+ line = f.readline()
+ return line.split()[2:]
+ if os.path.exists(self.hierarchy_filename[:-9] + "harrays"):
+ if self._parse_binary_hierarchy(): return
+ t1 = time.time()
+ pattern = r"Pointer: Grid\[(\d*)\]->NextGrid(Next|This)Level = (\d*)$"
+ patt = re.compile(pattern)
+ f = open(self.hierarchy_filename, "rb")
+ self.grids = [self.grid(1, self)]
+ self.grids[0].Level = 0
+ si, ei, LE, RE, fn, np = [], [], [], [], [], []
+ all = [si, ei, LE, RE, fn]
+ f.readline() # Blank at top
+ for grid_id in xrange(self.num_grids):
+ # We will unroll this list
+ si.append(_next_token_line("GridStartIndex", f))
+ ei.append(_next_token_line("GridEndIndex", f))
+ LE.append(_next_token_line("GridLeftEdge", f))
+ RE.append(_next_token_line("GridRightEdge", f))
+ nb = int(_next_token_line("NumberOfBaryonFields", f)[0])
+ fn.append(["-1"])
+ if nb > 0: fn[-1] = _next_token_line("BaryonFileName", f)
+ np.append(int(_next_token_line("NumberOfParticles", f)[0]))
+ if nb == 0 and np[-1] > 0: fn[-1] = _next_token_line("FileName", f)
+ line = f.readline()
+ while len(line) > 2:
+ if line.startswith("Pointer:"):
+ vv = patt.findall(line)[0]
+ self.__pointer_handler(vv)
+ line = f.readline()
+ continue
+ params = line.split()
+ line = f.readline()
+ self.grid_dimensions.flat[:] = ei
+ self.grid_dimensions -= na.array(si, self.float_type)
+ self.grid_dimensions += 1
+ self.grid_left_edge.flat[:] = LE
+ self.grid_right_edge.flat[:] = RE
+ self.grid_particle_count.flat[:] = np
+ self.grids = na.array(self.grids, dtype='object')
+ self.filenames = fn
+ self._store_binary_hierarchy()
+ t2 = time.time()
+
+ def __pointer_handler(self, m):
+ sgi = int(m[2])-1
+ if sgi == -1: return # if it's 0, then we're done with that lineage
+ # Okay, so, we have a pointer. We make a new grid, with an id of the length+1
+ # (recall, Enzo grids are 1-indexed)
+ self.grids.append(self.grid(len(self.grids)+1, self))
+ # We'll just go ahead and make a weakref to cache
+ second_grid = self.grids[sgi] # zero-indexed already
+ first_grid = self.grids[int(m[0])-1]
+ if m[1] == "Next":
+ first_grid._children_ids.append(second_grid.id)
+ second_grid._parent_id = first_grid.id
+ second_grid.Level = first_grid.Level + 1
+ elif m[1] == "This":
+ if first_grid.Parent is not None:
+ first_grid.Parent._children_ids.append(second_grid.id)
+ second_grid._parent_id = first_grid._parent_id
+ second_grid.Level = first_grid.Level
+ self.grid_levels[sgi] = second_grid.Level
+
+ def _parse_binary_hierarchy(self):
+ mylog.info("Getting the binary hierarchy")
+ f = h5py.File(self.hierarchy_filename[:-9] + "harrays")
+ self.grid_dimensions[:] = f["/ActiveDimensions"][:]
+ self.grid_left_edge[:] = f["/LeftEdges"][:]
+ self.grid_right_edge[:] = f["/RightEdges"][:]
+ self.grid_particle_count[:,0] = f["/NumberOfParticles"][:]
+ levels = f["/Level"][:]
+ parents = f["/ParentIDs"][:]
+ procs = f["/Processor"][:]
+ grids = []
+ self.filenames = []
+ grids = [self.grid(gi+1, self) for gi in xrange(self.num_grids)]
+ giter = izip(grids, levels, procs, parents)
+ bn = "%s.cpu%%04i" % (self.pf)
+ pmap = [(bn % P,) for P in xrange(procs.max()+1)]
+ for grid,L,P,Pid in giter:
+ grid.Level = L
+ grid._parent_id = Pid
+ if Pid > -1:
+ grids[Pid-1]._children_ids.append(grid.id)
+ self.filenames.append(pmap[P])
+ self.grids = na.array(grids, dtype='object')
+ f.close()
+ mylog.info("Finished with binary hierarchy reading")
+ return True
- def __del__(self):
- self._close_data_file()
- try:
- del self.eiTopGrid
- except:
- pass
- for gridI in xrange(self.num_grids):
- for g in self.gridTree[gridI]:
- del g
- del self.gridReverseTree
- del self.gridLeftEdge, self.gridRightEdge
- del self.gridLevels, self.gridStartIndices, self.gridEndIndices
- del self.gridTimes
- del self.gridTree
-
- def __set_all_filenames(self, fns):
- if self._strip_path:
- for fnI, fn in enumerate(fns):
- self.grids[fnI].filename = os.path.sep.join([self.directory,
- os.path.basename(fn)])
- elif fns[0][0] == os.path.sep:
- for fnI, fn in enumerate(fns):
- self.grids[fnI].filename = fn
- else:
- for fnI, fn in enumerate(fns):
- self.grids[fnI].filename = os.path.sep.join([self.directory,
- fn])
- mylog.debug("Done with baryon filenames")
- for g in self.grids:
- self.__setup_filemap(g)
- mylog.debug("Done with filemap")
+ def _store_binary_hierarchy(self):
+ # We don't do any of the logic here, we just check if the data file
+ # is open...
+ if self._data_file is None: return
+ if self.data_style != "enzo_packed_3d": return
+ mylog.info("Storing the binary hierarchy")
+ f = h5py.File(self.hierarchy_filename[:-9] + "harrays", "w")
+ f.create_dataset("/LeftEdges", data=self.grid_left_edge)
+ f.create_dataset("/RightEdges", data=self.grid_right_edge)
+ parents, procs, levels = [], [], []
+ for i,g in enumerate(self.grids):
+ if g.Parent is not None:
+ parents.append(g.Parent.id)
+ else:
+ parents.append(-1)
+ procs.append(int(self.filenames[i][0][-4:]))
+ levels.append(g.Level)
+
+ parents = na.array(parents, dtype='int64')
+ procs = na.array(procs, dtype='int64')
+ levels = na.array(levels, dtype='int64')
+ f.create_dataset("/ParentIDs", data=parents)
+ f.create_dataset("/Processor", data=procs)
+ f.create_dataset("/Level", data=levels)
- def __parse_hierarchy_file(self):
- def __split_convert(vals, func, toAdd, curGrid):
- """
- Quick function to split up a parameter and convert it and toss onto a grid
- """
- j = 0
- for v in vals.split():
- toAdd[curGrid-1,j] = func(v)
- j+=1
- 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", line_index)
- if len(line) < 2:
- continue
- param, vals = line.split("=")
- param = param.rstrip() # This slows things down considerably...
- # or so I used to think...
- if param == "Grid":
- curGrid = int(vals)
- self.grids[curGrid-1] = self.grid(curGrid)
- elif param == "GridDimension":
- __split_convert(vals, float, self.gridDimensions, curGrid)
- elif param == "GridStartIndex":
- __split_convert(vals, int, self.gridStartIndices, curGrid)
- elif param == "GridEndIndex":
- __split_convert(vals, int, self.gridEndIndices, curGrid)
- elif param == "GridLeftEdge":
- __split_convert(vals, float, self.gridLeftEdge, curGrid)
- elif param == "GridRightEdge":
- __split_convert(vals, float, self.gridRightEdge, curGrid)
- elif param == "Level":
- __split_convert(vals, int, self.gridLevels, curGrid)
- elif param == "Time":
- __split_convert(vals, float, self.gridTimes, curGrid)
- elif param == "NumberOfParticles":
- __split_convert(vals, int, self.gridNumberOfParticles, curGrid)
- elif param == "FileName":
- self.grids[curGrid-1].set_filename(vals[1:-1])
- elif param == "BaryonFileName":
- self.grids[curGrid-1].set_filename(vals[1:-1])
- mylog.info("Caching hierarchy information")
- allArrays = na.zeros((self.num_grids,18),'float64')
- allArrays[:,0:3] = self.gridDimensions[:]
- allArrays[:,3:6] = self.gridStartIndices[:]
- allArrays[:,6:9] = self.gridEndIndices[:]
- allArrays[:,9:12] = self.gridLeftEdge[:]
- allArrays[:,12:15] = self.gridRightEdge[:]
- allArrays[:,15:16] = self.gridLevels[:]
- allArrays[:,16:17] = self.gridTimes[:]
- allArrays[:,17:18] = self.gridNumberOfParticles[:]
- if self.num_grids > 1000:
- self.save_data(allArrays, "/","Hierarchy")
- del allArrays
-
- def __obtain_filenames(self):
- mylog.debug("Copied to local array.")
- # This needs to go elsewhere:
- # Now get the baryon filenames
- mylog.debug("Getting baryon filenames")
- for patt in ["BaryonFileName", "FileName", "ParticleFileName"]:
- re_FileName = constructRegularExpressions(patt,('s'))
- fn_results = re.findall(re_FileName, self.__hierarchy_string)
- if len(fn_results):
- self.__set_all_filenames(fn_results)
- return
+ f.create_dataset("/ActiveDimensions", data=self.grid_dimensions)
+ f.create_dataset("/NumberOfParticles", data=self.grid_particle_count[:,0])
- def __setup_grid_tree(self):
- mylog.debug("No cached tree found, creating")
- self.grids[0].Level = 0 # Bootstrap
- self.gridLevels[0] = 0 # Bootstrap
- p = re.compile(r"Pointer: Grid\[(\d*)\]->NextGrid(Next|This)Level = (\d*)$", re.M)
- # Now we assemble the grid tree
- # This is where all the time is spent.
- for m in p.finditer(self.__hierarchy_string):
- secondGrid = int(m.group(3))-1 # zero-index versus one-index
- if secondGrid == -1:
- continue
- firstGrid = int(m.group(1))-1
- if m.group(2) == "Next":
- self.gridTree[firstGrid].append(weakref.proxy(self.grids[secondGrid]))
- self.gridReverseTree[secondGrid] = firstGrid + 1
- self.grids[secondGrid].Level = self.grids[firstGrid].Level + 1
- self.gridLevels[secondGrid] = self.gridLevels[firstGrid] + 1
- elif m.group(2) == "This":
- parent = self.gridReverseTree[firstGrid]
- if parent and parent > -1:
- self.gridTree[parent-1].append(weakref.proxy(self.grids[secondGrid]))
- self.gridReverseTree[secondGrid] = parent
- self.grids[secondGrid].Level = self.grids[firstGrid].Level
- 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", force=True)
- self.save_data(na.array(self.gridReverseTree), "/", "ReverseTree", force=True)
- self.save_data(self.gridLevels, "/", "Levels", force=True)
+ f.close()
- @parallel_blocking_call
- def _populate_hierarchy(self):
- """
- Instantiates all of the grid objects, with their appropriate
- parameters. This is the work-horse.
- """
- if self.data_style == 'enzo_packed_3d':
- self.cpu_map = defaultdict(lambda: [][:])
- self.file_access = {}
- harray = self.get_data("/", "Hierarchy")
- if self.num_grids <= 1000:
- mylog.info("Skipping serialization!")
- if harray and self.num_grids > 1000:
- self._deserialize_hierarchy(harray)
- else:
- self.__parse_hierarchy_file()
- self.__obtain_filenames()
- treeArray = self.get_data("/", "Tree")
- if treeArray == None:
- self.__setup_grid_tree()
- else:
- mylog.debug("Grabbing serialized tree data")
- 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
- mylog.debug("Tree created")
- self.maxLevel = self.gridLevels.max()
- self.max_level = self.maxLevel
- # Now we do things that we need all the grids to do
- #self.fieldList = self.grids[0].getFields()
- # The rest of this can probably be done with list comprehensions, but
- # I think this way is clearer.
- mylog.debug("Preparing grids")
- for i, grid in enumerate(self.grids):
- if (i%1e4) == 0: mylog.debug("Prepared % 7i / % 7i grids", i, self.num_grids)
- grid._prepare_grid()
- self._setup_grid_dxs()
- mylog.debug("Prepared")
- self._setup_field_lists()
- self.levelIndices = {}
- self.levelNum = {}
- ad = self.gridEndIndices - self.gridStartIndices + 1
- for level in xrange(self.maxLevel+1):
- self.level_stats[level]['numgrids'] = na.where(self.gridLevels==level)[0].size
- li = na.where(self.gridLevels[:,0] == level)
- self.level_stats[level]['numcells'] = ad[li].prod(axis=1).sum()
- self.levelIndices[level] = self._select_level(level)
- self.levelNum[level] = len(self.levelIndices[level])
- mylog.debug("Hierarchy fully populated.")
+ def _rebuild_top_grids(self, level = 0):
+ #for level in xrange(self.max_level+1):
+ mylog.info("Rebuilding grids on level %s", level)
+ cmask = (self.grid_levels.flat == (level + 1))
+ cmsum = cmask.sum()
+ mask = na.zeros(self.num_grids, dtype='bool')
+ for grid in self.select_grids(level):
+ mask[:] = 0
+ LE = self.grid_left_edge[grid.id - grid._id_offset]
+ RE = self.grid_right_edge[grid.id - grid._id_offset]
+ grids, grid_i = self.get_box_grids(LE, RE)
+ mask[grid_i] = 1
+ grid._children_ids = []
+ cgrids = self.grids[ ( mask * cmask).astype('bool') ]
+ mylog.info("%s: %s / %s", grid, len(cgrids), cmsum)
+ for cgrid in cgrids:
+ grid._children_ids.append(cgrid.id)
+ cgrid._parent_id = grid.id
+ mylog.info("Finished rebuilding")
- def _setup_grid_dxs(self):
- mylog.debug("Setting up corners and dxs")
- self._setup_grid_corners()
- dx = (self.gridRightEdge[:,0] - self.gridLeftEdge[:,0]) / \
- (self.gridEndIndices[:,0]-self.gridStartIndices[:,0]+1)
- dy = (self.gridRightEdge[:,1] - self.gridLeftEdge[:,1]) / \
- (self.gridEndIndices[:,1]-self.gridStartIndices[:,1]+1)
- dz = (self.gridRightEdge[:,2] - self.gridLeftEdge[:,2]) / \
- (self.gridEndIndices[:,2]-self.gridStartIndices[:,2]+1)
- self.gridDxs[:,0] = dx[:]
- self.gridDys[:,0] = dy[:]
- self.gridDzs[:,0] = dz[:]
- mylog.debug("Flushing to grids")
- for grid in self.grids:
- grid._setup_dx()
- mylog.debug("Done flushing to grids")
- if ytcfg.getboolean("lagos","ReconstructHierarchy") == True:
- mylog.debug("Reconstructing hierarchy")
- for level in range(self.maxLevel+1):
- grids_to_recon = self.select_grids(level)
- pbar = None
- if len(self.grids) > 3e5:
- pbar = get_pbar('Reconstructing level % 2i / % 2i ' \
- % (level, self.maxLevel),
- len(grids_to_recon))
- for i,grid in enumerate(grids_to_recon):
- if pbar: pbar.update(i)
- if grid.Parent is not None: grid._guess_properties_from_parent()
- if pbar: pbar.finish()
-
- def _join_field_lists(self, field_list):
- return field_list
-
- def _setup_field_lists(self):
- field_list = self.get_data("/", "DataFields")
- if field_list is None:
- mylog.info("Gathering a field list (this may take a moment.)")
- field_list = set()
- random_sample = self._generate_random_grids()
- for grid in random_sample:
- if not hasattr(grid, 'filename'): continue
- try:
- gf = grid.getFields()
- except grid._read_exception:
- mylog.debug("Grid %s is a bit funky?", grid.id)
- continue
- mylog.debug("Grid %s has: %s", grid.id, gf)
- field_list = field_list.union(gf)
- field_list = self._join_field_lists(field_list)
- self.save_data(list(field_list),"/","DataFields")
+ def _populate_grid_objects(self):
+ for g,f in izip(self.grids, self.filenames):
+ g._prepare_grid()
+ g._setup_dx()
+ g.set_filename(f[0])
+ #if g.Parent is not None: g._guess_properties_from_parent()
+ del self.filenames # No longer needed.
+ self.max_level = self.grid_levels.max()
+
+ def _detect_fields(self):
+ # Do this only on the root processor to save disk work.
+ if self._mpi_get_rank() == 0 or self._mpi_get_rank() == None:
+ field_list = self.get_data("/", "DataFields")
+ if field_list is None:
+ mylog.info("Gathering a field list (this may take a moment.)")
+ field_list = set()
+ random_sample = self._generate_random_grids()
+ for grid in random_sample:
+ if not hasattr(grid, 'filename'): continue
+ try:
+ gf = self.io._read_field_names(grid)
+ except self.io._read_exception:
+ mylog.debug("Grid %s is a bit funky?", grid.id)
+ continue
+ mylog.debug("Grid %s has: %s", grid.id, gf)
+ field_list = field_list.union(gf)
+ else:
+ field_list = None
+ field_list = self._mpi_bcast_pickled(field_list)
+ self.save_data(list(field_list),"/","DataFields",passthrough=True)
self.field_list = list(field_list)
+ def _setup_unknown_fields(self):
+ for field in self.field_list:
+ if field in self.parameter_file.field_info: continue
+ mylog.info("Adding %s to list of fields", field)
+ cf = None
+ if self.parameter_file.has_key(field):
+ def external_wrapper(f):
+ def _convert_function(data):
+ return data.convert(f)
+ return _convert_function
+ cf = external_wrapper(field)
+ add_field(field, lambda a, b: None,
+ convert_function=cf, take_log=False)
+
+ def _setup_derived_fields(self):
+ self.derived_field_list = []
+ for field in self.parameter_file.field_info:
+ try:
+ fd = self.parameter_file.field_info[field].get_dependencies(
+ pf = self.parameter_file)
+ except:
+ continue
+ available = na.all([f in self.field_list for f in fd.requested])
+ if available: self.derived_field_list.append(field)
+ for field in self.field_list:
+ if field not in self.derived_field_list:
+ self.derived_field_list.append(field)
+
def _generate_random_grids(self):
if self.num_grids > 40:
starter = na.random.randint(0, 20)
@@ -992,96 +628,125 @@
random_sample = na.mgrid[0:max(len(self.grids)-1,1)].astype("int32")
return self.grids[(random_sample,)]
+ def find_particles_by_type(self, ptype, max_num=None, additional_fields=None):
+ """
+ Returns a structure of arrays with all of the particles'
+ positions, velocities, masses, types, IDs, and attributes for
+ a particle type **ptype** for a maximum of **max_num**
+ particles. If non-default particle fields are used, provide
+ them in **additional_fields**.
+ """
+ # Not sure whether this routine should be in the general HierarchyType.
+ if self.gridNumberOfParticles.sum() == 0:
+ mylog.info("Data contains no particles.");
+ return None
+ if additional_fields is None:
+ additional_fields = ['metallicity_fraction', 'creation_time',
+ 'dynamical_time']
+ pfields = [f for f in self.field_list if f.startswith('particle_')]
+ nattr = self.parameter_file['NumberOfParticleAttributes']
+ if nattr > 0:
+ pfields += additional_fields[:nattr]
+ # Find where the particles reside and count them
+ if max_num is None: max_num = 1e100
+ total = 0
+ pstore = []
+ for level in range(self.max_level, -1, -1):
+ for grid in self.select_grids(level):
+ index = na.where(grid['particle_type'] == ptype)[0]
+ total += len(index)
+ pstore.append(index)
+ if total >= max_num: break
+ if total >= max_num: break
+ result = None
+ if total > 0:
+ result = {}
+ for p in pfields:
+ result[p] = na.zeros(total, 'float64')
+ # Now we retrieve data for each field
+ ig = count = 0
+ for level in range(self.max_level, -1, -1):
+ for grid in self.select_grids(level):
+ nidx = len(pstore[ig])
+ if nidx > 0:
+ for p in pfields:
+ result[p][count:count+nidx] = grid[p][pstore[ig]]
+ count += nidx
+ ig += 1
+ if count >= total: break
+ if count >= total: break
+ # Crop data if retrieved more than max_num
+ if count > max_num:
+ for p in pfields:
+ result[p] = result[p][0:max_num]
+ return result
+
+
class EnzoHierarchyInMemory(EnzoHierarchy):
- _data_style = 'enzo_inline'
- def _obtain_enzo(self):
- import enzo; return enzo
+
+ grid = EnzoGridInMemory
+ _enzo = None
+
+ @property
+ def enzo(self):
+ if self._enzo is None:
+ import enzo
+ self._enzo = enzo
+ return self._enzo
def __init__(self, pf, data_style = None):
+ self.data_style = data_style
+ self.float_type = 'float64'
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 = enzo.hierarchy_information["GridLeftEdge"].dtype
- self.data_style = data_style # Mandated
+ self.float_type = self.enzo.hierarchy_information["GridLeftEdge"].dtype
self.directory = os.getcwd()
- self.num_grids = enzo.hierarchy_information["GridDimensions"].shape[0]
- self._setup_data_queue()
- AMRHierarchy.__init__(self, pf)
+ AMRHierarchy.__init__(self, pf, data_style)
- def _initialize_data_file(self):
+ def _initialize_data_storage(self):
pass
- def _join_field_lists(self, field_list):
- from mpi4py import MPI
- MPI.COMM_WORLD.Barrier()
- 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 = list(set(data))
- else:
- MPI.COMM_WORLD.send(data, dest=0, tag=0)
- MPI.COMM_WORLD.Barrier()
- return MPI.COMM_WORLD.bcast(data, root=0)
+ def _count_grids(self):
+ self.num_grids = self.enzo.hierarchy_information["GridDimensions"].shape[0]
- def _populate_hierarchy(self):
+ def _parse_hierarchy(self):
self._copy_hierarchy_structure()
- enzo = self._obtain_enzo()
mylog.debug("Copying reverse tree")
- self.gridReverseTree = enzo.hierarchy_information["GridParentIDs"].ravel().tolist()
+ reverse_tree = self.enzo.hierarchy_information["GridParentIDs"].ravel().tolist()
# Initial setup:
mylog.debug("Reconstructing parent-child relationships")
- #self.gridTree = [ [] for i in range(self.num_grids) ]
- for id,pid in enumerate(self.gridReverseTree):
+ self.grids = []
+ # We enumerate, so it's 0-indexed id and 1-indexed pid
+ self.filenames = ["-1"] * self.num_grids
+ for id,pid in enumerate(reverse_tree):
+ self.grids.append(self.grid(id+1, self))
+ self.grids[-1].Level = self.grid_levels[id]
if pid > 0:
- self.gridTree[pid-1].append(
- weakref.proxy(self.grids[id]))
- else:
- self.gridReverseTree[id] = None
- self.max_level = self.gridLevels.max()
- self.maxLevel = self.max_level
+ self.grids[-1]._parent_id = pid
+ self.grids[pid-1]._children_ids.append(self.grids[-1].id)
+ self.max_level = self.grid_levels.max()
mylog.debug("Preparing grids")
for i, grid in enumerate(self.grids):
if (i%1e4) == 0: mylog.debug("Prepared % 7i / % 7i grids", i, self.num_grids)
grid.filename = None
grid._prepare_grid()
- grid.proc_num = self.gridProcs[i,0]
- self._setup_grid_dxs()
+ grid.proc_num = self.grid_procs[i,0]
+ self.grids = na.array(self.grids, dtype='object')
mylog.debug("Prepared")
- self._setup_field_lists()
- self.levelIndices = {}
- self.levelNum = {}
- ad = self.gridEndIndices - self.gridStartIndices + 1
- for level in xrange(self.maxLevel+1):
- self.level_stats[level]['numgrids'] = na.where(self.gridLevels==level)[0].size
- li = na.where(self.gridLevels[:,0] == level)
- self.level_stats[level]['numcells'] = ad[li].prod(axis=1).sum()
- self.levelIndices[level] = self._select_level(level)
- self.levelNum[level] = len(self.levelIndices[level])
- mylog.debug("Hierarchy fully populated.")
- def _copy_hierarchy_structure(self):
- enzo = self._obtain_enzo()
- self.gridDimensions[:] = enzo.hierarchy_information["GridDimensions"][:]
- self.gridStartIndices[:] = enzo.hierarchy_information["GridStartIndices"][:]
- self.gridEndIndices[:] = enzo.hierarchy_information["GridEndIndices"][:]
- self.gridLeftEdge[:] = enzo.hierarchy_information["GridLeftEdge"][:]
- self.gridRightEdge[:] = enzo.hierarchy_information["GridRightEdge"][:]
- self.gridLevels[:] = enzo.hierarchy_information["GridLevels"][:]
- self.gridTimes[:] = enzo.hierarchy_information["GridTimes"][:]
- self.gridProcs = enzo.hierarchy_information["GridProcs"].copy()
- self.gridNumberOfParticles[:] = enzo.hierarchy_information["GridNumberOfParticles"][:]
+ def _initialize_grid_arrays(self):
+ EnzoHierarchy._initialize_grid_arrays(self)
+ self.grid_procs = na.zeros((self.num_grids,1),'int32')
- def _generate_random_grids(self):
- my_proc = ytcfg.getint("yt","__parallel_rank")
- gg = self.grids[self.gridProcs[:,0] == my_proc]
- if len(gg) > 40:
- starter = na.random.randint(0, 20)
- random_sample = na.mgrid[starter:len(gg)-1:20j].astype("int32")
- mylog.debug("Checking grids: %s", random_sample.tolist())
- else:
- random_sample = na.mgrid[0:max(len(gg)-1,1)].astype("int32")
- return gg[(random_sample,)]
+ def _copy_hierarchy_structure(self):
+ # Dimensions are important!
+ self.grid_dimensions[:] = self.enzo.hierarchy_information["GridEndIndices"][:]
+ self.grid_dimensions -= self.enzo.hierarchy_information["GridStartIndices"][:]
+ self.grid_dimensions += 1
+ self.grid_left_edge[:] = self.enzo.hierarchy_information["GridLeftEdge"][:]
+ self.grid_right_edge[:] = self.enzo.hierarchy_information["GridRightEdge"][:]
+ self.grid_levels[:] = self.enzo.hierarchy_information["GridLevels"][:]
+ self.grid_procs = self.enzo.hierarchy_information["GridProcs"].copy()
+ self.grid_particle_count[:] = self.enzo.hierarchy_information["GridNumberOfParticles"][:]
def save_data(self, *args, **kwargs):
pass
@@ -1185,7 +850,7 @@
self.readGlobalHeader(header_filename,self.parameter_file.paranoid_read) # also sets up the grid objects
self.__cache_endianness(self.levels[-1].grids[-1])
AMRHierarchy.__init__(self,pf)
- self._setup_data_queue()
+ self._setup_data_io()
self._setup_field_list()
def readGlobalHeader(self,filename,paranoid_read):
Added: trunk/yt/lagos/ObjectFindingMixin.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/ObjectFindingMixin.py Fri Nov 20 20:27:29 2009
@@ -0,0 +1,168 @@
+"""
+AMR hierarchy container class
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2007-2009 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from yt.arraytypes import *
+from yt.logger import lagosLogger as mylog
+from EnzoDefs import NUMTOCHECK
+
+class ObjectFindingMixin(object):
+
+ def find_ray_grids(self, coord, axis):
+ """
+ Returns the (objects, indices) of grids that an (x,y) ray intersects
+ along *axis*
+ """
+ # Let's figure out which grids are on the slice
+ mask=na.ones(self.num_grids)
+ # So if gRE > coord, we get a mask, if not, we get a zero
+ # if gLE > coord, we get a zero, if not, mask
+ # Thus, if the coordinate is between the two edges, we win!
+ na.choose(na.greater(self.grid_right_edge[:,x_dict[axis]],coord[0]),(0,mask),mask)
+ na.choose(na.greater(self.grid_left_edge[:,x_dict[axis]],coord[0]),(mask,0),mask)
+ na.choose(na.greater(self.grid_right_edge[:,y_dict[axis]],coord[1]),(0,mask),mask)
+ na.choose(na.greater(self.grid_left_edge[:,y_dict[axis]],coord[1]),(mask,0),mask)
+ ind = na.where(mask == 1)
+ return self.grids[ind], ind
+
+ def find_max(self, field, finest_levels = True):
+ """
+ Returns (value, center) of location of maximum for a given field.
+ """
+ if (field, finest_levels) in self._max_locations:
+ return self._max_locations[(field, finest_levels)]
+ mg, mc, mv, pos = self.find_max_cell_location(field, finest_levels)
+ self._max_locations[(field, finest_levels)] = (mv, pos)
+ return mv, pos
+
+ def find_max_cell_location(self, field, finest_levels = True):
+ if finest_levels is True:
+ gi = (self.grid_levels >= self.max_level - NUMTOCHECK).ravel()
+ source = self.grid_collection([0.0]*3, self.grids[gi])
+ else:
+ source = self.all_data()
+ mylog.debug("Searching %s grids for maximum value of %s",
+ len(source._grids), field)
+ max_val, maxi, mx, my, mz, mg = \
+ source.quantities["MaxLocation"]( field, lazy_reader=True)
+ max_grid = self.grids[mg]
+ mc = na.unravel_index(maxi, max_grid.ActiveDimensions)
+ mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s %s", \
+ max_val, mx, my, mz, max_grid, max_grid.Level, mc)
+ self.parameters["Max%sValue" % (field)] = max_val
+ self.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
+ return max_grid, mc, max_val, na.array((mx,my,mz), dtype='float64')
+
+ def find_min(self, field):
+ """
+ Returns (value, center) of location of minimum for a given field
+ """
+ gI = na.where(self.grid_levels >= 0) # Slow but pedantic
+ minVal = 1e100
+ for grid in self.grids[gI[0]]:
+ mylog.debug("Checking %s (level %s)", grid.id, grid.Level)
+ val, coord = grid.find_min(field)
+ if val < minVal:
+ minCoord = coord
+ minVal = val
+ minGrid = grid
+ mc = na.array(minCoord)
+ pos=minGrid.get_position(mc)
+ mylog.info("Min Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s", \
+ minVal, pos[0], pos[1], pos[2], minGrid, minGrid.Level)
+ self.center = pos
+ self.parameters["Min%sValue" % (field)] = minVal
+ self.parameters["Min%sPos" % (field)] = "%s" % (pos)
+ return minVal, pos
+
+ def find_point(self, coord):
+ """
+ Returns the (objects, indices) of grids containing an (x,y,z) point
+ """
+ mask=na.ones(self.num_grids)
+ for i in xrange(len(coord)):
+ na.choose(na.greater(self.grid_left_edge[:,i],coord[i]), (mask,0), mask)
+ na.choose(na.greater(self.grid_right_edge[:,i],coord[i]), (0,mask), mask)
+ ind = na.where(mask == 1)
+ return self.grids[ind], ind
+
+ def find_slice_grids(self, coord, axis):
+ """
+ Returns the (objects, indices) of grids that a slice intersects along
+ *axis*
+ """
+ # Let's figure out which grids are on the slice
+ mask=na.ones(self.num_grids)
+ # So if gRE > coord, we get a mask, if not, we get a zero
+ # if gLE > coord, we get a zero, if not, mask
+ # Thus, if the coordinate is between the edges, we win!
+ #ind = na.where( na.logical_and(self.grid_right_edge[:,axis] > coord, \
+ #self.grid_left_edge[:,axis] < coord))
+ na.choose(na.greater(self.grid_right_edge[:,axis],coord),(0,mask),mask)
+ na.choose(na.greater(self.grid_left_edge[:,axis],coord),(mask,0),mask)
+ ind = na.where(mask == 1)
+ return self.grids[ind], ind
+
+ def find_sphere_grids(self, center, radius):
+ """
+ Returns objects, indices of grids within a sphere
+ """
+ centers = (self.grid_right_edge + self.grid_left_edge)/2.0
+ long_axis = na.maximum.reduce(self.grid_right_edge - self.grid_left_edge, 1)
+ t = na.abs(centers - center)
+ DW = self.parameter_file["DomainRightEdge"] \
+ - self.parameter_file["DomainLeftEdge"]
+ na.minimum(t, na.abs(DW-t), t)
+ dist = na.sqrt(na.sum((t**2.0), axis=1))
+ gridI = na.where(dist < (radius + long_axis))
+ return self.grids[gridI], gridI
+
+ def get_box_grids(self, left_edge, right_edge):
+ """
+ Gets back all the grids between a left edge and right edge
+ """
+ grid_i = na.where((na.all(self.grid_right_edge > left_edge, axis=1)
+ & na.all(self.grid_left_edge < right_edge, axis=1)) == True)
+ return self.grids[grid_i], grid_i
+
+ def get_periodic_box_grids(self, left_edge, right_edge):
+ left_edge = na.array(left_edge)
+ right_edge = na.array(right_edge)
+ mask = na.zeros(self.grids.shape, dtype='bool')
+ dl = self.parameters["DomainLeftEdge"]
+ dr = self.parameters["DomainRightEdge"]
+ db = right_edge - left_edge
+ for off_x in [-1, 0, 1]:
+ nle = left_edge.copy()
+ nre = left_edge.copy()
+ nle[0] = dl[0] + (dr[0]-dl[0])*off_x + left_edge[0]
+ for off_y in [-1, 0, 1]:
+ nle[1] = dl[1] + (dr[1]-dl[1])*off_y + left_edge[1]
+ for off_z in [-1, 0, 1]:
+ nle[2] = dl[2] + (dr[2]-dl[2])*off_z + left_edge[2]
+ nre = nle + db
+ g, gi = self.get_box_grids(nle, nre)
+ mask[gi] = True
+ return self.grids[mask], na.where(mask)
+
Modified: trunk/yt/lagos/OutputTypes.py
==============================================================================
--- trunk/yt/lagos/OutputTypes.py (original)
+++ trunk/yt/lagos/OutputTypes.py Fri Nov 20 20:27:29 2009
@@ -44,6 +44,10 @@
mylog.debug("Registering: %s as %s", name, cls)
def __new__(cls, filename=None, *args, **kwargs):
+ if not isinstance(filename, types.StringTypes):
+ obj = object.__new__(cls)
+ obj.__init__(filename, *args, **kwargs)
+ return obj
apath = os.path.abspath(filename)
if not os.path.exists(apath): raise IOError(filename)
if apath not in _cached_pfs:
Modified: trunk/yt/lagos/ParallelTools.py
==============================================================================
--- trunk/yt/lagos/ParallelTools.py (original)
+++ trunk/yt/lagos/ParallelTools.py Fri Nov 20 20:27:29 2009
@@ -26,10 +26,11 @@
from yt.lagos import *
from yt.funcs import *
import yt.logger, logging
-import itertools, sys, cStringIO
+import itertools, sys, cStringIO, cPickle
+exe_name = os.path.basename(sys.executable)
# At import time, we determined whether or not we're being run in parallel.
-if os.path.basename(sys.executable) in \
+if exe_name in \
["mpi4py", "embed_enzo",
"python"+sys.version[:3]+"-mpi"] \
or "--parallel" in sys.argv or '_parallel' in dir(sys) \
@@ -42,6 +43,7 @@
ytcfg["yt","__parallel_rank"] = str(MPI.COMM_WORLD.rank)
ytcfg["yt","__parallel_size"] = str(MPI.COMM_WORLD.size)
ytcfg["yt","__parallel"] = "True"
+ if exe_name == "embed_enzo": ytcfg["yt","inline"] = "True"
# I believe we do not need to turn this off manually
#ytcfg["yt","StoreParameterFiles"] = "False"
# Now let's make sure we have the right options set.
@@ -263,7 +265,7 @@
x = na.mgrid[0:1:(cc[0]+1)*1j][cx:cx+2]
y = na.mgrid[0:1:(cc[1]+1)*1j][cy:cy+2]
- DLE, DRE = self.pf["DomainLeftEdge"], self.pf["DomainRightEdge"]
+ DLE, DRE = self.pf["DomainLeftEdge"].copy(), self.pf["DomainRightEdge"].copy()
LE = na.ones(3, dtype='float64') * DLE
RE = na.ones(3, dtype='float64') * DRE
LE[xax] = x[0] * (DRE[xax]-DLE[xax]) + DLE[xax]
@@ -276,9 +278,22 @@
return True, reg
def _partition_hierarchy_3d(self, padding=0.0):
- LE, RE = self.pf["DomainLeftEdge"], self.pf["DomainRightEdge"]
+ LE, RE = self.pf["DomainLeftEdge"].copy(), self.pf["DomainRightEdge"].copy()
if not self._distributed:
return False, LE, RE, self.hierarchy.grid_collection(self.center, self.hierarchy.grids)
+ elif ytcfg.getboolean("yt", "inline"):
+ # At this point, we want to identify the root grid tile to which
+ # this processor is assigned.
+ # The only way I really know how to do this is to get the level-0
+ # grid that belongs to this processor.
+ grids = self.pf.h.select_grids(0)
+ root_grids = [g for g in grids
+ if g.proc_num == MPI.COMM_WORLD.rank]
+ if len(root_grids) != 1: raise RuntimeError
+ #raise KeyError
+ LE = root_grids[0].LeftEdge
+ RE = root_grids[0].RightEdge
+ return True, LE, RE, self.hierarchy.region(self.center, LE, RE)
cc = MPI.Compute_dims(MPI.COMM_WORLD.size, 3)
mi = MPI.COMM_WORLD.rank
@@ -295,6 +310,468 @@
LE, RE, self.hierarchy.periodic_region_strict(self.center, LE-padding, RE+padding)
return False, LE, RE, self.hierarchy.region_strict(self.center, LE, RE)
+
+ def _partition_hierarchy_3d_bisection_list(self):
+ """
+ Returns an array that is used to drive _partition_hierarchy_3d_bisection,
+ below.
+ """
+
+ def factor(n):
+ if n == 1: return [1]
+ i = 2
+ limit = n**0.5
+ while i <= limit:
+ if n % i == 0:
+ ret = factor(n/i)
+ ret.append(i)
+ return ret
+ i += 1
+ return [n]
+
+ cc = MPI.Compute_dims(MPI.COMM_WORLD.size, 3)
+ si = MPI.COMM_WORLD.size
+
+ factors = factor(si)
+ xyzfactors = [factor(cc[0]), factor(cc[1]), factor(cc[2])]
+
+ # Each entry of cuts is a two element list, that is:
+ # [cut dim, number of cuts]
+ cuts = []
+ # The higher cuts are in the beginning.
+ # We're going to do our best to make the cuts cyclic, i.e. x, then y,
+ # then z, etc...
+ lastdim = 0
+ for f in factors:
+ nextdim = (lastdim + 1) % 3
+ while True:
+ if f in xyzfactors[nextdim]:
+ cuts.append([nextdim, f])
+ topop = xyzfactors[nextdim].index(f)
+ temp = xyzfactors[nextdim].pop(topop)
+ lastdim = nextdim
+ break
+ nextdim = (nextdim + 1) % 3
+ return cuts
+
+
+ def _partition_hierarchy_3d_bisection(self, axis, bins, counts, top_bounds = None,\
+ old_group = None, old_comm = None, cut=None, old_cc=None):
+ """
+ Partition the volume into evenly weighted subvolumes using the distribution
+ in counts. The bisection happens in the MPI communicator group old_group.
+ You may need to set "MPI_COMM_MAX" and "MPI_GROUP_MAX" environment
+ variables.
+ """
+ counts = counts.astype('int64')
+ if not self._distributed:
+ LE, RE = self.pf["DomainLeftEdge"].copy(), self.pf["DomainRightEdge"].copy()
+ return False, LE, RE, self.hierarchy.grid_collection(self.center, self.hierarchy.grids)
+
+ # First time through the world is the current group.
+ if old_group == None or old_comm == None:
+ old_group = MPI.COMM_WORLD.Get_group()
+ old_comm = MPI.COMM_WORLD
+
+ # Figure out the gridding based on the deepness of cuts.
+ if old_cc is None:
+ cc = MPI.Compute_dims(MPI.COMM_WORLD.size, 3)
+ else:
+ cc = old_cc
+ cc[cut[0]] /= cut[1]
+ # Set the boundaries of the full bounding box for this group.
+ if top_bounds == None:
+ LE, RE = self.pf["DomainLeftEdge"].copy(), self.pf["DomainRightEdge"].copy()
+ else:
+ LE, RE = top_bounds
+
+ ra = old_group.Get_rank() # In this group, not WORLD, unless it's the first time.
+
+ # First find the total number of particles in my group.
+ parts = old_comm.allreduce(int(counts.sum()), op=MPI.SUM)
+ # Now the full sum in the bins along this axis in this group.
+ full_counts = na.empty(counts.size, dtype='int64')
+ old_comm.Allreduce([counts, MPI.LONG], [full_counts, MPI.LONG], op=MPI.SUM)
+ # Find the bin that passes the cut points.
+ midpoints = [LE[axis]]
+ sum = 0
+ bin = 0
+ for step in xrange(1,cut[1]):
+ while sum < ((parts*step)/cut[1]):
+ lastsum = sum
+ sum += full_counts[bin]
+ bin += 1
+ # Bin edges
+ left_edge = bins[bin-1]
+ right_edge = bins[bin]
+ # Find a better approx of the midpoint cut line using a linear approx.
+ a = float(sum - lastsum) / (right_edge - left_edge)
+ midpoints.append(left_edge + (0.5 - (float(lastsum) / parts / 2)) / a)
+ #midpoint = (left_edge + right_edge) / 2.
+ midpoints.append(RE[axis])
+ # Now we need to split the members of this group into chunks.
+ # The values that go into the _ranks are the ranks of the tasks
+ # in *this* communicator group, which go zero to size - 1. They are not
+ # the same as the global ranks!
+ groups = {}
+ ranks = {}
+ old_group_size = old_group.Get_size()
+ for step in xrange(cut[1]):
+ groups[step] = na.arange(step*old_group_size/cut[1], (step+1)*old_group_size/cut[1])
+ # [ (start, stop, step), ]
+ ranks[step] = [ (groups[step][0], groups[step][-1], 1), ]
+
+ # Based on where we are, adjust our LE or RE, depending on axis. At the
+ # same time assign the new MPI group membership.
+ for step in xrange(cut[1]):
+ if ra in groups[step]:
+ LE[axis] = midpoints[step]
+ RE[axis] = midpoints[step+1]
+ new_group = old_group.Range_incl(ranks[step])
+ new_comm = old_comm.Create(new_group)
+
+ if old_cc is not None:
+ old_group.Free()
+ old_comm.Free()
+
+ new_top_bounds = (LE,RE)
+
+ # Using the new boundaries, regrid.
+ mi = new_comm.rank
+ cx, cy, cz = na.unravel_index(mi, cc)
+ x = na.mgrid[LE[0]:RE[0]:(cc[0]+1)*1j][cx:cx+2]
+ y = na.mgrid[LE[1]:RE[1]:(cc[1]+1)*1j][cy:cy+2]
+ z = na.mgrid[LE[2]:RE[2]:(cc[2]+1)*1j][cz:cz+2]
+
+ my_LE = na.array([x[0], y[0], z[0]], dtype='float64')
+ my_RE = na.array([x[1], y[1], z[1]], dtype='float64')
+
+ # Return a new subvolume and associated stuff.
+ return new_group, new_comm, my_LE, my_RE, new_top_bounds, cc,\
+ self.hierarchy.region_strict(self.center, my_LE, my_RE)
+
+ def _partition_hierarchy_3d_weighted_1d(self, weight=None, bins=None, padding=0.0, axis=0, min_sep=.1):
+ LE, RE = self.pf["DomainLeftEdge"].copy(), self.pf["DomainRightEdge"].copy()
+ if not self._distributed:
+ return False, LE, RE, self.hierarchy.grid_collection(self.center, self.hierarchy.grids)
+
+ cc = MPI.Compute_dims(MPI.COMM_WORLD.size, 3)
+ mi = MPI.COMM_WORLD.rank
+ si = MPI.COMM_WORLD.size
+ cx, cy, cz = na.unravel_index(mi, cc)
+
+ gridx = na.mgrid[LE[0]:RE[0]:(cc[0]+1)*1j]
+ gridy = na.mgrid[LE[1]:RE[1]:(cc[1]+1)*1j]
+ gridz = na.mgrid[LE[2]:RE[2]:(cc[2]+1)*1j]
+
+ x = gridx[cx:cx+2]
+ y = gridy[cy:cy+2]
+ z = gridz[cz:cz+2]
+
+ LE = na.array([x[0], y[0], z[0]], dtype='float64')
+ RE = na.array([x[1], y[1], z[1]], dtype='float64')
+
+ # Default to normal if we don't have a weight, or our subdivisions are
+ # not enough to warrant this procedure.
+ if weight is None or cc[axis] < 1:
+ if padding > 0:
+ return True, \
+ LE, RE, self.hierarchy.periodic_region_strict(self.center, LE-padding, RE+padding)
+
+ return False, LE, RE, self.hierarchy.region_strict(self.center, LE, RE)
+
+ # Find the densest subvolumes globally
+ local_weight = na.zeros((si, weight.size),dtype='float64')
+ local_weight[mi,:] = weight
+ weights = self._mpi_allsum(local_weight)
+ avg_weight = weights.mean()
+ weights = weights.max(axis=0)
+
+ moved_count = 0
+ moved = {}
+ w_copy = weights.copy()
+
+ if mi == 0:
+ print 'w_copy',w_copy,'gridx',gridx
+
+ while moved_count < (cc[axis]-1):
+ con = False
+ # Find the current peak
+ hi_mark = na.argmax(w_copy)
+ # If this peak isn't high enough, we're done
+ height = w_copy[hi_mark]
+ if height < 10.*avg_weight:
+ if mi == 0:
+ print 'breaking',moved_count, height, avg_weight
+ break
+ # If this mark is too close to a previous one, avg this one out
+ # and restart a search.
+ new_cen = (bins[hi_mark] + bins[hi_mark+1])/2.
+ if mi==0:
+ print 'moved',moved
+ for source in moved:
+ if mi == 0:
+ print 'moved',abs(moved[source] - new_cen)
+ if abs(moved[source] - new_cen) < min_sep:
+ w_copy[hi_mark] = avg_weight
+ if mi == 0:
+ print 'continued'
+ con = True
+ if con:
+ continue
+ # Find the lowest value entry
+ lo_mark = na.argmin(w_copy)
+ # Record this as a new mapping.
+ moved[(bins[lo_mark] + bins[lo_mark+1])/2.] = (bins[hi_mark] + bins[hi_mark+1])/2.
+ # Fix the values so they're not pulled again.
+ w_copy[hi_mark] = avg_weight
+ w_copy[lo_mark] = avg_weight
+ moved_count += 1
+
+ # Now for each key in moved, we move the axis closest to that value to
+ # the value in the dict.
+ temp_gridx = []
+ for source in moved:
+ tomove = na.argmin(abs(gridx - source))
+ temp_gridx.append(moved[source])
+ gridx[tomove] = -1.
+
+ for g in gridx:
+ if g >= 0.:
+ temp_gridx.append(g)
+
+ temp_gridx.sort()
+ gridx = na.array(temp_gridx)
+ if mi == 0:
+ print 'gridx',gridx,'len=',len(gridx)
+ x = gridx[cx:cx+2]
+ y = gridy[cy:cy+2]
+ z = gridz[cz:cz+2]
+
+ LE = na.array([x[0], y[0], z[0]], dtype='float64')
+ RE = na.array([x[1], y[1], z[1]], dtype='float64')
+
+ if padding > 0:
+ return True, \
+ LE, RE, self.hierarchy.periodic_region_strict(self.center, LE-padding, RE+padding)
+
+ return False, LE, RE, self.hierarchy.region_strict(self.center, LE, RE)
+
+
+
+ def _partition_hierarchy_3d_weighted(self, weight=None, padding=0.0, agg=8.):
+ LE, RE = self.pf["DomainLeftEdge"].copy(), self.pf["DomainRightEdge"].copy()
+ if not self._distributed:
+ return False, LE, RE, self.hierarchy.grid_collection(self.center, self.hierarchy.grids)
+
+ cc = MPI.Compute_dims(MPI.COMM_WORLD.size, 3)
+ mi = MPI.COMM_WORLD.rank
+ cx, cy, cz = na.unravel_index(mi, cc)
+
+ gridx = na.mgrid[LE[0]:RE[0]:(cc[0]+1)*1j]
+ gridy = na.mgrid[LE[1]:RE[1]:(cc[1]+1)*1j]
+ gridz = na.mgrid[LE[2]:RE[2]:(cc[2]+1)*1j]
+
+ x = gridx[cx:cx+2]
+ y = gridy[cy:cy+2]
+ z = gridz[cz:cz+2]
+
+ LE = na.array([x[0], y[0], z[0]], dtype='float64')
+ RE = na.array([x[1], y[1], z[1]], dtype='float64')
+
+ old_vol = ((RE - LE)**2).sum()
+
+ # Default to normal if we don't have a weight, or our subdivisions are
+ # not enough to warrant this procedure.
+ if weight is None or cc[0] < 2 or cc[1] < 2 or cc[2] < 2:
+ if padding > 0:
+ return True, \
+ LE, RE, self.hierarchy.periodic_region_strict(self.center, LE-padding, RE+padding)
+
+ return False, LE, RE, self.hierarchy.region_strict(self.center, LE, RE)
+
+ # Build the matrix of weights.
+ weights = na.zeros(cc, dtype='float64')
+ weights[cx,cy,cz] = weight
+ weights = self._mpi_allsum(weights)
+ weights = weights / weights.sum()
+
+ # Figure out the sums of weights along the axes
+ xface = weights.sum(axis=0)
+ yface = weights.sum(axis=1)
+ zface = weights.sum(axis=2)
+
+ xedge = yface.sum(axis=1)
+ yedge = xface.sum(axis=1)
+ zedge = xface.sum(axis=0)
+
+ # Get a polynomial fit to each axis weight distribution
+ xcen = gridx[:-1]
+ xcen += xcen[1]/2.
+ ycen = gridy[:-1]
+ ycen += ycen[1]/2.
+ zcen = gridz[:-1]
+ zcen += zcen[1]/2.
+
+ xfit = na.polyfit(xcen, xedge, 3)
+ yfit = na.polyfit(ycen, yedge, 3)
+ zfit = na.polyfit(zcen, zedge, 3)
+
+ # Find the normalized weights with trapizoidal integration
+ # We also apply an aggression factor to the values to make the
+ # boundaries shift more.
+ div_count = int(1. / padding)
+ divs = na.arange(div_count+1, dtype='float64') / div_count
+ xvals = na.polyval(xfit, divs)
+ for i,xv in enumerate(xvals):
+ if xv > xedge.mean(): xvals[i] *= agg
+ yvals = na.polyval(yfit, divs)
+ for i,yv in enumerate(yvals):
+ if yv > yedge.mean(): yvals[i] *= agg
+ zvals = na.polyval(zfit, divs)
+ for i,zv in enumerate(zvals):
+ if zv > zedge.mean(): zvals[i] *= agg
+ xnorm = na.trapz(xvals, x=divs)
+ ynorm = na.trapz(yvals, x=divs)
+ znorm = na.trapz(zvals, x=divs)
+
+ # We want to start the integration from the side of the axis where
+ # the highest density is, so that it gets small regions.
+ xstart = float(na.argmax(xedge))/2.
+ if xstart > 0.5: xstart = div_count
+ else: xstart = 0
+
+ ystart = float(na.argmax(yedge))/2.
+ if ystart > 0.5: ystart = div_count
+ else: ystart = 0
+
+ zstart = float(na.argmax(zedge))/2.
+ if zstart > 0.5: zstart = div_count
+ else: zstart = 0
+
+
+ # Find the boundaries. We are assured that none of the boundaries are
+ # too small because each step of div is big enough because it's set
+ # by the padding.
+ nextx = 1./xedge.size
+ nexty = 1./yedge.size
+ nextz = 1./zedge.size
+ boundx = [0.]
+ boundy = [0.]
+ boundz = [0.]
+ donex, doney, donez = False, False, False
+ for i in xrange(div_count):
+ if xstart == 0:
+ xi = 0
+ xv = i
+ xa = i
+ else:
+ xi = div_count - i
+ xf = div_count
+ xa = xi
+ if (na.trapz(xvals[xi:xf], x=divs[xi:xf])/xnorm) >= nextx and not donex:
+ boundx.append(divs[xa])
+ if len(boundx) == cc[0]:
+ donex = True
+ nextx += 1./xedge.size
+ if ystart == 0:
+ yi = 0
+ yf = i
+ ya = i
+ else:
+ yi = div_count - i
+ yf = div_count
+ ya = yi
+ if (na.trapz(yvals[yi:yf], x=divs[yi:yf])/ynorm) >= nexty and not doney:
+ boundy.append(divs[ya])
+ if len(boundy) == cc[1]:
+ doney = True
+ nexty += 1./yedge.size
+ if zstart == 0:
+ zi = 0
+ zf = i
+ za = i
+ else:
+ zi = div_count - i
+ zf = div_count
+ za = zi
+ if (na.trapz(zvals[zi:zf], x=divs[zi:zf])/znorm) >= nextz and not donez:
+ boundz.append(divs[za])
+ if len(boundz) == cc[2]:
+ donez = True
+ nextz += 1./zedge.size
+
+ boundx.sort()
+ boundy.sort()
+ boundz.sort()
+
+ # Check for problems, fatally for now because I'm the only one using this
+ # and I don't mind that, it will help me fix things.
+ if len(boundx) < cc[0] or len(boundy) < cc[1] or len(boundz) < cc[2]:
+ print 'weighted stuff broken.'
+ print 'cc', cc
+ print len(boundx), len(boundy), len(boundz)
+ sys.exit()
+
+ boundx.append(1.)
+ boundy.append(1.)
+ boundz.append(1.)
+
+ if mi == 0:
+ print 'x',boundx
+ print 'y',boundy
+ print 'z',boundz
+
+ # Update the boundaries
+ new_LE = na.array([boundx[cx], boundy[cy], boundz[cz]], dtype='float64')
+ new_RE = na.array([boundx[cx+1], boundy[cy+1], boundz[cz+1]], dtype='float64')
+
+ new_vol = ((new_RE - new_LE) **2).sum()
+ print 'P%04d weight %f old_vol %f new_vol %f ratio %f' % \
+ (mi, weight, old_vol, new_vol, new_vol/old_vol)
+
+ if padding > 0:
+ return True, \
+ new_LE, new_RE, self.hierarchy.periodic_region_strict(self.center, new_LE-padding, new_RE+padding)
+
+ return False, new_LE, new_RE, self.hierarchy.region_strict(self.center, new_LE, new_RE)
+
+
+ def _mpi_find_neighbor_3d(self, shift):
+ """ Given a shift array, 1x3 long, find the task ID
+ of that neighbor. For example, shift=[1,0,0] finds the neighbor
+ immediately to the right in the positive x direction. Each task
+ has 26 neighbors, of which some may be itself depending on the number
+ and arrangement of tasks.
+ """
+ if not self._distributed: return 0
+ shift = na.array(shift)
+ cc = na.array(MPI.Compute_dims(MPI.COMM_WORLD.size, 3))
+ mi = MPI.COMM_WORLD.rank
+ si = MPI.COMM_WORLD.size
+ # store some facts about myself
+ mi_cx,mi_cy,mi_cz = na.unravel_index(mi,cc)
+ mi_ar = na.array([mi_cx,mi_cy,mi_cz])
+ # these are identical on all tasks
+ # should these be calculated once and stored?
+ #dLE = na.empty((si,3), dtype='float64') # positions not needed yet...
+ #dRE = na.empty((si,3), dtype='float64')
+ tasks = na.empty((cc[0],cc[1],cc[2]), dtype='int64')
+
+ for i in range(si):
+ cx,cy,cz = na.unravel_index(i,cc)
+ tasks[cx,cy,cz] = i
+ #x = na.mgrid[LE[0]:RE[0]:(cc[0]+1)*1j][cx:cx+2]
+ #y = na.mgrid[LE[1]:RE[1]:(cc[1]+1)*1j][cy:cy+2]
+ #z = na.mgrid[LE[2]:RE[2]:(cc[2]+1)*1j][cz:cz+2]
+ #dLE[i, :] = na.array([x[0], y[0], z[0]], dtype='float64')
+ #dRE[i, :] = na.array([x[1], y[1], z[1]], dtype='float64')
+
+ # find the neighbor
+ ne = (mi_ar + shift) % cc
+ ne = tasks[ne[0],ne[1],ne[2]]
+ return ne
+
def _barrier(self):
if not self._distributed: return
@@ -303,25 +780,17 @@
@parallel_passthrough
def _mpi_catdict(self, data):
- self._barrier()
field_keys = data.keys()
field_keys.sort()
- np = MPI.COMM_WORLD.size
+ size = data[field_keys[0]].shape[-1]
+ # MPI_Scan is an inclusive scan
+ sizes = MPI.COMM_WORLD.alltoall( [size]*MPI.COMM_WORLD.size )
+ offsets = na.add.accumulate([0] + sizes)[:-1]
+ arr_size = MPI.COMM_WORLD.allreduce(size, op=MPI.SUM)
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:
- temp_data = []
- if data[key] is not None: temp_data.append(data[key])
- for i in range(1,np):
- 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()
- data[key] = _bcast_array(data[key])
- self._barrier()
+ dd = data[key]
+ rv = _alltoallv_array(dd, arr_size, offsets, sizes)
+ data[key] = rv
return data
@parallel_passthrough
@@ -337,6 +806,324 @@
return data
@parallel_passthrough
+ def _mpi_joindict_unpickled_double(self, data):
+ self._barrier()
+ size = 0
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1,MPI.COMM_WORLD.size):
+ size = MPI.COMM_WORLD.recv(source=i, tag=0)
+ keys = na.empty(size, dtype='int64')
+ values = na.empty(size, dtype='float64')
+ MPI.COMM_WORLD.Recv([keys, MPI.LONG], i, 0)
+ MPI.COMM_WORLD.Recv([values, MPI.DOUBLE], i, 0)
+ for i,key in enumerate(keys):
+ data[key] = values[i]
+ # Now convert root's data to arrays.
+ size = len(data)
+ root_keys = na.empty(size, dtype='int64')
+ root_values = na.empty(size, dtype='float64')
+ count = 0
+ for key in data:
+ root_keys[count] = key
+ root_values[count] = data[key]
+ count += 1
+ else:
+ MPI.COMM_WORLD.send(len(data), 0, 0)
+ keys = na.empty(len(data), dtype='int64')
+ values = na.empty(len(data), dtype='float64')
+ count = 0
+ for key in data:
+ keys[count] = key
+ values[count] = data[key]
+ count += 1
+ MPI.COMM_WORLD.Send([keys, MPI.LONG], 0, 0)
+ MPI.COMM_WORLD.Send([values, MPI.DOUBLE], 0, 0)
+ # Now send it back as arrays.
+ size = MPI.COMM_WORLD.bcast(size, root=0)
+ if MPI.COMM_WORLD.rank != 0:
+ del keys, values
+ root_keys = na.empty(size, dtype='int64')
+ root_values = na.empty(size, dtype='float64')
+ MPI.COMM_WORLD.Bcast([root_keys, MPI.LONG], root=0)
+ MPI.COMM_WORLD.Bcast([root_values, MPI.DOUBLE], root=0)
+ # Convert back to a dict.
+ del data
+ data = dict(itertools.izip(root_keys, root_values))
+ return data
+
+ @parallel_passthrough
+ def _mpi_joindict_unpickled_long(self, data):
+ self._barrier()
+ size = 0
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1,MPI.COMM_WORLD.size):
+ size = MPI.COMM_WORLD.recv(source=i, tag=0)
+ keys = na.empty(size, dtype='int64')
+ values = na.empty(size, dtype='int64')
+ MPI.COMM_WORLD.Recv([keys, MPI.LONG], i, 0)
+ MPI.COMM_WORLD.Recv([values, MPI.LONG], i, 0)
+ for i,key in enumerate(keys):
+ data[key] = values[i]
+ # Now convert root's data to arrays.
+ size = len(data)
+ root_keys = na.empty(size, dtype='int64')
+ root_values = na.empty(size, dtype='int64')
+ count = 0
+ for key in data:
+ root_keys[count] = key
+ root_values[count] = data[key]
+ count += 1
+ else:
+ MPI.COMM_WORLD.send(len(data), 0, 0)
+ keys = na.empty(len(data), dtype='int64')
+ values = na.empty(len(data), dtype='int64')
+ count = 0
+ for key in data:
+ keys[count] = key
+ values[count] = data[key]
+ count += 1
+ MPI.COMM_WORLD.Send([keys, MPI.LONG], 0, 0)
+ MPI.COMM_WORLD.Send([values, MPI.LONG], 0, 0)
+ # Now send it back as arrays.
+ size = MPI.COMM_WORLD.bcast(size, root=0)
+ if MPI.COMM_WORLD.rank != 0:
+ del keys, values
+ root_keys = na.empty(size, dtype='int64')
+ root_values = na.empty(size, dtype='int64')
+ MPI.COMM_WORLD.Bcast([root_keys, MPI.LONG], root=0)
+ MPI.COMM_WORLD.Bcast([root_values, MPI.LONG], root=0)
+ # Convert back to a dict.
+ del data
+ data = dict(itertools.izip(root_keys,root_values))
+ return data
+
+ @parallel_passthrough
+ def _mpi_concatenate_array_long(self, data):
+ self._barrier()
+ size = 0
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1, MPI.COMM_WORLD.size):
+ size = MPI.COMM_WORLD.recv(source=i, tag=0)
+ new_data = na.empty(size, dtype='int64')
+ MPI.COMM_WORLD.Recv([new_data, MPI.LONG], i, 0)
+ data = na.concatenate((data, new_data))
+ size = data.size
+ del new_data
+ else:
+ MPI.COMM_WORLD.send(data.size, 0, 0)
+ MPI.COMM_WORLD.Send([data, MPI.LONG], 0, 0)
+ # Now we distribute the full array.
+ size = MPI.COMM_WORLD.bcast(size, root=0)
+ if MPI.COMM_WORLD.rank != 0:
+ del data
+ data = na.empty(size, dtype='int64')
+ MPI.COMM_WORLD.Bcast([data, MPI.LONG], root=0)
+ return data
+
+ @parallel_passthrough
+ def _mpi_concatenate_array_double(self, data):
+ self._barrier()
+ size = 0
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1, MPI.COMM_WORLD.size):
+ size = MPI.COMM_WORLD.recv(source=i, tag=0)
+ new_data = na.empty(size, dtype='float64')
+ MPI.COMM_WORLD.Recv([new_data, MPI.DOUBLE], i, 0)
+ data = na.concatenate((data, new_data))
+ size = data.size
+ del new_data
+ else:
+ MPI.COMM_WORLD.send(data.size, 0, 0)
+ MPI.COMM_WORLD.Send([data, MPI.DOUBLE], 0, 0)
+ # Now we distribute the full array.
+ size = MPI.COMM_WORLD.bcast(size, root=0)
+ if MPI.COMM_WORLD.rank != 0:
+ del data
+ data = na.empty(size, dtype='float64')
+ MPI.COMM_WORLD.Bcast([data, MPI.DOUBLE], root=0)
+ return data
+
+ @parallel_passthrough
+ def _mpi_minimum_array_long(self, data):
+ """
+ Specifically for parallelHOP. For the identical array on each task,
+ it merges the arrays together, taking the lower value at each index.
+ """
+ self._barrier()
+ size = data.size # They're all the same size, of course
+ if MPI.COMM_WORLD.rank == 0:
+ new_data = na.empty(size, dtype='int64')
+ for i in range(1, MPI.COMM_WORLD.size):
+ MPI.COMM_WORLD.Recv([new_data, MPI.LONG], i, 0)
+ data = na.minimum(data, new_data)
+ del new_data
+ else:
+ MPI.COMM_WORLD.Send([data, MPI.LONG], 0, 0)
+ # Redistribute from root
+ MPI.COMM_WORLD.Bcast([data, MPI.LONG], root=0)
+ return data
+
+ @parallel_passthrough
+ def _mpi_bcast_long_dict_unpickled(self, data):
+ self._barrier()
+ size = 0
+ if MPI.COMM_WORLD.rank == 0:
+ size = len(data)
+ size = MPI.COMM_WORLD.bcast(size, root=0)
+ root_keys = na.empty(size, dtype='int64')
+ root_values = na.empty(size, dtype='int64')
+ if MPI.COMM_WORLD.rank == 0:
+ count = 0
+ for key in data:
+ root_keys[count] = key
+ root_values[count] = data[key]
+ count += 1
+ MPI.COMM_WORLD.Bcast([root_keys, MPI.LONG], root=0)
+ MPI.COMM_WORLD.Bcast([root_values, MPI.LONG], root=0)
+ if MPI.COMM_WORLD.rank != 0:
+ data = {}
+ for i,key in enumerate(root_keys):
+ data[key] = root_values[i]
+ return data
+
+ @parallel_passthrough
+ def _mpi_maxdict(self, data):
+ """
+ For each key in data, find the maximum value across all tasks, and
+ then broadcast it back.
+ """
+ self._barrier()
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1,MPI.COMM_WORLD.size):
+ temp_data = MPI.COMM_WORLD.recv(source=i, tag=0)
+ for key in temp_data:
+ try:
+ old_value = data[key]
+ except KeyError:
+ # This guarantees the new value gets added.
+ old_value = None
+ if old_value < temp_data[key]:
+ data[key] = temp_data[key]
+ else:
+ MPI.COMM_WORLD.send(data, dest=0, tag=0)
+ data = MPI.COMM_WORLD.bcast(data, root=0)
+ self._barrier()
+ return data
+
+ def _mpi_maxdict_dict(self, data):
+ """
+ Similar to above, but finds maximums for dicts of dicts. This is
+ specificaly for a part of chainHOP.
+ """
+ if not self._distributed:
+ top_keys = []
+ bot_keys = []
+ vals = []
+ for top_key in data:
+ for bot_key in data[top_key]:
+ top_keys.append(top_key)
+ bot_keys.append(bot_key)
+ vals.append(data[top_key][bot_key])
+ top_keys = na.array(top_keys, dtype='int64')
+ bot_keys = na.array(bot_keys, dtype='int64')
+ vals = na.array(vals, dtype='float64')
+ return (top_keys, bot_keys, vals)
+ self._barrier()
+ size = 0
+ top_keys = []
+ bot_keys = []
+ vals = []
+ for top_key in data:
+ for bot_key in data[top_key]:
+ top_keys.append(top_key)
+ bot_keys.append(bot_key)
+ vals.append(data[top_key][bot_key])
+ top_keys = na.array(top_keys, dtype='int64')
+ bot_keys = na.array(bot_keys, dtype='int64')
+ vals = na.array(vals, dtype='float64')
+ del data
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1,MPI.COMM_WORLD.size):
+ size = MPI.COMM_WORLD.recv(source=i, tag=0)
+ mylog.info('Global Hash Table Merge %d of %d size %d' % \
+ (i,MPI.COMM_WORLD.size, size))
+ recv_top_keys = na.empty(size, dtype='int64')
+ recv_bot_keys = na.empty(size, dtype='int64')
+ recv_vals = na.empty(size, dtype='float64')
+ MPI.COMM_WORLD.Recv([recv_top_keys, MPI.LONG], source=i, tag=0)
+ MPI.COMM_WORLD.Recv([recv_bot_keys, MPI.LONG], source=i, tag=0)
+ MPI.COMM_WORLD.Recv([recv_vals, MPI.DOUBLE], source=i, tag=0)
+ top_keys = na.concatenate([top_keys, recv_top_keys])
+ bot_keys = na.concatenate([bot_keys, recv_bot_keys])
+ vals = na.concatenate([vals, recv_vals])
+# for j, top_key in enumerate(top_keys):
+# if j%1000 == 0: mylog.info(j)
+# # Make sure there's an entry for top_key in data
+# try:
+# test = data[top_key]
+# except KeyError:
+# data[top_key] = {}
+# try:
+# old_value = data[top_key][bot_keys[j]]
+# except KeyError:
+# # This guarantees the new value gets added.
+# old_value = None
+# if old_value < vals[j]:
+# data[top_key][bot_keys[j]] = vals[j]
+ else:
+# top_keys = []
+# bot_keys = []
+# vals = []
+# for top_key in data:
+# for bot_key in data[top_key]:
+# top_keys.append(top_key)
+# bot_keys.append(bot_key)
+# vals.append(data[top_key][bot_key])
+# top_keys = na.array(top_keys, dtype='int64')
+# bot_keys = na.array(bot_keys, dtype='int64')
+# vals = na.array(vals, dtype='float64')
+ size = top_keys.size
+ MPI.COMM_WORLD.send(size, dest=0, tag=0)
+ MPI.COMM_WORLD.Send([top_keys, MPI.LONG], dest=0, tag=0)
+ MPI.COMM_WORLD.Send([bot_keys, MPI.LONG], dest=0, tag=0)
+ MPI.COMM_WORLD.Send([vals, MPI.DOUBLE], dest=0, tag=0)
+ # Getting ghetto here, we're going to decompose the dict into arrays,
+ # send that, and then reconstruct it. When data is too big the pickling
+ # of the dict fails.
+ if MPI.COMM_WORLD.rank == 0:
+# data = defaultdict(dict)
+# for i,top_key in enumerate(top_keys):
+# try:
+# old = data[top_key][bot_keys[i]]
+# except KeyError:
+# old = None
+# if old < vals[i]:
+# data[top_key][bot_keys[i]] = vals[i]
+# top_keys = []
+# bot_keys = []
+# vals = []
+# for top_key in data:
+# for bot_key in data[top_key]:
+# top_keys.append(top_key)
+# bot_keys.append(bot_key)
+# vals.append(data[top_key][bot_key])
+# del data
+# top_keys = na.array(top_keys, dtype='int64')
+# bot_keys = na.array(bot_keys, dtype='int64')
+# vals = na.array(vals, dtype='float64')
+ size = top_keys.size
+ # Broadcast them using array methods
+ size = MPI.COMM_WORLD.bcast(size, root=0)
+ if MPI.COMM_WORLD.rank != 0:
+ top_keys = na.empty(size, dtype='int64')
+ bot_keys = na.empty(size, dtype='int64')
+ vals = na.empty(size, dtype='float64')
+ MPI.COMM_WORLD.Bcast([top_keys,MPI.LONG], root=0)
+ MPI.COMM_WORLD.Bcast([bot_keys,MPI.LONG], root=0)
+ MPI.COMM_WORLD.Bcast([vals, MPI.DOUBLE], root=0)
+ return (top_keys, bot_keys, vals)
+
+ @parallel_passthrough
def __mpi_recvlist(self, data):
# First we receive, then we make a new list.
data = ensure_list(data)
@@ -379,15 +1166,39 @@
self._barrier()
return data
+ @parallel_passthrough
+ def _mpi_bcast_pickled(self, data):
+ self._barrier()
+ data = MPI.COMM_WORLD.bcast(data, root=0)
+ return data
+
def _should_i_write(self):
if not self._distributed: return True
return (MPI.COMM_WORLD == 0)
- def _preload(self, grids, fields, queue):
+ def _preload(self, grids, fields, io_handler):
# This will preload if it detects we are parallel capable and
# if so, we load *everything* that we need. Use with some care.
+ mylog.debug("Preloading %s from %s grids", fields, len(grids))
if not self._distributed: return
- queue.preload(grids, fields)
+ io_handler.preload(grids, fields)
+
+ @parallel_passthrough
+ def _mpi_double_array_max(self,data):
+ """
+ Finds the na.maximum of a distributed array and returns the result
+ back to all. The array should be the same length on all tasks!
+ """
+ self._barrier()
+ if MPI.COMM_WORLD.rank == 0:
+ recv_data = na.empty(data.size, dtype='float64')
+ for i in xrange(1, MPI.COMM_WORLD.size):
+ MPI.COMM_WORLD.Recv([recv_data, MPI.DOUBLE], source=i, tag=0)
+ data = na.maximum(data, recv_data)
+ else:
+ MPI.COMM_WORLD.Send([data, MPI.DOUBLE], dest=0, tag=0)
+ MPI.COMM_WORLD.Bcast([data, MPI.DOUBLE], root=0)
+ return data
@parallel_passthrough
def _mpi_allsum(self, data):
@@ -396,6 +1207,70 @@
# relatively small ( < 1e7 elements )
return MPI.COMM_WORLD.allreduce(data, op=MPI.SUM)
+ @parallel_passthrough
+ def _mpi_Allsum_double(self, data):
+ self._barrier()
+ # Non-pickling float allsum of a float array, data.
+ temp = data.copy()
+ MPI.COMM_WORLD.Allreduce([temp, MPI.DOUBLE], [data, MPI.DOUBLE], op=MPI.SUM)
+ del temp
+ return data
+
+ @parallel_passthrough
+ def _mpi_Allsum_long(self, data):
+ self._barrier()
+ # Non-pickling float allsum of an int array, data.
+ temp = data.copy()
+ MPI.COMM_WORLD.Allreduce([temp, MPI.LONG], [data, MPI.LONG], op=MPI.SUM)
+ del temp
+ return data
+
+ @parallel_passthrough
+ def _mpi_allmax(self, data):
+ self._barrier()
+ return MPI.COMM_WORLD.allreduce(data, op=MPI.MAX)
+
+ @parallel_passthrough
+ def _mpi_allmin(self, data):
+ self._barrier()
+ return MPI.COMM_WORLD.allreduce(data, op=MPI.MIN)
+
+ ###
+ # Non-blocking stuff.
+ ###
+
+ def _mpi_Irecv_long(self, data, source):
+ if not self._distributed: return -1
+ return MPI.COMM_WORLD.Irecv([data, MPI.LONG], source, 0)
+
+ def _mpi_Irecv_double(self, data, source):
+ if not self._distributed: return -1
+ return MPI.COMM_WORLD.Irecv([data, MPI.DOUBLE], source, 0)
+
+ def _mpi_Isend_long(self, data, dest):
+ if not self._distributed: return -1
+ return MPI.COMM_WORLD.Isend([data, MPI.LONG], dest, 0)
+
+ def _mpi_Isend_double(self, data, dest):
+ if not self._distributed: return -1
+ return MPI.COMM_WORLD.Isend([data, MPI.DOUBLE], dest, 0)
+
+ def _mpi_Request_Waitall(self, hooks):
+ if not self._distributed: return
+ MPI.Request.Waitall(hooks)
+
+ ###
+ # End non-blocking stuff.
+ ###
+
+ def _mpi_get_size(self):
+ if not self._distributed: return None
+ return MPI.COMM_WORLD.size
+
+ def _mpi_get_rank(self):
+ if not self._distributed: return None
+ return MPI.COMM_WORLD.rank
+
def _mpi_info_dict(self, info):
if not self._distributed: return 0, {0:info}
self._barrier()
@@ -423,6 +1298,11 @@
obj._owner = MPI.COMM_WORLD.rank
obj._distributed = True
+ def _do_not_claim_object(self, obj):
+ if not self._distributed: return
+ obj._owner = -1
+ obj._distributed = True
+
def _write_on_root(self, fn):
if not self._distributed: return open(fn, "w")
if MPI.COMM_WORLD.rank == 0:
@@ -430,9 +1310,12 @@
else:
return cStringIO.StringIO()
- def _get_filename(self, prefix):
+ def _get_filename(self, prefix, rank=None):
if not self._distributed: return prefix
- return "%s_%03i" % (prefix, MPI.COMM_WORLD.rank)
+ if rank == None:
+ return "%s_%04i" % (prefix, MPI.COMM_WORLD.rank)
+ else:
+ return "%s_%04i" % (prefix, rank)
def _is_mine(self, obj):
if not obj._distributed: return True
@@ -470,3 +1353,25 @@
tmp = arr.view(__tocast)
MPI.COMM_WORLD.Bcast([tmp, MPI.CHAR], root=root)
return arr
+
+def _alltoallv_array(send, total_size, offsets, sizes):
+ if len(send.shape) > 1:
+ recv = []
+ for i in range(send.shape[0]):
+ recv.append(_alltoallv_array(send[i,:].copy(), total_size, offsets, sizes))
+ recv = na.array(recv)
+ return recv
+ offset = offsets[MPI.COMM_WORLD.rank]
+ tmp_send = send.view(__tocast)
+ recv = na.empty(total_size, dtype=send.dtype)
+ recv[offset:offset+send.size] = send[:]
+ dtr = send.dtype.itemsize / tmp_send.dtype.itemsize # > 1
+ soff = [0] * MPI.COMM_WORLD.size
+ ssize = [tmp_send.size] * MPI.COMM_WORLD.size
+ roff = [off * dtr for off in offsets]
+ rsize = [siz * dtr for siz in sizes]
+ tmp_recv = recv.view(__tocast)
+ MPI.COMM_WORLD.Alltoallv((tmp_send, (ssize, soff), MPI.CHAR),
+ (tmp_recv, (rsize, roff), MPI.CHAR))
+ return recv
+
Added: trunk/yt/lagos/ParticleIO.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/ParticleIO.py Fri Nov 20 20:27:29 2009
@@ -0,0 +1,123 @@
+"""
+The particle-IO handler
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2009 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from yt.funcs import *
+from yt.lagos import *
+
+particle_handler_registry = defaultdict()
+
+class ParticleIOHandler(object):
+ class __metaclass__(type):
+ def __init__(cls, name, b, d):
+ type.__init__(cls, name, b, d)
+ if hasattr(cls, "_source_type"):
+ particle_handler_registry[cls._source_type] = cls
+
+ _source_type = None
+ def __init__(self, pf, source):
+ self.pf = pf
+ self.data = {}
+ self.source = source
+
+ def __getitem__(self, key):
+ if key not in self.data:
+ self.get_data(key)
+ return self.data[key]
+
+ def __setitem__(self, key, val):
+ self.data[key] = val
+
+ def __delitem__(self, key):
+ del self.data[key]
+
+ def iter(self):
+ for val in self.data.keys(): yield val
+
+ def get_data(self, fields):
+ fields = ensure_list(fields)
+ self.source.get_data(fields)
+ for field in fields:
+ self[field] = self.source[field]
+
+particle_handler_registry.default_factory = lambda: ParticleIOHandler
+
+class ParticleIOHandlerRegion(ParticleIOHandler):
+ periodic = False
+ _source_type = "region"
+ def __init__(self, pf, source):
+ self.left_edge = source.left_edge
+ self.right_edge = source.right_edge
+ ParticleIOHandler.__init__(self, pf, source)
+
+ def get_data(self, fields):
+ mylog.info("Getting %s using ParticleIO" % str(fields))
+ fields = ensure_list(fields)
+ if not self.pf.h.io._particle_reader:
+ mylog.info("not self.pf.h.io._particle_reader")
+ return self.source.get_data(fields)
+ rtype = 0
+ DLE = na.array(self.pf["DomainLeftEdge"], dtype='float64')
+ DRE = na.array(self.pf["DomainRightEdge"], dtype='float64')
+ args = (na.array(self.left_edge), na.array(self.right_edge),
+ int(self.periodic), DLE, DRE)
+ count_list, grid_list = [], []
+ for grid in self.source._grids:
+ if grid.NumberOfParticles == 0: continue
+ grid_list.append(grid)
+ if self.source._is_fully_enclosed(grid):
+ count_list.append(grid.NumberOfParticles)
+ else:
+ count_list.append(-1)
+ # region type, left_edge, right_edge, periodic, grid_list
+ fields_to_read = []
+ conv_factors = []
+ for field in fields:
+ f = self.pf.field_info[field]
+ if f._particle_convert_function is None:
+ fields_to_read.append(field)
+ conv_factors.append(na.ones(len(grid_list), dtype='float64'))
+ else:
+ to_add = f.get_dependencies(pf = self.pf).requested
+ if len(to_add) != 1: raise KeyError
+ fields_to_read += to_add
+ conv_factors.append(
+ na.fromiter((f.particle_convert(g) for g in grid_list),
+ count=len(grid_list), dtype='float64'))
+ conv_factors = na.array(conv_factors).transpose()
+ self.conv_factors = conv_factors
+ rv = self.pf.h.io._read_particles(
+ fields_to_read, rtype, args, grid_list, count_list,
+ conv_factors)
+ for field, v in zip(fields, rv): self[field] = v
+
+class ParticleIOHandlerRegionStrict(ParticleIOHandlerRegion):
+ _source_type = "region_strict"
+
+class ParticleIOHandlerPeriodicRegion(ParticleIOHandlerRegion):
+ periodic = True
+ _source_type = "periodic_region"
+
+class ParticleIOHandlerPeriodicRegionStrict(ParticleIOHandlerPeriodicRegion):
+ _source_type = "periodic_region_strict"
Modified: trunk/yt/lagos/Profiles.py
==============================================================================
--- trunk/yt/lagos/Profiles.py (original)
+++ trunk/yt/lagos/Profiles.py Fri Nov 20 20:27:29 2009
@@ -72,7 +72,7 @@
def _initialize_parallel(self, fields):
g_objs = [g for g in self._get_grid_objs()]
self._preload(g_objs, self._get_dependencies(fields),
- self._data_source.hierarchy.queue)
+ self._data_source.hierarchy.io)
def _lazy_add_fields(self, fields, weight, accumulation):
self._ngrids = 0
Modified: trunk/yt/lagos/UniversalFields.py
==============================================================================
--- trunk/yt/lagos/UniversalFields.py (original)
+++ trunk/yt/lagos/UniversalFields.py Fri Nov 20 20:27:29 2009
@@ -53,20 +53,20 @@
# I violate it here in order to keep the name/func_name relationship
def _dx(field, data):
- return data.dx
- return na.ones(data.ActiveDimensions, dtype='float64') * data.dx
+ return data.dds[0]
+ return na.ones(data.ActiveDimensions, dtype='float64') * data.dds[0]
add_field('dx', function=_dx, display_field=False,
validators=[ValidateSpatial(0)])
def _dy(field, data):
- return data.dy
- return na.ones(data.ActiveDimensions, dtype='float64') * data.dy
+ return data.dds[1]
+ return na.ones(data.ActiveDimensions, dtype='float64') * data.dds[1]
add_field('dy', function=_dy, display_field=False,
validators=[ValidateSpatial(0)])
def _dz(field, data):
- return data.dz
- return na.ones(data.ActiveDimensions, dtype='float64') * data.dz
+ return data.dds[2]
+ return na.ones(data.ActiveDimensions, dtype='float64') * data.dds[2]
add_field('dz', function=_dz,
display_field=False, validators=[ValidateSpatial(0)])
@@ -130,29 +130,24 @@
add_field("SoundSpeed", function=_SoundSpeed,
units=r"\rm{cm}/\rm{s}")
-def particle_func(p_field, dtype='float64', ptype=None):
- dspace = None
- if ptype is not None:
- dspace = "AddressParticleType%02i" % ptype
+def particle_func(p_field, dtype='float64'):
def _Particles(field, data):
+ io = data.hierarchy.io
if not data.NumberOfParticles > 0:
return na.array([], dtype=dtype)
try:
- return data._read_data(p_field, dspace=dspace).astype(dtype)
- except data._read_exception:
+ return io._read_data_set(data, p_field).astype(dtype)
+ except io._read_exception:
pass
# This is bad. But it's the best idea I have right now.
- return data._read_data(p_field.replace("_"," "), dspace=dspace).astype(dtype)
+ return data._read_data(p_field.replace("_"," ")).astype(dtype)
return _Particles
-for ptype in [None] + range(10):
- for pf in ["type", "mass"] + \
- ["position_%s" % ax for ax in 'xyz']:
- pname = "particle_%s" % pf
- if ptype is not None: pname = "PT%02i_%s" % (ptype, pname)
- pfunc = particle_func("particle_%s" % (pf), ptype=ptype)
- add_field(pname, function=pfunc,
- validators = [ValidateSpatial(0)],
- particle_type=True)
+for pf in ["type", "mass"] + \
+ ["position_%s" % ax for ax in 'xyz']:
+ pfunc = particle_func("particle_%s" % (pf))
+ add_field("particle_%s" % pf, function=pfunc,
+ validators = [ValidateSpatial(0)],
+ particle_type=True)
def _convRetainInt(data):
return 1
@@ -178,7 +173,7 @@
validators = [ValidateSpatial(0),
ValidateDataField(pf)],
particle_type=True)
-add_field("particle mass", function=particle_func("particle_mass"),
+add_field("particle_mass", function=particle_func("particle_mass"),
validators=[ValidateSpatial(0)], particle_type=True)
add_field("Dark matter density", function=lambda a,b: None,
@@ -204,12 +199,17 @@
return data.convert("Density")*(data.convert("cm")**3.0)
def _convertParticleMassMsun(data):
return data.convert("Density")*((data.convert("cm")**3.0)/1.989e33)
+def _IOLevelParticleMassMsun(grid):
+ dd = dict(particle_mass = na.ones(1), CellVolumeCode=grid["CellVolumeCode"])
+ cf = (_ParticleMass(None, dd) * _convertParticleMassMsun(grid))[0]
+ return cf
add_field("ParticleMass",
function=_ParticleMass, validators=[ValidateSpatial(0)],
particle_type=True, convert_function=_convertParticleMass)
add_field("ParticleMassMsun",
function=_ParticleMass, validators=[ValidateSpatial(0)],
- particle_type=True, convert_function=_convertParticleMassMsun)
+ particle_type=True, convert_function=_convertParticleMassMsun,
+ particle_convert_function=_IOLevelParticleMassMsun)
def _RadialMachNumber(field, data):
"""M{|v|/t_sound}"""
@@ -372,7 +372,7 @@
convert_function=_convertCellMassMsun)
def _StarMass(field,data):
- return data["star_density"] * data["CellVolume"]
+ return data["star_density_pyx"] * data["CellVolume"]
add_field("StarMassMsun", units=r"M_{\odot}",
function=_StarMass,
convert_function=_convertCellMassMsun)
Modified: trunk/yt/lagos/__init__.py
==============================================================================
--- trunk/yt/lagos/__init__.py (original)
+++ trunk/yt/lagos/__init__.py Fri Nov 20 20:27:29 2009
@@ -88,9 +88,11 @@
from ClusterFiles import *
from ContourFinder import *
from Clump import *
+from ParticleIO import *
from BaseDataTypes import *
from BaseGridType import *
from EnzoRateData import *
+from ObjectFindingMixin import *
from HierarchyType import *
from OutputTypes import *
from Profiles import *
Modified: trunk/yt/lagos/hop/EnzoHop.c
==============================================================================
--- trunk/yt/lagos/hop/EnzoHop.c (original)
+++ trunk/yt/lagos/hop/EnzoHop.c Fri Nov 20 20:27:29 2009
@@ -24,6 +24,7 @@
//
#include "Python.h"
+#include "structmember.h"
#include <stdio.h>
#include <math.h>
#include <signal.h>
@@ -39,65 +40,81 @@
void regroup_main(float dens_outer, HC *my_comm);
static PyObject *_HOPerror;
-static PyObject *
-Py_EnzoHop(PyObject *obj, PyObject *args)
+int convert_particle_arrays(
+ PyObject *oxpos, PyObject *oypos, PyObject *ozpos, PyObject *omass,
+ PyArrayObject **xpos, PyArrayObject **ypos, PyArrayObject **zpos,
+ PyArrayObject **mass)
{
- PyObject *oxpos, *oypos, *ozpos,
- *omass;
-
- PyArrayObject *xpos, *ypos, *zpos,
- *mass;
- xpos=ypos=zpos=mass=NULL;
- npy_float64 totalmass = 0.0;
- float normalize_to = 1.0;
- float thresh = 160.0;
-
- int i;
-
- if (!PyArg_ParseTuple(args, "OOOO|ff",
- &oxpos, &oypos, &ozpos, &omass, &thresh, &normalize_to))
- return PyErr_Format(_HOPerror,
- "EnzoHop: Invalid parameters.");
/* First the regular source arrays */
- xpos = (PyArrayObject *) PyArray_FromAny(oxpos,
+ *xpos = (PyArrayObject *) PyArray_FromAny(oxpos,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
- if(!xpos){
+ if(!*xpos){
PyErr_Format(_HOPerror,
"EnzoHop: xpos didn't work.");
- goto _fail;
+ return -1;
}
- int num_particles = PyArray_SIZE(xpos);
+ int num_particles = PyArray_SIZE(*xpos);
- ypos = (PyArrayObject *) PyArray_FromAny(oypos,
+ *ypos = (PyArrayObject *) PyArray_FromAny(oypos,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
- if((!ypos)||(PyArray_SIZE(ypos) != num_particles)) {
+ if((!*ypos)||(PyArray_SIZE(*ypos) != num_particles)) {
PyErr_Format(_HOPerror,
"EnzoHop: xpos and ypos must be the same length.");
- goto _fail;
+ return -1;
}
- zpos = (PyArrayObject *) PyArray_FromAny(ozpos,
+ *zpos = (PyArrayObject *) PyArray_FromAny(ozpos,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
- if((!zpos)||(PyArray_SIZE(zpos) != num_particles)) {
+ if((!*zpos)||(PyArray_SIZE(*zpos) != num_particles)) {
PyErr_Format(_HOPerror,
"EnzoHop: xpos and zpos must be the same length.");
- goto _fail;
+ return -1;
}
- mass = (PyArrayObject *) PyArray_FromAny(omass,
+ *mass = (PyArrayObject *) PyArray_FromAny(omass,
PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
- if((!mass)||(PyArray_SIZE(mass) != num_particles)) {
+ if((!*mass)||(PyArray_SIZE(*mass) != num_particles)) {
PyErr_Format(_HOPerror,
"EnzoHop: xpos and mass must be the same length.");
- goto _fail;
+ return -1;
}
+ return num_particles;
+
+}
+
+
+static PyObject *
+Py_EnzoHop(PyObject *obj, PyObject *args)
+{
+ PyObject *oxpos, *oypos, *ozpos,
+ *omass;
+
+ PyArrayObject *xpos, *ypos, *zpos,
+ *mass;
+ xpos=ypos=zpos=mass=NULL;
+ npy_float64 totalmass = 0.0;
+ float normalize_to = 1.0;
+ float thresh = 160.0;
+
+ int i;
+
+ if (!PyArg_ParseTuple(args, "OOOO|ff",
+ &oxpos, &oypos, &ozpos, &omass, &thresh, &normalize_to))
+ return PyErr_Format(_HOPerror,
+ "EnzoHop: Invalid parameters.");
+
+ int num_particles = convert_particle_arrays(
+ oxpos, oypos, ozpos, omass,
+ &xpos, &ypos, &zpos, &mass);
+ if (num_particles < 0) goto _fail;
+
for(i = 0; i < num_particles; i++)
totalmass+=*(npy_float64*)PyArray_GETPTR1(mass,i);
totalmass /= normalize_to;
@@ -201,6 +218,209 @@
__declspec(dllexport)
#endif
+//
+// Now a fun wrapper class for the kD-tree
+//
+
+typedef struct {
+ PyObject_HEAD
+ KD kd;
+ PyArrayObject *xpos, *ypos, *zpos;
+ PyArrayObject *mass, *densities;
+ int num_particles;
+} kDTreeType;
+
+static int
+kDTreeType_init(kDTreeType *self, PyObject *args, PyObject *kwds)
+{
+ int nBuckets = 16, i;
+ float normalize_to = 1.0;
+ static char *kwlist[] = {"xpos", "ypos", "zpos", "mass",
+ "nbuckets", "norm", NULL};
+ PyObject *oxpos, *oypos, *ozpos,
+ *omass;
+ self->xpos=self->ypos=self->zpos=self->mass=NULL;
+
+
+ if (! PyArg_ParseTupleAndKeywords(args, kwds, "OOOO|if", kwlist,
+ &oxpos, &oypos, &ozpos, &omass,
+ &nBuckets, &normalize_to))
+ return -1; /* Should this give an error? */
+
+ kdInit(&self->kd, nBuckets);
+
+ self->num_particles = convert_particle_arrays(
+ oxpos, oypos, ozpos, omass,
+ &self->xpos, &self->ypos, &self->zpos, &self->mass);
+
+ self->kd->nActive = self->num_particles;
+ self->kd->p = malloc(sizeof(PARTICLE)*self->num_particles);
+ if (self->kd->p == NULL) {
+ fprintf(stderr, "failed allocating particles.\n");
+ goto _fail;
+ }
+
+ /* Now we set up our Density array */
+ self->densities = (PyArrayObject *)
+ PyArray_SimpleNewFromDescr(1, PyArray_DIMS(self->xpos),
+ PyArray_DescrFromType(NPY_FLOAT64));
+
+ npy_float64 totalmass = 0.0;
+ for(i= 0; i < self->num_particles; i++) {
+ self->kd->p[i].np_index = i;
+ *(npy_float64*)(PyArray_GETPTR1(self->densities, i)) = 0.0;
+ totalmass+=*(npy_float64*)PyArray_GETPTR1(self->mass,i);
+ }
+ totalmass /= normalize_to;
+
+
+ self->kd->np_masses = self->mass;
+ self->kd->np_pos[0] = self->xpos;
+ self->kd->np_pos[1] = self->ypos;
+ self->kd->np_pos[2] = self->zpos;
+ self->kd->np_densities = self->densities;
+ self->kd->totalmass = totalmass;
+
+ PrepareKD(self->kd);
+ kdBuildTree(self->kd);
+
+ return 0;
+
+ _fail:
+ Py_XDECREF(self->xpos);
+ Py_XDECREF(self->ypos);
+ Py_XDECREF(self->zpos);
+ Py_XDECREF(self->mass);
+
+ if(self->kd->p!=NULL)free(self->kd->p);
+
+ return -1;
+}
+
+static void
+kDTreeType_dealloc(kDTreeType *self)
+{
+ kdFinish(self->kd);
+ Py_XDECREF(self->xpos);
+ Py_XDECREF(self->ypos);
+ Py_XDECREF(self->zpos);
+ Py_XDECREF(self->mass);
+
+ self->ob_type->tp_free((PyObject*)self);
+}
+
+static PyObject *
+kDTreeType_up_pass(kDTreeType *self, PyObject *args) {
+ int iCell;
+
+ if (!PyArg_ParseTuple(args, "i", &iCell))
+ return PyErr_Format(_HOPerror,
+ "kDTree.up_pass: invalid parameters.");
+
+ if(iCell >= self->num_particles)
+ return PyErr_Format(_HOPerror,
+ "kDTree.up_pass: iCell cannot be >= num_particles!");
+
+ kdUpPass(self->kd, iCell);
+ return Py_None;
+}
+
+static PyObject *
+kDTreeType_median_jst(kDTreeType *self, PyObject *args) {
+ int d, l, u;
+
+ if (!PyArg_ParseTuple(args, "iii", &d, &l, &u))
+ return PyErr_Format(_HOPerror,
+ "kDTree.median_jst: invalid parameters.");
+
+ if(d >= 3)
+ return PyErr_Format(_HOPerror,
+ "kDTree.median_jst: d cannot be >= 3!");
+
+ if(l >= self->num_particles)
+ return PyErr_Format(_HOPerror,
+ "kDTree.median_jst: l cannot be >= num_particles!");
+
+ if(u >= self->num_particles)
+ return PyErr_Format(_HOPerror,
+ "kDTree.median_jst: u cannot be >= num_particles!");
+
+ int median = kdMedianJst(self->kd, d, l, u);
+
+ PyObject *omedian = PyInt_FromLong((long)median);
+ return omedian;
+}
+
+static PyMemberDef kDTreeType_members[] = {
+ { "xpos", T_OBJECT, offsetof(kDTreeType, xpos), 0,
+ "The xposition array."},
+ { "ypos", T_OBJECT, offsetof(kDTreeType, ypos), 0,
+ "The yposition array."},
+ { "zpos", T_OBJECT, offsetof(kDTreeType, zpos), 0,
+ "The zposition array."},
+ { "mass", T_OBJECT, offsetof(kDTreeType, mass), 0,
+ "The mass array."},
+ { "densities", T_OBJECT, offsetof(kDTreeType, densities), 0,
+ "The density array."},
+ { "num_particles", T_INT, offsetof(kDTreeType, num_particles), 0,
+ "The number of particles"},
+ { NULL }
+};
+
+static PyMethodDef
+kDTreeType_methods[] = {
+ { "up_pass", (PyCFunction) kDTreeType_up_pass, METH_VARARGS,
+ "Pass up something or another, I'm not really sure."},
+ { "median_jst", (PyCFunction) kDTreeType_median_jst, METH_VARARGS,
+ "Use the JST Median algorithm on two points along a dimension."},
+ // typically there would be more here...
+
+ { NULL }
+};
+
+static PyTypeObject
+kDTreeTypeDict = {
+ PyObject_HEAD_INIT(NULL)
+ 0, /* ob_size */
+ "kDTree", /* tp_name */
+ sizeof(kDTreeType), /* tp_basicsize */
+ 0, /* tp_itemsize */
+ (destructor)kDTreeType_dealloc, /* tp_dealloc */
+ 0, /* tp_print */
+ 0, /* tp_getattr */
+ 0, /* tp_setattr */
+ 0, /* tp_compare */
+ 0, /* tp_repr */
+ 0, /* tp_as_number */
+ 0, /* tp_as_sequence */
+ 0, /* tp_as_mapping */
+ 0, /* tp_hash */
+ 0, /* tp_call */
+ 0, /* tp_str */
+ 0, /* tp_getattro */
+ 0, /* tp_setattro */
+ 0, /* tp_as_buffer */
+ Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /* tp_flags*/
+ "kDTree object", /* tp_doc */
+ 0, /* tp_traverse */
+ 0, /* tp_clear */
+ 0, /* tp_richcompare */
+ 0, /* tp_weaklistoffset */
+ 0, /* tp_iter */
+ 0, /* tp_iternext */
+ kDTreeType_methods, /* tp_methods */
+ kDTreeType_members, /* tp_members */
+ 0, /* tp_getset */
+ 0, /* tp_base */
+ 0, /* tp_dict */
+ 0, /* tp_descr_get */
+ 0, /* tp_descr_set */
+ 0, /* tp_dictoffset */
+ (initproc)kDTreeType_init, /* tp_init */
+ 0, /* tp_alloc */
+ 0, /* tp_new */
+};
+
void initEnzoHop(void)
{
PyObject *m, *d;
@@ -208,7 +428,16 @@
d = PyModule_GetDict(m);
_HOPerror = PyErr_NewException("EnzoHop.HOPerror", NULL, NULL);
PyDict_SetItemString(d, "error", _HOPerror);
- import_array();
+
+ kDTreeTypeDict.tp_new = PyType_GenericNew;
+ if (PyType_Ready(&kDTreeTypeDict) < 0) {
+ return;
+ }
+
+ Py_INCREF(&kDTreeTypeDict);
+ PyModule_AddObject(m, "kDTree", (PyObject*)&kDTreeTypeDict);
+
+ import_array();
}
/*
Modified: trunk/yt/lagos/hop/hop_hop.c
==============================================================================
--- trunk/yt/lagos/hop/hop_hop.c (original)
+++ trunk/yt/lagos/hop/hop_hop.c Fri Nov 20 20:27:29 2009
@@ -62,7 +62,7 @@
int bDensity,bGroup,bSym,bMerge,nDens,nHop,nMerge,bTopHat;
float fDensThresh;
- nBucket = 16;
+ nBucket = 16;
nSmooth = 64;
nDens = 64;
nHop = -1;
@@ -117,7 +117,7 @@
if (bTopHat) smSmooth(smx,smDensityTH);
else if (bSym) smSmooth(smx,smDensitySym);
else smSmooth(smx,smDensity);
- } /* Else, we've read them */
+ } /* Else, we've read them */
if (bGroup) {
INFORM("Finding Densest Neighbors...\n");
if (bDensity && nHop<nSmooth) smReSmooth(smx,smHop);
Modified: trunk/yt/lagos/hop/hop_kd.c
==============================================================================
--- trunk/yt/lagos/hop/hop_kd.c (original)
+++ trunk/yt/lagos/hop/hop_kd.c Fri Nov 20 20:27:29 2009
@@ -47,6 +47,7 @@
kd = (KD)malloc(sizeof(struct kdContext));
assert(kd != NULL);
kd->nBucket = nBucket;
+ kd->kdNodes = NULL;
*pkd = kd;
return(1);
}
@@ -217,8 +218,8 @@
void kdFinish(KD kd)
{
- free(kd->p);
- free(kd->kdNodes);
+ if(kd->p!=NULL)free(kd->p);
+ if(kd->kdNodes!=NULL)free(kd->kdNodes);
free(kd);
}
Modified: trunk/yt/performance_counters.py
==============================================================================
--- trunk/yt/performance_counters.py (original)
+++ trunk/yt/performance_counters.py Fri Nov 20 20:27:29 2009
@@ -66,41 +66,33 @@
return func_wrapper
def print_stats(self):
- print "Current counter status:\n"
+ mylog.info("Current counter status:\n")
times = []
for i in self.counters:
insort(times, [self.starttime[i], i, 1]) # 1 for 'on'
if not self.counting[i]:
insort(times, [self.endtime[i], i, 0]) # 0 for 'off'
- #print times
- shift = -1
+ shifts = {}
+ order = []
+ endtimes = {}
+ shift = 0
multi = 5
- max = 20
- endline = ""
for i in times:
- # A starting entry
+ # a starting entry
if i[2] == 1:
+ shifts[i[1]] = shift
+ order.append(i[1])
shift += 1
- # An ending entry
if i[2] == 0:
- # if shift > 1, this is a nested entry, so we want to record
- # this line to be printed later when the top level finish entry
- # is encountered.
- if shift > 0:
- if self.counting[i[1]]:
- endline = "%s%i : %s : still running\n%s" % (" "*shift*multi,shift, i[1],endline)
- else:
- endline = "%s%i : %s : %0.3e\n%s" % (" "*shift*multi,shift, i[1], self.counters[i[1]], endline)
- shift -= 1
- # A top level entry.
- else:
- if self.counting[i[1]]:
- line = "%i : %s : still running\n%s" % (shift, i[1],endline)
- else:
- line = "%i : %s : %0.3e\n%s" % (shift, i[1], self.counters[i[1]],endline)
- shift -= 1
- endline = ""
- print line
+ shift -= 1
+ endtimes[i[1]] = self.counters[i[1]]
+ line = ''
+ for i in order:
+ if self.counting[i]:
+ line = "%s%s%i : %s : still running\n" % (line, " "*shifts[i]*multi, shifts[i], i)
+ else:
+ line = "%s%s%i : %s : %0.3e\n" % (line, " "*shifts[i]*multi, shifts[i], i, self.counters[i])
+ mylog.info("\n" + line)
def exit(self):
if self._on:
Modified: trunk/yt/raven/Callbacks.py
==============================================================================
--- trunk/yt/raven/Callbacks.py (original)
+++ trunk/yt/raven/Callbacks.py Fri Nov 20 20:27:29 2009
@@ -215,8 +215,8 @@
pxs, pys = na.mgrid[-1:1:3j,-1:1:3j]
else:
pxs, pys = na.mgrid[0:0:1j,0:0:1j]
- GLE = plot.data.gridLeftEdge
- GRE = plot.data.gridRightEdge
+ GLE = plot.data.grid_left_edge
+ GRE = plot.data.grid_right_edge
for px_off, py_off in zip(pxs.ravel(), pys.ravel()):
pxo = px_off * dom[px_index]
pyo = py_off * dom[py_index]
Modified: trunk/yt/raven/PlotTypes.py
==============================================================================
--- trunk/yt/raven/PlotTypes.py (original)
+++ trunk/yt/raven/PlotTypes.py Fri Nov 20 20:27:29 2009
@@ -403,7 +403,7 @@
self["Unit"] = str(unit)
self["Width"] = float(width)
if isinstance(unit, types.StringTypes):
- unit = self.data.hierarchy[str(unit)]
+ unit = self.data.pf[str(unit)]
self.width = width / unit
self._refresh_display_width()
More information about the yt-svn
mailing list