Merge branch 'develop' into URR_ptable_LCG_approach

Conflicts:
	tests/test_mgxs_library_distribcell/results_true.dat

openmc/arithmetic.py

@@ -1,5 +1,7 @@
 import sys
+import copy
 from numbers import Integral
+from collections import Iterable
 
 import numpy as np
 
@@ -14,7 +16,7 @@ if sys.version_info[0] >= 3:
 _TALLY_ARITHMETIC_OPS = ['+', '-', '*', '/', '^']
 
 # Acceptable tally aggregation operations
-_TALLY_AGGREGATE_OPS = ['sum', 'mean']
+_TALLY_AGGREGATE_OPS = ['sum', 'avg']
 
 
 class CrossScore(object):
@@ -186,6 +188,23 @@ class CrossNuclide(object):
             return existing
 
     def __repr__(self):
+        return self.name
+
+
+    @property
+    def left_nuclide(self):
+        return self._left_nuclide
+
+    @property
+    def right_nuclide(self):
+        return self._right_nuclide
+
+    @property
+    def binary_op(self):
+        return self._binary_op
+
+    @property
+    def name(self):
 
         string = ''
 
@@ -207,18 +226,6 @@ class CrossNuclide(object):
 
         return string
 
-    @property
-    def left_nuclide(self):
-        return self._left_nuclide
-
-    @property
-    def right_nuclide(self):
-        return self._right_nuclide
-
-    @property
-    def binary_op(self):
-        return self._binary_op
-
     @left_nuclide.setter
     def left_nuclide(self, left_nuclide):
         cv.check_type('left_nuclide', left_nuclide,
@@ -430,7 +437,7 @@ class CrossFilter(object):
         filter_index = left_index * self.right_filter.num_bins + right_index
         return filter_index
 
-    def get_pandas_dataframe(self, datasize, summary=None):
+    def get_pandas_dataframe(self, data_size, summary=None):
         """Builds a Pandas DataFrame for the CrossFilter's bins.
 
         This method constructs a Pandas DataFrame object for the CrossFilter
@@ -445,7 +452,7 @@ class CrossFilter(object):
 
         Parameters
         ----------
-        datasize : Integral
+        data_size : Integral
             The total number of bins in the tally corresponding to this filter
         summary : None or Summary
             An optional Summary object to be used to construct columns for
@@ -472,19 +479,18 @@ class CrossFilter(object):
 
         # If left and right filters are identical, do not combine bins
         if self.left_filter == self.right_filter:
-            df = self.left_filter.get_pandas_dataframe(datasize, summary)
+            df = self.left_filter.get_pandas_dataframe(data_size, summary)
 
         # If left and right filters are different, combine their bins
         else:
-            left_df = self.left_filter.get_pandas_dataframe(datasize, summary)
-            right_df = self.right_filter.get_pandas_dataframe(datasize, summary)
+            left_df = self.left_filter.get_pandas_dataframe(data_size, summary)
+            right_df = self.right_filter.get_pandas_dataframe(data_size, summary)
             left_df = left_df.astype(str)
             right_df = right_df.astype(str)
             df = '(' + left_df + ' ' + self.binary_op + ' ' + right_df + ')'
 
         return df
 
-
 class AggregateScore(object):
     """A special-purpose tally score used to encapsulate an aggregate of a
     subset or all of tally's scores for tally aggregation.
@@ -494,7 +500,7 @@ class AggregateScore(object):
     scores : Iterable of str or CrossScore
         The scores included in the aggregation
     aggregate_op : str
-        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        The tally aggregation operator (e.g., 'sum', 'avg', etc.) used
         to aggregate across a tally's scores with this AggregateScore
 
     Attributes
@@ -502,7 +508,7 @@ class AggregateScore(object):
     scores : Iterable of str or CrossScore
         The scores included in the aggregation
     aggregate_op : str
-        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        The tally aggregation operator (e.g., 'sum', 'avg', etc.) used
         to aggregate across a tally's scores with this AggregateScore
 
     """
@@ -556,10 +562,16 @@ class AggregateScore(object):
     def aggregate_op(self):
         return self._aggregate_op
 
+    @property
+    def name(self):
+
+        # Append each score in the aggregate to the string
+        string = '(' + ', '.join(self.scores) + ')'
+        return string
+
     @scores.setter
     def scores(self, scores):
-        cv.check_iterable_type('scores', scores,
-            (basestring, CrossScore, AggregateScore))
+        cv.check_iterable_type('scores', scores, basestring)
         self._scores = scores
 
     @aggregate_op.setter
@@ -578,7 +590,7 @@ class AggregateNuclide(object):
     nuclides : Iterable of str or Nuclide or CrossNuclide
         The nuclides included in the aggregation
     aggregate_op : str
-        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        The tally aggregation operator (e.g., 'sum', 'avg', etc.) used
         to aggregate across a tally's nuclides with this AggregateNuclide
 
     Attributes
@@ -586,7 +598,7 @@ class AggregateNuclide(object):
     nuclides : Iterable of str or Nuclide or CrossNuclide
         The nuclides included in the aggregation
     aggregate_op : str
-        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        The tally aggregation operator (e.g., 'sum', 'avg', etc.) used
         to aggregate across a tally's nuclides with this AggregateNuclide
 
     """
@@ -644,10 +656,19 @@ class AggregateNuclide(object):
     def aggregate_op(self):
         return self._aggregate_op
 
+    @property
+    def name(self):
+
+        # Append each nuclide in the aggregate to the string
+        names = [nuclide.name if isinstance(nuclide, Nuclide) else str(nuclide)
+                 for nuclide in self.nuclides]
+        string = '(' + ', '.join(map(str, names)) + ')'
+        return string
+
     @nuclides.setter
     def nuclides(self, nuclides):
         cv.check_iterable_type('nuclides', nuclides,
-           (basestring, Nuclide, CrossNuclide, AggregateNuclide))
+            (basestring, Nuclide, CrossNuclide))
         self._nuclides = nuclides
 
     @aggregate_op.setter
@@ -668,7 +689,7 @@ class AggregateFilter(object):
     bins : Iterable of tuple
         The filter bins included in the aggregation
     aggregate_op : str
-        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        The tally aggregation operator (e.g., 'sum', 'avg', etc.) used
         to aggregate across a tally filter's bins with this AggregateFilter
 
     Attributes
@@ -678,7 +699,7 @@ class AggregateFilter(object):
     aggregate_filter : filter
         The filter included in the aggregation
     aggregate_op : str
-        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        The tally aggregation operator (e.g., 'sum', 'avg', etc.) used
         to aggregate across a tally filter's bins with this AggregateFilter
     bins : Iterable of tuple
         The filter bins included in the aggregation
@@ -715,6 +736,21 @@ class AggregateFilter(object):
     def __ne__(self, other):
         return not self == other
 
+    def __gt__(self, other):
+        if self.type != other.type:
+            if self.aggregate_filter.type in _FILTER_TYPES and \
+              other.aggregate_filter.type in _FILTER_TYPES:
+                delta = _FILTER_TYPES.index(self.aggregate_filter.type) - \
+                        _FILTER_TYPES.index(other.aggregate_filter.type)
+                return delta > 0
+            else:
+                return False
+        else:
+            return False
+
+    def __lt__(self, other):
+        return not self > other
+
     def __repr__(self):
         string = 'AggregateFilter\n'
         string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type)
@@ -759,7 +795,7 @@ class AggregateFilter(object):
 
     @property
     def num_bins(self):
-        return 1 if self.aggregate_filter else 0
+        return len(self.bins) if self.aggregate_filter else 0
 
     @property
     def stride(self):
@@ -776,14 +812,13 @@ class AggregateFilter(object):
 
     @aggregate_filter.setter
     def aggregate_filter(self, aggregate_filter):
-        cv.check_type('aggregate_filter', aggregate_filter,
-            (Filter, CrossFilter, AggregateFilter))
+        cv.check_type('aggregate_filter', aggregate_filter, (Filter, CrossFilter))
         self._aggregate_filter = aggregate_filter
 
     @bins.setter
     def bins(self, bins):
-        cv.check_iterable_type('bins', bins, (Integral, tuple))
-        self._bins = bins
+        cv.check_iterable_type('bins', bins, Iterable)
+        self._bins = list(map(tuple, bins))
 
     @aggregate_op.setter
     def aggregate_op(self, aggregate_op):
@@ -823,15 +858,14 @@ class AggregateFilter(object):
 
         """
 
-        if filter_bin not in self.bins and \
-           filter_bin != self._aggregate_filter.bins:
+        if filter_bin not in self.bins:
             msg = 'Unable to get the bin index for AggregateFilter since ' \
                   '"{0}" is not one of the bins'.format(filter_bin)
             raise ValueError(msg)
         else:
-            return 0
+            return self.bins.index(filter_bin)
 
-    def get_pandas_dataframe(self, datasize, summary=None):
+    def get_pandas_dataframe(self, data_size, summary=None):
         """Builds a Pandas DataFrame for the AggregateFilter's bins.
 
         This method constructs a Pandas DataFrame object for the AggregateFilter
@@ -840,7 +874,7 @@ class AggregateFilter(object):
 
         Parameters
         ----------
-        datasize : Integral
+        data_size : Integral
             The total number of bins in the tally corresponding to this filter
         summary : None or Summary
             An optional Summary object to be used to construct columns for
@@ -868,14 +902,80 @@ class AggregateFilter(object):
 
         import pandas as pd
 
-        # Construct a sring representing the filter aggregation
-        aggregate_bin = '{0}('.format(self.aggregate_op)
-        aggregate_bin += ', '.join(map(str, self.bins)) + ')'
+        # Create NumPy array of the bin tuples for repeating / tiling
+        filter_bins = np.empty(self.num_bins, dtype=tuple)
+        for i, bin in enumerate(self.bins):
+            filter_bins[i] = bin
 
-        # Construct NumPy array of bin repeated for each element in dataframe
-        aggregate_bin_array = np.array([aggregate_bin])
-        aggregate_bin_array = np.repeat(aggregate_bin_array, datasize)
+        # Repeat and tile bins as needed for DataFrame
+        filter_bins = np.repeat(filter_bins, self.stride)
+        tile_factor = data_size / len(filter_bins)
+        filter_bins = np.tile(filter_bins, tile_factor)
 
-        # Construct Pandas DataFrame for the AggregateFilter
-        df = pd.DataFrame({self.type: aggregate_bin_array})
+        # Create DataFrame with aggregated bins
+        df = pd.DataFrame({self.type: filter_bins})
         return df
+
+    def can_merge(self, other):
+        """Determine if AggregateFilter can be merged with another.
+
+        Parameters
+        ----------
+        other : AggregateFilter
+            Filter to compare with
+
+        Returns
+        -------
+        bool
+            Whether the filter can be merged
+
+        """
+
+        if not isinstance(other, AggregateFilter):
+            return False
+
+        # Filters must be of the same type
+        elif self.type != other.type:
+            return False
+
+        # None of the bins in this filter should match in the other filter
+        for bin in self.bins:
+            if bin in other.bins:
+                return False
+
+        # If all conditional checks passed then filters are mergeable
+        return True
+
+    def merge(self, other):
+        """Merge this aggregatefilter with another.
+
+        Parameters
+        ----------
+        other : AggregateFilter
+            Filter to merge with
+
+        Returns
+        -------
+        merged_filter : AggregateFilter
+            Filter resulting from the merge
+
+        """
+
+        if not self.can_merge(other):
+            msg = 'Unable to merge "{0}" with "{1}" ' \
+                  'filters'.format(self.type, other.type)
+            raise ValueError(msg)
+
+        # Create deep copy of filter to return as merged filter
+        merged_filter = copy.deepcopy(self)
+
+        # Merge unique filter bins
+        merged_bins = self.bins + other.bins
+
+        # Sort energy bin edges
+        if 'energy' in self.type:
+            merged_bins = sorted(merged_bins)
+
+        # Assign merged bins to merged filter
+        merged_filter.bins = list(merged_bins)
+        return merged_filter

openmc/element.py

@@ -57,6 +57,15 @@ class Element(object):
     def __ne__(self, other):
         return not self == other
 
+    def __gt__(self, other):
+        return repr(self) > repr(other)
+
+    def __lt__(self, other):
+        return not self > other
+
+    def __hash__(self):
+        return hash(repr(self))
+
     def __hash__(self):
         return hash(repr(self))
 

openmc/filter.py

@@ -77,6 +77,25 @@ class Filter(object):
     def __ne__(self, other):
         return not self == other
 
+    def __gt__(self, other):
+        if self.type != other.type:
+            if self.type in _FILTER_TYPES and other.type in _FILTER_TYPES:
+                delta = _FILTER_TYPES.index(self.type) - \
+                        _FILTER_TYPES.index(other.type)
+                return delta > 0
+            else:
+                return False
+        else:
+            # Compare largest/smallest energy bin edges in energy filters
+            # This logic is used when merging tallies with energy filters
+            if 'energy' in self.type and 'energy' in other.type:
+                return self.bins[0] >= other.bins[-1]
+            else:
+                return max(self.bins) > max(other.bins)
+
+    def __lt__(self, other):
+        return not self > other
+
     def __hash__(self):
         return hash(repr(self))
 
@@ -246,20 +265,28 @@ class Filter(object):
             return False
 
         # Filters must be of the same type
-        elif self.type != other.type:
+        if self.type != other.type:
             return False
 
         # Distribcell filters cannot have more than one bin
-        elif self.type == 'distribcell':
+        if self.type == 'distribcell':
             return False
 
         # Mesh filters cannot have more than one bin
         elif self.type == 'mesh':
             return False
 
-        # Different energy bins are not mergeable
+        # Different energy bins structures must be mutually exclusive and
+        # share only one shared bin edge at the minimum or maximum energy
         elif 'energy' in self.type:
-            return False
+            # This low energy edge coincides with other's high energy edge
+            if self.bins[0] == other.bins[-1]:
+                return True
+            # This high energy edge coincides with other's low energy edge
+            elif self.bins[-1] == other.bins[0]:
+                return True
+            else:
+                return False
 
         else:
             return True
@@ -288,9 +315,21 @@ class Filter(object):
         merged_filter = copy.deepcopy(self)
 
         # Merge unique filter bins
-        merged_bins = list(set(np.concatenate((self.bins, other.bins))))
-        merged_filter.bins = merged_bins
-        merged_filter.num_bins = len(merged_bins)
+        merged_bins = np.concatenate((self.bins, other.bins))
+        merged_bins = np.unique(merged_bins)
+
+        # Sort energy bin edges
+        if 'energy' in self.type:
+            merged_bins = sorted(merged_bins)
+
+        # Assign merged bins to merged filter
+        merged_filter.bins = list(merged_bins)
+
+        # Count bins in the merged filter
+        if 'energy' in merged_filter.type:
+            merged_filter.num_bins = len(merged_bins) - 1
+        else:
+            merged_filter.num_bins = len(merged_bins)
 
         return merged_filter
 
@@ -521,14 +560,8 @@ class Filter(object):
 
         """
 
-        # Attempt to import Pandas
-        try:
-            import pandas as pd
-        except ImportError:
-            msg = 'The Pandas Python package must be installed on your system'
-            raise ImportError(msg)
-
         # Initialize Pandas DataFrame
+        import pandas as pd
         df = pd.DataFrame()
 
         # mesh filters
@@ -707,7 +740,6 @@ class Filter(object):
             filter_bins = np.repeat(filter_bins, self.stride)
             tile_factor = data_size / len(filter_bins)
             filter_bins = np.tile(filter_bins, tile_factor)
-            filter_bins = filter_bins
             df = pd.DataFrame({self.type : filter_bins})
 
             # If OpenCG level info DataFrame was created, concatenate

openmc/geometry.py

@@ -65,8 +65,10 @@ class Geometry(object):
 
         # Find the distribcell index of the cell.
         cells = self.get_all_cells()
-        if path[-1] in cells:
-            distribcell_index = cells[path[-1]].distribcell_index
+        for cell in cells:
+            if cell.id == path[-1]:
+                distribcell_index = cell.distribcell_index
+                break
         else:
             raise RuntimeError('Could not find cell {} specified in a \
                                 distribcell filter'.format(path[-1]))
@@ -94,7 +96,16 @@ class Geometry(object):
 
         """
 
-        return self._root_universe.get_all_cells()
+        all_cells = self._root_universe.get_all_cells()
+        cells = set()
+
+        for cell in all_cells.values():
+            if cell._type == 'normal':
+                cells.add(cell)
+
+        cells = list(cells)
+        cells.sort(key=lambda x: x.id)
+        return cells
 
     def get_all_universes(self):
         """Return all universes defined
@@ -106,7 +117,15 @@ class Geometry(object):
 
         """
 
-        return self._root_universe.get_all_universes()
+        all_universes = self._root_universe.get_all_universes()
+        universes = set()
+
+        for universe in all_universes.values():
+            universes.add(universe)
+
+        universes = list(universes)
+        universes.sort(key=lambda x: x.id)
+        return universes
 
     def get_all_nuclides(self):
         """Return all nuclides assigned to a material in the geometry
@@ -150,10 +169,19 @@ class Geometry(object):
         return materials
 
     def get_all_material_cells(self):
+        """Return all cells filled by a material
+
+        Returns
+        -------
+        list of openmc.universe.Cell
+            Cells filled by Materials in the geometry
+
+        """
+
         all_cells = self.get_all_cells()
         material_cells = set()
 
-        for cell_id, cell in all_cells.items():
+        for cell in all_cells:
             if cell._type == 'normal':
                 material_cells.add(cell)
 
@@ -174,9 +202,9 @@ class Geometry(object):
         all_universes = self.get_all_universes()
         material_universes = set()
 
-        for universe_id, universe in all_universes.items():
-            cells = universe._cells
-            for cell_id, cell in cells.items():
+        for universe in all_universes:
+            cells = universe.cells
+            for cell in cells:
                 if cell._type == 'normal':
                     material_universes.add(universe)
 
@@ -184,6 +212,227 @@ class Geometry(object):
         material_universes.sort(key=lambda x: x.id)
         return material_universes
 
+    def get_all_lattices(self):
+        """Return all lattices defined
+
+        Returns
+        -------
+        list of openmc.universe.Lattice
+            Lattices in the geometry
+
+        """
+
+        cells = self.get_all_cells()
+        lattices = set()
+
+        for cell in cells:
+            if isinstance(cell.fill, openmc.Lattice):
+                lattices.add(cell.fill)
+
+        lattices = list(lattices)
+        lattices.sort(key=lambda x: x.id)
+        return lattices
+
+    def get_materials_by_name(self, name, case_sensitive=False, matching=False):
+        """Return a list of materials with matching names.
+
+        Parameters
+        ----------
+        name : str
+            The name to match
+        case_sensitive : bool
+            Whether to distinguish upper and lower case letters in each
+            material's name (default is True)
+        matching : bool
+            Whether the names must match completely (default is True)
+
+        Returns
+        -------
+        list of openmc.material.Material
+            Materials matching the queried name
+
+        """
+
+        if not case_sensitive:
+            name = name.lower()
+
+        all_materials = self.get_all_materials()
+        materials = set()
+
+        for material in all_materials:
+            material_name = material.name
+            if not case_sensitive:
+                material_name = material_name.lower()
+
+            if material_name == name:
+                materials.add(material)
+            elif not matching and name in material_name:
+                materials.add(material)
+
+        materials = list(materials)
+        materials.sort(key=lambda x: x.id)
+        return materials
+
+    def get_cells_by_name(self, name, case_sensitive=False, matching=False):
+        """Return a list of cells with matching names.
+
+        Parameters
+        ----------
+        name : str
+            The name to search match
+        case_sensitive : bool
+            Whether to distinguish upper and lower case letters in each
+            cell's name (default is True)
+        matching : bool
+            Whether the names must match completely (default is True)
+
+        Returns
+        -------
+        list of openmc.universe.Cell
+            Cells matching the queried name
+
+        """
+
+        if not case_sensitive:
+            name = name.lower()
+
+        all_cells = self.get_all_cells()
+        cells = set()
+
+        for cell in all_cells:
+            cell_name = cell.name
+            if not case_sensitive:
+                cell_name = cell_name.lower()
+
+            if cell_name == name:
+                cells.add(cell)
+            elif not matching and name in cell_name:
+                cells.add(cell)
+
+        cells = list(cells)
+        cells.sort(key=lambda x: x.id)
+        return cells
+
+    def get_cells_by_fill_name(self, name, case_sensitive=False, matching=False):
+        """Return a list of cells with fills with matching names.
+
+        Parameters
+        ----------
+        name : str
+            The name to match
+        case_sensitive : bool
+            Whether to distinguish upper and lower case letters in each
+            cell's name (default is True)
+        matching : bool
+            Whether the names must match completely (default is True)
+
+        Returns
+        -------
+        list of openmc.universe.Cell
+            Cells with fills matching the queried name
+
+        """
+
+        if not case_sensitive:
+            name = name.lower()
+
+        all_cells = self.get_all_cells()
+        cells = set()
+
+        for cell in all_cells:
+            cell_fill_name = cell.fill.name
+            if not case_sensitive:
+                cell_fill_name = cell_fill_name.lower()
+
+            if cell_fill_name == name:
+                cells.add(cell)
+            elif not matching and name in cell_fill_name:
+                cells.add(cell)
+
+        cells = list(cells)
+        cells.sort(key=lambda x: x.id)
+        return cells
+
+    def get_universes_by_name(self, name, case_sensitive=False, matching=False):
+        """Return a list of universes with matching names.
+
+        Parameters
+        ----------
+        name : str
+            The name to match
+        case_sensitive : bool
+            Whether to distinguish upper and lower case letters in each
+            universe's name (default is True)
+        matching : bool
+            Whether the names must match completely (default is True)
+
+        Returns
+        -------
+        list of openmc.universe.Universe
+            Universes matching the queried name
+
+        """
+
+        if not case_sensitive:
+            name = name.lower()
+
+        all_universes = self.get_all_universes()
+        universes = set()
+
+        for universe in all_universes:
+            universe_name = universe.name
+            if not case_sensitive:
+                universe_name = universe_name.lower()
+
+            if universe_name == name:
+                universes.add(universe)
+            elif not matching and name in universe_name:
+                universes.add(universe)
+
+        universes = list(universes)
+        universes.sort(key=lambda x: x.id)
+        return universes
+
+    def get_lattices_by_name(self, name, case_sensitive=False, matching=False):
+        """Return a list of lattices with matching names.
+
+        Parameters
+        ----------
+        name : str
+            The name to match
+        case_sensitive : bool
+            Whether to distinguish upper and lower case letters in each
+            lattice's name (default is True)
+        matching : bool
+            Whether the names must match completely (default is True)
+
+        Returns
+        -------
+        list of openmc.universe.Lattice
+            Lattices matching the queried name
+
+        """
+
+        if not case_sensitive:
+            name = name.lower()
+
+        all_lattices = self.get_all_lattices()
+        lattices = set()
+
+        for lattice in all_lattices:
+            lattice_name = lattice.name
+            if not case_sensitive:
+                lattice_name = lattice_name.lower()
+
+            if lattice_name == name:
+                lattices.add(lattice)
+            elif not matching and name in lattice_name:
+                lattices.add(lattice)
+
+        lattices = list(lattices)
+        lattices.sort(key=lambda x: x.id)
+        return lattices
+
 
 class GeometryFile(object):
     """Geometry file used for an OpenMC simulation. Corresponds directly to the

openmc/mgxs/groups.py

@@ -54,10 +54,12 @@ class EnergyGroups(object):
     def __eq__(self, other):
         if not isinstance(other, EnergyGroups):
             return False
-        elif (self.group_edges != other.group_edges).all():
+        elif self.num_groups != other.num_groups:
             return False
-        else:
+        elif np.allclose(self.group_edges, other.group_edges):
             return True
+        else:
+            return False
 
     def __ne__(self, other):
         return not self == other
@@ -236,3 +238,64 @@ class EnergyGroups(object):
         condensed_groups.group_edges = group_edges
 
         return condensed_groups
+
+    def can_merge(self, other):
+        """Determine if energy groups can be merged with another.
+
+        Parameters
+        ----------
+        other : EnergyGroups
+            EnergyGroups to compare with
+
+        Returns
+        -------
+        bool
+            Whether the energy groups can be merged
+
+        """
+
+        if not isinstance(other, EnergyGroups):
+            return False
+
+        # If the energy group structures match then groups are mergeable
+        if self == other:
+            return True
+
+        # This low energy edge coincides with other's high energy edge
+        if self.group_edges[0] == other.group_edges[-1]:
+            return True
+        # This high energy edge coincides with other's low energy edge
+        elif self.group_edges[-1] == other.group_edges[0]:
+            return True
+        else:
+            return False
+
+    def merge(self, other):
+        """Merge this energy groups with another.
+
+        Parameters
+        ----------
+        other : EnergyGroups
+            EnergyGroups to merge with
+
+        Returns
+        -------
+        merged_groups : EnergyGroups
+            EnergyGroups resulting from the merge
+
+        """
+
+        if not self.can_merge(other):
+            raise ValueError('Unable to merge energy groups')
+
+        # Create deep copy to return as merged energy groups
+        merged_groups = copy.deepcopy(self)
+
+        # Merge unique filter bins
+        merged_edges = np.concatenate((self.group_edges, other.group_edges))
+        merged_edges = np.unique(merged_edges)
+        merged_edges = sorted(merged_edges)
+
+        # Assign merged edges to merged groups
+        merged_groups.group_edges = list(merged_edges)
+        return merged_groups

openmc/mgxs/library.py

@@ -116,11 +116,11 @@ class Library(object):
             clone._by_nuclide = self.by_nuclide
             clone._mgxs_types = self.mgxs_types
             clone._domain_type = self.domain_type
-            clone._domains = self.domains
+            clone._domains = copy.deepcopy(self.domains)
             clone._correction = self.correction
             clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
             clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
-            clone._all_mgxs = self.all_mgxs
+            clone._all_mgxs = copy.deepcopy(self.all_mgxs)
             clone._sp_filename = self._sp_filename
             clone._keff = self._keff
             clone._sparse = self.sparse

openmc/mgxs/mgxs.py

@@ -67,6 +67,10 @@ class MGXS(object):
         The energy group structure for energy condensation
     by_nuclide : bool
         If true, computes cross sections for each nuclide in domain
+    nuclides : Iterable of basestring
+        The user-specified nuclides to compute cross sections. If by_nuclide
+        is True but nuclides are not specified by the user, all nuclides in the
+        spatial domain will be used.
     name : str, optional
         Name of the multi-group cross section. Used as a label to identify
         tallies in OpenMC 'tallies.xml' file.
@@ -111,6 +115,8 @@ class MGXS(object):
     sparse : bool
         Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format
         for compressed data storage
+    derived : bool
+        Whether or not the MGXS is merged from one or more other MGXS
 
     """
 
@@ -123,6 +129,7 @@ class MGXS(object):
         self._name = ''
         self._rxn_type = None
         self._by_nuclide = None
+        self._nuclides = None
         self._domain = None
         self._domain_type = None
         self._energy_groups = None
@@ -131,6 +138,7 @@ class MGXS(object):
         self._rxn_rate_tally = None
         self._xs_tally = None
         self._sparse = False
+        self._derived = False
 
         self.name = name
         self.by_nuclide = by_nuclide
@@ -151,6 +159,7 @@ class MGXS(object):
             clone._name = self.name
             clone._rxn_type = self.rxn_type
             clone._by_nuclide = self.by_nuclide
+            clone._nuclides = copy.deepcopy(self._nuclides)
             clone._domain = self.domain
             clone._domain_type = self.domain_type
             clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
@@ -158,6 +167,7 @@ class MGXS(object):
             clone._rxn_rate_tally = copy.deepcopy(self._rxn_rate_tally, memo)
             clone._xs_tally = copy.deepcopy(self._xs_tally, memo)
             clone._sparse = self.sparse
+            clone._derived = self.derived
 
             clone._tallies = OrderedDict()
             for tally_type, tally in self.tallies.items():
@@ -232,8 +242,7 @@ class MGXS(object):
 
     @property
     def num_subdomains(self):
-        tally = list(self.tallies.values())[0]
-        domain_filter = tally.find_filter(self.domain_type)
+        domain_filter = self.xs_tally.find_filter(self.domain_type)
         return domain_filter.num_bins
 
     @property
@@ -250,6 +259,10 @@ class MGXS(object):
         else:
             return 'sum'
 
+    @property
+    def derived(self):
+        return self._derived
+
     @name.setter
     def name(self, name):
         cv.check_type('name', name, basestring)
@@ -260,6 +273,11 @@ class MGXS(object):
         cv.check_type('by_nuclide', by_nuclide, bool)
         self._by_nuclide = by_nuclide
 
+    @nuclides.setter
+    def nuclides(self, nuclides):
+        cv.check_iterable_type('nuclides', nuclides, basestring)
+        self._nuclides = nuclides
+
     @domain.setter
     def domain(self, domain):
         cv.check_type('domain', domain, tuple(_DOMAINS))
@@ -389,8 +407,14 @@ class MGXS(object):
         if self.domain is None:
             raise ValueError('Unable to get all nuclides without a domain')
 
-        nuclides = self.domain.get_all_nuclides()
-        return nuclides.keys()
+        # If the user defined nuclides, return them
+        if self._nuclides:
+            return self._nuclides
+
+        # Otherwise, return all nuclides in the spatial domain
+        else:
+            nuclides = self.domain.get_all_nuclides()
+            return nuclides.keys()
 
     def get_nuclide_density(self, nuclide):
         """Get the atomic number density in units of atoms/b-cm for a nuclide
@@ -553,7 +577,7 @@ class MGXS(object):
         # If computing xs for each nuclide, replace CrossNuclides with originals
         if self.by_nuclide:
             self.xs_tally._nuclides = []
-            nuclides = self.domain.get_all_nuclides()
+            nuclides = self.get_all_nuclides()
             for nuclide in nuclides:
                 self.xs_tally.add_nuclide(openmc.Nuclide(nuclide))
 
@@ -682,7 +706,7 @@ class MGXS(object):
 
         # Construct a collection of the domain filter bins
         if not isinstance(subdomains, basestring):
-            cv.check_iterable_type('subdomains', subdomains, Integral)
+            cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=2)
             for subdomain in subdomains:
                 filters.append(self.domain_type)
                 filter_bins.append((subdomain,))
@@ -855,19 +879,181 @@ class MGXS(object):
 
         # Clone this MGXS to initialize the subdomain-averaged version
         avg_xs = copy.deepcopy(self)
-        avg_xs._rxn_rate_tally = None
-        avg_xs._xs_tally = None
 
-        # Average each of the tallies across subdomains
-        for tally_type, tally in avg_xs.tallies.items():
-            tally_avg = tally.summation(filter_type=self.domain_type,
-                                        filter_bins=subdomains)
-            avg_xs.tallies[tally_type] = tally_avg
+        if self.derived:
+            avg_xs._rxn_rate_tally = avg_xs.rxn_rate_tally.average(
+                filter_type=self.domain_type, filter_bins=subdomains)
+        else:
+            avg_xs._rxn_rate_tally = None
+            avg_xs._xs_tally = None
 
-        avg_xs._domain_type = 'sum({0})'.format(self.domain_type)
+            # Average each of the tallies across subdomains
+            for tally_type, tally in avg_xs.tallies.items():
+                tally_avg = tally.average(filter_type=self.domain_type,
+                                          filter_bins=subdomains)
+                avg_xs.tallies[tally_type] = tally_avg
+
+        avg_xs._domain_type = 'avg({0})'.format(self.domain_type)
         avg_xs.sparse = self.sparse
         return avg_xs
 
+    def get_slice(self, nuclides=[], groups=[]):
+        """Build a sliced MGXS for the specified nuclides and energy groups.
+
+        This method constructs a new MGXS to encapsulate a subset of the data
+        represented by this MGXS. The subset of data to include in the tally
+        slice is determined by the nuclides and energy groups specified in
+        the input parameters.
+
+        Parameters
+        ----------
+        nuclides : list of str
+            A list of nuclide name strings
+            (e.g., ['U-235', 'U-238']; default is [])
+        groups : list of Integral
+            A list of energy group indices starting at 1 for the high energies
+            (e.g., [1, 2, 3]; default is [])
+
+        Returns
+        -------
+        MGXS
+            A new tally which encapsulates the subset of data requested for the
+            nuclide(s) and/or energy group(s) requested in the parameters.
+
+        """
+
+        cv.check_iterable_type('nuclides', nuclides, basestring)
+        cv.check_iterable_type('energy_groups', groups, Integral)
+
+        # Build lists of filters and filter bins to slice
+        if len(groups) == 0:
+            filters = []
+            filter_bins = []
+        else:
+            filter_bins = []
+            for group in groups:
+                group_bounds = self.energy_groups.get_group_bounds(group)
+                filter_bins.append(group_bounds)
+            filter_bins = [tuple(filter_bins)]
+            filters = ['energy']
+
+        # Clone this MGXS to initialize the sliced version
+        slice_xs = copy.deepcopy(self)
+        slice_xs._rxn_rate_tally = None
+        slice_xs._xs_tally = None
+
+        # Slice each of the tallies across nuclides and energy groups
+        for tally_type, tally in slice_xs.tallies.items():
+            slice_nuclides = [nuc for nuc in nuclides if nuc in tally.nuclides]
+            if len(groups) != 0 and tally.contains_filter('energy'):
+                tally_slice = tally.get_slice(filters=filters,
+                    filter_bins=filter_bins, nuclides=slice_nuclides)
+            else:
+                tally_slice = tally.get_slice(nuclides=slice_nuclides)
+            slice_xs.tallies[tally_type] = tally_slice
+
+        # Assign sliced energy group structure to sliced MGXS
+        if groups:
+            new_group_edges = []
+            for group in groups:
+                group_edges = self.energy_groups.get_group_bounds(group)
+                new_group_edges.extend(group_edges)
+            new_group_edges = np.unique(new_group_edges)
+            slice_xs.energy_groups.group_edges = sorted(new_group_edges)
+
+        # Assign sliced nuclides to sliced MGXS
+        if nuclides:
+            slice_xs.nuclides = nuclides
+
+        slice_xs.sparse = self.sparse
+        return slice_xs
+
+    def can_merge(self, other):
+        """Determine if another MGXS can be merged with this one
+
+        If results have been loaded from a statepoint, then MGXS are only
+        mergeable along one and only one of enegy groups or nuclides.
+
+        Parameters
+        ----------
+        other : MGXS
+            MGXS to check for merging
+
+        """
+
+        if not isinstance(other, type(self)):
+            return False
+
+        # Compare reaction type, energy groups, nuclides, domain type
+        if self.rxn_type != other.rxn_type:
+            return False
+        elif not self.energy_groups.can_merge(other.energy_groups):
+            return False
+        elif self.by_nuclide != other.by_nuclide:
+            return False
+        elif self.domain_type != other.domain_type:
+            return False
+        elif 'distribcell' not in self.domain_type and self.domain != other.domain:
+            return False
+        elif not self.xs_tally.can_merge(other.xs_tally):
+            return False
+        elif not self.rxn_rate_tally.can_merge(other.rxn_rate_tally):
+            return False
+
+        # If all conditionals pass then MGXS are mergeable
+        return True
+
+    def merge(self, other):
+        """Merge another MGXS with this one
+
+        MGXS are only mergeable if their energy groups and nuclides are either
+        identical or mutually exclusive. If results have been loaded from a
+        statepoint, then MGXS are only mergeable along one and only one of
+        energy groups or nuclides.
+
+        Parameters
+        ----------
+        other : MGXS
+            MGXS to merge with this one
+
+        Returns
+        -------
+        merged_mgxs : MGXS
+            Merged MGXS
+
+        """
+
+        if not self.can_merge(other):
+            raise ValueError('Unable to merge MGXS')
+
+        # Create deep copy of tally to return as merged tally
+        merged_mgxs = copy.deepcopy(self)
+        merged_mgxs._derived = True
+
+        # Merge energy groups
+        if self.energy_groups != other.energy_groups:
+            merged_groups = self.energy_groups.merge(other.energy_groups)
+            merged_mgxs.energy_groups = merged_groups
+
+        # Merge nuclides
+        if self.nuclides != other.nuclides:
+
+            # The nuclides must be mutually exclusive
+            for nuclide in self.nuclides:
+                if nuclide in other.nuclides:
+                    msg = 'Unable to merge MGXS with shared nuclides'
+                    raise ValueError(msg)
+
+            # Concatenate lists of nuclides for the merged MGXS
+            merged_mgxs.nuclides = self.nuclides + other.nuclides
+
+        # Null base tallies but merge reaction rate and cross section tallies
+        merged_mgxs._tallies = OrderedDict()
+        merged_mgxs._rxn_rate_tally = self.rxn_rate_tally.merge(other.rxn_rate_tally)
+        merged_mgxs._xs_tally = self.xs_tally.merge(other.xs_tally)
+
+        return merged_mgxs
+
     def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'):
         """Print a string representation for the multi-group cross section.
 
@@ -1022,6 +1208,9 @@ class MGXS(object):
             cv.check_iterable_type('subdomains', subdomains, Integral)
         elif self.domain_type == 'distribcell':
             subdomains = np.arange(self.num_subdomains, dtype=np.int)
+        elif self.domain_type == 'avg(distribcell)':
+            domain_filter = self.xs_tally.find_filter('avg(distribcell)')
+            subdomains = domain_filter.bins
         else:
             subdomains = [self.domain.id]
 
@@ -1239,14 +1428,13 @@ class MGXS(object):
             df.rename(columns={'energy low [MeV]': 'group in'},
                       inplace=True)
             in_groups = np.tile(all_groups, self.num_subdomains)
-            in_groups = np.repeat(in_groups, self.num_groups)
+            in_groups = np.repeat(in_groups, df.shape[0] / in_groups.size)
             df['group in'] = in_groups
             del df['energy high [MeV]']
 
             df.rename(columns={'energyout low [MeV]': 'group out'},
                       inplace=True)
-            out_groups = \
-                np.tile(all_groups, self.num_subdomains * self.num_groups)
+            out_groups = np.tile(all_groups, df.shape[0] / all_groups.size)
             df['group out'] = out_groups
             del df['energyout high [MeV]']
             columns = ['group in', 'group out']
@@ -1285,8 +1473,7 @@ class MGXS(object):
 
         # Sort the dataframe by domain type id (e.g., distribcell id) and
         # energy groups such that data is from fast to thermal
-        df.sort([self.domain_type] + columns, inplace=True)
-
+        df.sort_values(by=[self.domain_type] + columns, inplace=True)
         return df
 
 
@@ -1329,7 +1516,7 @@ class TotalXS(MGXS):
 
     @property
     def rxn_rate_tally(self):
-        if self._rxn_rate_tally is None:
+        if self._rxn_rate_tally is None :
             self._rxn_rate_tally = self.tallies['total']
             self._rxn_rate_tally.sparse = self.sparse
         return self._rxn_rate_tally
@@ -1735,6 +1922,58 @@ class ScatterMatrixXS(MGXS):
         cv.check_value('correction', correction, ('P0', None))
         self._correction = correction
 
+    def get_slice(self, nuclides=[], in_groups=[], out_groups=[]):
+        """Build a sliced ScatterMatrix for the specified nuclides and
+        energy groups.
+
+        This method constructs a new MGXS to encapsulate a subset of the data
+        represented by this MGXS. The subset of data to include in the tally
+        slice is determined by the nuclides and energy groups specified in
+        the input parameters.
+
+        Parameters
+        ----------
+        nuclides : list of str
+            A list of nuclide name strings
+            (e.g., ['U-235', 'U-238']; default is [])
+        in_groups : list of Integral
+            A list of incoming energy group indices starting at 1 for the high
+            energies (e.g., [1, 2, 3]; default is [])
+        out_groups : list of Integral
+            A list of outgoing energy group indices starting at 1 for the high
+            energies (e.g., [1, 2, 3]; default is [])
+
+        Returns
+        -------
+        MGXS
+            A new tally which encapsulates the subset of data requested for the
+            nuclide(s) and/or energy group(s) requested in the parameters.
+
+        """
+
+        # Call super class method and null out derived tallies
+        slice_xs = super(ScatterMatrixXS, self).get_slice(nuclides, in_groups)
+        slice_xs._rxn_rate_tally = None
+        slice_xs._xs_tally = None
+
+        # Slice outgoing energy groups if needed
+        if len(out_groups) != 0:
+            filter_bins = []
+            for group in out_groups:
+                group_bounds = self.energy_groups.get_group_bounds(group)
+                filter_bins.append(group_bounds)
+            filter_bins = [tuple(filter_bins)]
+
+            # Slice each of the tallies across energyout groups
+            for tally_type, tally in slice_xs.tallies.items():
+                if tally.contains_filter('energyout'):
+                    tally_slice = tally.get_slice(filters=['energyout'],
+                                                  filter_bins=filter_bins)
+                    slice_xs.tallies[tally_type] = tally_slice
+
+        slice_xs.sparse = self.sparse
+        return slice_xs
+
     def get_xs(self, in_groups='all', out_groups='all',
                subdomains='all', nuclides='all', xs_type='macro',
                order_groups='increasing', value='mean'):
@@ -1789,7 +2028,7 @@ class ScatterMatrixXS(MGXS):
 
         # Construct a collection of the domain filter bins
         if not isinstance(subdomains, basestring):
-            cv.check_iterable_type('subdomains', subdomains, Integral)
+            cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=2)
             for subdomain in subdomains:
                 filters.append(self.domain_type)
                 filter_bins.append((subdomain,))
@@ -2078,6 +2317,113 @@ class Chi(MGXS):
 
         return self._xs_tally
 
+    def get_slice(self, nuclides=[], groups=[]):
+        """Build a sliced Chi for the specified nuclides and energy groups.
+
+        This method constructs a new MGXS to encapsulate a subset of the data
+        represented by this MGXS. The subset of data to include in the tally
+        slice is determined by the nuclides and energy groups specified in
+        the input parameters.
+
+        Parameters
+        ----------
+        nuclides : list of str
+            A list of nuclide name strings
+            (e.g., ['U-235', 'U-238']; default is [])
+        groups : list of Integral
+            A list of energy group indices starting at 1 for the high energies
+            (e.g., [1, 2, 3]; default is [])
+
+        Returns
+        -------
+        MGXS
+            A new tally which encapsulates the subset of data requested for the
+            nuclide(s) and/or energy group(s) requested in the parameters.
+
+        """
+
+        # Temporarily remove energy filter from nu-fission-in since its
+        # group structure will work in super MGXS.get_slice(...) method
+        nu_fission_in = self.tallies['nu-fission-in']
+        energy_filter = nu_fission_in.find_filter('energy')
+        nu_fission_in.remove_filter(energy_filter)
+
+        # Call super class method and null out derived tallies
+        slice_xs = super(Chi, self).get_slice(nuclides, groups)
+        slice_xs._rxn_rate_tally = None
+        slice_xs._xs_tally = None
+
+        # Slice energy groups if needed
+        if len(groups) != 0:
+            filter_bins = []
+            for group in groups:
+                group_bounds = self.energy_groups.get_group_bounds(group)
+                filter_bins.append(group_bounds)
+            filter_bins = [tuple(filter_bins)]
+
+            # Slice nu-fission-out tally along energyout filter
+            nu_fission_out = slice_xs.tallies['nu-fission-out']
+            tally_slice = nu_fission_out.get_slice(filters=['energyout'],
+                                                   filter_bins=filter_bins)
+            slice_xs._tallies['nu-fission-out'] = tally_slice
+
+        # Add energy filter back to nu-fission-in tallies
+        self.tallies['nu-fission-in'].add_filter(energy_filter)
+        slice_xs._tallies['nu-fission-in'].add_filter(energy_filter)
+
+        slice_xs.sparse = self.sparse
+        return slice_xs
+
+    def merge(self, other):
+        """Merge another Chi with this one
+
+        If results have been loaded from a statepoint, then Chi are only
+        mergeable along one and only one of energy groups or nuclides.
+
+        Parameters
+        ----------
+        other : MGXS
+            MGXS to merge with this one
+
+        Returns
+        -------
+        merged_mgxs : MGXS
+            Merged MGXS
+        """
+
+        if not self.can_merge(other):
+            raise ValueError('Unable to merge Chi')
+
+        # Create deep copy of tally to return as merged tally
+        merged_mgxs = copy.deepcopy(self)
+        merged_mgxs._derived = True
+        merged_mgxs._rxn_rate_tally = None
+        merged_mgxs._xs_tally = None
+
+        # Merge energy groups
+        if self.energy_groups != other.energy_groups:
+            merged_groups = self.energy_groups.merge(other.energy_groups)
+            merged_mgxs.energy_groups = merged_groups
+
+        # Merge nuclides
+        if self.nuclides != other.nuclides:
+
+            # The nuclides must be mutually exclusive
+            for nuclide in self.nuclides:
+                if nuclide in other.nuclides:
+                    msg = 'Unable to merge Chi with shared nuclides'
+                    raise ValueError(msg)
+
+            # Concatenate lists of nuclides for the merged MGXS
+            merged_mgxs.nuclides = self.nuclides + other.nuclides
+
+        # Merge tallies
+        for tally_key in self.tallies:
+            merged_tally = self.tallies[tally_key].merge(other.tallies[tally_key])
+            merged_mgxs.tallies[tally_key] = merged_tally
+
+        return merged_mgxs
+
     def get_xs(self, groups='all', subdomains='all', nuclides='all',
                xs_type='macro', order_groups='increasing', value='mean'):
         """Returns an array of the fission spectrum.
@@ -2129,7 +2475,7 @@ class Chi(MGXS):
 
         # Construct a collection of the domain filter bins
         if not isinstance(subdomains, basestring):
-            cv.check_iterable_type('subdomains', subdomains, Integral)
+            cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=2)
             for subdomain in subdomains:
                 filters.append(self.domain_type)
                 filter_bins.append((subdomain,))

openmc/nuclide.py

@@ -60,6 +60,12 @@ class Nuclide(object):
     def __ne__(self, other):
         return not self == other
 
+    def __gt__(self, other):
+        return repr(self) > repr(other)
+
+    def __lt__(self, other):
+        return not self > other
+
     def __hash__(self):
         return hash(repr(self))
 

openmc/tallies.py

@@ -302,7 +302,7 @@ class Tally(object):
 
     @property
     def sum(self):
-        if not self._sp_filename:
+        if not self._sp_filename or self.derived:
             return None
 
         if not self._results_read:
@@ -674,72 +674,193 @@ class Tally(object):
 
         self._nuclides.remove(nuclide)
 
-    def can_merge(self, tally):
-        """Determine if another tally can be merged with this one
+    def _can_merge_filters(self, other):
+        """Determine if another tally's filters can be merged with this one's
+
+        The types of filters between the two tallies must match identically.
+        The bins in all of the filters must match identically, or be mergeable
+        in only one filter. This is a helper method for the can_merge(...)
+        and merge(...) methods.
 
         Parameters
         ----------
-        tally : Tally
-            Tally to check for merging
+        other : Tally
+            Tally to check for mergeable filters
 
         """
 
-        if not isinstance(tally, Tally):
+        # Two tallys must have the same number of filters
+        if len(self.filters) != len(other.filters):
             return False
 
-        # Must have same estimator
-        if self.estimator != tally.estimator:
-            return False
-
-        # Must have same nuclides
-        if len(self.nuclides) != len(tally.nuclides):
-            return False
-
-        for nuclide in self.nuclides:
-            if nuclide not in tally.nuclides:
-                return False
-
-        # Must have same or mergeable filters
-        if len(self.filters) != len(tally.filters):
-            return False
-
-        # Check if only one tally contains a delayed group filter
-        tally1_dg = False
-        for filter1 in self.filters:
-            if filter1.type == 'delayedgroup':
-                tally1_dg = True
-
-        tally2_dg = False
-        for filter2 in tally.filters:
-            if filter2.type == 'delayedgroup':
-                tally2_dg = True
-
         # Return False if only one tally has a delayed group filter
-        if (tally1_dg or tally2_dg) and not (tally1_dg and tally2_dg):
+        tally1_dg = self.contains_filter('delayedgroup')
+        tally2_dg = other.contains_filter('delayedgroup')
+        if sum([tally1_dg, tally2_dg]) == 1:
             return False
 
         # Look to see if all filters are the same, or one or more can be merged
         for filter1 in self.filters:
+            merge_filters = False
             mergeable_filter = False
 
-            for filter2 in tally.filters:
-                if filter1 == filter2 or filter1.can_merge(filter2):
+            for filter2 in other.filters:
+
+                # If filters match, they are mergeable
+                if filter1 == filter2:
                     mergeable_filter = True
                     break
 
+                # If filters are first mergeable filters encountered
+                elif filter1.can_merge(filter2) and not merge_filters:
+                    merge_filters = True
+                    mergeable_filter = True
+                    break
+
+                # If filters are the second mergeable filters encountered
+                elif filter1.can_merge(filter2) and merge_filters:
+                    return False
+
             # If no mergeable filter was found, the tallies are not mergeable
             if not mergeable_filter:
                 return False
 
-        # Tallies are mergeable if all conditional checks passed
+        # Tally filters are mergeable if all conditional checks passed
         return True
 
-    def merge(self, tally):
-        """Merge another tally with this one
+    def _can_merge_nuclides(self, other):
+        """Determine if another tally's nuclides can be merged with this one's
+
+        The nuclides between the two tallies must be mutually exclusive or
+        identically matching. This is a helper method for the can_merge(...)
+        and merge(...) methods.
 
         Parameters
         ----------
-        tally : Tally
+        other : Tally
+            Tally to check for mergeable nuclides
+
+        """
+
+        no_nuclides_match = True
+        all_nuclides_match = True
+
+        # Search for each of this tally's nuclides in the other tally
+        for nuclide in self.nuclides:
+            if nuclide not in other.nuclides:
+                all_nuclides_match = False
+            else:
+                no_nuclides_match = False
+
+        # Search for each of the other tally's nuclides in this tally
+        for nuclide in other.nuclides:
+            if nuclide not in self.nuclides:
+                all_nuclides_match = False
+            else:
+                no_nuclides_match = False
+
+        # Either all nuclides should match, or none should
+        if no_nuclides_match or all_nuclides_match:
+            return True
+        else:
+            return False
+
+    def _can_merge_scores(self, other):
+        """Determine if another tally's scores can be merged with this one's
+
+        The scores between the two tallies must be mutually exclusive or
+        identically matching. This is a helper method for the can_merge(...)
+        and merge(...) methods.
+
+        Parameters
+        ----------
+        other : Tally
+            Tally to check for mergeable scores
+
+        """
+
+        no_scores_match = True
+        all_scores_match = True
+
+        # Search for each of this tally's scores in the other tally
+        for score in self.scores:
+            if score not in other.scores:
+                all_scores_match = False
+            else:
+                no_scores_match = False
+
+        # Search for each of the other tally's scores in this tally
+        for score in other.scores:
+            if score not in self.scores:
+                all_scores_match = False
+            else:
+                no_scores_match = False
+
+        # Nuclides cannot be specified on 'flux' scores
+        if 'flux' in self.scores or 'flux' in other.scores:
+            if self.nuclides != other.nuclides:
+                return False
+
+        # Either all scores should match, or none should
+        if no_scores_match or all_scores_match:
+            return True
+        else:
+            return False
+
+    def can_merge(self, other):
+        """Determine if another tally can be merged with this one
+
+        If results have been loaded from a statepoint, then tallies are only
+        mergeable along one and only one of filter bins, nuclides or scores.
+
+        Parameters
+        ----------
+        other : Tally
+            Tally to check for merging
+
+        """
+
+        if not isinstance(other, Tally):
+            return False
+
+        # Must have same estimator
+        if self.estimator != other.estimator:
+            return False
+
+        equal_filters = sorted(self.filters) == sorted(other.filters)
+        equal_nuclides = sorted(self.nuclides) == sorted(other.nuclides)
+        equal_scores = sorted(self.scores) == sorted(other.scores)
+        equality = [equal_filters, equal_nuclides, equal_scores]
+
+        # If all filters, nuclides and scores match then tallies are mergeable
+        if equal_filters and equal_nuclides and equal_scores:
+            return True
+
+        # Variables to indicate matching filter bins, nuclides and scores
+        merge_filters = self._can_merge_filters(other)
+        merge_nuclides = self._can_merge_nuclides(other)
+        merge_scores = self._can_merge_scores(other)
+        mergeability = [merge_filters, merge_nuclides, merge_scores]
+
+        if not all(mergeability):
+            return False
+
+        # If the tally results have been read from the statepoint, we can only
+        # at least two of filters, nuclides and scores must match
+        elif self._results_read and sum(equality) < 2:
+            return False
+        else:
+            return True
+
+    def merge(self, other):
+        """Merge another tally with this one
+
+        If results have been loaded from a statepoint, then tallies are only
+        mergeable along one and only one of filter bins, nuclides or scores.
+
+        Parameters
+        ----------
+        other : Tally
             Tally to merge with this one
 
         Returns
@@ -749,9 +870,9 @@ class Tally(object):
 
         """
 
-        if not self.can_merge(tally):
-            msg = 'Unable to merge tally ID="{0}" with ' + \
-                   '"{1}"'.format(tally.id, self.id)
+        if not self.can_merge(other):
+            msg = 'Unable to merge tally ID="{0}" with ' \
+                   '"{1}"'.format(other.id, self.id)
             raise ValueError(msg)
 
         # Create deep copy of tally to return as merged tally
@@ -760,23 +881,127 @@ class Tally(object):
         # Differentiate Tally with a new auto-generated Tally ID
         merged_tally.id = None
 
-        # Merge filters
-        for i, filter1 in enumerate(merged_tally.filters):
-            for filter2 in tally.filters:
-                if filter1 != filter2 and filter1.can_merge(filter2):
-                    merged_filter = filter1.merge(filter2)
-                    merged_tally.filters[i] = merged_filter
-                    break
+        # If the two tallies are equal, simply return copy
+        if self == other:
+            return merged_tally
 
-        # Add unique scores from second tally to merged tally
-        for score in tally.scores:
-            if score not in merged_tally.scores:
-                merged_tally.add_score(score)
+        # Create deep copy of other tally to use for array concatenation
+        other_copy = copy.deepcopy(other)
 
-        # Add triggers from second tally to merged tally
-        for trigger in tally.triggers:
+        # Identify if filters, nuclides and scores are mergeable and/or equal
+        merge_filters = self._can_merge_filters(other)
+        merge_nuclides = self._can_merge_nuclides(other)
+        merge_scores = self._can_merge_scores(other)
+        equal_filters = sorted(self.filters) == sorted(other.filters)
+        equal_nuclides = sorted(self.nuclides) == sorted(other.nuclides)
+        equal_scores = sorted(self.scores) == sorted(other.scores)
+
+        # If two tallies can be merged along a filter's bins
+        if merge_filters and not equal_filters:
+
+            # Search for mergeable filters
+            for i, filter1 in enumerate(self.filters):
+                for j, filter2 in enumerate(other.filters):
+                    if filter1 != filter2 and filter1.can_merge(filter2):
+                        other_copy._swap_filters(other_copy.filters[i], filter2)
+                        merged_tally.filters[i] = filter1.merge(filter2)
+                        join_right = filter1 < filter2
+                        merge_axis = i
+                        break
+
+        # If two tallies can be merged along nuclide bins
+        if merge_nuclides and not equal_nuclides:
+            merge_axis = self.num_filters
+            join_right = True
+
+            # Add unique nuclides from other tally to merged tally
+            for nuclide in other.nuclides:
+                if nuclide not in merged_tally.nuclides:
+                    merged_tally.add_nuclide(nuclide)
+
+        # If two tallies can be merged along score bins
+        if merge_scores and not equal_scores:
+            merge_axis = self.num_filters + 1
+            join_right = True
+
+            # Add unique scores from other tally to merged tally
+            for score in other.scores:
+                if score not in merged_tally.scores:
+                    merged_tally.add_score(score)
+
+        # Add triggers from other tally to merged tally
+        for trigger in other.triggers:
             merged_tally.add_trigger(trigger)
 
+        # If results have not been read, then return tally for input generation
+        if self._results_read is None:
+            return merged_tally
+        # Otherwise, this is a derived tally which needs merged results arrays
+        else:
+            self._derived = True
+
+        # Update filter strides in merged tally
+        merged_tally._update_filter_strides()
+
+        # Concatenate sum arrays if present in both tallies
+        if self.sum is not None and other_copy.sum is not None:
+            self_sum = self.get_reshaped_data(value='sum')
+            other_sum = other_copy.get_reshaped_data(value='sum')
+
+            if join_right:
+                merged_sum = \
+                    np.concatenate((self_sum, other_sum), axis=merge_axis)
+            else:
+                merged_sum = \
+                    np.concatenate((other_sum, self_sum), axis=merge_axis)
+
+            merged_tally._sum = np.reshape(merged_sum, merged_tally.shape)
+
+        # Concatenate sum_sq arrays if present in both tallies
+        if self.sum_sq is not None and other.sum_sq is not None:
+            self_sum_sq = self.get_reshaped_data(value='sum_sq')
+            other_sum_sq = other_copy.get_reshaped_data(value='sum_sq')
+
+            if join_right:
+                merged_sum_sq = \
+                    np.concatenate((self_sum_sq, other_sum_sq), axis=merge_axis)
+            else:
+                merged_sum_sq = \
+                    np.concatenate((other_sum_sq, self_sum_sq), axis=merge_axis)
+
+            merged_tally._sum_sq = np.reshape(merged_sum_sq, merged_tally.shape)
+
+        # Concatenate mean arrays if present in both tallies
+        if self.mean is not None and other.mean is not None:
+            self_mean = self.get_reshaped_data(value='mean')
+            other_mean = other_copy.get_reshaped_data(value='mean')
+
+            if join_right:
+                merged_mean = \
+                    np.concatenate((self_mean, other_mean), axis=merge_axis)
+            else:
+                merged_mean = \
+                    np.concatenate((other_mean, self_mean), axis=merge_axis)
+
+            merged_tally._mean = np.reshape(merged_mean, merged_tally.shape)
+
+        # Concatenate std. dev. arrays if present in both tallies
+        if self.std_dev is not None and other.std_dev is not None:
+            self_std_dev = self.get_reshaped_data(value='std_dev')
+            other_std_dev = other_copy.get_reshaped_data(value='std_dev')
+
+            if join_right:
+                merged_std_dev = \
+                    np.concatenate((self_std_dev, other_std_dev), axis=merge_axis)
+            else:
+                merged_std_dev = \
+                    np.concatenate((other_std_dev, self_std_dev), axis=merge_axis)
+
+            merged_tally._std_dev = np.reshape(merged_std_dev, merged_tally.shape)
+
+        # Sparsify merged tally if both tallies are sparse
+        merged_tally.sparse = self.sparse and other.sparse
+
         return merged_tally
 
     def get_tally_xml(self):
@@ -847,6 +1072,32 @@ class Tally(object):
 
         return element
 
+    def contains_filter(self, filter_type):
+        """Looks for a filter in the tally that matches a specified type
+
+        Parameters
+        ----------
+        filter_type : str
+            Type of the filter, e.g. 'mesh'
+
+        Returns
+        -------
+        filter_found : bool
+            True if the tally contains a filter of the requested type;
+            otherwise false
+
+        """
+
+        filter_found = False
+
+        # Look through all of this Tally's Filters for the type requested
+        for test_filter in self.filters:
+            if test_filter.type == filter_type:
+                filter_found = True
+                break
+
+        return filter_found
+
     def find_filter(self, filter_type):
         """Return a filter in the tally that matches a specified type
 
@@ -1302,14 +1553,8 @@ class Tally(object):
                   'Summary info'.format(self.id)
             raise KeyError(msg)
 
-        # Attempt to import Pandas
-        try:
-            import pandas as pd
-        except ImportError:
-            msg = 'The Pandas Python package must be installed on your system'
-            raise ImportError(msg)
-
         # Initialize a pandas dataframe for the tally data
+        import pandas as pd
         df = pd.DataFrame()
 
         # Find the total length of the tally data array
@@ -1326,23 +1571,36 @@ class Tally(object):
         # Include DataFrame column for nuclides if user requested it
         if nuclides:
             nuclides = []
+            column_name = 'nuclide'
 
             for nuclide in self.nuclides:
-                # Write Nuclide name if Summary info was linked with StatePoint
                 if isinstance(nuclide, Nuclide):
                     nuclides.append(nuclide.name)
+                elif isinstance(nuclide, AggregateNuclide):
+                    nuclides.append(nuclide.name)
+                    column_name = '{0}(nuclide)'.format(nuclide.aggregate_op)
                 else:
                     nuclides.append(nuclide)
 
             # Tile the nuclide bins into a DataFrame column
             nuclides = np.repeat(nuclides, len(self.scores))
             tile_factor = data_size / len(nuclides)
-            df['nuclide'] = np.tile(nuclides, int(tile_factor))
+            df[column_name] = np.tile(nuclides, int(tile_factor))
 
         # Include column for scores if user requested it
         if scores:
+            scores = []
+            column_name = 'score'
+
+            for score in self.scores:
+                if isinstance(score, (basestring, CrossScore)):
+                    scores.append(score)
+                elif isinstance(score, AggregateScore):
+                    scores.append(score.name)
+                    column_name = '{0}(score)'.format(score.aggregate_op)
+
             tile_factor = data_size / len(self.scores)
-            df['score'] = np.tile(self.scores, int(tile_factor))
+            df[column_name] = np.tile(scores, int(tile_factor))
 
         # Append columns with mean, std. dev. for each tally bin
         df['mean'] = self.mean.ravel()
@@ -1951,19 +2209,12 @@ class Tally(object):
 
         """
 
-        # Check that results have been read
-        if not self.derived and self.sum is None:
-            msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \
-                  'since it does not contain any results.'.format(self.id)
-            raise ValueError(msg)
-
-        cv.check_type('filter1', filter1, Filter)
-        cv.check_type('filter2', filter2, Filter)
+        cv.check_type('filter1', filter1, (Filter, CrossFilter, AggregateFilter))
+        cv.check_type('filter2', filter2, (Filter, CrossFilter, AggregateFilter))
 
         # Check that the filters exist in the tally and are not the same
         if filter1 == filter2:
-            msg = 'Unable to swap a filter with itself'
-            raise ValueError(msg)
+            return
         elif filter1 not in self.filters:
             msg = 'Unable to swap "{0}" filter1 in Tally ID="{1}" since it ' \
                   'does not contain such a filter'.format(filter1.type, self.id)
@@ -2684,7 +2935,13 @@ class Tally(object):
                   'since it does not contain any results.'.format(self.id)
             raise ValueError(msg)
 
+        # Create deep copy of tally to return as sliced tally
         new_tally = copy.deepcopy(self)
+        new_tally._derived = True
+
+        # Differentiate Tally with a new auto-generated Tally ID
+        new_tally.id = None
+
         new_tally.sparse = False
 
         if not self.derived and self.sum is not None:
@@ -2769,7 +3026,7 @@ class Tally(object):
     def summation(self, scores=[], filter_type=None,
                   filter_bins=[], nuclides=[], remove_filter=False):
         """Vectorized sum of tally data across scores, filter bins and/or
-        nuclides using tally addition.
+        nuclides using tally aggregation.
 
         This method constructs a new tally to encapsulate the sum of the data
         represented by the summation of the data in this tally. The tally data
@@ -2811,7 +3068,7 @@ class Tally(object):
         tally_sum._derived = True
         tally_sum._estimator = self.estimator
         tally_sum._num_realizations = self.num_realizations
-        tally_sum.with_batch_statistics = self.with_batch_statistics
+        tally_sum._with_batch_statistics = self.with_batch_statistics
         tally_sum._with_summary = self.with_summary
         tally_sum._sp_filename = self._sp_filename
         tally_sum._results_read = self._results_read
@@ -2852,7 +3109,7 @@ class Tally(object):
                     # Add AggregateFilter to the tally sum
                     if not remove_filter:
                         filter_sum = \
-                            AggregateFilter(self_filter, filter_bins, 'sum')
+                            AggregateFilter(self_filter, [tuple(filter_bins)], 'sum')
                         tally_sum.add_filter(filter_sum)
 
                 # Add a copy of each filter not summed across to the tally sum
@@ -2914,6 +3171,157 @@ class Tally(object):
         tally_sum.sparse = self.sparse
         return tally_sum
 
+    def average(self, scores=[], filter_type=None,
+                  filter_bins=[], nuclides=[], remove_filter=False):
+        """Vectorized average of tally data across scores, filter bins and/or
+        nuclides using tally aggregation.
+
+        This method constructs a new tally to encapsulate the average of the
+        data represented by the average of the data in this tally. The tally
+        data average is determined by the scores, filter bins and nuclides
+        specified in the input parameters.
+
+        Parameters
+        ----------
+        scores : list of str
+            A list of one or more score strings to average across
+            (e.g., ['absorption', 'nu-fission']; default is [])
+        filter_type : str
+            A filter type string (e.g., 'cell', 'energy') corresponding to the
+            filter bins to average across
+        filter_bins : Iterable of Integral or tuple
+            A list of the filter bins corresponding to the filter_type parameter
+            Each bin in the list is the integer ID for 'material', 'surface',
+            'cell', 'cellborn', and 'universe' Filters. Each bin is an integer
+            for the cell instance ID for 'distribcell' Filters. Each bin is a
+            2-tuple of floats for 'energy' and 'energyout' filters corresponding
+            to the energy boundaries of the bin of interest. Each bin is an
+            (x,y,z) 3-tuple for 'mesh' filters corresponding to the mesh cell of
+            interest.
+        nuclides : list of str
+            A list of nuclide name strings to average across
+            (e.g., ['U-235', 'U-238']; default is [])
+        remove_filter : bool
+            If a filter is being averaged over, this bool indicates whether to
+            remove that filter in the returned tally. Default is False.
+
+        Returns
+        -------
+        Tally
+            A new tally which encapsulates the average of data requested.
+        """
+
+        # Create new derived Tally for average
+        tally_avg = Tally()
+        tally_avg._derived = True
+        tally_avg._estimator = self.estimator
+        tally_avg._num_realizations = self.num_realizations
+        tally_avg._with_batch_statistics = self.with_batch_statistics
+        tally_avg._with_summary = self.with_summary
+        tally_avg._sp_filename = self._sp_filename
+        tally_avg._results_read = self._results_read
+
+        # Get tally data arrays reshaped with one dimension per filter
+        mean = self.get_reshaped_data(value='mean')
+        std_dev = self.get_reshaped_data(value='std_dev')
+
+        # Average across any filter bins specified by the user
+        if filter_type in _FILTER_TYPES:
+            find_filter = self.find_filter(filter_type)
+
+            # If user did not specify filter bins, average across all bins
+            if len(filter_bins) == 0:
+                bin_indices = np.arange(find_filter.num_bins)
+
+                if filter_type == 'distribcell':
+                    filter_bins = np.arange(find_filter.num_bins)
+                else:
+                    num_bins = find_filter.num_bins
+                    filter_bins = \
+                        [(find_filter.get_bin(i)) for i in range(num_bins)]
+
+            # Only average across bins specified by the user
+            else:
+                bin_indices = \
+                    [find_filter.get_bin_index(bin) for bin in filter_bins]
+
+            # Average across the bins in the user-specified filter
+            for i, self_filter in enumerate(self.filters):
+                if self_filter.type == filter_type:
+                    mean = np.take(mean, indices=bin_indices, axis=i)
+                    std_dev = np.take(std_dev, indices=bin_indices, axis=i)
+                    mean = np.mean(mean, axis=i, keepdims=True)
+                    std_dev = np.mean(std_dev**2, axis=i, keepdims=True)
+                    std_dev /= len(bin_indices)
+                    std_dev = np.sqrt(std_dev)
+
+                    # Add AggregateFilter to the tally avg
+                    if not remove_filter:
+                        filter_sum = \
+                            AggregateFilter(self_filter, [tuple(filter_bins)], 'avg')
+                        tally_avg.add_filter(filter_sum)
+
+                # Add a copy of each filter not averaged across to the tally avg
+                else:
+                    tally_avg.add_filter(copy.deepcopy(self_filter))
+
+        # Add a copy of this tally's filters to the tally avg
+        else:
+            tally_avg._filters = copy.deepcopy(self.filters)
+
+        # Sum across any nuclides specified by the user
+        if len(nuclides) != 0:
+            nuclide_bins = [self.get_nuclide_index(nuclide) for nuclide in nuclides]
+            axis_index = self.num_filters
+            mean = np.take(mean, indices=nuclide_bins, axis=axis_index)
+            std_dev = np.take(std_dev, indices=nuclide_bins, axis=axis_index)
+            mean = np.mean(mean, axis=axis_index, keepdims=True)
+            std_dev = np.mean(std_dev**2, axis=axis_index, keepdims=True)
+            std_dev /= len(nuclide_bins)
+            std_dev = np.sqrt(std_dev)
+
+            # Add AggregateNuclide to the tally avg
+            nuclide_avg = AggregateNuclide(nuclides, 'avg')
+            tally_avg.add_nuclide(nuclide_avg)
+
+        # Add a copy of this tally's nuclides to the tally avg
+        else:
+            tally_avg._nuclides = copy.deepcopy(self.nuclides)
+
+        # Sum across any scores specified by the user
+        if len(scores) != 0:
+            score_bins = [self.get_score_index(score) for score in scores]
+            axis_index = self.num_filters + 1
+            mean = np.take(mean, indices=score_bins, axis=axis_index)
+            std_dev = np.take(std_dev, indices=score_bins, axis=axis_index)
+            mean = np.sum(mean, axis=axis_index, keepdims=True)
+            std_dev = np.sum(std_dev**2, axis=axis_index, keepdims=True)
+            std_dev /= len(score_bins)
+            std_dev = np.sqrt(std_dev)
+
+            # Add AggregateScore to the tally avg
+            score_sum = AggregateScore(scores, 'avg')
+            tally_avg.add_score(score_sum)
+
+        # Add a copy of this tally's scores to the tally avg
+        else:
+            tally_avg._scores = copy.deepcopy(self.scores)
+
+        # Update the tally avg's filter strides
+        tally_avg._update_filter_strides()
+
+        # Reshape condensed data arrays with one dimension for all filters
+        mean = np.reshape(mean, tally_avg.shape)
+        std_dev = np.reshape(std_dev, tally_avg.shape)
+
+        # Assign tally avg's data with the new arrays
+        tally_avg._mean = mean
+        tally_avg._std_dev = std_dev
+
+        # If original tally was sparse, sparsify the tally average
+        tally_avg.sparse = self.sparse
+        return tally_avg
+
     def diagonalize_filter(self, new_filter):
         """Diagonalize the tally data array along a new axis of filter bins.
 

src/tracking.F90

@@ -202,7 +202,7 @@ contains
 
         ! Advance xs_seed N times ahead to avoid re-using prn
         if (p % E /= p % last_E) &
-          xs_seed = prn_skip_ahead(n_nuc_zaid_total, xs_seed)
+             xs_seed = prn_skip_ahead(n_nuc_zaid_total, xs_seed)
 
         ! Set all uvws to base level -- right now, after a collision, only the
         ! base level uvws are changed

tests/test_asymmetric_lattice/test_asymmetric_lattice.py

@@ -19,13 +19,10 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness):
         # Build full core geometry from underlying input set
         self._input_set.build_default_materials_and_geometry()
 
-        # Extract all universes from the full core geometry
-        geometry = self._input_set.geometry.geometry
-        all_univs = geometry.get_all_universes()
-
         # Extract universes encapsulating fuel and water assemblies
-        water = all_univs[7]
-        fuel = all_univs[8]
+        geometry = self._input_set.geometry.geometry
+        water = geometry.get_universes_by_name('water assembly (hot)')[0]
+        fuel = geometry.get_universes_by_name('fuel assembly (hot)')[0]
 
         # Construct a 3x3 lattice of fuel assemblies
         core_lat = openmc.RectLattice(name='3x3 Core Lattice', lattice_id=202)
@@ -102,9 +99,10 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness):
         outstr += ', '.join(map(str, tally.std_dev.flatten())) + '\n'
 
         # Extract fuel assembly lattices from the summary
-        all_cells = su.openmc_geometry.get_all_cells()
-        fuel = all_cells[80].fill
-        core = all_cells[1].fill
+        core = su.get_cell_by_id(1)
+        fuel = su.get_cell_by_id(80)
+        fuel = fuel.fill
+        core = core.fill
 
         # Append a string of lattice distribcell offsets to the string
         outstr += ', '.join(map(str, fuel.offsets.flatten())) + '\n'

tests/test_mgxs_library_distribcell/results_true.dat

@@ -1,5 +1,5 @@
-                                    sum(distribcell)  group in nuclide      mean  std. dev.
-0  sum(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ...         1   total  0.651951   1.469284                                    sum(distribcell)  group in nuclide  mean  std. dev.
-0  sum(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ...         1   total     0          0                                    sum(distribcell)  group in  group out nuclide     mean  std. dev.
-0  sum(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ...         1          1   total  0.53214   1.320678                                    sum(distribcell)  group out nuclide  mean  std. dev.
-0  sum(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ...          1   total     0          0
+                                    avg(distribcell)  group in nuclide      mean  std. dev.
+0  (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...         1   total  0.651951   1.469284                                    avg(distribcell)  group in nuclide  mean  std. dev.
+0  (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...         1   total     0          0                                    avg(distribcell)  group in  group out nuclide     mean  std. dev.
+0  (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...         1          1   total  0.53214   1.320678                                    avg(distribcell)  group out nuclide  mean  std. dev.
+0  (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...          1   total     0          0

tests/test_tally_slice_merge/inputs_true.dat

@@ -0,0 +1 @@
+8d1ab9e4add51b99045e990ac9c3dad9447e9720d811bc430d4bfdd7c2c035424bcb7750e4a4d0ec0460ea1ef4be46ac58372ed01d55f5d8cfeebbce75559066

tests/test_tally_slice_merge/results_true.dat

@@ -0,0 +1,49 @@
+   energy low [MeV]  energy high [MeV]  cell nuclide    score     mean  std. dev.
+0          0.00e+00           6.25e-07    21   U-235  fission 7.75e-02   6.68e-03   energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0          0.00e+00           6.25e-07    21   U-235  nu-fission 1.89e-01   1.63e-02   energy low [MeV]  energy high [MeV]  cell nuclide    score     mean  std. dev.
+0          0.00e+00           6.25e-07    21   U-238  fission 1.09e-07   9.57e-09   energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0          0.00e+00           6.25e-07    21   U-238  nu-fission 2.71e-07   2.39e-08   energy low [MeV]  energy high [MeV]  cell nuclide    score     mean  std. dev.
+0          6.25e-07           2.00e+01    21   U-235  fission 1.92e-02   1.23e-03   energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0          6.25e-07           2.00e+01    21   U-235  nu-fission 4.69e-02   2.99e-03   energy low [MeV]  energy high [MeV]  cell nuclide    score     mean  std. dev.
+0          6.25e-07           2.00e+01    21   U-238  fission 1.22e-02   1.16e-03   energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0          6.25e-07           2.00e+01    21   U-238  nu-fission 3.41e-02   3.39e-03   energy low [MeV]  energy high [MeV]  cell nuclide    score     mean  std. dev.
+0          0.00e+00           6.25e-07    27   U-235  fission 8.32e-02   1.82e-03   energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0          0.00e+00           6.25e-07    27   U-235  nu-fission 2.03e-01   4.44e-03   energy low [MeV]  energy high [MeV]  cell nuclide    score     mean  std. dev.
+0          0.00e+00           6.25e-07    27   U-238  fission 1.17e-07   2.95e-09   energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0          0.00e+00           6.25e-07    27   U-238  nu-fission 2.93e-07   7.36e-09   energy low [MeV]  energy high [MeV]  cell nuclide    score     mean  std. dev.
+0          6.25e-07           2.00e+01    27   U-235  fission 2.60e-02   1.70e-03   energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0          6.25e-07           2.00e+01    27   U-235  nu-fission 6.38e-02   4.14e-03   energy low [MeV]  energy high [MeV]  cell nuclide    score     mean  std. dev.
+0          6.25e-07           2.00e+01    27   U-238  fission 1.47e-02   6.22e-04   energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0          6.25e-07           2.00e+01    27   U-238  nu-fission 4.12e-02   1.97e-03    energy low [MeV]  energy high [MeV]  cell nuclide       score     mean  std. dev.
+0           0.00e+00           6.25e-07    21   U-235     fission 7.75e-02   6.68e-03
+1           0.00e+00           6.25e-07    21   U-235  nu-fission 1.89e-01   1.63e-02
+2           0.00e+00           6.25e-07    21   U-238     fission 1.09e-07   9.57e-09
+3           0.00e+00           6.25e-07    21   U-238  nu-fission 2.71e-07   2.39e-08
+4           0.00e+00           6.25e-07    27   U-235     fission 8.32e-02   1.82e-03
+5           0.00e+00           6.25e-07    27   U-235  nu-fission 2.03e-01   4.44e-03
+6           0.00e+00           6.25e-07    27   U-238     fission 1.17e-07   2.95e-09
+7           0.00e+00           6.25e-07    27   U-238  nu-fission 2.93e-07   7.36e-09
+8           6.25e-07           2.00e+01    21   U-235     fission 1.92e-02   1.23e-03
+9           6.25e-07           2.00e+01    21   U-235  nu-fission 4.69e-02   2.99e-03
+10          6.25e-07           2.00e+01    21   U-238     fission 1.22e-02   1.16e-03
+11          6.25e-07           2.00e+01    21   U-238  nu-fission 3.41e-02   3.39e-03
+12          6.25e-07           2.00e+01    27   U-235     fission 2.60e-02   1.70e-03
+13          6.25e-07           2.00e+01    27   U-235  nu-fission 6.38e-02   4.14e-03
+14          6.25e-07           2.00e+01    27   U-238     fission 1.47e-02   6.22e-04
+15          6.25e-07           2.00e+01    27   U-238  nu-fission 4.12e-02   1.97e-03         sum(distribcell)  energy low [MeV]  energy high [MeV] nuclide       score     mean  std. dev.
+0   (0, 100, 2000, 30000)          0.00e+00           6.25e-07   U-235     fission 0.00e+00   0.00e+00
+1   (0, 100, 2000, 30000)          0.00e+00           6.25e-07   U-235  nu-fission 0.00e+00   0.00e+00
+2   (0, 100, 2000, 30000)          0.00e+00           6.25e-07   U-238     fission 0.00e+00   0.00e+00
+3   (0, 100, 2000, 30000)          0.00e+00           6.25e-07   U-238  nu-fission 0.00e+00   0.00e+00
+4   (0, 100, 2000, 30000)          6.25e-07           2.00e+01   U-235     fission 0.00e+00   0.00e+00
+5   (0, 100, 2000, 30000)          6.25e-07           2.00e+01   U-235  nu-fission 0.00e+00   0.00e+00
+6   (0, 100, 2000, 30000)          6.25e-07           2.00e+01   U-238     fission 0.00e+00   0.00e+00
+7   (0, 100, 2000, 30000)          6.25e-07           2.00e+01   U-238  nu-fission 0.00e+00   0.00e+00
+8      (500, 5000, 50000)          0.00e+00           6.25e-07   U-235     fission 0.00e+00   0.00e+00
+9      (500, 5000, 50000)          0.00e+00           6.25e-07   U-235  nu-fission 0.00e+00   0.00e+00
+10     (500, 5000, 50000)          0.00e+00           6.25e-07   U-238     fission 0.00e+00   0.00e+00
+11     (500, 5000, 50000)          0.00e+00           6.25e-07   U-238  nu-fission 0.00e+00   0.00e+00
+12     (500, 5000, 50000)          6.25e-07           2.00e+01   U-235     fission 0.00e+00   0.00e+00
+13     (500, 5000, 50000)          6.25e-07           2.00e+01   U-235  nu-fission 0.00e+00   0.00e+00
+14     (500, 5000, 50000)          6.25e-07           2.00e+01   U-238     fission 0.00e+00   0.00e+00
+15     (500, 5000, 50000)          6.25e-07           2.00e+01   U-238  nu-fission 0.00e+00   0.00e+00

tests/test_tally_slice_merge/test_tally_slice_merge.py

@@ -0,0 +1,165 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import glob
+import hashlib
+import itertools
+sys.path.insert(0, os.pardir)
+from testing_harness import PyAPITestHarness
+import openmc
+
+
+class TallySliceMergeTestHarness(PyAPITestHarness):
+    def _build_inputs(self):
+
+        # The summary.h5 file needs to be created to read in the tallies
+        self._input_set.settings.output = {'summary': True}
+
+        # Initialize the tallies file
+        tallies_file = openmc.TalliesFile()
+
+        # Define nuclides and scores to add to both tallies
+        self.nuclides = ['U-235', 'U-238']
+        self.scores = ['fission', 'nu-fission']
+
+        # Define filters for energy and spatial domain
+
+        low_energy = openmc.Filter(type='energy', bins=[0., 0.625e-6])
+        high_energy = openmc.Filter(type='energy', bins=[0.625e-6, 20.])
+        merged_energies = low_energy.merge(high_energy)
+
+        cell_21 = openmc.Filter(type='cell', bins=[21])
+        cell_27 = openmc.Filter(type='cell', bins=[27])
+        distribcell_filter = openmc.Filter(type='distribcell', bins=[21])
+
+        self.cell_filters = [cell_21, cell_27]
+        self.energy_filters = [low_energy, high_energy]
+
+        # Initialize cell tallies with filters, nuclides and scores
+        tallies = []
+        for cell_filter in self.energy_filters:
+            for energy_filter in self.cell_filters:
+                for nuclide in self.nuclides:
+                    for score in self.scores:
+                        tally = openmc.Tally()
+                        tally.estimator = 'tracklength'
+                        tally.add_score(score)
+                        tally.add_nuclide(nuclide)
+                        tally.add_filter(cell_filter)
+                        tally.add_filter(energy_filter)
+                        tallies.append(tally)
+
+        # Merge all cell tallies together
+        while len(tallies) != 1:
+            halfway = int(len(tallies) / 2)
+            zip_split = zip(tallies[:halfway], tallies[halfway:])
+            tallies = list(map(lambda xy: xy[0].merge(xy[1]), zip_split))
+
+        # Specify a name for the tally
+        tallies[0].name = 'cell tally'
+
+        # Initialize a distribcell tally
+        distribcell_tally = openmc.Tally(name='distribcell tally')
+        distribcell_tally.estimator = 'tracklength'
+        distribcell_tally.add_filter(distribcell_filter)
+        distribcell_tally.add_filter(merged_energies)
+        for score in self.scores:
+            distribcell_tally.add_score(score)
+        for nuclide in self.nuclides:
+            distribcell_tally.add_nuclide(nuclide)
+
+        # Add tallies to a TalliesFile
+        tallies_file = openmc.TalliesFile()
+        tallies_file.add_tally(tallies[0])
+        tallies_file.add_tally(distribcell_tally)
+
+        # Export tallies to file
+        self._input_set.tallies = tallies_file
+        super(TallySliceMergeTestHarness, self)._build_inputs()
+
+    def _get_results(self, hash_output=False):
+        """Digest info in the statepoint and return as a string."""
+
+        # Read the statepoint file.
+        statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
+        sp = openmc.StatePoint(statepoint)
+
+        # Read the summary file.
+        summary = glob.glob(os.path.join(os.getcwd(), 'summary.h5'))[0]
+        su = openmc.Summary(summary)
+        sp.link_with_summary(su)
+
+        # Extract the cell tally
+        tallies = [sp.get_tally(name='cell tally')]
+
+        # Slice the tallies by cell filter bins
+        cell_filter_prod = itertools.product(tallies, self.cell_filters)
+        tallies = map(lambda tf: tf[0].get_slice(filters=[tf[1].type], 
+                      filter_bins=[tf[1].get_bin(0)]), cell_filter_prod)
+
+        # Slice the tallies by energy filter bins
+        energy_filter_prod = itertools.product(tallies, self.energy_filters)
+        tallies = map(lambda tf: tf[0].get_slice(filters=[tf[1].type],
+                      filter_bins=[(tf[1].get_bin(0),)]), energy_filter_prod)
+
+        # Slice the tallies by nuclide
+        nuclide_prod = itertools.product(tallies, self.nuclides)
+        tallies = map(lambda tn: tn[0].get_slice(nuclides=[tn[1]]), nuclide_prod)
+
+        # Slice the tallies by score
+        score_prod = itertools.product(tallies, self.scores)
+        tallies = map(lambda ts: ts[0].get_slice(scores=[ts[1]]), score_prod)
+        tallies = list(tallies)
+
+        # Initialize an output string
+        outstr = ''
+
+        # Append sliced Tally Pandas DataFrames to output string
+        for tally in tallies:
+            df = tally.get_pandas_dataframe()
+            outstr += df.to_string()
+
+        # Merge all tallies together
+        while len(tallies) != 1:
+            halfway = int(len(tallies) / 2)
+            zip_split = zip(tallies[:halfway], tallies[halfway:])
+            tallies = list(map(lambda xy: xy[0].merge(xy[1]), zip_split))
+
+        # Append merged Tally Pandas DataFrame to output string
+        df = tallies[0].get_pandas_dataframe()
+        outstr += df.to_string()
+
+        # Extract the distribcell tally
+        distribcell_tally = sp.get_tally(name='distribcell tally')
+
+        # Sum up a few subdomains from the distribcell tally 
+        sum1 = distribcell_tally.summation(filter_type='distribcell', 
+                                           filter_bins=[0,100,2000,30000])
+        # Sum up a few subdomains from the distribcell tally
+        sum2 = distribcell_tally.summation(filter_type='distribcell', 
+                                           filter_bins=[500,5000,50000])
+
+        # Merge the distribcell tally slices
+        merge_tally = sum1.merge(sum2)
+
+        # Append merged Tally Pandas DataFrame to output string
+        df = merge_tally.get_pandas_dataframe()
+        outstr += df.to_string()
+
+        # Hash the results if necessary
+        if hash_output:
+            sha512 = hashlib.sha512()
+            sha512.update(outstr.encode('utf-8'))
+            outstr = sha512.hexdigest()
+
+        return outstr
+
+    def _cleanup(self):
+        super(TallySliceMergeTestHarness, self)._cleanup()
+        f = os.path.join(os.getcwd(), 'tallies.xml')
+        if os.path.exists(f): os.remove(f)
+
+if __name__ == '__main__':
+    harness = TallySliceMergeTestHarness('statepoint.10.h5', True)
+    harness.main()