Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Fixes #360 #361

Merged
merged 2 commits into from

2 participants

@scopatz
Owner

This issue fixes #360 and will hopefully put a lot of our IMesh Tag woes behind us. This also adds an iteration methods for Mesh that only yields the volume element (iter_ve()). This is useful for abstracting whether it is a structured mesh or not.

@scopatz scopatz added this to the v0.4 milestone
@scopatz
Owner

Mostly for @elliottbiondo

@elliottbiondo

This might just be par for the course, but I ran this issue when I merged this branch against my #359 branch:

In [1]: from pyne.mesh import Mesh, IMeshTag

In [2]: import numpy as np

In [3]: m = Mesh(structured=True, structured_coords=[[0,1,2],[0,1,2],[0,1,2]])

In [4]: m.catfish = IMeshTag(1, float)

In [5]: m.catfish[:] = np.ones(shape=(8,1))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-5-2d5a5a9f27cd> in <module>()
----> 1 m.catfish[:] = np.ones(shape=(8,1))

/home/elliott/.local/lib/python2.7/site-packages/pyne/mesh.pyc in __setitem__(self, key, value)
    379             v = np.empty(len(key), self.tag.type) if tsize == 1 else \
    380                 np.empty((len(key), tsize), self.tag.type)
--> 381             v[...] = value
    382             mtag[key] = v
    383         elif isinstance(key, np.ndarray) and key.dtype == np.bool:

ValueError: could not broadcast input array from shape (8,1) into shape (8)

In [6]: m.catfish[:] = [1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8]
@scopatz
Owner

Try using the correct shape: m.catfish[:] = np.ones(shape=(8,))

@elliottbiondo

Yes that works fine:

In [7]: m.catfish[:] = np.ones(shape=(8,)

But what about situations like:

result = np.empty(shape=(num_e_groups, num_vol_elements))

where num_e_groups might be 1, 42, or 175?

@elliottbiondo

Note that this also works:

In [8]: m.catfish[:] = np.ones(shape=(8,1))[0]

You can probably just add a conditional to treat this one corner case.

@elliottbiondo

Luckily all iMesh tag data is 0D or 1D, so it is not like this can happen for N dimensions.

@scopatz
Owner

Then the tag has a different size. This is actually an issue with iMesh / PyTAPS. For tag size == 1, the shape of the array has to be (len(key),). For tag size > 1, the shape of the array has to be (len(key), tag_size). We ensuring that without collapsing trivial dimensions. We could try to add some logic to do that, but I think it would trigger a whole bunch of bugs in the short term and is really just a convenience for people using the wrong shapes in the first place.

@scopatz
Owner

Luckily all iMesh tag data is 0D or 1D, so it is not like this can happen for N dimensions.

If it happened for N-dimensions it would be easier, because then they would have just used the normal numpy broadcasting rules :)

@scopatz
Owner

np.ones(shape=(8,1)) is not supported by MOAB for size 1 tags, so I am not sure that we should support it either...

@elliottbiondo

np.ones(shape=(8,1)) is not supported by MOAB for size 1 tags, so I am not sure that we should support it either..

Agreed.

@elliottbiondo

Everything else looks good and tests pass. Merged.

@elliottbiondo elliottbiondo merged commit 6ed1949 into pyne:staging
@scopatz
Owner

Sweet Thanks!

@scopatz scopatz deleted the scopatz:i360 branch
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Mar 8, 2014
  1. @scopatz
  2. @scopatz
This page is out of date. Refresh to see the latest.
Showing with 65 additions and 29 deletions.
  1. +43 −26 pyne/mesh.py
  2. +22 −3 tests/test_mesh.py
View
69 pyne/mesh.py
@@ -293,7 +293,7 @@ class IMeshTag(Tag):
features to this process.
"""
- def __init__(self, size=1, dtype='f8', mesh=None, name=None, doc=None):
+ def __init__(self, size=1, dtype='f8', default=0.0, mesh=None, name=None, doc=None):
"""Parameters
----------
size : int, optional
@@ -301,6 +301,9 @@ def __init__(self, size=1, dtype='f8', mesh=None, name=None, doc=None):
dtype : np.dtype or similar, optional
The data type of this tag from int, float, and byte. See PyTAPS
tags for more details.
+ default : dtype or None, optional
+ The default value to fill this tag with upon creation. If None, then
+ the tag is created empty.
mesh : Mesh, optional
The PyNE mesh to tag.
name : str, optional
@@ -313,11 +316,17 @@ def __init__(self, size=1, dtype='f8', mesh=None, name=None, doc=None):
if mesh is None or name is None:
self._lazy_args['size'] = size
self._lazy_args['dtype'] = dtype
- return
+ self._lazy_args['default'] = default
+ return
+ self.size = size
+ self.dtype = dtype
+ self.default = default
try:
self.tag = self.mesh.mesh.getTagHandle(self.name)
except iBase.TagNotFoundError:
self.tag = self.mesh.mesh.createTag(self.name, size, dtype)
+ if default is not None:
+ self[:] = default
def __delete__(self, mesh):
super(IMeshTag, self).__delete__(mesh)
@@ -327,7 +336,7 @@ def __getitem__(self, key):
m = self.mesh.mesh
size = len(self.mesh)
mtag = self.tag
- miter = m.iterate(iBase.Type.region, iMesh.Topology.all)
+ miter = self.mesh.iter_ve()
if isinstance(key, _INTEGRAL_TYPES):
if key >= size:
raise IndexError("key index {0} greater than the size of the "
@@ -349,32 +358,37 @@ def __getitem__(self, key):
"or fancy index.".format(key))
def __setitem__(self, key, value):
+ # get value into canonical form
+ tsize = self.size
+ value = np.asarray(value, self.tag.type)
+ value = np.atleast_1d(value) if tsize == 1 else np.atleast_2d(value)
+ # set up mesh to be iterated over
m = self.mesh.mesh
- size = len(self.mesh)
+ msize = len(self.mesh)
mtag = self.tag
- miter = m.iterate(iBase.Type.region, iMesh.Topology.all)
+ miter = self.mesh.iter_ve()
if isinstance(key, _INTEGRAL_TYPES):
- if key >= size:
+ if key >= msize:
raise IndexError("key index {0} greater than the size of the "
- "mesh {1}".format(key, size))
+ "mesh {1}".format(key, msize))
for i_ve in zip(range(key+1), miter):
pass
mtag[i_ve[1]] = value
elif isinstance(key, slice):
- idx = range(*key.indices(size))
- if not (isinstance(value, Sequence) and len(value) == len(idx)):
- value = [value] * len(idx)
- mtag[list(miter)[key]] = value
+ key = list(miter)[key]
+ v = np.empty(len(key), self.tag.type) if tsize == 1 else \
+ np.empty((len(key), tsize), self.tag.type)
+ v[...] = value
+ mtag[key] = v
elif isinstance(key, np.ndarray) and key.dtype == np.bool:
- if len(key) != size:
+ if len(key) != msize:
raise KeyError("boolean mask must match the length of the mesh.")
- ntrues = key.sum()
- if not (isinstance(value, Sequence) and len(value) == ntrues):
- value = [value] * ntrues
- mtag[[ve for b, ve in zip(key, miter) if b]] = value
+ key = [ve for b, ve in zip(key, miter) if b]
+ v = np.empty(len(key), self.tag.type) if tsize == 1 else \
+ np.empty((len(key), tsize), self.tag.type)
+ v[...] = value
+ mtag[key] = v
elif isinstance(key, Iterable):
- if not (isinstance(value, Sequence) and len(value) == len(key)):
- value = [value] * len(key)
ves = list(miter)
mtag[[ves[i] for i in key]] = value
else:
@@ -385,7 +399,7 @@ def __delitem__(self, key):
m = self.mesh.mesh
size = len(self.mesh)
mtag = self.tag
- miter = m.iterate(iBase.Type.region, iMesh.Topology.all)
+ miter = self.mesh.iter_ve()
if isinstance(key, _INTEGRAL_TYPES):
if key >= size:
raise IndexError("key index {0} greater than the size of the "
@@ -510,8 +524,7 @@ class Mesh(object):
def __init__(self, mesh=None, mesh_file=None, structured=False, \
structured_coords=None, structured_set=None,
structured_ordering='xyz', mats=None):
- """
- Parameters
+ """Parameters
----------
mesh : iMesh instance, optional
mesh_file : str, optional
@@ -649,11 +662,7 @@ def __init__(self, mesh=None, mesh_file=None, structured=False, \
self.mats = mats
# tag with volume id and ensure mats exist.
- if self.structured:
- ves = list(self.structured_iterate_hex(self.structured_ordering))
- else:
- ves = list(self.mesh.iterate(iBase.Type.region, iMesh.Topology.all))
-
+ ves = list(self.iter_ve())
tags = self.mesh.getAllTags(ves[0])
tags = set(tag.name for tag in tags)
if 'idx' in tags:
@@ -713,6 +722,14 @@ def __iter__(self):
iMesh.Topology.all)):
yield i, mats[i], ve
+ def iter_ve(self):
+ """Returns an iterator that yields on the volume elements.
+ """
+ if self.structured:
+ return self.structured_iterate_hex(self.structured_ordering)
+ else:
+ return self.mesh.iterate(iBase.Type.region, iMesh.Topology.all)
+
def __contains__(self, i):
return i < len(self)
View
25 tests/test_mesh.py
@@ -661,9 +661,7 @@ def test_imeshtag():
}
m = gen_mesh(mats=mats)
m.f = IMeshTag(mesh=m, name='f')
- ftag = m.mesh.getTagHandle('f')
- ftag[list(m.mesh.iterate(iBase.Type.region, iMesh.Topology.all))] = \
- [1.0, 2.0, 3.0, 4.0]
+ m.f[:] = [1.0, 2.0, 3.0, 4.0]
# Getting tags
assert_equal(m.f[0], 1.0)
@@ -733,6 +731,16 @@ def test_lazytaginit():
assert_in('cactus', m.tags)
assert_array_equal(m.cactus[0], [42, 43, 44])
+ x = np.arange(len(m))[:,np.newaxis] * np.array([42, 43, 44])
+ m.cactus[:] = x
+ assert_array_equal(m.cactus[2], x[2])
+
+def test_issue360():
+ a = Mesh(structured=True, structured_coords=[[0,1,2],[0,1],[0,1]])
+ a.cat = IMeshTag(3, float)
+ a.cat[:] = [[0.11, 0.22, 0.33],[0.44, 0.55, 0.66]]
+ a.cat[:] = np.array([[0.11, 0.22, 0.33],[0.44, 0.55, 0.66]])
+
def test_iter():
mats = {
0: Material({'H1': 1.0, 'K39': 1.0}, density=42.0),
@@ -748,6 +756,17 @@ def test_iter():
assert_is(mats[i], mat)
assert_equal(j, idx_tag[ve])
j += 1
+
+def test_iter_ve():
+ mats = {
+ 0: Material({'H1': 1.0, 'K39': 1.0}, density=42.0),
+ 1: Material({'H1': 0.1, 'O16': 1.0}, density=43.0),
+ 2: Material({'He4': 42.0}, density=44.0),
+ 3: Material({'Tm171': 171.0}, density=45.0),
+ }
+ m = gen_mesh(mats=mats)
+ ves1 = set(ve for _, _, ve in m)
+ ves2 = set(m.iter_ve())
def test_contains():
Something went wrong with that request. Please try again.