From 68b131522047c6023260be4043b9367bfdac2f17 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sat, 12 Dec 2015 23:44:15 -0500 Subject: [PATCH 01/12] Initial implementation of tally sparsification in Python API --- openmc/tallies.py | 65 +++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 60 insertions(+), 5 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 13a219dedd..2c45d981d4 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -10,6 +10,7 @@ from xml.etree import ElementTree as ET import sys import numpy as np +import scipy.sparse as sps from openmc import Mesh, Filter, Trigger, Nuclide from openmc.cross import CrossScore, CrossNuclide, CrossFilter @@ -104,6 +105,7 @@ class Tally(object): self._std_dev = None self._with_batch_statistics = False self._derived = False + self._sparse = False self._sp_filename = None self._results_read = False @@ -126,6 +128,7 @@ class Tally(object): clone._with_summary = self.with_summary clone._with_batch_statistics = self.with_batch_statistics clone._derived = self.derived + clone._sparse = self.sparse clone._sp_filename = self._sp_filename clone._results_read = self._results_read @@ -315,13 +318,21 @@ class Tally(object): self._sum = sum self._sum_sq = sum_sq + # Convert NumPy arrays to SciPy sparse matrices + if self.sparse: + self._sum = sps.lil_matrix(self._sum) + self._sum_sq = sps.lil_matrix(self._sum_sq) + # Indicate that Tally results have been read self._results_read = True # Close the HDF5 statepoint file f.close() - return self._sum + if self.sparse: + return self._sum.toarray() + else: + return self._sum @property def sum_sq(self): @@ -332,7 +343,10 @@ class Tally(object): # Force reading of sum and sum_sq self.sum - return self._sum_sq + if self.sparse: + return self._sum_sq.toarray() + else: + return self._sum_sq @property def mean(self): @@ -341,7 +355,15 @@ class Tally(object): return None self._mean = self.sum / self.num_realizations - return self._mean + + # Convert NumPy arrays to SciPy sparse matrices + if self.sparse: + self._mean = sps.lil_matrix(self._mean) + + if self.sparse: + return self._mean.toarray() + else: + return self._mean @property def std_dev(self): @@ -354,8 +376,17 @@ class Tally(object): self._std_dev = np.zeros_like(self.mean) self._std_dev[nonzero] = np.sqrt((self.sum_sq[nonzero]/n - self.mean[nonzero]**2)/(n - 1)) + + # Convert NumPy arrays to SciPy sparse matrices + if self.sparse: + self._std_dev = sps.lil_matrix(self._std_dev) + self.with_batch_statistics = True - return self._std_dev + + if self.sparse: + return self._std_dev.toarray() + else: + return self._std_dev @property def with_batch_statistics(self): @@ -365,6 +396,10 @@ class Tally(object): def derived(self): return self._derived + @property + def sparse(self): + return self._sparse + @estimator.setter def estimator(self, estimator): cv.check_value('estimator', estimator, @@ -619,7 +654,8 @@ class Tally(object): """ if not self.can_merge(tally): - msg = 'Unable to merge tally ID="{0}" with "{1}"'.format(tally.id, self.id) + msg = 'Unable to merge tally ID="{0}" with ' + \ + '"{1}"'.format(tally.id, self.id) raise ValueError(msg) # Create deep copy of tally to return as merged tally @@ -1107,6 +1143,25 @@ class Tally(object): return data + def sparsify(self): + """ + """ + + self._sparse = True + + # FIXME: get_values + # FIXME: each of the properties which load in the data + # FIXME: pandas dataframes + # summation, slicing + if self._sum: + self._sum = sps.lil_matrix(self._sum) + if self._sum_sq: + self._sum_sq = sps.lil_matrix(self._sum_sq) + if self._mean: + self._mean = sps.lil_matrix(self._mean) + if self._std_dev: + self._std_dev = sps.lil_matrix(self._std_dev) + def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True, summary=None): """Build a Pandas DataFrame for the Tally data. From 963238194b7b4a264c357e147709810a0e4dbea8 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sun, 13 Dec 2015 01:30:10 -0500 Subject: [PATCH 02/12] Fixed some bugs in tally sparsification - tally arithmetic ipython notebook now working --- openmc/tallies.py | 197 ++++++++++++++++++++++++++++++---------------- 1 file changed, 127 insertions(+), 70 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 67e9a08388..9cac3a8482 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -83,6 +83,11 @@ class Tally(object): An array containing the sample mean for each bin std_dev : ndarray An array containing the sample standard deviation for each bin + derived : bool + Whether or not the tally is derived from one or more other tallies + sparse : bool + Whether or not the tally uses SciPy's LIL sparse matrix format for + compressed data storage """ @@ -246,6 +251,10 @@ class Tally(object): def scores(self): return self._scores + @property + def shape(self): + return (self.num_filter_bins, self.num_nuclides, self.num_score_bins) + @property def num_scores(self): return len(self._scores) @@ -308,21 +317,18 @@ class Tally(object): return 1 if not val else val # Reshape the results arrays - new_shape = (nonzero(self.num_filter_bins), - nonzero(self.num_nuclides), - nonzero(self.num_score_bins)) - - sum = np.reshape(sum, new_shape) - sum_sq = np.reshape(sum_sq, new_shape) + sum = np.reshape(sum, self.shape) + sum_sq = np.reshape(sum_sq, self.shape) # Set the data for this Tally self._sum = sum self._sum_sq = sum_sq - # Convert NumPy arrays to SciPy sparse matrices + # Convert NumPy arrays to SciPy sparse LIL matrices if self.sparse: - self._sum = sps.lil_matrix(self._sum) - self._sum_sq = sps.lil_matrix(self._sum_sq) + self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape) + self._sum_sq = \ + sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) # Indicate that Tally results have been read self._results_read = True @@ -331,7 +337,7 @@ class Tally(object): f.close() if self.sparse: - return self._sum.toarray() + return np.reshape(self._sum.toarray(), self.shape) else: return self._sum @@ -345,7 +351,7 @@ class Tally(object): self.sum if self.sparse: - return self._sum_sq.toarray() + return np.reshape(self._sum_sq.toarray(), self.shape) else: return self._sum_sq @@ -357,12 +363,13 @@ class Tally(object): self._mean = self.sum / self.num_realizations - # Convert NumPy arrays to SciPy sparse matrices + # Convert NumPy array to SciPy sparse LIL matrix if self.sparse: - self._mean = sps.lil_matrix(self._mean) + self._mean = \ + sps.lil_matrix(self._mean.flatten(), self._mean.shape) if self.sparse: - return self._mean.toarray() + return np.reshape(self._mean.toarray(), self.shape) else: return self._mean @@ -378,14 +385,15 @@ class Tally(object): self._std_dev[nonzero] = np.sqrt((self.sum_sq[nonzero]/n - self.mean[nonzero]**2)/(n - 1)) - # Convert NumPy arrays to SciPy sparse matrices + # Convert NumPy array to SciPy sparse LIL matrix if self.sparse: - self._std_dev = sps.lil_matrix(self._std_dev) + self._std_dev = \ + sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape) self.with_batch_statistics = True if self.sparse: - return self._std_dev.toarray() + return np.reshape(self._std_dev.toarray(), self.shape) else: return self._std_dev @@ -1099,22 +1107,19 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. ValueError is also thrown - if the input parameters do not correspond to the Tally's attributes, + When this method is called before the Tally is populated with data, + or the input parameters do not correspond to the Tally's attributes, e.g., if the score(s) do not match those in the Tally. """ - # Ensure that StatePoint.read_results() was called first + # Ensure that the tally has data if (value == 'mean' and self.mean is None) or \ (value == 'std_dev' and self.std_dev is None) or \ (value == 'rel_err' and self.mean is None) or \ (value == 'sum' and self.sum is None) or \ (value == 'sum_sq' and self.sum_sq is None): - msg = 'The Tally ID="{0}" has no data to return. Call the ' \ - 'StatePoint.read_results() method before using ' \ - 'Tally.get_values(...)'.format(self.id) + msg = 'The Tally ID="{0}" has no data to return'.format(self.id) raise ValueError(msg) # Get filter, nuclide and score indices @@ -1148,20 +1153,41 @@ class Tally(object): """ """ + # FIXME: pandas dataframes, summation + if self._sum is not None: + self._sum = \ + sps.lil_matrix(self._sum.flatten(), self._sum.shape) + if self._sum_sq is not None: + self._sum_sq = \ + sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) + if self._mean is not None: + self._mean = \ + sps.lil_matrix(self._mean.flatten(), self._mean.shape) + if self._std_dev is not None: + self._std_dev = \ + sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape) + self._sparse = True - # FIXME: get_values - # FIXME: each of the properties which load in the data - # FIXME: pandas dataframes - # summation, slicing - if self._sum: - self._sum = sps.lil_matrix(self._sum) - if self._sum_sq: - self._sum_sq = sps.lil_matrix(self._sum_sq) - if self._mean: - self._mean = sps.lil_matrix(self._mean) - if self._std_dev: - self._std_dev = sps.lil_matrix(self._std_dev) + def densify(self): + """ + """ + + # If the tally is already dense, simply return + if not self._sparse: + return + + # FIXME: pandas dataframes, summation + if self._sum is not None: + self._sum = np.reshape(self._sum.toarray(), self.shape) + if self._sum_sq is not None: + self._sum_sq = np.reshape(self._sum_sq.toarray(), self.shape) + if self._mean is not None: + self._mean = np.reshape(self._mean.toarray(), self.shape) + if self._std_dev is not None: + self._std_dev = np.reshape(self._std_dev.toarray(), self.shape) + + self._sparse = False def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True, summary=None): @@ -1200,17 +1226,14 @@ class Tally(object): ------ KeyError When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. ImportError When Pandas can not be found on the caller's system """ - # Ensure that StatePoint.read_results() was called first + # Ensure that the tally has data if self.mean is None or self.std_dev is None: - msg = 'The Tally ID="{0}" has no data to return. Call the ' \ - 'StatePoint.read_results() method before using ' \ - 'Tally.get_pandas_dataframe(...)'.format(self.id) + msg = 'The Tally ID="{0}" has no data to return'.format(self.id) raise KeyError(msg) # If using Summary, ensure StatePoint.link_with_summary(...) was called @@ -1356,16 +1379,13 @@ class Tally(object): Raises ------ KeyError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ - # Ensure that StatePoint.read_results() was called first + # Ensure that the tally has data if self._sum is None or self._sum_sq is None and not self.derived: - msg = 'The Tally ID="{0}" has no data to export. Call the ' \ - 'StatePoint.read_results() method before using ' \ - 'Tally.export_results(...)'.format(self.id) + msg = 'The Tally ID="{0}" has no data to export'.format(self.id) raise KeyError(msg) if not isinstance(filename, basestring): @@ -1527,7 +1547,7 @@ class Tally(object): ------ ValueError When this method is called before the other tally is populated - with data by the StatePoint.read_results() method. + with data. """ @@ -1708,6 +1728,12 @@ class Tally(object): """ + # Use dense NumPy arrays for data alignment operations + self_sparse = self.sparse + other_sparse = other.sparse + self.densify() + other.densify() + # Get the set of filters that each tally is missing other_missing_filters = set(self.filters).difference(set(other.filters)) self_missing_filters = set(other.filters).difference(set(self.filters)) @@ -1769,7 +1795,7 @@ class Tally(object): # If necessary, swap other nuclide if other_index != i: - other.swap_nuclides(nuclide, other.nuclides[i]) + other.swap_nuclides(nuclide, other.nuclides[i]) # Repeat and tile the data by score in preparation for performing # the tensor product across scores. @@ -1820,6 +1846,12 @@ class Tally(object): filter.stride = stride stride *= filter.num_bins + # Restore tally operands to sparse storage + if self.sparse: + self.sparsify() + if other.sparse: + other.sparsify() + # Deep copy the mean and std dev data self_mean = copy.deepcopy(self.mean) self_std_dev = copy.deepcopy(self.std_dev) @@ -1864,7 +1896,7 @@ class Tally(object): ------ ValueError If this is a derived tally or this method is called before the tally - is populated with data by the StatePoint.read_results() method. + is populated with data. """ @@ -1981,7 +2013,7 @@ class Tally(object): ------ ValueError If this is a derived tally or this method is called before the tally - is populated with data by the StatePoint.read_results() method. + is populated with data. """ @@ -2036,7 +2068,7 @@ class Tally(object): ------ ValueError If this is a derived tally or this method is called before the tally - is populated with data by the StatePoint.read_results() method. + is populated with data. """ @@ -2103,8 +2135,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -2136,6 +2167,10 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If original tally operands were sparse, sparsify the sliced tally + if self.sparse and other.sparse: + new_tally.sparsify() + else: msg = 'Unable to add "{0}" to Tally ID="{1}"'.format(other, self.id) raise ValueError(msg) @@ -2173,8 +2208,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -2205,6 +2239,10 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If original tally operands were sparse, sparsify the sliced tally + if self.sparse and other.sparse: + new_tally.sparsify() + else: msg = 'Unable to subtract "{0}" from Tally ' \ 'ID="{1}"'.format(other, self.id) @@ -2243,8 +2281,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -2275,6 +2312,10 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If original tally operands were sparse, sparsify the sliced tally + if self.sparse and other.sparse: + new_tally.sparsify() + else: msg = 'Unable to multiply Tally ID="{0}" ' \ 'by "{1}"'.format(self.id, other) @@ -2313,8 +2354,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -2345,6 +2385,10 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If original tally operands were sparse, sparsify the sliced tally + if self.sparse and other.sparse: + new_tally.sparsify() + else: msg = 'Unable to divide Tally ID="{0}" ' \ 'by "{1}"'.format(self.id, other) @@ -2386,8 +2430,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -2419,6 +2462,10 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If original tally was sparse, sparsify the exponentiated tally + if self.sparse: + new_tally.sparsify() + else: msg = 'Unable to raise Tally ID="{0}" to ' \ 'power "{1}"'.format(self.id, power) @@ -2569,12 +2616,11 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ - # Ensure that StatePoint.read_results() was called first + # Ensure that the tally has data 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) @@ -2657,6 +2703,10 @@ class Tally(object): filter.stride = stride stride *= filter.num_bins + # If original tally was sparse, sparsify the sliced tally + if self.sparse: + new_tally.sparsify() + return new_tally def summation(self, scores=[], filter_type=None, @@ -2759,6 +2809,10 @@ class Tally(object): filters[i] = CrossFilter(filters[i-1], filters[i], '+') tally_sum.add_filter(filters[-1]) + # If original tally was sparse, sparsify the tally summation + if self.sparse: + tally_sum.sparsify() + return tally_sum def diagonalize_filter(self, new_filter): @@ -2799,7 +2853,6 @@ class Tally(object): num_filter_bins = new_tally.num_filter_bins num_nuclides = new_tally.num_nuclides num_score_bins = new_tally.num_score_bins - new_shape = (num_filter_bins, num_nuclides, num_score_bins) # Determine "base" indices along the new "diagonal", and the factor # by which the "base" indices should be repeated to account for all @@ -2816,16 +2869,16 @@ class Tally(object): # Inject this Tally's data along the diagonal of the diagonalized Tally if self.sum is not None: - new_tally._sum = np.zeros(new_shape, dtype=np.float64) + new_tally._sum = np.zeros(self.shape, dtype=np.float64) new_tally._sum[diag_indices, :, :] = self.sum if self.sum_sq is not None: - new_tally._sum_sq = np.zeros(new_shape, dtype=np.float64) + new_tally._sum_sq = np.zeros(self.shape, dtype=np.float64) new_tally._sum_sq[diag_indices, :, :] = self.sum_sq if self.mean is not None: - new_tally._mean = np.zeros(new_shape, dtype=np.float64) + new_tally._mean = np.zeros(self.shape, dtype=np.float64) new_tally._mean[diag_indices, :, :] = self.mean if self.std_dev is not None: - new_tally._std_dev = np.zeros(new_shape, dtype=np.float64) + new_tally._std_dev = np.zeros(self.shape, dtype=np.float64) new_tally._std_dev[diag_indices, :, :] = self.std_dev # Correct each Filter's stride @@ -2834,6 +2887,10 @@ class Tally(object): filter.stride = stride stride *= filter.num_bins + # If original tally was sparse, sparsify the diagonalized tally + if self.sparse: + new_tally.sparsify + return new_tally From f2deefcf472cd607277208ee8c5f8f447179c70b Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sun, 13 Dec 2015 01:39:21 -0500 Subject: [PATCH 03/12] Fixed bugs in sparse tallies with filter diagonalization --- openmc/tallies.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 9cac3a8482..f9d3599bda 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -2869,16 +2869,16 @@ class Tally(object): # Inject this Tally's data along the diagonal of the diagonalized Tally if self.sum is not None: - new_tally._sum = np.zeros(self.shape, dtype=np.float64) + new_tally._sum = np.zeros(new_tally.shape, dtype=np.float64) new_tally._sum[diag_indices, :, :] = self.sum if self.sum_sq is not None: - new_tally._sum_sq = np.zeros(self.shape, dtype=np.float64) + new_tally._sum_sq = np.zeros(new_tally.shape, dtype=np.float64) new_tally._sum_sq[diag_indices, :, :] = self.sum_sq if self.mean is not None: - new_tally._mean = np.zeros(self.shape, dtype=np.float64) + new_tally._mean = np.zeros(new_tally.shape, dtype=np.float64) new_tally._mean[diag_indices, :, :] = self.mean if self.std_dev is not None: - new_tally._std_dev = np.zeros(self.shape, dtype=np.float64) + new_tally._std_dev = np.zeros(new_tally.shape, dtype=np.float64) new_tally._std_dev[diag_indices, :, :] = self.std_dev # Correct each Filter's stride From 61aed1df1095cd750964e7289a3ee23eb8e433ca Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sun, 13 Dec 2015 01:47:38 -0500 Subject: [PATCH 04/12] Added docstrings for Tally sparsify and densify methods --- openmc/tallies.py | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index f9d3599bda..25bf2acb6a 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1150,10 +1150,21 @@ class Tally(object): return data def sparsify(self): - """ + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices. + + This method may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within the Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + See also + -------- + Tally.densify() + """ - # FIXME: pandas dataframes, summation + # Convert NumPy arrays to SciPy sparse LIL matrices if self._sum is not None: self._sum = \ sps.lil_matrix(self._sum.flatten(), self._sum.shape) @@ -1170,14 +1181,24 @@ class Tally(object): self._sparse = True def densify(self): - """ + """Convert tally data from SciPy list of lists (LIL) sparse matrices + to NumPy arrrays. + + This method may be used to restore a sparse tally to its original + state. This method will have no effect on the state of the tally + unless the Tally.sparse() method has been called. + + See also + -------- + Tally.sparsify() + """ # If the tally is already dense, simply return if not self._sparse: return - # FIXME: pandas dataframes, summation + # Convert SciPy sparse LIL matrices to NumPy arrays if self._sum is not None: self._sum = np.reshape(self._sum.toarray(), self.shape) if self._sum_sq is not None: From e85b659977fd85005783ae9a6e9cd64297e5fd83 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sun, 13 Dec 2015 02:04:19 -0500 Subject: [PATCH 05/12] Added tally sparsification to MGXS class --- openmc/mgxs/mgxs.py | 78 +++++++++++++++++++++++++++++++++++++++++++++ openmc/tallies.py | 2 +- 2 files changed, 79 insertions(+), 1 deletion(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 597e28e292..3804199cc9 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -103,6 +103,9 @@ class MGXS(object): nuclides : list of str or 'sum' A list of nuclide string names (e.g., 'U-238', 'O-16') when by_nuclide is True and 'sum' when by_nuclide is False. + sparse : bool + Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format + for compressed data storage """ @@ -121,6 +124,7 @@ class MGXS(object): self._tally_trigger = None self._tallies = None self._xs_tally = None + self._sparse = False self.name = name self.by_nuclide = by_nuclide @@ -146,6 +150,7 @@ class MGXS(object): clone._energy_groups = copy.deepcopy(self.energy_groups, memo) clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo) clone._xs_tally = copy.deepcopy(self.xs_tally, memo) + clone._sparse = self.sparse clone._tallies = OrderedDict() for tally_type, tally in self.tallies.items(): @@ -199,6 +204,10 @@ class MGXS(object): def xs_tally(self): return self._xs_tally + @property + def sparse(self): + return self._sparse + @property def num_subdomains(self): tally = list(self.tallies.values())[0] @@ -504,6 +513,10 @@ class MGXS(object): self._xs_tally._mean = np.nan_to_num(self.xs_tally.mean) self._xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev) + # Sparsify the MGXS derived tally + if self.sparse: + self.xs_tally.sparsify() + def load_from_statepoint(self, statepoint): """Extracts tallies in an OpenMC StatePoint with the data needed to compute multi-group cross sections. @@ -567,6 +580,9 @@ class MGXS(object): filter_bins, tally.nuclides) self.tallies[tally_type] = sp_tally + if self.sparse: + sp_tally.sparsify() + # Compute the cross section from the tallies self.compute_xs() @@ -765,6 +781,11 @@ class MGXS(object): # Compute the energy condensed multi-group cross section condensed_xs.compute_xs() + + # Sparsify the condensed MGXS + if self.sparse: + condensed_xs.sparsify() + return condensed_xs def get_subdomain_avg_xs(self, subdomains='all'): @@ -848,6 +869,10 @@ class MGXS(object): # Compute the subdomain-averaged multi-group cross section avg_xs.compute_xs() + # Sparsify the subdomain-averaged MGXS + if self.sparse: + avg_xs.sparsify() + return avg_xs def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'): @@ -2078,6 +2103,59 @@ class Chi(MGXS): super(Chi, self).compute_xs() + def sparsify(self): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices. + + This method may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within the Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + See also + -------- + MGXS.densify(), Tally.sparsify(), Tally.densify() + + """ + + # Sparsify the derived MGXS tally + if self.xs_tally: + self.xs_tally.sparsify() + + # Sparsify the MGXS' base tallies + for tally_name in self.tallies: + self.tallies[tally_name].sparsify() + + self._sparse = True + + def densify(self): + """Convert tally data from SciPy list of lists (LIL) sparse matrices + to NumPy arrrays. + + This method may be used to restore a sparse tally to its original + state. This method will have no effect on the state of the tally + unless the Tally.sparse() method has been called. + + See also + -------- + MGXS.sparsify(), Tally.sparsify(), Tally.densify() + + """ + + # If the MGXS is already dense, simply return + if not self.sparse: + return + + # Densify the derived MGXS tally + if self.xs_tally: + self.xs_tally.densify() + + # Densify the MGXS' base tallies + for tally_name in self.tallies: + self.tallies[tally_name].densify() + + self._sparse = False + 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. diff --git a/openmc/tallies.py b/openmc/tallies.py index 25bf2acb6a..1ee1d88907 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1195,7 +1195,7 @@ class Tally(object): """ # If the tally is already dense, simply return - if not self._sparse: + if not self.sparse: return # Convert SciPy sparse LIL matrices to NumPy arrays From beb3189d8d988e66346de247cc667a1aa7c83915 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sun, 13 Dec 2015 11:14:01 -0500 Subject: [PATCH 06/12] Fixed bugs in tally sparsification which broke MGXS library --- openmc/mgxs/library.py | 32 +++++ openmc/mgxs/mgxs.py | 101 +++++----------- openmc/tallies.py | 266 +++++++++++++++++++++-------------------- 3 files changed, 197 insertions(+), 202 deletions(-) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index ef3e90fc41..ca9fa5d669 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -73,6 +73,9 @@ class Library(object): name : str, optional Name of the multi-group cross section library. Used as a label to identify tallies in OpenMC 'tallies.xml' file. + sparse : bool + Whether or not the Library's tallies use SciPy's LIL sparse matrix + format for compressed data storage """ @@ -92,6 +95,7 @@ class Library(object): self._all_mgxs = OrderedDict() self._sp_filename = None self._keff = None + self._sparse = False self.name = name self.openmc_geometry = openmc_geometry @@ -119,6 +123,7 @@ class Library(object): clone._all_mgxs = self.all_mgxs clone._sp_filename = self._sp_filename clone._keff = self._keff + clone._sparse = self.sparse clone._all_mgxs = OrderedDict() for domain in self.domains: @@ -208,6 +213,10 @@ class Library(object): def keff(self): return self._keff + @property + def sparse(self): + return self._sparse + @openmc_geometry.setter def openmc_geometry(self, openmc_geometry): cv.check_type('openmc_geometry', openmc_geometry, openmc.Geometry) @@ -286,6 +295,28 @@ class Library(object): cv.check_type('tally trigger', tally_trigger, openmc.Trigger) self._tally_trigger = tally_trigger + @sparse.setter + def sparse(self, sparse): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices, and vice versa. + + This property may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within the Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + """ + + cv.check_type('sparse', sparse, bool) + + # Sparsify or densify each MGXS in the Library + for domain in self.domains: + for mgxs_type in self.mgxs_types: + mgxs = self.get_mgxs(domain, mgxs_type) + mgxs.sparse = self.sparse + + self._sparse = sparse + def build_library(self): """Initialize MGXS objects in each domain and for each reaction type in the library. @@ -381,6 +412,7 @@ class Library(object): for mgxs_type in self.mgxs_types: mgxs = self.get_mgxs(domain, mgxs_type) mgxs.load_from_statepoint(statepoint) + mgxs.sparse = self.sparse def get_mgxs(self, domain, mgxs_type): """Return the MGXS object for some domain and reaction rate type. diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 3804199cc9..1b451affc9 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -258,6 +258,29 @@ class MGXS(object): cv.check_type('tally trigger', tally_trigger, openmc.Trigger) self._tally_trigger = tally_trigger + @sparse.setter + def sparse(self, sparse): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices, and vice versa. + + This property may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within the Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + """ + + cv.check_type('sparse', sparse, bool) + + # Sparsify or densify the derived MGXS tally and its base tallies + if self.xs_tally: + self.xs_tally.sparse = sparse + + for tally_name in self.tallies: + self.tallies[tally_name].sparse = sparse + + self._sparse = sparse + @staticmethod def get_mgxs(mgxs_type, domain=None, domain_type=None, energy_groups=None, by_nuclide=False, name=''): @@ -512,10 +535,7 @@ class MGXS(object): # Remove NaNs which may have resulted from divide-by-zero operations self._xs_tally._mean = np.nan_to_num(self.xs_tally.mean) self._xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev) - - # Sparsify the MGXS derived tally - if self.sparse: - self.xs_tally.sparsify() + self.xs_tally.sparse = self.sparse def load_from_statepoint(self, statepoint): """Extracts tallies in an OpenMC StatePoint with the data needed to @@ -578,11 +598,9 @@ class MGXS(object): estimator=tally.estimator) sp_tally = sp_tally.get_slice(tally.scores, filters, filter_bins, tally.nuclides) + sp_tally.sparse = self.sparse self.tallies[tally_type] = sp_tally - if self.sparse: - sp_tally.sparsify() - # Compute the cross section from the tallies self.compute_xs() @@ -732,6 +750,7 @@ class MGXS(object): # Clone this MGXS to initialize the condensed version condensed_xs = copy.deepcopy(self) + condensed_xs.sparse = False condensed_xs.energy_groups = coarse_groups # Build energy indices to sum across @@ -781,11 +800,7 @@ class MGXS(object): # Compute the energy condensed multi-group cross section condensed_xs.compute_xs() - - # Sparsify the condensed MGXS - if self.sparse: - condensed_xs.sparsify() - + condensed_xs.sparse = self.sparse return condensed_xs def get_subdomain_avg_xs(self, subdomains='all'): @@ -828,8 +843,9 @@ class MGXS(object): # Clone this MGXS to initialize the subdomain-averaged version avg_xs = copy.deepcopy(self) + avg_xs.sparse = False - # If domain is distribcell, make the new domain 'cell' + # If domain is distribcell, make the new domain 'cell' if self.domain_type == 'distribcell': avg_xs.domain_type = 'cell' @@ -868,11 +884,7 @@ class MGXS(object): # Compute the subdomain-averaged multi-group cross section avg_xs.compute_xs() - - # Sparsify the subdomain-averaged MGXS - if self.sparse: - avg_xs.sparsify() - + avg_xs.sparse = self.sparse return avg_xs def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'): @@ -2103,59 +2115,6 @@ class Chi(MGXS): super(Chi, self).compute_xs() - def sparsify(self): - """Convert tally data from NumPy arrays to SciPy list of lists (LIL) - sparse matrices. - - This method may be used to reduce the amount of data in memory during - tally data processing. The tally data will be stored as SciPy LIL - matrices internally within the Tally object. All tally data access - properties and methods will return data as a dense NumPy array. - - See also - -------- - MGXS.densify(), Tally.sparsify(), Tally.densify() - - """ - - # Sparsify the derived MGXS tally - if self.xs_tally: - self.xs_tally.sparsify() - - # Sparsify the MGXS' base tallies - for tally_name in self.tallies: - self.tallies[tally_name].sparsify() - - self._sparse = True - - def densify(self): - """Convert tally data from SciPy list of lists (LIL) sparse matrices - to NumPy arrrays. - - This method may be used to restore a sparse tally to its original - state. This method will have no effect on the state of the tally - unless the Tally.sparse() method has been called. - - See also - -------- - MGXS.sparsify(), Tally.sparsify(), Tally.densify() - - """ - - # If the MGXS is already dense, simply return - if not self.sparse: - return - - # Densify the derived MGXS tally - if self.xs_tally: - self.xs_tally.densify() - - # Densify the MGXS' base tallies - for tally_name in self.tallies: - self.tallies[tally_name].densify() - - self._sparse = False - 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. diff --git a/openmc/tallies.py b/openmc/tallies.py index 1ee1d88907..b57a0d228f 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -326,7 +326,8 @@ class Tally(object): # Convert NumPy arrays to SciPy sparse LIL matrices if self.sparse: - self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape) + self._sum = \ + sps.lil_matrix(self._sum.flatten(), self._sum.shape) self._sum_sq = \ sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) @@ -536,6 +537,48 @@ class Tally(object): cv.check_type('sum_sq', sum_sq, Iterable) self._sum_sq = sum_sq + @sparse.setter + def sparse(self, sparse): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices, and vice versa. + + This property may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within the Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + """ + + cv.check_type('sparse', sparse, bool) + + # Convert NumPy arrays to SciPy sparse LIL matrices + if sparse and not self.sparse: + if self._sum is not None: + self._sum = \ + sps.lil_matrix(self._sum.flatten(), self._sum.shape) + if self._sum_sq is not None: + self._sum_sq = \ + sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) + if self._mean is not None: + self._mean = \ + sps.lil_matrix(self._mean.flatten(), self._mean.shape) + if self._std_dev is not None: + self._std_dev = \ + sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape) + self._sparse = True + + # Convert SciPy sparse LIL matrices to NumPy arrays + elif not sparse and self.sparse: + if self._sum is not None: + self._sum = np.reshape(self._sum.toarray(), self.shape) + if self._sum_sq is not None: + self._sum_sq = np.reshape(self._sum_sq.toarray(), self.shape) + if self._mean is not None: + self._mean = np.reshape(self._mean.toarray(), self.shape) + if self._std_dev is not None: + self._std_dev = np.reshape(self._std_dev.toarray(), self.shape) + self._sparse = False + def remove_score(self, score): """Remove a score from the tally @@ -1149,67 +1192,6 @@ class Tally(object): return data - def sparsify(self): - """Convert tally data from NumPy arrays to SciPy list of lists (LIL) - sparse matrices. - - This method may be used to reduce the amount of data in memory during - tally data processing. The tally data will be stored as SciPy LIL - matrices internally within the Tally object. All tally data access - properties and methods will return data as a dense NumPy array. - - See also - -------- - Tally.densify() - - """ - - # Convert NumPy arrays to SciPy sparse LIL matrices - if self._sum is not None: - self._sum = \ - sps.lil_matrix(self._sum.flatten(), self._sum.shape) - if self._sum_sq is not None: - self._sum_sq = \ - sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) - if self._mean is not None: - self._mean = \ - sps.lil_matrix(self._mean.flatten(), self._mean.shape) - if self._std_dev is not None: - self._std_dev = \ - sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape) - - self._sparse = True - - def densify(self): - """Convert tally data from SciPy list of lists (LIL) sparse matrices - to NumPy arrrays. - - This method may be used to restore a sparse tally to its original - state. This method will have no effect on the state of the tally - unless the Tally.sparse() method has been called. - - See also - -------- - Tally.sparsify() - - """ - - # If the tally is already dense, simply return - if not self.sparse: - return - - # Convert SciPy sparse LIL matrices to NumPy arrays - if self._sum is not None: - self._sum = np.reshape(self._sum.toarray(), self.shape) - if self._sum_sq is not None: - self._sum_sq = np.reshape(self._sum_sq.toarray(), self.shape) - if self._mean is not None: - self._mean = np.reshape(self._mean.toarray(), self.shape) - if self._std_dev is not None: - self._std_dev = np.reshape(self._std_dev.toarray(), self.shape) - - self._sparse = False - def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True, summary=None): """Build a Pandas DataFrame for the Tally data. @@ -1626,6 +1608,9 @@ class Tally(object): self_copy = copy.deepcopy(self) other_copy = copy.deepcopy(other) + self_copy.sparse = False + other_copy.sparse = False + # Align the tally data based on desired hybrid product data = self_copy._align_tally_data(other_copy, filter_product, nuclide_product, score_product) @@ -1691,7 +1676,8 @@ class Tally(object): else: all_nuclides = [self_copy.nuclides, other_copy.nuclides] for self_nuclide, other_nuclide in itertools.product(*all_nuclides): - new_nuclide = CrossNuclide(self_nuclide, other_nuclide, binary_op) + new_nuclide = \ + CrossNuclide(self_nuclide, other_nuclide, binary_op) new_tally.add_nuclide(new_nuclide) # Add scores to the new tally @@ -1700,7 +1686,8 @@ class Tally(object): for self_score in self_copy.scores: new_tally.add_score(self_score) else: - new_tally.num_score_bins = self_copy.num_score_bins * other_copy.num_score_bins + new_tally.num_score_bins = \ + self_copy.num_score_bins * other_copy.num_score_bins all_scores = [self_copy.scores, other_copy.scores] for self_score, other_score in itertools.product(*all_scores): new_score = CrossScore(self_score, other_score, binary_op) @@ -1717,7 +1704,6 @@ class Tally(object): def _align_tally_data(self, other, filter_product, nuclide_product, score_product): """Aligns data from two tallies for tally arithmetic. - This is a helper method to construct a dict of dicts of the "aligned" data arrays from each tally for tally arithmetic. The method analyzes the filters, scores and nuclides in both tallies and determines how to @@ -1726,7 +1712,6 @@ class Tally(object): 'tile' and 'repeat' operations to the new data arrays such that all possible combinations of the data in each tally's bins will be made when the arithmetic operation is applied to the arrays. - Parameters ---------- other : Tally @@ -1740,24 +1725,18 @@ class Tally(object): score_product : str The type of product (tensor or entrywise) to be performed between score data. - Returns ------- dict A dictionary of dictionaries to "aligned" 'mean' and 'std. dev' NumPy arrays for each tally's data. - """ - # Use dense NumPy arrays for data alignment operations - self_sparse = self.sparse - other_sparse = other.sparse - self.densify() - other.densify() - # Get the set of filters that each tally is missing - other_missing_filters = set(self.filters).difference(set(other.filters)) - self_missing_filters = set(other.filters).difference(set(self.filters)) + other_missing_filters = \ + set(self.filters).difference(set(other.filters)) + self_missing_filters = \ + set(other.filters).difference(set(self.filters)) # Add filters present in self but not in other to other for filter in other_missing_filters: @@ -1784,30 +1763,40 @@ class Tally(object): # Repeat and tile the data by nuclide in preparation for performing # the tensor product across nuclides. if nuclide_product == 'tensor': - self._mean = np.repeat(self.mean, other.num_nuclides, axis=1) - self._std_dev = np.repeat(self.std_dev, other.num_nuclides, axis=1) - other._mean = np.tile(other.mean, (1, self.num_nuclides, 1)) - other._std_dev = np.tile(other.std_dev, (1, self.num_nuclides, 1)) + self._mean = \ + np.repeat(self.mean, other.num_nuclides, axis=1) + self._std_dev = \ + np.repeat(self.std_dev, other.num_nuclides, axis=1) + other._mean = \ + np.tile(other.mean, (1, self.num_nuclides, 1)) + other._std_dev = \ + np.tile(other.std_dev, (1, self.num_nuclides, 1)) # Add nuclides to each tally such that each tally contains the complete - # set of nuclides necessary to perform an entrywise product. New nuclides - # added to a tally will have all their scores set to zero. + # set of nuclides necessary to perform an entrywise product. New + # nuclides added to a tally will have all their scores set to zero. else: # Get the set of nuclides that each tally is missing - other_missing_nuclides = set(self.nuclides).difference(set(other.nuclides)) - self_missing_nuclides = set(other.nuclides).difference(set(self.nuclides)) + other_missing_nuclides = \ + set(self.nuclides).difference(set(other.nuclides)) + self_missing_nuclides = \ + set(other.nuclides).difference(set(self.nuclides)) # Add nuclides present in self but not in other to other for nuclide in other_missing_nuclides: - other._mean = np.insert(other.mean, other.num_nuclides, 0, axis=1) - other._std_dev = np.insert(other.std_dev, other.num_nuclides, 0, axis=1) + other._mean = \ + np.insert(other.mean, other.num_nuclides, 0, axis=1) + other._std_dev = \ + np.insert(other.std_dev, other.num_nuclides, 0, axis=1) other.add_nuclide(nuclide) # Add nuclides present in other but not in self to self for nuclide in self_missing_nuclides: - self._mean = np.insert(self.mean, self.num_nuclides, 0, axis=1) - self._std_dev = np.insert(self.std_dev, self.num_nuclides, 0, axis=1) + self._mean = \ + np.insert(self.mean, self.num_nuclides, 0, axis=1) + self._std_dev = \ + np.insert(self.std_dev, self.num_nuclides, 0, axis=1) self.add_nuclide(nuclide) # Align other nuclides with self nuclides @@ -1821,30 +1810,40 @@ class Tally(object): # Repeat and tile the data by score in preparation for performing # the tensor product across scores. if score_product == 'tensor': - self._mean = np.repeat(self.mean, other.num_score_bins, axis=2) - self._std_dev = np.repeat(self.std_dev, other.num_score_bins, axis=2) - other._mean = np.tile(other.mean, (1, 1, self.num_score_bins)) - other._std_dev = np.tile(other.std_dev, (1, 1, self.num_score_bins)) + self._mean = \ + np.repeat(self.mean, other.num_score_bins, axis=2) + self._std_dev = \ + np.repeat(self.std_dev, other.num_score_bins, axis=2) + other._mean = \ + np.tile(other.mean, (1, 1, self.num_score_bins)) + other._std_dev = \ + np.tile(other.std_dev, (1, 1, self.num_score_bins)) - # Add scores to each tally such that each tally contains the complete set - # of scores necessary to perform an entrywise product. New scores added - # to a tally will be set to zero. + # Add scores to each tally such that each tally contains the complete + # set of scores necessary to perform an entrywise product. New scores + # added to a tally will be set to zero. else: # Get the set of scores that each tally is missing - other_missing_scores = set(self.scores).difference(set(other.scores)) - self_missing_scores = set(other.scores).difference(set(self.scores)) + other_missing_scores = \ + set(self.scores).difference(set(other.scores)) + self_missing_scores = \ + set(other.scores).difference(set(self.scores)) # Add scores present in self but not in other to other for score in other_missing_scores: - other._mean = np.insert(other.mean, other.num_score_bins, 0, axis=2) - other._std_dev = np.insert(other.std_dev, other.num_score_bins, 0, axis=2) + other._mean = \ + np.insert(other.mean, other.num_score_bins, 0, axis=2) + other._std_dev = \ + np.insert(other.std_dev, other.num_score_bins, 0, axis=2) other.add_score(score) # Add scores present in other but not in self to self for score in self_missing_scores: - self._mean = np.insert(self.mean, self.num_score_bins, 0, axis=2) - self._std_dev = np.insert(self.std_dev, self.num_score_bins, 0, axis=2) + self._mean = \ + np.insert(self.mean, self.num_score_bins, 0, axis=2) + self._std_dev = \ + np.insert(self.std_dev, self.num_score_bins, 0, axis=2) self.add_score(score) # Align other scores with self scores @@ -1867,12 +1866,6 @@ class Tally(object): filter.stride = stride stride *= filter.num_bins - # Restore tally operands to sparse storage - if self.sparse: - self.sparsify() - if other.sparse: - other.sparsify() - # Deep copy the mean and std dev data self_mean = copy.deepcopy(self.mean) self_std_dev = copy.deepcopy(self.std_dev) @@ -1886,6 +1879,7 @@ class Tally(object): data['other']['mean'] = other_mean data['self']['std. dev.'] = self_std_dev data['other']['std. dev.'] = other_std_dev + return data def swap_filters(self, filter1, filter2, inplace=False): @@ -2169,6 +2163,10 @@ class Tally(object): if isinstance(other, Tally): new_tally = self.hybrid_product(other, binary_op='+') + # If both tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True + elif isinstance(other, Real): new_tally = Tally(name='derived') new_tally._derived = True @@ -2188,9 +2186,8 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) - # If original tally operands were sparse, sparsify the sliced tally - if self.sparse and other.sparse: - new_tally.sparsify() + # If this tally operand is sparse, sparsify the new tally + new_tally.sparse = self.sparse else: msg = 'Unable to add "{0}" to Tally ID="{1}"'.format(other, self.id) @@ -2242,6 +2239,10 @@ class Tally(object): if isinstance(other, Tally): new_tally = self.hybrid_product(other, binary_op='-') + # If both tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True + elif isinstance(other, Real): new_tally = Tally(name='derived') new_tally._derived = True @@ -2260,9 +2261,8 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) - # If original tally operands were sparse, sparsify the sliced tally - if self.sparse and other.sparse: - new_tally.sparsify() + # If this tally operand is sparse, sparsify the new tally + new_tally.sparse = self.sparse else: msg = 'Unable to subtract "{0}" from Tally ' \ @@ -2315,6 +2315,10 @@ class Tally(object): if isinstance(other, Tally): new_tally = self.hybrid_product(other, binary_op='*') + # If original tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True + elif isinstance(other, Real): new_tally = Tally(name='derived') new_tally._derived = True @@ -2333,9 +2337,8 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) - # If original tally operands were sparse, sparsify the sliced tally - if self.sparse and other.sparse: - new_tally.sparsify() + # If this tally operand is sparse, sparsify the new tally + new_tally.sparse = self.sparse else: msg = 'Unable to multiply Tally ID="{0}" ' \ @@ -2388,6 +2391,10 @@ class Tally(object): if isinstance(other, Tally): new_tally = self.hybrid_product(other, binary_op='/') + # If original tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True + elif isinstance(other, Real): new_tally = Tally(name='derived') new_tally._derived = True @@ -2406,9 +2413,8 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) - # If original tally operands were sparse, sparsify the sliced tally - if self.sparse and other.sparse: - new_tally.sparsify() + # If this tally operand is sparse, sparsify the new tally + new_tally.sparse = self.sparse else: msg = 'Unable to divide Tally ID="{0}" ' \ @@ -2464,6 +2470,10 @@ class Tally(object): if isinstance(power, Tally): new_tally = self.hybrid_product(power, binary_op='^') + # If original tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True + elif isinstance(power, Real): new_tally = Tally(name='derived') new_tally._derived = True @@ -2484,8 +2494,7 @@ class Tally(object): new_tally.add_score(score) # If original tally was sparse, sparsify the exponentiated tally - if self.sparse: - new_tally.sparsify() + new_tally.sparse = self.sparse else: msg = 'Unable to raise Tally ID="{0}" to ' \ @@ -2648,6 +2657,7 @@ class Tally(object): raise ValueError(msg) new_tally = copy.deepcopy(self) + new_tally.sparse = False if self.sum is not None: new_sum = self.get_values(scores, filters, filter_bins, @@ -2663,7 +2673,7 @@ class Tally(object): new_tally._mean = new_mean if self.std_dev is not None: new_std_dev = self.get_values(scores, filters, filter_bins, - nuclides, 'std_dev') + nuclides, 'std_dev') new_tally._std_dev = new_std_dev # SCORES @@ -2725,9 +2735,7 @@ class Tally(object): stride *= filter.num_bins # If original tally was sparse, sparsify the sliced tally - if self.sparse: - new_tally.sparsify() - + new_tally.sparse = self.sparse return new_tally def summation(self, scores=[], filter_type=None, @@ -2831,9 +2839,7 @@ class Tally(object): tally_sum.add_filter(filters[-1]) # If original tally was sparse, sparsify the tally summation - if self.sparse: - tally_sum.sparsify() - + tally_sum.sparse = self.sparse return tally_sum def diagonalize_filter(self, new_filter): @@ -2909,9 +2915,7 @@ class Tally(object): stride *= filter.num_bins # If original tally was sparse, sparsify the diagonalized tally - if self.sparse: - new_tally.sparsify - + new_tally.sparse = self.sparse return new_tally From 16fcc9300caf8e5bb0a52ca648ecacf61a8b3aee Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Mon, 14 Dec 2015 09:32:10 -0500 Subject: [PATCH 07/12] Remove unused nonzero method and added docstring for shape property to Tally Python API class --- openmc/tallies.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index b57a0d228f..0954122d80 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -70,6 +70,9 @@ class Tally(object): Total number of filter bins accounting for all filters num_bins : Integral Total number of bins for the tally + shape : 3-tuple of Integral + The shape of the tally data array ordered as the number of filter bins, + nuclide bins and score bins num_realizations : Integral Total number of realizations with_summary : bool @@ -251,10 +254,6 @@ class Tally(object): def scores(self): return self._scores - @property - def shape(self): - return (self.num_filter_bins, self.num_nuclides, self.num_score_bins) - @property def num_scores(self): return len(self._scores) @@ -279,6 +278,10 @@ class Tally(object): num_bins *= self.num_score_bins return num_bins + @property + def shape(self): + return (self.num_filter_bins, self.num_nuclides, self.num_score_bins) + @property def estimator(self): return self._estimator @@ -312,10 +315,6 @@ class Tally(object): sum = data['sum'] sum_sq = data['sum_sq'] - # Define a routine to convert 0 to 1 - def nonzero(val): - return 1 if not val else val - # Reshape the results arrays sum = np.reshape(sum, self.shape) sum_sq = np.reshape(sum_sq, self.shape) From 38676476c39d806b513f167df554337792b69049 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Mon, 4 Jan 2016 19:40:49 -0500 Subject: [PATCH 08/12] Removed manual tuple construction for tally shape in place of new shape property --- openmc/mgxs/mgxs.py | 12 ++++-------- openmc/tallies.py | 6 ------ 2 files changed, 4 insertions(+), 14 deletions(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index e9ce00254d..968504a06d 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -789,10 +789,8 @@ class MGXS(object): std_dev = np.sqrt(std_dev) # Reshape condensed data arrays with one dimension for all filters - new_shape = \ - (tally.num_filter_bins, tally.num_nuclides, tally.num_scores,) - mean = np.reshape(mean, new_shape) - std_dev = np.reshape(std_dev, new_shape) + mean = np.reshape(mean, tally.shape) + std_dev = np.reshape(std_dev, tally.shape) # Override tally's data with the new condensed data tally._mean = mean @@ -873,10 +871,8 @@ class MGXS(object): domain_filter.num_bins = 1 # Reshape averaged data arrays with one dimension for all filters - new_shape = \ - (tally.num_filter_bins, tally.num_nuclides, tally.num_scores,) - mean = np.reshape(mean, new_shape) - std_dev = np.reshape(std_dev, new_shape) + mean = np.reshape(mean, tally.shape) + std_dev = np.reshape(std_dev, tally.shape) # Override tally's data with the new condensed data tally._mean = mean diff --git a/openmc/tallies.py b/openmc/tallies.py index 8764e2e17c..3a8f102ce3 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -2841,12 +2841,6 @@ class Tally(object): new_tally = copy.deepcopy(self) new_tally.add_filter(new_filter) - # Determine the shape of data in the new diagonalized Tally - num_filter_bins = new_tally.num_filter_bins - num_nuclides = new_tally.num_nuclides - num_scores = new_tally.num_scores - new_shape = (num_filter_bins, num_nuclides, num_scores) - # Determine "base" indices along the new "diagonal", and the factor # by which the "base" indices should be repeated to account for all # other filter bins in the diagonalized tally From 1af9c8b6a290764abc86b256dbcdd5fd1e1f80af Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Mon, 4 Jan 2016 20:46:26 -0500 Subject: [PATCH 09/12] Fixed bug when getting Pandas DataFrame from distribcell tally with 2D Lattice --- openmc/opencg_compatible.py | 4 +++- openmc/statepoint.py | 41 +++++++++++++++++++++++++++++++++---- openmc/universe.py | 2 +- 3 files changed, 41 insertions(+), 6 deletions(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index 61485b5099..0bda48c160 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -901,7 +901,9 @@ def get_opencg_lattice(openmc_lattice): # Convert 2D universes array to 3D for OpenCG if len(universes.shape) == 2: - universes.shape = (1,) + universes.shape + new_universes = universes.copy() + new_universes.shape = (1,) + universes.shape + universes = new_universes # Initialize an empty array for the OpenCG nested Universes in this Lattice universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 2b1e57dfa4..3c0759f124 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -75,6 +75,9 @@ class StatePoint(object): energy of the source site. source_present : bool Indicate whether source sites are present + sparse : bool + Whether or not the tallies uses SciPy's LIL sparse matrix format for + compressed data storage tallies : dict Dictionary whose keys are tally IDs and whose values are Tally objects tallies_present : bool @@ -110,6 +113,7 @@ class StatePoint(object): self._tallies_read = False self._summary = False self._global_tallies = None + self._sparse = False def close(self): self._f.close() @@ -318,6 +322,10 @@ class StatePoint(object): def source_present(self): return self._f['source_present'].value > 0 + @property + def sparse(self): + return self._sparse + @property def tallies(self): if not self._tallies_read: @@ -340,7 +348,8 @@ class StatePoint(object): for tally_key in tally_keys: # Read the Tally size specifications - n_realizations = self._f['{0}{1}/n_realizations'.format(base, tally_key)].value + n_realizations = \ + self._f['{0}{1}/n_realizations'.format(base, tally_key)].value # Create Tally object and assign basic properties tally = openmc.Tally(tally_id=tally_key) @@ -350,7 +359,8 @@ class StatePoint(object): tally.num_realizations = n_realizations # Read the number of Filters - n_filters = self._f['{0}{1}/n_filters'.format(base, tally_key)].value + n_filters = \ + self._f['{0}{1}/n_filters'.format(base, tally_key)].value subbase = '{0}{1}/filter '.format(base, tally_key) @@ -358,7 +368,8 @@ class StatePoint(object): for j in range(1, n_filters+1): # Read the Filter type - filter_type = self._f['{0}{1}/type'.format(subbase, j)].value.decode() + filter_type = \ + self._f['{0}{1}/type'.format(subbase, j)].value.decode() n_bins = self._f['{0}{1}/n_bins'.format(subbase, j)].value @@ -380,7 +391,8 @@ class StatePoint(object): tally.add_filter(filter) # Read Nuclide bins - nuclide_names = self._f['{0}{1}/nuclides'.format(base, tally_key)].value + nuclide_names = \ + self._f['{0}{1}/nuclides'.format(base, tally_key)].value # Add all Nuclides to the Tally for name in nuclide_names: @@ -415,6 +427,7 @@ class StatePoint(object): tally.add_score(score) # Add Tally to the global dictionary of all Tallies + tally.sparse = self.sparse self._tallies[tally_key] = tally self._tallies_read = True @@ -439,6 +452,26 @@ class StatePoint(object): def with_summary(self): return False if self.summary is None else True + @sparse.setter + def sparse(self, sparse): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices, and vice versa. + + This property may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within each Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + """ + + cv.check_type('sparse', sparse, bool) + self._sparse = sparse + + # Update tally sparsities + if self._tallies_read: + for tally_id in self.tallies: + self.tallies[tally_id].sparse = self.sparse + def get_tally(self, scores=[], filters=[], nuclides=[], name=None, id=None, estimator=None): """Finds and returns a Tally object with certain properties. diff --git a/openmc/universe.py b/openmc/universe.py index b997e8845c..8417f27c71 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -1081,7 +1081,7 @@ class RectLattice(Lattice): # For 2D Lattices if len(self._dimension) == 2: offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_index-1] - offset += self._universes[i[1]][i[2]].get_cell_instance(path, + offset += self._universes[i[1]-1][i[2]-1].get_cell_instance(path, distribcell_index) # For 3D Lattices From b516520e6c79d63204bee107ee89e70dad2f9616 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Tue, 5 Jan 2016 09:24:35 -0500 Subject: [PATCH 10/12] Reinserted blank lines into docstrings in tallies.py --- openmc/tallies.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index 3a8f102ce3..8bfdeb3fc9 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1701,6 +1701,7 @@ class Tally(object): def _align_tally_data(self, other, filter_product, nuclide_product, score_product): """Aligns data from two tallies for tally arithmetic. + This is a helper method to construct a dict of dicts of the "aligned" data arrays from each tally for tally arithmetic. The method analyzes the filters, scores and nuclides in both tallies and determines how to @@ -1709,6 +1710,7 @@ class Tally(object): 'tile' and 'repeat' operations to the new data arrays such that all possible combinations of the data in each tally's bins will be made when the arithmetic operation is applied to the arrays. + Parameters ---------- other : Tally @@ -1722,11 +1724,13 @@ class Tally(object): score_product : {'tensor', 'entrywise'} The type of product (tensor or entrywise) to be performed between score data. + Returns ------- dict A dictionary of dictionaries to "aligned" 'mean' and 'std. dev' NumPy arrays for each tally's data. + """ # Get the set of filters that each tally is missing From d2c9e1808270f72827a403c562ef5df12925dbe7 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Tue, 5 Jan 2016 09:37:38 -0500 Subject: [PATCH 11/12] Added scipy to install_requires in setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 907d80e31b..935e8a3658 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,7 @@ kwargs = {'name': 'openmc', if have_setuptools: kwargs.update({ # Required dependencies - 'install_requires': ['numpy', 'h5py', 'matplotlib'], + 'install_requires': ['numpy', 'h5py', 'matplotlib', 'scipy'], # Optional dependencies 'extras_require': { From 3eb10519ac566dc92f0c6775f5a465f7e1b7f5a4 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Wed, 6 Jan 2016 08:45:25 -0500 Subject: [PATCH 12/12] Made SciPy imports for sparse tallies optional --- openmc/tallies.py | 10 +++++++++- setup.py | 3 ++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 8bfdeb3fc9..7ec1ca50cb 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -10,7 +10,6 @@ from xml.etree import ElementTree as ET import sys import numpy as np -import scipy.sparse as sps from openmc import Mesh, Filter, Trigger, Nuclide from openmc.cross import CrossScore, CrossNuclide, CrossFilter @@ -323,6 +322,8 @@ class Tally(object): # Convert NumPy arrays to SciPy sparse LIL matrices if self.sparse: + import scipy.sparse as sps + self._sum = \ sps.lil_matrix(self._sum.flatten(), self._sum.shape) self._sum_sq = \ @@ -363,6 +364,8 @@ class Tally(object): # Convert NumPy array to SciPy sparse LIL matrix if self.sparse: + import scipy.sparse as sps + self._mean = \ sps.lil_matrix(self._mean.flatten(), self._mean.shape) @@ -385,6 +388,8 @@ class Tally(object): # Convert NumPy array to SciPy sparse LIL matrix if self.sparse: + import scipy.sparse as sps + self._std_dev = \ sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape) @@ -550,6 +555,8 @@ class Tally(object): # Convert NumPy arrays to SciPy sparse LIL matrices if sparse and not self.sparse: + import scipy.sparse as sps + if self._sum is not None: self._sum = \ sps.lil_matrix(self._sum.flatten(), self._sum.shape) @@ -562,6 +569,7 @@ class Tally(object): if self._std_dev is not None: self._std_dev = \ sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape) + self._sparse = True # Convert SciPy sparse LIL matrices to NumPy arrays diff --git a/setup.py b/setup.py index 935e8a3658..d401f5fbaf 100644 --- a/setup.py +++ b/setup.py @@ -32,11 +32,12 @@ kwargs = {'name': 'openmc', if have_setuptools: kwargs.update({ # Required dependencies - 'install_requires': ['numpy', 'h5py', 'matplotlib', 'scipy'], + 'install_requires': ['numpy', 'h5py', 'matplotlib'], # Optional dependencies 'extras_require': { 'pandas': ['pandas'], + 'sparse' : ['scipy'], 'vtk': ['vtk', 'silomesh'], 'validate': ['lxml'] }})