From d79524b9ce99d6b6eed66d9e02880d37142e37a7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 3 Apr 2019 08:06:01 -0500 Subject: [PATCH 1/3] Replace float_endf with a fast C/Cython version --- openmc/data/_endf.pyx | 8 ++++++++ openmc/data/endf.c | 37 +++++++++++++++++++++++++++++++++++ openmc/data/endf.py | 11 ++++++++++- setup.py | 4 ++-- tests/unit_tests/test_endf.py | 2 ++ 5 files changed, 59 insertions(+), 3 deletions(-) create mode 100644 openmc/data/_endf.pyx create mode 100644 openmc/data/endf.c diff --git a/openmc/data/_endf.pyx b/openmc/data/_endf.pyx new file mode 100644 index 0000000000..991ee015b2 --- /dev/null +++ b/openmc/data/_endf.pyx @@ -0,0 +1,8 @@ +# cython: c_string_type=str, c_string_encoding=ascii + +cdef extern from "endf.c": + double cfloat_endf(const char* buffer, int n) + +def float_endf(s): + cdef const char* c_string = s + return cfloat_endf(c_string, len(s)) diff --git a/openmc/data/endf.c b/openmc/data/endf.c new file mode 100644 index 0000000000..b83e130c73 --- /dev/null +++ b/openmc/data/endf.c @@ -0,0 +1,37 @@ +#include + +double cfloat_endf(const char* buffer, int n) +{ + char arr[12]; // 11 characters plus a null terminator + int j = 0; // current position in arr + int found_decimal = 0; + int found_exponent = 0; + for (int i = 0; i < n; ++i) { + // Skip whitespace characters + char c = buffer[i]; + if (c == ' ') continue; + + if (found_decimal) { + if (!found_exponent) { + if (c == '+' || c == '-') { + // In the case that we encounter +/- and we haven't yet encountered + // e/E, we manually add it + arr[j++] = 'e'; + found_exponent = 1; + + } else if (c == 'e' || c == 'E') { + found_exponent = 1; + } + } + } else if (c == '.') { + found_decimal = 1; + } + + // Copy character + arr[j++] = c; + } + + // Done copying. Add null terminator and convert to double + arr[j] = '\0'; + return atof(arr); +} diff --git a/openmc/data/endf.py b/openmc/data/endf.py index 6131f3c6b9..92a558e217 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -20,6 +20,11 @@ from numpy.polynomial.polynomial import Polynomial from .data import ATOMIC_SYMBOL, gnd_name from .function import Tabulated1D, INTERPOLATION_SCHEME from openmc.stats.univariate import Uniform, Tabular, Legendre +try: + from ._endf import float_endf + _CYTHON = True +except ImportError: + _CYTHON = False _LIBRARY = {0: 'ENDF/B', 1: 'ENDF/A', 2: 'JEFF', 3: 'EFF', @@ -69,7 +74,7 @@ SUM_RULES = {1: [2, 3], ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]) ?(\d+)') -def float_endf(s): +def py_float_endf(s): """Convert string of floating point number in ENDF to float. The ENDF-6 format uses an 'e-less' floating point number format, @@ -92,6 +97,10 @@ def float_endf(s): return float(ENDF_FLOAT_RE.sub(r'\1e\2\3', s)) +if not _CYTHON: + float_endf = py_float_endf + + def int_endf(s): """Convert string of integer number in ENDF to int. diff --git a/setup.py b/setup.py index 778df040c0..38ee57f716 100755 --- a/setup.py +++ b/setup.py @@ -68,10 +68,10 @@ kwargs = { }, } -# If Cython is present, add resonance reconstruction capability +# If Cython is present, add resonance reconstruction and fast float_endf if have_cython: kwargs.update({ - 'ext_modules': cythonize('openmc/data/reconstruct.pyx'), + 'ext_modules': cythonize('openmc/data/*.pyx'), 'include_dirs': [np.get_include()] }) diff --git a/tests/unit_tests/test_endf.py b/tests/unit_tests/test_endf.py index 147e2db873..b6f79291de 100644 --- a/tests/unit_tests/test_endf.py +++ b/tests/unit_tests/test_endf.py @@ -9,6 +9,8 @@ def test_float_endf(): assert endf.float_endf('6.022-23') == approx(6.022e-23) assert endf.float_endf(' +1.01+ 2') == approx(101.0) assert endf.float_endf(' -1.01- 2') == approx(-0.0101) + assert endf.float_endf('+ 2 . 3+ 1') == approx(23.0) + assert endf.float_endf('-7 .8 -1') == approx(-0.78) def test_int_endf(): From ae4b0106dc8171c563d51ec76b1d5d71e2289267 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 3 Apr 2019 08:06:40 -0500 Subject: [PATCH 2/3] Don't create Legendre polynomial object unless actually needed --- openmc/stats/univariate.py | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index c172c7c64b..e61f216ef0 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -312,7 +312,7 @@ class Normal(Univariate): r"""Normally distributed sampling. The Normal Distribution is characterized by two parameters - :math:`\mu` and :math:`\sigma` and has density function + :math:`\mu` and :math:`\sigma` and has density function :math:`p(X) dX = 1/(\sqrt{2\pi}\sigma) e^{-(X-\mu)^2/(2\sigma^2)}` Parameters @@ -325,9 +325,9 @@ class Normal(Univariate): Attributes ---------- mean_value : float - Mean of the Normal distribution + Mean of the Normal distribution std_dev : float - Standard deviation of the Normal distribution + Standard deviation of the Normal distribution """ def __init__(self, mean_value, std_dev): @@ -380,17 +380,17 @@ class Normal(Univariate): class Muir(Univariate): """Muir energy spectrum. - The Muir energy spectrum is a Gaussian spectrum, but for + The Muir energy spectrum is a Gaussian spectrum, but for convenience reasons allows the user 3 parameters to define the distribution, e0 the mean energy of particles, the mass - of reactants m_rat, and the ion temperature kt. + of reactants m_rat, and the ion temperature kt. Parameters ---------- e0 : float Mean of the Muir distribution in units of eV m_rat : float - Ratio of the sum of the masses of the reaction inputs to an + Ratio of the sum of the masses of the reaction inputs to an AMU kt : float Ion temperature for the Muir distribution in units of eV @@ -400,7 +400,7 @@ class Muir(Univariate): e0 : float Mean of the Muir distribution in units of eV m_rat : float - Ratio of the sum of the masses of the reaction inputs to an + Ratio of the sum of the masses of the reaction inputs to an AMU kt : float Ion temperature for the Muir distribution in units of eV @@ -582,26 +582,27 @@ class Legendre(Univariate): def __init__(self, coefficients): self.coefficients = coefficients + self._legendre_poly = None def __call__(self, x): - return self._legendre_polynomial(x) + # Create Legendre polynomial if we haven't yet + if self._legendre_poly is None: + l = np.arange(len(self._coefficients)) + coeffs = (2.*l + 1.)/2. * self._coefficients + self._legendre_poly = np.polynomial.Legendre(coeffs) + + return self._legendre_poly(x) def __len__(self): - return len(self._legendre_polynomial.coef) + return len(self._coefficients) @property def coefficients(self): - poly = self._legendre_polynomial - l = np.arange(poly.degree() + 1) - return 2./(2.*l + 1.) * poly.coef + return self._coefficients @coefficients.setter def coefficients(self, coefficients): - cv.check_type('Legendre expansion coefficients', coefficients, - Iterable, Real) - l = np.arange(len(coefficients)) - coeffs = (2.*l + 1.)/2. * np.array(coefficients) - self._legendre_polynomial = np.polynomial.Legendre(coeffs) + self._coefficients = np.asarray(coefficients) def to_xml_element(self, element_name): raise NotImplementedError From ef0d216f662327316c26072412f591f342a7312e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 4 Apr 2019 16:40:28 -0500 Subject: [PATCH 3/3] Fix float_endf for decimal-less numbers. Also secretly support d/D --- openmc/data/endf.c | 12 +++++++----- tests/unit_tests/test_endf.py | 11 +++++++++++ 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/openmc/data/endf.c b/openmc/data/endf.c index b83e130c73..fee6af179e 100644 --- a/openmc/data/endf.c +++ b/openmc/data/endf.c @@ -4,14 +4,14 @@ double cfloat_endf(const char* buffer, int n) { char arr[12]; // 11 characters plus a null terminator int j = 0; // current position in arr - int found_decimal = 0; + int found_significand = 0; int found_exponent = 0; for (int i = 0; i < n; ++i) { // Skip whitespace characters char c = buffer[i]; if (c == ' ') continue; - if (found_decimal) { + if (found_significand) { if (!found_exponent) { if (c == '+' || c == '-') { // In the case that we encounter +/- and we haven't yet encountered @@ -19,12 +19,14 @@ double cfloat_endf(const char* buffer, int n) arr[j++] = 'e'; found_exponent = 1; - } else if (c == 'e' || c == 'E') { + } else if (c == 'e' || c == 'E' || c == 'd' || c == 'D') { + arr[j++] = 'e'; found_exponent = 1; + continue; } } - } else if (c == '.') { - found_decimal = 1; + } else if (c == '.' || (c >= '0' && c <= '9')) { + found_significand = 1; } // Copy character diff --git a/tests/unit_tests/test_endf.py b/tests/unit_tests/test_endf.py index b6f79291de..11a8db4432 100644 --- a/tests/unit_tests/test_endf.py +++ b/tests/unit_tests/test_endf.py @@ -11,6 +11,17 @@ def test_float_endf(): assert endf.float_endf(' -1.01- 2') == approx(-0.0101) assert endf.float_endf('+ 2 . 3+ 1') == approx(23.0) assert endf.float_endf('-7 .8 -1') == approx(-0.78) + assert endf.float_endf('3.14e0') == approx(3.14) + assert endf.float_endf('3.14E0') == approx(3.14) + assert endf.float_endf('3.14e-1') == approx(0.314) + assert endf.float_endf('3.14d0') == approx(3.14) + assert endf.float_endf('3.14D0') == approx(3.14) + assert endf.float_endf('3.14d-1') == approx(0.314) + assert endf.float_endf('1+2') == approx(100.0) + assert endf.float_endf('-1+2') == approx(-100.0) + assert endf.float_endf('1.+2') == approx(100.0) + assert endf.float_endf('-1.+2') == approx(-100.0) + assert endf.float_endf(' ') == 0.0 def test_int_endf():