More interpolation types in Tabular.

include/openmc/math_functions.h

@@ -201,6 +201,18 @@ std::complex<double> faddeeva(std::complex<double> z);
 //! \return Derivative of Faddeeva function evaluated at z
 std::complex<double> w_derivative(std::complex<double> 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
 //!

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')

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();
   }
 }
 

src/math_functions.cpp

@@ -919,6 +919,24 @@ std::complex<double> w_derivative(std::complex<double> 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(

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,17 +274,20 @@ 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
     samples, weights = d.sample(n_samples)
     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')