[Yt-svn] yt-commit r545 - in trunk: tests yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Sun Jun 8 23:16:57 PDT 2008
Author: mturk
Date: Sun Jun 8 23:16:55 2008
New Revision: 545
URL: http://yt.spacepope.org/changeset/545
Log:
Added an interpolated covering grid. It operates via a cascading interpolation
scheme, so the first level gets filled then interpolated to the second, which
freely overwrites the data from the previous level, and so on. It seems to
give good results. Some aliasing occurs, but to my visual inspection of the
data it is unavoidable in a first-order scheme. I added some unit tests for
it.
This is possibly not the final form the smoothing grid will take, but it's good
enough for now. The S2PLOT interface is going to move to this as soon as I can
test it on my work machine.
Modified:
trunk/tests/test_lagos.py
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/HelperFunctions.py
trunk/yt/lagos/HierarchyType.py
Modified: trunk/tests/test_lagos.py
==============================================================================
--- trunk/tests/test_lagos.py (original)
+++ trunk/tests/test_lagos.py Sun Jun 8 23:16:55 2008
@@ -90,7 +90,7 @@
self.assertEqual(ms,mr) # Asserting equality between the two
def testProjectionCorrectnessMultipleFields(self):
- p = self.hierarchy.proj(0,["Density","Ones"]) # Unweighted
+ p = self.hierarchy.proj(0,["Density","Ones"], weight=None) # Unweighted
self.assertTrue(na.all(p["Ones"] == 1.0))
def testProjectionMakingMultipleFields(self):
@@ -98,16 +98,24 @@
# One for each field, pdx, pdy, px, py, and one for the weight
self.assertEqual(len(p.data.keys()), 8)
- def testProjectionMaking(self):
+ def testProjectionSuccess(self):
p = self.hierarchy.proj(0,"Density") # Unweighted
p = self.hierarchy.proj(1,"Temperature","Density") # Weighted
p = self.hierarchy.proj(2,"Entropy") # Derived field
- def testProjectionCorrectness(self):
+ def testUnweightedProjectionCorrectness(self):
# Now we test that we get good answers
for axis in range(3):
p = self.hierarchy.proj(axis, "Ones") # Derived field
self.assertTrue(na.all(p["Ones"] == 1.0))
+ # Regardless of weighting, we want ones back
+
+ def testWeightedProjectionCorrectness(self):
+ # Now we test that we get good answers
+ for axis in range(3):
+ # Regardless of weighting, we want ones back
+ p = self.hierarchy.proj(axis, "Ones", "Density")
+ self.assertTrue(na.all(p["Ones"] == 1.0))
# Now we test each datatype in turn
@@ -207,6 +215,43 @@
accumulation_y)
setattr(Data3DBase, name, func)
+class TestSmoothedCoveringGrid(LagosTestingBase, unittest.TestCase):
+ def setUp(self):
+ LagosTestingBase.setUp(self)
+
+ def testAllCover(self):
+ DIMS = 32
+ for i in range(self.hierarchy.max_level+1):
+ dx = (DIMS*2**i)**-1
+ LE = na.array([0.5,0.5,0.5])-(dx*DIMS/2.0)
+ RE = na.array([0.5,0.5,0.5])+(dx*DIMS/2.0)
+ cg = self.hierarchy.smoothed_covering_grid(
+ level=i, left_edge=LE, right_edge=RE,
+ dims=[DIMS]*3, fields=["Density"])
+ self.assertFalse(na.any(na.isnan(cg["Density"])))
+ self.assertFalse(na.any(cg["Density"]==-999))
+
+ def testWrongDims(self):
+ self.failUnlessRaises(ValueError,
+ self.hierarchy.smoothed_covering_grid,
+ level=2, left_edge=[0.25, 0.25, 0.25],
+ right_edge=[0.75, 0.75, 0.75],
+ dims=[123, 124, 125])
+
+ def testAddField(self):
+ DIMS = 64
+ i = 5
+ dx = (DIMS*2**i)**-1
+ LE = na.array([0.5,0.5,0.5])-(dx*DIMS/2.0)
+ RE = na.array([0.5,0.5,0.5])+(dx*DIMS/2.0)
+ cg = self.hierarchy.smoothed_covering_grid(
+ level=i, left_edge=LE, right_edge=RE,
+ dims=[DIMS]*3, fields=["Density"])
+ self.assertFalse(na.any(na.isnan(cg["Temperature"])))
+ self.assertFalse(na.any(cg["Temperature"]==-999))
+
+
+
class TestDataCube(LagosTestingBase, unittest.TestCase):
def setUp(self):
LagosTestingBase.setUp(self)
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Sun Jun 8 23:16:55 2008
@@ -1516,3 +1516,101 @@
grid.LeftEdge, g_dx, g_fields, grid.child_mask,
self.left_edge, self.right_edge, c_dx, c_fields,
ll)
+
+class EnzoSmoothedCoveringGrid(EnzoCoveringGrid):
+ def __init__(self, *args, **kwargs):
+ dlog2 = na.log10(kwargs['dims'])/na.log10(2)
+ if not na.all(na.floor(dlog2) == na.ceil(dlog2)):
+ mylog.warning("Must be power of two dimensions")
+ raise ValueError
+ kwargs['num_ghost_zones'] = 0
+ EnzoCoveringGrid.__init__(self, *args, **kwargs)
+ if na.any(self.left_edge == 0):
+ self.left_edge += self.dx
+ self.ActiveDimensions -= 1
+ if na.any(self.right_edge == 0):
+ self.right_edge -= self.dx
+ self.ActiveDimensions -= 1
+
+ def _get_list_of_grids(self):
+ grids, ind = self.pf.hierarchy.get_box_grids(self.left_edge-self.dx,
+ self.right_edge+self.dx)
+ level_ind = na.where(self.pf.hierarchy.gridLevels.ravel()[ind] <= self.level)
+ sort_ind = na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind])
+ self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)]
+
+ def _get_level_array(self, level, fields):
+ fields = ensure_list(fields)
+ # We assume refinement by a factor of two
+ rf = 2**(self.level - level)
+ dims = na.maximum(1,self.ActiveDimensions/rf) + 2
+ dx = (self.right_edge-self.left_edge)/(dims-2)
+ x,y,z = (na.mgrid[0:dims[0],0:dims[1],0:dims[2]].astype('float64')+0.5)\
+ * dx[0]
+ x += self.left_edge[0] - dx[0]
+ y += self.left_edge[1] - dx[1]
+ z += self.left_edge[2] - dx[2]
+ bounds = [self.left_edge[0]-self['cdx'], self.right_edge[0]+self['cdx'],
+ self.left_edge[1]-self['cdy'], self.right_edge[1]+self['cdy'],
+ self.left_edge[2]-self['cdz'], self.right_edge[2]+self['cdz']]
+ fake_grid = {'x':x,'y':y,'z':z,'dx':dx[0],'dy':dx[1],'dz':dx[2]}
+ for ax in 'xyz': self['cd%s'%ax] = fake_grid['d%s'%ax]
+ for field in fields:
+ # Generate the new grid field
+ if field in fieldInfo and fieldInfo[field].take_log:
+ interpolator = TrilinearFieldInterpolator(
+ na.log10(self[field]), bounds, ['x','y','z'])
+ self[field] = 10**interpolator(fake_grid)
+ else:
+ interpolator = TrilinearFieldInterpolator(
+ self[field], bounds, ['x','y','z'])
+ self[field] = interpolator(fake_grid)
+ return fake_grid
+
+ def get_data(self, field=None):
+ self._get_list_of_grids()
+ # We don't generate coordinates here.
+ if field == None:
+ fields_to_get = self.fields
+ else:
+ fields_to_get = ensure_list(field)
+ for field in fields_to_get:
+ grid_count = 0
+ if self.data.has_key(field):
+ continue
+ mylog.debug("Getting field %s from %s possible grids",
+ field, len(self._grids))
+ self[field] = na.zeros(self.ActiveDimensions, dtype='float64') -999
+ if self._use_pbar: pbar = \
+ get_pbar('Searching grids for values ', len(self._grids))
+ # 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)
+ 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
+ for level in range(self.level+1):
+ for grid in self.select_grids(level):
+ if self._use_pbar: pbar.update(grid_count)
+ self._get_data_from_grid(grid, field)
+ grid_count += 1
+ if level < self.level: self._get_level_array(level+1, field)
+ self[field] = self[field][1:-1,1:-1,1:-1]
+ if self._use_pbar: pbar.finish()
+
+ @restore_grid_state
+ def _get_data_from_grid(self, grid, fields):
+ fields = ensure_list(fields)
+ g_dx = na.array([grid.dx, grid.dy, grid.dz])
+ c_dx = na.array([self['cdx'],self['cdy'],self['cdz']])
+ g_fields = [grid[field] for field in fields]
+ c_fields = [self[field] for field in fields]
+ total = PointCombine.DataCubeRefine(
+ grid.LeftEdge, g_dx, g_fields, grid.child_mask,
+ self.left_edge-c_dx, self.right_edge+c_dx,
+ c_dx, c_fields,
+ 1)
+
+ def flush_data(self, *args, **kwargs):
+ raise KeyError("Can't do this")
Modified: trunk/yt/lagos/HelperFunctions.py
==============================================================================
--- trunk/yt/lagos/HelperFunctions.py (original)
+++ trunk/yt/lagos/HelperFunctions.py Sun Jun 8 23:16:55 2008
@@ -37,7 +37,7 @@
orig_shape = data_object[self.x_name].shape
x_vals = data_object[self.x_name].ravel()
- x_i = na.digitize(data_object[self.x_name], self.x_bins) - 1
+ x_i = na.digitize(x_vals, self.x_bins) - 1
if na.any((x_i == -1) | (x_i == len(self.x_bins)-1)):
mylog.error("Sorry, but your values are outside" + \
" the table! Dunno what to do, so dying.")
@@ -63,8 +63,8 @@
x_vals = data_object[self.x_name].ravel()
y_vals = data_object[self.y_name].ravel()
- x_i = na.digitize(data_object[self.x_name], self.x_bins) - 1
- y_i = na.digitize(data_object[self.y_name], self.y_bins) - 1
+ x_i = na.digitize(x_vals, self.x_bins) - 1
+ y_i = na.digitize(y_vals, self.y_bins) - 1
if na.any((x_i == -1) | (x_i == len(self.x_bins)-1)) \
or na.any((y_i == -1) | (y_i == len(self.y_bins)-1)):
mylog.error("Sorry, but your values are outside" + \
@@ -98,9 +98,9 @@
y_vals = data_object[self.y_name].ravel()
z_vals = data_object[self.z_name].ravel()
- x_i = na.digitize(data_object[self.x_name], self.x_bins) - 1
- y_i = na.digitize(data_object[self.y_name], self.y_bins) - 1
- z_i = na.digitize(data_object[self.z_name], self.z_bins) - 1
+ x_i = na.digitize(x_vals, self.x_bins) - 1
+ y_i = na.digitize(y_vals, self.y_bins) - 1
+ z_i = na.digitize(z_vals, self.z_bins) - 1
if na.any((x_i == -1) | (x_i == len(self.x_bins)-1)) \
or na.any((y_i == -1) | (y_i == len(self.y_bins)-1)) \
or na.any((z_i == -1) | (z_i == len(self.z_bins)-1)):
@@ -117,6 +117,12 @@
xm = (self.x_bins[x_i+1] - x_vals) / (self.x_bins[x_i+1] - self.x_bins[x_i])
ym = (self.y_bins[y_i+1] - y_vals) / (self.y_bins[y_i+1] - self.y_bins[y_i])
zm = (self.z_bins[z_i+1] - z_vals) / (self.z_bins[z_i+1] - self.z_bins[z_i])
+ if na.any(na.isnan(self.table)):
+ raise ValueError
+ if na.any(na.isnan(x) | na.isnan(y) | na.isnan(z)):
+ raise ValueError
+ if na.any(na.isnan(xm) | na.isnan(ym) | na.isnan(zm)):
+ raise ValueError
my_vals = \
self.table[x_i ,y_i ,z_i ] * (xm*ym*zm) \
+ self.table[x_i+1,y_i ,z_i ] * (x *ym*zm) \
@@ -126,4 +132,4 @@
+ self.table[x_i ,y_i+1,z_i+1] * (xm*y *z ) \
+ self.table[x_i+1,y_i+1,z_i ] * (x *y *zm) \
+ self.table[x_i+1,y_i+1,z_i+1] * (x *y *z )
- return my_vals.reshape(orig_shape)
\ No newline at end of file
+ return my_vals.reshape(orig_shape)
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Sun Jun 8 23:16:55 2008
@@ -176,6 +176,7 @@
self.slice = classobj("EnzoSlice",(EnzoSliceBase,), dd)
self.region = classobj("EnzoRegion",(EnzoRegionBase,), dd)
self.covering_grid = classobj("EnzoCoveringGrid",(EnzoCoveringGrid,), dd)
+ self.smoothed_covering_grid = classobj("EnzoSmoothedCoveringGrid",(EnzoSmoothedCoveringGrid,), dd)
self.sphere = classobj("EnzoSphere",(EnzoSphereBase,), dd)
self.cutting = classobj("EnzoCuttingPlane",(EnzoCuttingPlaneBase,), dd)
self.ray = classobj("EnzoOrthoRay",(EnzoOrthoRayBase,), dd)
More information about the yt-svn
mailing list