diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index 0d960c33db7..d67b8305ef5 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -201,6 +201,18 @@ std::complex faddeeva(std::complex z); //! \return Derivative of Faddeeva function evaluated at z std::complex w_derivative(std::complex z, int order); +//! Evaluate relative exponential function +//! +//! \param x Real argument +//! \return (exp(x)-1)/x without loss of precision near 0 +double exprel(double x); + +//! Evaluate relative logarithm function +//! +//! \param x Real argument +//! \return log(1+x)/x without loss of precision near 0 +double log1prel(double x); + //! Helper function to get index and interpolation function on an incident //! energy grid //! diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 272367f6f27..edd0469076e 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -10,6 +10,7 @@ import lxml.etree as ET import numpy as np from scipy.integrate import trapezoid +from scipy.special import exprel, hyp1f1, lambertw import scipy import openmc.checkvalue as cv @@ -25,6 +26,16 @@ } +def exprel2(x): + """Evaluate 2*(exp(x)-1-x)/x^2 without loss of precision near 0""" + return hyp1f1(1, 3, x) + + +def log1prel(x): + """Evaluate log(1+x)/x without loss of precision near 0""" + return np.where(np.abs(x) < 1e-16, 1.0, np.log1p(x) / x) + + class Univariate(EqualityMixin, ABC): """Probability distribution of a single random variable. @@ -1299,44 +1310,76 @@ def cdf(self): c[1:] = p[:x.size-1] * np.diff(x) elif self.interpolation == 'linear-linear': c[1:] = 0.5 * (p[:-1] + p[1:]) * np.diff(x) + elif self.interpolation == "linear-log": + m = np.diff(p) / np.diff(np.log(x)) + c[1:] = p[:-1] * np.diff(x) + m * ( + x[1:] * (np.diff(np.log(x)) - 1.0) + x[:-1] + ) + elif self.interpolation == "log-linear": + m = np.diff(np.log(p)) / np.diff(x) + c[1:] = p[:-1] * np.diff(x) * exprel(m * np.diff(x)) + elif self.interpolation == "log-log": + m = np.diff(np.log(x * p)) / np.diff(np.log(x)) + c[1:] = (x * p)[:-1] * np.diff(np.log(x)) * exprel(m * np.diff(np.log(x))) else: - raise NotImplementedError('Can only generate CDFs for tabular ' - 'distributions using histogram or ' - 'linear-linear interpolation') - + raise NotImplementedError( + f"Cannot generate CDFs for tabular " + f"distributions using {self.interpolation} interpolation" + ) return np.cumsum(c) def mean(self): """Compute the mean of the tabular distribution""" - if self.interpolation == 'linear-linear': - mean = 0.0 - for i in range(1, len(self.x)): - y_min = self.p[i-1] - y_max = self.p[i] - x_min = self.x[i-1] - x_max = self.x[i] - - m = (y_max - y_min) / (x_max - x_min) - exp_val = (1./3.) * m * (x_max**3 - x_min**3) - exp_val += 0.5 * m * x_min * (x_min**2 - x_max**2) - exp_val += 0.5 * y_min * (x_max**2 - x_min**2) - mean += exp_val - - elif self.interpolation == 'histogram': - x_l = self.x[:-1] - x_r = self.x[1:] - p_l = self.p[:self.x.size-1] - mean = (0.5 * (x_l + x_r) * (x_r - x_l) * p_l).sum() + # use normalized probabilities when computing mean + p = self.p / self.cdf().max() + x = self.x + x_min = x[:-1] + x_max = x[1:] + p_min = p[: x.size - 1] + p_max = p[1:] + + if self.interpolation == "linear-linear": + m = np.diff(p) / np.diff(x) + mean = ( + (1.0 / 3.0) * m * np.diff(x**3) + + 0.5 * (p_min - m * x_min) * np.diff(x**2) + ).sum() + elif self.interpolation == "linear-log": + m = np.diff(p) / np.diff(np.log(x)) + mean = ( + (1.0 / 4.0) + * m + * x_min**2 + * ((x_max / x_min) ** 2 * (2 * np.diff(np.log(x)) - 1) + 1) + + +0.5 * p_min * np.diff(x**2) + ).sum() + elif self.interpolation == "log-linear": + m = np.diff(np.log(p)) / np.diff(x) + mean = ( + p_min + * ( + np.diff(x) ** 2 + * ((0.5 * exprel2(m * np.diff(x)) * (m * np.diff(x) - 1) + 1)) + + np.diff(x) * x_min * exprel(m * np.diff(x)) + ) + ).sum() + elif self.interpolation == "log-log": + m = np.diff(np.log(p)) / np.diff(np.log(x)) + mean = ( + p_min + * x_min**2 + * np.diff(np.log(x)) + * exprel((m + 2) * np.diff(np.log(x))) + ).sum() + elif self.interpolation == "histogram": + mean = (0.5 * (x_min + x_max) * np.diff(x) * p_min).sum() else: - raise NotImplementedError('Can only compute mean for tabular ' - 'distributions using histogram ' - 'or linear-linear interpolation.') - - # Normalize for when integral of distribution is not 1 - mean /= self.integral() - + raise NotImplementedError( + f"Cannot compute mean for tabular " + f"distributions using {self.interpolation} interpolation" + ) return mean def normalize(self): @@ -1395,11 +1438,56 @@ def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None): quad[quad < 0.0] = 0.0 m[non_zero] = x_i[non_zero] + (np.sqrt(quad) - p_i[non_zero]) / m[non_zero] samples_out = m + elif self.interpolation == "linear-log": + # get variable and probability values for the + # next entry + x_i1 = self.x[cdf_idx + 1] + p_i1 = p[cdf_idx + 1] + # compute slope between entries + m = (p_i1 - p_i) / np.log(x_i1 / x_i) + # set values for zero slope + zero = m == 0.0 + m[zero] = x_i[zero] + (xi[zero] - c_i[zero]) / p_i[zero] + + positive = m > 0 + negative = m < 0 + a = p_i / m - 1 + m[positive] = ( + x_i + * ((xi - c_i) / (m * x_i) + a) + / np.real(lambertw((((xi - c_i) / (m * x_i) + a)) * np.exp(a))) + )[positive] + m[negative] = ( + x_i + * ((xi - c_i) / (m * x_i) + a) + / np.real(lambertw((((xi - c_i) / (m * x_i) + a)) * np.exp(a), -1.0)) + )[negative] + samples_out = m + elif self.interpolation == "log-linear": + # get variable and probability values for the + # next entry + x_i1 = self.x[cdf_idx + 1] + p_i1 = p[cdf_idx + 1] + # compute slope between entries + m = np.log(p_i1 / p_i) / (x_i1 - x_i) + f = (xi - c_i) / p_i + samples_out = x_i + f * log1prel(m * f) + elif self.interpolation == "log-log": + # get variable and probability values for the + # next entry + x_i1 = self.x[cdf_idx + 1] + p_i1 = p[cdf_idx + 1] + # compute slope between entries + m = np.log((x_i1 * p_i1) / (x_i * p_i)) / np.log(x_i1 / x_i) + f = (xi - c_i) / (x_i * p_i) + + samples_out = x_i * np.exp(f * log1prel(m * f)) else: - raise NotImplementedError('Can only sample tabular distributions ' - 'using histogram or ' - 'linear-linear interpolation') + raise NotImplementedError( + f"Cannot sample tabular distributions " + f"for {self.inteprolation} interpolation " + ) assert all(samples_out < self.x[-1]) return samples_out @@ -1497,6 +1585,23 @@ def integral(self): return np.sum(np.diff(self.x) * self.p[:self.x.size-1]) elif self.interpolation == 'linear-linear': return trapezoid(self.p, self.x) + elif self.interpolation == "linear-log": + m = np.diff(self.p) / np.diff(np.log(self.x)) + return np.sum( + self.p[:-1] * np.diff(self.x) + + m * (self.x[1:] * (np.diff(np.log(self.x)) - 1.0) + self.x[:-1]) + ) + elif self.interpolation == "log-linear": + m = np.diff(np.log(self.p)) / np.diff(self.x) + return np.sum(self.p[:-1] * np.diff(self.x) * exprel(m * np.diff(self.x))) + elif self.interpolation == "log-log": + m = np.diff(np.log(self.p)) / np.diff(np.log(self.x)) + return np.sum( + self.p[:-1] + * self.x[:-1] + * np.diff(np.log(self.x)) + * exprel((m + 1) * np.diff(np.log(self.x))) + ) else: raise NotImplementedError( f'integral() not supported for {self.inteprolation} interpolation') diff --git a/src/distribution.cpp b/src/distribution.cpp index 537a56171d9..756b098ed91 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -401,6 +401,10 @@ Tabular::Tabular(pugi::xml_node node) interp_ = Interpolation::histogram; } else if (temp == "linear-linear") { interp_ = Interpolation::lin_lin; + } else if (temp == "log-linear") { + interp_ = Interpolation::log_lin; + } else if (temp == "log-log") { + interp_ = Interpolation::log_log; } else { openmc::fatal_error( "Unsupported interpolation type for distribution: " + temp); @@ -437,13 +441,6 @@ void Tabular::init( std::copy(x, x + n, std::back_inserter(x_)); std::copy(p, p + n, std::back_inserter(p_)); - // Check interpolation parameter - if (interp_ != Interpolation::histogram && - interp_ != Interpolation::lin_lin) { - openmc::fatal_error("Only histogram and linear-linear interpolation " - "for tabular distribution is supported."); - } - // Calculate cumulative distribution function if (c) { std::copy(c, c + n, std::back_inserter(c_)); @@ -455,6 +452,18 @@ void Tabular::init( c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]); } else if (interp_ == Interpolation::lin_lin) { c_[i] = c_[i - 1] + 0.5 * (p_[i - 1] + p_[i]) * (x_[i] - x_[i - 1]); + } else if (interp_ == Interpolation::log_lin) { + double m = std::log(p_[i] / p_[i - 1]) / (x_[i] - x_[i - 1]); + c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]) * + exprel(m * (x_[i] - x_[i - 1])); + } else if (interp_ == Interpolation::log_log) { + double m = std::log((x_[i] * p_[i]) / (x_[i - 1] * p_[i - 1])) / + std::log(x_[i] / x_[i - 1]); + c_[i] = c_[i - 1] + x_[i - 1] * p_[i - 1] * + std::log(x_[i] / x_[i - 1]) * + exprel(m * std::log(x_[i] / x_[i - 1])); + } else { + UNREACHABLE(); } } } @@ -495,7 +504,7 @@ double Tabular::sample_unbiased(uint64_t* seed) const } else { return x_i; } - } else { + } else if (interp_ == Interpolation::lin_lin) { // Linear-linear interpolation double x_i1 = x_[i + 1]; double p_i1 = p_[i + 1]; @@ -508,6 +517,24 @@ double Tabular::sample_unbiased(uint64_t* seed) const (std::sqrt(std::max(0.0, p_i * p_i + 2 * m * (c - c_i))) - p_i) / m; } + } else if (interp_ == Interpolation::log_lin) { + // Log-linear interpolation + double x_i1 = x_[i + 1]; + double p_i1 = p_[i + 1]; + + double m = std::log(p_i1 / p_i) / (x_i1 - x_i); + double f = (c - c_i) / p_i; + return x_i + f * log1prel(m * f); + } else if (interp_ == Interpolation::log_log) { + // Log-Log interpolation + double x_i1 = x_[i + 1]; + double p_i1 = p_[i + 1]; + + double m = std::log((x_i1 * p_i1) / (x_i * p_i)) / std::log(x_i1 / x_i); + double f = (c - c_i) / (p_i * x_i); + return x_i * std::exp(f * log1prel(m * f)); + } else { + UNREACHABLE(); } } diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 9473f928b2c..f45165414db 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -919,6 +919,24 @@ std::complex w_derivative(std::complex z, int order) } } +double exprel(double x) +{ + if (std::abs(x) < 1e-16) + return 1.0; + else { + return std::expm1(x) / x; + } +} + +double log1prel(double x) +{ + if (std::abs(x) < 1e-16) + return 1.0; + else { + return std::log1p(x) / x; + } +} + // Helper function to get index and interpolation function on an incident energy // grid void get_energy_index( diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index e92e9e13520..97da09f424b 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -4,6 +4,7 @@ import pytest import openmc import openmc.stats +from openmc.stats.univariate import _INTERPOLATION_SCHEMES from scipy.integrate import trapezoid from tests.unit_tests import assert_sample_mean @@ -273,7 +274,7 @@ def test_watt(): @pytest.mark.flaky(reruns=1) def test_tabular(): # test linear-linear sampling - x = np.array([0.0, 5.0, 7.0, 10.0]) + x = np.array([0.001, 5.0, 7.0, 10.0]) p = np.array([10.0, 20.0, 5.0, 6.0]) d = openmc.stats.Tabular(x, p, 'linear-linear') n_samples = 100_000 @@ -281,9 +282,12 @@ def test_tabular(): assert_sample_mean(samples, d.mean()) assert np.all(weights == 1.0) - # test linear-linear normalization - d.normalize() - assert d.integral() == pytest.approx(1.0) + for scheme in _INTERPOLATION_SCHEMES: + # test sampling + d = openmc.stats.Tabular(x, p, scheme) + n_samples = 100_000 + samples = d.sample(n_samples)[0] + assert_sample_mean(samples, d.mean()) # test histogram sampling d = openmc.stats.Tabular(x, p, interpolation='histogram')