From 1c219600a1459374bcd9d2fa6203c36d61075214 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 21 Sep 2021 18:34:13 +0000 Subject: [PATCH 01/15] Implementation of openmc.stats.Mixture distribution Added to_xml_element python method create nodes with probability and distribution definition in order to have the data bundled together. A different approach might be to save just two lists of "values". Implemented c++ Mixture class Parsing xml-data Sampling the distribution from the cummulative probability table --- include/openmc/distribution.h | 34 ++++++++++++++++++++++++++----- openmc/stats/univariate.py | 23 ++++++++++++++++++++- src/distribution.cpp | 38 +++++++++++++++++++++++++++++++++++ 3 files changed, 89 insertions(+), 6 deletions(-) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 44fd985d2..dcd31f429 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -24,6 +24,14 @@ public: virtual double sample(uint64_t* seed) const = 0; }; +using UPtrDist = unique_ptr; + +//! Return univariate probability distribution specified in XML file +//! \param[in] node XML node representing distribution +//! \return Unique pointer to distribution +UPtrDist distribution_from_xml(pugi::xml_node node); + + //============================================================================== //! A discrete distribution (probability mass function) //============================================================================== @@ -222,12 +230,28 @@ private: vector x_; //! Possible outcomes }; -using UPtrDist = unique_ptr; +//============================================================================== +//! Mixture distribution +//============================================================================== + +class Mixture : public Distribution { +public: + explicit Mixture(pugi::xml_node node); + + //! Sample a value from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return Sampled value + double sample(uint64_t* seed) const; + + // d and p property + const vector& d() const { return d_; } + const vector& p() const { return p_; } +private: + vector d_; //!< Pointer to sub-distributions + vector p_; //!< tabulated probability density + vector c_; //!< cumulative distribution +}; -//! Return univariate probability distribution specified in XML file -//! \param[in] node XML node representing distribution -//! \return Unique pointer to distribution -UPtrDist distribution_from_xml(pugi::xml_node node); } // namespace openmc diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 4bbc50927..c7b6dadbe 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -812,7 +812,28 @@ class Mixture(Univariate): self._distribution = distribution def to_xml_element(self, element_name): - raise NotImplementedError + """Return XML representation of the mixture distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing mixture distribution data + + """ + element = ET.Element(element_name) + element.set("type", "mixture") + + for p,d in zip(self.probability, self.distribution): + data = ET.SubElement(element, "pair") + data.set("probability", str(p)) + data.append(d.to_xml_element("dist")) + + return element @classmethod def from_xml_element(cls, elem): diff --git a/src/distribution.cpp b/src/distribution.cpp index 4713110ab..f473513ba 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -287,6 +287,42 @@ double Equiprobable::sample(uint64_t* seed) const return xl + ((n - 1) * r - i) * (xr - xl); } +//============================================================================== +// Mixture implementation +//============================================================================== + +Mixture::Mixture(pugi::xml_node node) +{ + // Read and initialize tabular distribution + //auto params = get_node_array(node, "parameters"); + + double cumsum = 0.0; + for (pugi::xml_node pair = node.child("pair"); pair; pair = pair.next_sibling("pair")) { + if (!pair.attribute("probability")) openmc::fatal_error("Mixture pair element does not have probability."); + if (!pair.child("dist")) openmc::fatal_error("Mixture pair element does not have a distribution."); + double p = std::stod(pair.attribute("probability").value()); + p_.push_back(p); + cumsum += p; + c_.push_back(cumsum); + d_.push_back(distribution_from_xml(pair.child("dist"))); + } + std::transform(c_.begin(), c_.end(), c_.begin(), [cumsum](double c){ return c/cumsum;}); +} + +double Mixture::sample(uint64_t* seed) const { + // + // Sample value of CDF + double c = prn(seed); + + // Find first CDF bin which is above the sampled value and sample corresponding sub-distribution + for (size_t i = 0; i < c_.size(); i++) { + if (c<=c_[i]) return d_[i]->sample(seed); + } + + // This should not happen. Catch it + openmc::fatal_error("Bad Sampling in Mixture Distribution."); +} + //============================================================================== // Helper function //============================================================================== @@ -315,6 +351,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node) dist = UPtrDist {new Discrete(node)}; } else if (type == "tabular") { dist = UPtrDist {new Tabular(node)}; + } else if (type == "mixture") { + dist = UPtrDist{new Mixture(node)}; } else { openmc::fatal_error("Invalid distribution type: " + type); } From e7fa663b3d7bab2dedc0b2327ac0ffa55588554e Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Thu, 23 Sep 2021 09:34:11 +0000 Subject: [PATCH 02/15] Update mixture distribution update storeage of c++ class update sampling -> uses std::lower_bound for O(log(N)) sampling of the distribution implemented from_xml_element updated settings.rst include a mixture in the regression tests --- docs/source/io_formats/settings.rst | 18 +++++--- include/openmc/distribution.h | 16 ++++--- openmc/stats/univariate.py | 22 +++++++++- src/distribution.cpp | 43 ++++++++++++------- tests/regression_tests/source/inputs_true.dat | 23 ++++++++++ .../regression_tests/source/results_true.dat | 2 +- tests/regression_tests/source/test.py | 4 +- 7 files changed, 98 insertions(+), 30 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 6ea8cfc56..b76b23e63 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -619,16 +619,17 @@ variable and whose sub-elements/attributes are as follows: :type: The type of the distribution. Valid options are "uniform", "discrete", - "tabular", "maxwell", and "watt". The "uniform" option produces variates - sampled from a uniform distribution over a finite interval. The "discrete" - option produces random variates that can assume a finite number of values - (i.e., a distribution characterized by a probability mass function). The - "tabular" option produces random variates sampled from a tabulated + "tabular", "maxwell", "watt", and "mixture". The "uniform" option producess + variates sampled from a uniform distribution over a finite interval. The + "discrete" option produces random variates that can assume a finite number + of values (i.e., a distribution characterized by a probability mass function). + The "tabular" option produces random variates sampled from a tabulated distribution where the density function is either a histogram or linearly-interpolated between tabulated points. The "watt" option produces random variates is sampled from a Watt fission spectrum (only used for energies). The "maxwell" option produce variates sampled from a Maxwell - fission spectrum (only used for energies). + fission spectrum (only used for energies). The "mixture" option produces samples + from univariate sub-distributions with given probabilities. *Default*: None @@ -649,6 +650,11 @@ variable and whose sub-elements/attributes are as follows: number :math:`a` that parameterizes the distribution :math:`p(x) dx = c x e^{-x/a} dx`. + For a "mixture" distribution, ``parameters`` provide the :math:'(p,d)' pairs + connecting the probabilites :math:'p' with the different sub-distributions + :math:'d'. All probabilities :math:'p' are given first followed by the corresponding + distributions :math:'d'. + .. note:: The above format should be used even when using the multi-group :ref:`energy_mode`. :interpolation: diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index dcd31f429..30640d07b 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -243,13 +243,17 @@ public: //! \return Sampled value double sample(uint64_t* seed) const; - // d and p property - const vector& d() const { return d_; } - const vector& p() const { return p_; } private: - vector d_; //!< Pointer to sub-distributions - vector p_; //!< tabulated probability density - vector c_; //!< cumulative distribution + struct Pair { + double cummulative_probability_ {0.0}; + UPtrDist distribution_; + // bool operator<(const Pair& o) const { + // return cummulative_probability_ < o.cummulative_probability_; + //} + bool operator<(double p) const { return cummulative_probability_ < p; } + }; + + vector distribution_; //!< sub-distributions + cummulative probabilities }; diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index c7b6dadbe..b5997b1cd 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -837,4 +837,24 @@ class Mixture(Univariate): @classmethod def from_xml_element(cls, elem): - raise NotImplementedError + """Generate mixture distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Mixture + Mixture distribution generated from XML element + + """ + P = [] + D = [] + for pair in elem: + if pair.tag == "pair": + P.append( float(get_text(pair, 'probability')) ) + D.append( Univariate.from_xml_element(pair.find("dist")) ) + + return cls(P,D) diff --git a/src/distribution.cpp b/src/distribution.cpp index f473513ba..4b52a35e0 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -297,30 +297,43 @@ Mixture::Mixture(pugi::xml_node node) //auto params = get_node_array(node, "parameters"); double cumsum = 0.0; - for (pugi::xml_node pair = node.child("pair"); pair; pair = pair.next_sibling("pair")) { + for (pugi::xml_node pair = node.child("pair"); pair; + pair = pair.next_sibling("pair")) { + // Check that required data exists if (!pair.attribute("probability")) openmc::fatal_error("Mixture pair element does not have probability."); if (!pair.child("dist")) openmc::fatal_error("Mixture pair element does not have a distribution."); - double p = std::stod(pair.attribute("probability").value()); - p_.push_back(p); - cumsum += p; - c_.push_back(cumsum); - d_.push_back(distribution_from_xml(pair.child("dist"))); + + // cummulative sum of probybilities + cumsum += std::stod(pair.attribute("probability").value()); + + // Save cummulative probybility and distrubution + distribution_.push_back( + {cumsum, distribution_from_xml(pair.child("dist"))}); + } + + // Normalize cummulative probabilities to 1 + for (auto& pair : distribution_) { + pair.cummulative_probability_ /= cumsum; } - std::transform(c_.begin(), c_.end(), c_.begin(), [cumsum](double c){ return c/cumsum;}); } -double Mixture::sample(uint64_t* seed) const { - // +double Mixture::sample(uint64_t* seed) const +{ // Sample value of CDF - double c = prn(seed); + const double p = prn(seed); - // Find first CDF bin which is above the sampled value and sample corresponding sub-distribution - for (size_t i = 0; i < c_.size(); i++) { - if (c<=c_[i]) return d_[i]->sample(seed); - } + // find matching distribution + const auto it = + std::lower_bound(distribution_.cbegin(), distribution_.cend(), p); // This should not happen. Catch it - openmc::fatal_error("Bad Sampling in Mixture Distribution."); + // TODO: Remove this check, when the code is trusted to always operate + if (it == distribution_.cend()) { + fatal_error("Bad Sampling in Mixture Distribution."); + } + + // Sample the choosen distribution + return it->distribution_->sample(seed); } //============================================================================== diff --git a/tests/regression_tests/source/inputs_true.dat b/tests/regression_tests/source/inputs_true.dat index de81a4179..eb60821c6 100644 --- a/tests/regression_tests/source/inputs_true.dat +++ b/tests/regression_tests/source/inputs_true.dat @@ -76,4 +76,27 @@ 1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23 + + + + + + -2.0 0.0 2.0 0.2 0.3 0.2 + + + + + + + + + + + + + 1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23 + + + + diff --git a/tests/regression_tests/source/results_true.dat b/tests/regression_tests/source/results_true.dat index 66b610b3a..1815324f6 100644 --- a/tests/regression_tests/source/results_true.dat +++ b/tests/regression_tests/source/results_true.dat @@ -1,2 +1,2 @@ k-combined: -3.033600E-01 1.676108E-03 +2.980096E-01 9.632798E-04 diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index 91316acc6..229ac8383 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -55,18 +55,20 @@ class SourceTestHarness(PyAPITestHarness): energy1 = openmc.stats.Maxwell(1.2895e6) energy2 = openmc.stats.Watt(0.988e6, 2.249e-6) energy3 = openmc.stats.Tabular(E, p, interpolation='histogram') + energy4 = openmc.stats.Mixture([1,2,3], [energy1, energy2, energy3]) source1 = openmc.Source(spatial1, angle1, energy1, strength=0.5) source2 = openmc.Source(spatial2, angle2, energy2, strength=0.3) source3 = openmc.Source(spatial3, angle3, energy3, strength=0.1) source4 = openmc.Source(spatial4, angle3, energy3, strength=0.1) source5 = openmc.Source(spatial5, angle3, energy3, strength=0.1) + source6 = openmc.Source(spatial5, angle3, energy4, strength=0.1) settings = openmc.Settings() settings.batches = 10 settings.inactive = 5 settings.particles = 1000 - settings.source = [source1, source2, source3, source4, source5] + settings.source = [source1, source2, source3, source4, source5, source6] settings.export_to_xml() From 60b5985746cce01c2542a54a0d2c936c9959981d Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Thu, 23 Sep 2021 09:54:39 +0000 Subject: [PATCH 03/15] update Mixture distribution use std::pair instead of custom type lambda for comparison in std::lower_bound instead of operator< --- include/openmc/distribution.h | 13 ++++--------- src/distribution.cpp | 10 +++++----- 2 files changed, 9 insertions(+), 14 deletions(-) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 30640d07b..8cdd78f78 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -244,16 +244,11 @@ public: double sample(uint64_t* seed) const; private: - struct Pair { - double cummulative_probability_ {0.0}; - UPtrDist distribution_; - // bool operator<(const Pair& o) const { - // return cummulative_probability_ < o.cummulative_probability_; - //} - bool operator<(double p) const { return cummulative_probability_ < p; } - }; + // Storrage for probability + distribution + using DistPair = std::pair; - vector distribution_; //!< sub-distributions + cummulative probabilities + vector + distribution_; //!< sub-distributions + cummulative probabilities }; diff --git a/src/distribution.cpp b/src/distribution.cpp index 4b52a35e0..e4775a5b7 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -308,12 +308,12 @@ Mixture::Mixture(pugi::xml_node node) // Save cummulative probybility and distrubution distribution_.push_back( - {cumsum, distribution_from_xml(pair.child("dist"))}); + std::make_pair(cumsum, distribution_from_xml(pair.child("dist")))); } // Normalize cummulative probabilities to 1 for (auto& pair : distribution_) { - pair.cummulative_probability_ /= cumsum; + pair.first /= cumsum; } } @@ -323,8 +323,8 @@ double Mixture::sample(uint64_t* seed) const const double p = prn(seed); // find matching distribution - const auto it = - std::lower_bound(distribution_.cbegin(), distribution_.cend(), p); + const auto it = std::lower_bound(distribution_.cbegin(), distribution_.cend(), + p, [](const DistPair& pair, double p) { return pair.first < p; }); // This should not happen. Catch it // TODO: Remove this check, when the code is trusted to always operate @@ -333,7 +333,7 @@ double Mixture::sample(uint64_t* seed) const } // Sample the choosen distribution - return it->distribution_->sample(seed); + return it->second->sample(seed); } //============================================================================== From 6612b4ba0479cfad4aea711cd640fba314ed0412 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 08:59:43 +0000 Subject: [PATCH 04/15] updated unit_test for Mixture distribution --- tests/unit_tests/test_stats.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index e15958478..9e6237ba1 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -101,8 +101,12 @@ def test_mixture(): assert mix.distribution == [d1, d2] assert len(mix) == 4 - with pytest.raises(NotImplementedError): - mix.to_xml_element('distribution') + elem = mix.to_xml_element('distribution') + + d = openmc.stats.Mixture.from_xml_element(elem) + assert d.probability == p + assert d.distribution == [d1, d2] + assert len(d) == 4 def test_polar_azimuthal(): From f93887aa21c8231e807bcd830f133dd1ab0de9be Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 17:07:31 +0200 Subject: [PATCH 05/15] Update docs/source/io_formats/settings.rst Co-authored-by: Paul Romano --- docs/source/io_formats/settings.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index b76b23e63..e73bcff28 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -619,7 +619,7 @@ variable and whose sub-elements/attributes are as follows: :type: The type of the distribution. Valid options are "uniform", "discrete", - "tabular", "maxwell", "watt", and "mixture". The "uniform" option producess + "tabular", "maxwell", "watt", and "mixture". The "uniform" option produces variates sampled from a uniform distribution over a finite interval. The "discrete" option produces random variates that can assume a finite number of values (i.e., a distribution characterized by a probability mass function). From 9123ee9239761cbb56167e8dac8876b263ff74a6 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 17:08:49 +0200 Subject: [PATCH 06/15] Update src/distribution.cpp Co-authored-by: Paul Romano --- src/distribution.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/distribution.cpp b/src/distribution.cpp index e4775a5b7..f5b413e0d 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -293,9 +293,6 @@ double Equiprobable::sample(uint64_t* seed) const Mixture::Mixture(pugi::xml_node node) { - // Read and initialize tabular distribution - //auto params = get_node_array(node, "parameters"); - double cumsum = 0.0; for (pugi::xml_node pair = node.child("pair"); pair; pair = pair.next_sibling("pair")) { From cb27271afcc4872c7d81241da05abb97084862e1 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 17:09:17 +0200 Subject: [PATCH 07/15] Update src/distribution.cpp replace with range based for loop Co-authored-by: Paul Romano --- src/distribution.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/distribution.cpp b/src/distribution.cpp index f5b413e0d..a95f6bef8 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -294,8 +294,7 @@ double Equiprobable::sample(uint64_t* seed) const Mixture::Mixture(pugi::xml_node node) { double cumsum = 0.0; - for (pugi::xml_node pair = node.child("pair"); pair; - pair = pair.next_sibling("pair")) { + for (pugi::xml_node pair : node.children("pair")) { // Check that required data exists if (!pair.attribute("probability")) openmc::fatal_error("Mixture pair element does not have probability."); if (!pair.child("dist")) openmc::fatal_error("Mixture pair element does not have a distribution."); From d7622c1b5db6395ad7f094d8f47ce2768f591210 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 17:09:42 +0200 Subject: [PATCH 08/15] Update src/distribution.cpp Co-authored-by: Paul Romano --- src/distribution.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/distribution.cpp b/src/distribution.cpp index a95f6bef8..2573695bf 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -296,8 +296,8 @@ Mixture::Mixture(pugi::xml_node node) double cumsum = 0.0; for (pugi::xml_node pair : node.children("pair")) { // Check that required data exists - if (!pair.attribute("probability")) openmc::fatal_error("Mixture pair element does not have probability."); - if (!pair.child("dist")) openmc::fatal_error("Mixture pair element does not have a distribution."); + if (!pair.attribute("probability")) fatal_error("Mixture pair element does not have probability."); + if (!pair.child("dist")) fatal_error("Mixture pair element does not have a distribution."); // cummulative sum of probybilities cumsum += std::stod(pair.attribute("probability").value()); From 3fde9412f4f28494866624cbef40a7d2fe8a911d Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 17:10:11 +0200 Subject: [PATCH 09/15] Update src/distribution.cpp Co-authored-by: Paul Romano --- src/distribution.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/distribution.cpp b/src/distribution.cpp index 2573695bf..d9505d3b2 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -323,10 +323,7 @@ double Mixture::sample(uint64_t* seed) const p, [](const DistPair& pair, double p) { return pair.first < p; }); // This should not happen. Catch it - // TODO: Remove this check, when the code is trusted to always operate - if (it == distribution_.cend()) { - fatal_error("Bad Sampling in Mixture Distribution."); - } + Ensures(it != distribution_.cend()); // Sample the choosen distribution return it->second->sample(seed); From 1c47fb205917c28bd67da200a925ebc466d70ba3 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 17:10:34 +0200 Subject: [PATCH 10/15] Update src/distribution.cpp Co-authored-by: Paul Romano --- src/distribution.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/distribution.cpp b/src/distribution.cpp index d9505d3b2..af3aac5af 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -325,7 +325,7 @@ double Mixture::sample(uint64_t* seed) const // This should not happen. Catch it Ensures(it != distribution_.cend()); - // Sample the choosen distribution + // Sample the chosen distribution return it->second->sample(seed); } From 9dc9b2af458f305c251d4572e1c1e1a335b3ea87 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 17:10:49 +0200 Subject: [PATCH 11/15] Update tests/regression_tests/source/test.py Co-authored-by: Paul Romano --- tests/regression_tests/source/test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index 229ac8383..006358ca7 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -55,7 +55,7 @@ class SourceTestHarness(PyAPITestHarness): energy1 = openmc.stats.Maxwell(1.2895e6) energy2 = openmc.stats.Watt(0.988e6, 2.249e-6) energy3 = openmc.stats.Tabular(E, p, interpolation='histogram') - energy4 = openmc.stats.Mixture([1,2,3], [energy1, energy2, energy3]) + energy4 = openmc.stats.Mixture([1, 2, 3], [energy1, energy2, energy3]) source1 = openmc.Source(spatial1, angle1, energy1, strength=0.5) source2 = openmc.Source(spatial2, angle2, energy2, strength=0.3) From 0b3cd2544baa0046f9dbf2bcffca5aca42887bf4 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 17:11:30 +0200 Subject: [PATCH 12/15] Update openmc/stats/univariate.py better naming Co-authored-by: Paul Romano --- openmc/stats/univariate.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index b5997b1cd..576b2247c 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -850,11 +850,10 @@ class Mixture(Univariate): Mixture distribution generated from XML element """ - P = [] - D = [] - for pair in elem: - if pair.tag == "pair": - P.append( float(get_text(pair, 'probability')) ) - D.append( Univariate.from_xml_element(pair.find("dist")) ) + probability = [] + distribution = [] + for pair in elem.findall('pair'): + probability.append(float(get_text(pair, 'probability'))) + distribution.append(Univariate.from_xml_element(pair.find("dist"))) - return cls(P,D) + return cls(probability, distribution) From fef3c2dc9d2e5d46c3c84f23b514c9d93fbc95cb Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 29 Sep 2021 19:54:18 +0000 Subject: [PATCH 13/15] update settings.rst --- docs/source/io_formats/settings.rst | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index e73bcff28..8a83de7ac 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -650,19 +650,24 @@ variable and whose sub-elements/attributes are as follows: number :math:`a` that parameterizes the distribution :math:`p(x) dx = c x e^{-x/a} dx`. - For a "mixture" distribution, ``parameters`` provide the :math:'(p,d)' pairs - connecting the probabilites :math:'p' with the different sub-distributions - :math:'d'. All probabilities :math:'p' are given first followed by the corresponding - distributions :math:'d'. - .. note:: The above format should be used even when using the multi-group :ref:`energy_mode`. + :interpolation: For a "tabular" distribution, ``interpolation`` can be set to "histogram" or "linear-linear" thereby specifying how tabular points are to be interpolated. *Default*: histogram +:pair: + For a "mixture" distribution, this element provides a distribution and its corresponding probability.a + + :probability: + An attribute or ``pair`` that provides the probability of a univatiate distribution within a "mixture" distribution. + + :dist: + This sub-element of a ``pair`` element provides information on the corresponding univariate distribution. + ------------------------- ```` Element ------------------------- From fe9d929ed0a40cc191faa6d5210b02e13681f2f9 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Fri, 1 Oct 2021 12:17:47 +0200 Subject: [PATCH 14/15] Apply suggestions from code review Co-authored-by: Paul Romano --- docs/source/io_formats/settings.rst | 4 ++-- openmc/stats/univariate.py | 2 +- src/distribution.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 8a83de7ac..c797212c1 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -660,10 +660,10 @@ variable and whose sub-elements/attributes are as follows: *Default*: histogram :pair: - For a "mixture" distribution, this element provides a distribution and its corresponding probability.a + For a "mixture" distribution, this element provides a distribution and its corresponding probability. :probability: - An attribute or ``pair`` that provides the probability of a univatiate distribution within a "mixture" distribution. + An attribute or ``pair`` that provides the probability of a univariate distribution within a "mixture" distribution. :dist: This sub-element of a ``pair`` element provides information on the corresponding univariate distribution. diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 576b2247c..af6c9a6d1 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -828,7 +828,7 @@ class Mixture(Univariate): element = ET.Element(element_name) element.set("type", "mixture") - for p,d in zip(self.probability, self.distribution): + for p, d in zip(self.probability, self.distribution): data = ET.SubElement(element, "pair") data.set("probability", str(p)) data.append(d.to_xml_element("dist")) diff --git a/src/distribution.cpp b/src/distribution.cpp index af3aac5af..c7952aba3 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -358,7 +358,7 @@ UPtrDist distribution_from_xml(pugi::xml_node node) } else if (type == "tabular") { dist = UPtrDist {new Tabular(node)}; } else if (type == "mixture") { - dist = UPtrDist{new Mixture(node)}; + dist = UPtrDist {new Mixture(node)}; } else { openmc::fatal_error("Invalid distribution type: " + type); } From 8991a104322554104006d769abce87499cc19f32 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Fri, 1 Oct 2021 10:23:39 +0000 Subject: [PATCH 15/15] need to include gsl-lite.hpp for Ensures() --- src/distribution.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/distribution.cpp b/src/distribution.cpp index c7952aba3..62296cb25 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -7,6 +7,8 @@ #include // for runtime_error #include // for string, stod +#include + #include "openmc/error.h" #include "openmc/math_functions.h" #include "openmc/random_dist.h"