mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
More interpolation types in Tabular. (#3413)
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
GuySten
GitHub
parent
360ec24b41
commit
96383fcb2b
5 changed files with 201 additions and 45 deletions
|
|
@ -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
|
||||
//!
|
||||
|
|
|
|||
|
|
@ -10,6 +10,7 @@ from warnings import warn
|
|||
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 @@ _INTERPOLATION_SCHEMES = {
|
|||
}
|
||||
|
||||
|
||||
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,64 @@ class Tabular(Univariate):
|
|||
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)
|
||||
# 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]
|
||||
|
||||
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()
|
||||
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 +1426,56 @@ class Tabular(Univariate):
|
|||
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,9 +1573,22 @@ class Tabular(Univariate):
|
|||
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')
|
||||
f'integral() not supported for {self.interpolation} interpolation')
|
||||
|
||||
|
||||
class Legendre(Univariate):
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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(
|
||||
|
|
|
|||
|
|
@ -4,6 +4,7 @@ import numpy as np
|
|||
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')
|
||||
|
|
@ -291,6 +295,12 @@ def test_tabular():
|
|||
assert_sample_mean(samples, d.mean())
|
||||
assert np.all(weights == 1.0)
|
||||
|
||||
# Multiplying the probabilities should preserve the mean but change the integral
|
||||
d2 = openmc.stats.Tabular(x, p*2, interpolation='histogram')
|
||||
assert d2.mean() == pytest.approx(d.mean())
|
||||
assert d2.integral() == pytest.approx(2.0*d.integral())
|
||||
|
||||
# Normalizing should result in an integral of 1
|
||||
d.normalize()
|
||||
assert d.integral() == pytest.approx(1.0)
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue