15568603301

Committed 10 Jun 2025 07:33PM UTC coverage: 85.1% (-0.06%) from 85.158%

Build # 15568603301

Build Type

Pull #3413

github

Committed by

web-flow

Commit Message

Merge 981cb1b97 into f796fa04e

Pull Request Pull Request #3413: More interpolation types in Tabular.

Run Details

69 of 121 new or added lines in 3 files covered. (57.02%)

34 existing lines in 2 files now uncovered.

52396 of 61570 relevant lines covered (85.1%)

36601278.62 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

70.94

/src/distribution.cpp

#include "openmc/distribution.h"

#include <algorithm> // for copy
#include <array>
#include <cmath>     // for sqrt, floor, max
#include <iterator>  // for back_inserter
#include <numeric>   // for accumulate
#include <stdexcept> // for runtime_error
#include <string>    // for string, stod

#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_dist.h"
#include "openmc/random_lcg.h"
#include "openmc/xml_interface.h"

namespace openmc {

//==============================================================================
// DiscreteIndex implementation
//==============================================================================

DiscreteIndex::DiscreteIndex(pugi::xml_node node)
{
  auto params = get_node_array<double>(node, "parameters");
  std::size_t n = params.size() / 2;

  assign({params.data() + n, n});
}

DiscreteIndex::DiscreteIndex(span<const double> p)
{
  assign(p);
}

void DiscreteIndex::assign(span<const double> p)
{
  prob_.assign(p.begin(), p.end());

  this->init_alias();
}

void DiscreteIndex::init_alias()
{
  normalize();

  // The initialization and sampling method is based on Vose
  // (DOI: 10.1109/32.92917)
  // Vectors for large and small probabilities based on 1/n
  vector<size_t> large;
  vector<size_t> small;

  size_t n = prob_.size();

  // Set and allocate memory
  alias_.assign(n, 0);

  // Fill large and small vectors based on 1/n
  for (size_t i = 0; i < n; i++) {
    prob_[i] *= n;
    if (prob_[i] > 1.0) {
      large.push_back(i);
    } else {
      small.push_back(i);
    }
  }

  while (!large.empty() && !small.empty()) {
    int j = small.back();
    int k = large.back();

    // Remove last element of small
    small.pop_back();

    // Update probability and alias based on Vose's algorithm
    prob_[k] += prob_[j] - 1.0;
    alias_[j] = k;

    // Move large index to small vector, if it is no longer large
    if (prob_[k] < 1.0) {
      small.push_back(k);
      large.pop_back();
    }
  }
}

size_t DiscreteIndex::sample(uint64_t* seed) const
{
  // Alias sampling of discrete distribution
  size_t n = prob_.size();
  if (n > 1) {
    size_t u = prn(seed) * n;
    if (prn(seed) < prob_[u]) {
      return u;
    } else {
      return alias_[u];
    }
  } else {
    return 0;
  }
}

void DiscreteIndex::normalize()
{
  // Renormalize density function so that it sums to unity. Note that we save
  // the integral of the distribution so that if it is used as part of another
  // distribution (e.g., Mixture), we know its relative strength.
  integral_ = std::accumulate(prob_.begin(), prob_.end(), 0.0);
  for (auto& p_i : prob_) {
    p_i /= integral_;
  }
}

//==============================================================================
// Discrete implementation
//==============================================================================

Discrete::Discrete(pugi::xml_node node) : di_(node)
{
  auto params = get_node_array<double>(node, "parameters");

  std::size_t n = params.size() / 2;

  x_.assign(params.begin(), params.begin() + n);
}

Discrete::Discrete(const double* x, const double* p, size_t n) : di_({p, n})
{

  x_.assign(x, x + n);
}

double Discrete::sample(uint64_t* seed) const
{
  return x_[di_.sample(seed)];
}

//==============================================================================
// Uniform implementation
//==============================================================================

Uniform::Uniform(pugi::xml_node node)
{
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 2) {
    fatal_error("Uniform distribution must have two "
                "parameters specified.");
  }

  a_ = params.at(0);
  b_ = params.at(1);
}

double Uniform::sample(uint64_t* seed) const
{
  return a_ + prn(seed) * (b_ - a_);
}

//==============================================================================
// PowerLaw implementation
//==============================================================================

PowerLaw::PowerLaw(pugi::xml_node node)
{
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 3) {
    fatal_error("PowerLaw distribution must have three "
                "parameters specified.");
  }

  const double a = params.at(0);
  const double b = params.at(1);
  const double n = params.at(2);

  offset_ = std::pow(a, n + 1);
  span_ = std::pow(b, n + 1) - offset_;
  ninv_ = 1 / (n + 1);
}

double PowerLaw::sample(uint64_t* seed) const
{
  return std::pow(offset_ + prn(seed) * span_, ninv_);
}

//==============================================================================
// Maxwell implementation
//==============================================================================

Maxwell::Maxwell(pugi::xml_node node)
{
  theta_ = std::stod(get_node_value(node, "parameters"));
}

double Maxwell::sample(uint64_t* seed) const
{
  return maxwell_spectrum(theta_, seed);
}

//==============================================================================
// Watt implementation
//==============================================================================

Watt::Watt(pugi::xml_node node)
{
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 2)
    openmc::fatal_error("Watt energy distribution must have two "
                        "parameters specified.");

  a_ = params.at(0);
  b_ = params.at(1);
}

double Watt::sample(uint64_t* seed) const
{
  return watt_spectrum(a_, b_, seed);
}

//==============================================================================
// Normal implementation
//==============================================================================
Normal::Normal(pugi::xml_node node)
{
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 2) {
    openmc::fatal_error("Normal energy distribution must have two "
                        "parameters specified.");
  }

  mean_value_ = params.at(0);
  std_dev_ = params.at(1);
}

double Normal::sample(uint64_t* seed) const
{
  return normal_variate(mean_value_, std_dev_, seed);
}

//==============================================================================
// Tabular implementation
//==============================================================================

Tabular::Tabular(pugi::xml_node node)
{
  if (check_for_node(node, "interpolation")) {
    std::string temp = get_node_value(node, "interpolation");
    if (temp == "histogram") {
      interp_ = Interpolation::histogram;
    } else if (temp == "linear-linear") {
      interp_ = Interpolation::lin_lin;
    } else if (temp == "log-linear") {
      interp_ = Interpolation::log_lin;
    } else if (temp == "linear-log") {
      interp_ = Interpolation::lin_log;
    } else if (temp == "log-log") {
      interp_ = Interpolation::log_log;
    } else {
      openmc::fatal_error(
        "Unsupported interpolation type for distribution: " + temp);
    }
  } else {
    interp_ = Interpolation::histogram;
  }

  // Read and initialize tabular distribution. If number of parameters is odd,
  // add an extra zero for the 'p' array.
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() % 2 != 0) {
    params.push_back(0.0);
  }
  std::size_t n = params.size() / 2;
  const double* x = params.data();
  const double* p = x + n;
  init(x, p, n);
}

Tabular::Tabular(const double* x, const double* p, int n, Interpolation interp,
  const double* c)
  : interp_ {interp}
{
  init(x, p, n, c);
}

void Tabular::init(
  const double* x, const double* p, std::size_t n, const double* c)
{
  // Copy x/p arrays into vectors
  std::copy(x, x + n, std::back_inserter(x_));
  std::copy(p, p + n, std::back_inserter(p_));

  // Calculate cumulative distribution function
  if (c) {
    std::copy(c, c + n, std::back_inserter(c_));
  } else {
    c_.resize(n);
    c_[0] = 0.0;
    for (int i = 1; i < n; ++i) {
      if (interp_ == Interpolation::histogram) {
        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::lin_log) {
        double m = (p_[i] - p_[i - 1]) / std::log(x_[i] / x_[i - 1]);
        c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]) +
                m * (x_[i] * (std::log(x_[i] / x_[i - 1]) - 1.0) + 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] = p_[i - 1] * std::log(x_[i] / x_[i - 1]) *
                exprel(m * std::log(x_[i] / x_[i - 1]));
      } else {
        UNREACHABLE();
      }
    }
  }

  // Normalize density and distribution functions. Note that we save the
  // integral of the distribution so that if it is used as part of another
  // distribution (e.g., Mixture), we know its relative strength.
  integral_ = c_[n - 1];
  for (int i = 0; i < n; ++i) {
    p_[i] = p_[i] / integral_;
    c_[i] = c_[i] / integral_;
  }
}

double Tabular::sample(uint64_t* seed) const
{
  // Sample value of CDF
  double c = prn(seed);

  // Find first CDF bin which is above the sampled value
  double c_i = c_[0];
  int i;
  std::size_t n = c_.size();
  for (i = 0; i < n - 1; ++i) {
    if (c <= c_[i + 1])
      break;
    c_i = c_[i + 1];
  }

  // Determine bounding PDF values
  double x_i = x_[i];
  double p_i = p_[i];

  if (interp_ == Interpolation::histogram) {
    // Histogram interpolation
    if (p_i > 0.0) {
      return x_i + (c - c_i) / p_i;
    } else {
      return x_i;
    }
  } else if (interp_ == Interpolation::lin_lin) {
    // Linear-linear interpolation
    double x_i1 = x_[i + 1];
    double p_i1 = p_[i + 1];

    double m = (p_i1 - p_i) / (x_i1 - x_i);
    if (m == 0.0) {
      return x_i + (c - c_i) / p_i;
    } else {
      return x_i +
             (std::sqrt(std::max(0.0, p_i * p_i + 2 * m * (c - c_i))) - p_i) /
               m;
    }
  } else if (interp_ == Interpolation::lin_log) {
    // Linear-Log interpolation
    double x_i1 = x_[i + 1];
    double p_i1 = p_[i + 1];

    double m = (p_i1 - p_i) / std::log(x_i1 / x_i);
    double a = p_i / m - 1;
    if (m > 0.0) {
      return x_i * ((c - c_i) / (m * x_i) + a) /
             lambert_w0(((c - c_i) / (m * x_i) + a) * std::exp(a));
    } else if (m < 0.0) {
      return x_i * ((c - c_i) / (m * x_i) + a) /
             lambert_wm1(((c - c_i) / (m * x_i) + a) * std::exp(a));
    } else {
      return x_i + (c - c_i) / p_i;
    }
  } 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);
    if (m == 0.0) {
      return x_i + (c - c_i) / p_i;
    } else {
      return x_i + 1.0 / m * std::log1p(m * (c - c_i) / p_i);
    }
  } 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);
    if (m == 0.0) {
      return x_i * std::exp((c - c_i) / (p_i * x_i));
    } else {
      return x_i * std::pow(1.0 + m * (c - c_i) / (p_i * x_i), 1.0 / m);
    }
  } else {
    UNREACHABLE();
  }
}

//==============================================================================
// Equiprobable implementation
//==============================================================================

double Equiprobable::sample(uint64_t* seed) const
{
  std::size_t n = x_.size();

  double r = prn(seed);
  int i = std::floor((n - 1) * r);

  double xl = x_[i];
  double xr = x_[i + i];
  return xl + ((n - 1) * r - i) * (xr - xl);
}

//==============================================================================
// Mixture implementation
//==============================================================================

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"))
      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 probabilities
    double p = std::stod(pair.attribute("probability").value());

    // Save cummulative probability and distribution
    auto dist = distribution_from_xml(pair.child("dist"));
    cumsum += p * dist->integral();

    distribution_.push_back(std::make_pair(cumsum, std::move(dist)));
  }

  // Save integral of distribution
  integral_ = cumsum;

  // Normalize cummulative probabilities to 1
  for (auto& pair : distribution_) {
    pair.first /= cumsum;
  }
}

double Mixture::sample(uint64_t* seed) const
{
  // Sample value of CDF
  const double p = prn(seed);

  // find matching distribution
  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
  assert(it != distribution_.cend());

  // Sample the chosen distribution
  return it->second->sample(seed);
}

//==============================================================================
// Helper function
//==============================================================================

UPtrDist distribution_from_xml(pugi::xml_node node)
{
  if (!check_for_node(node, "type"))
    openmc::fatal_error("Distribution type must be specified.");

  // Determine type of distribution
  std::string type = get_node_value(node, "type", true, true);

  // Allocate extension of Distribution
  UPtrDist dist;
  if (type == "uniform") {
    dist = UPtrDist {new Uniform(node)};
  } else if (type == "powerlaw") {
    dist = UPtrDist {new PowerLaw(node)};
  } else if (type == "maxwell") {
    dist = UPtrDist {new Maxwell(node)};
  } else if (type == "watt") {
    dist = UPtrDist {new Watt(node)};
  } else if (type == "normal") {
    dist = UPtrDist {new Normal(node)};
  } else if (type == "discrete") {
    dist = UPtrDist {new Discrete(node)};
  } else if (type == "tabular") {
    dist = UPtrDist {new Tabular(node)};
  } else if (type == "mixture") {
    dist = UPtrDist {new Mixture(node)};
  } else if (type == "muir") {
    openmc::fatal_error(
      "'muir' distributions are now specified using the openmc.stats.muir() "
      "function in Python. Please regenerate your XML files.");
  } else {
    openmc::fatal_error("Invalid distribution type: " + type);
  }
  return dist;
}

} // namespace openmc

1	#include "openmc/distribution.h"
2
3	#include <algorithm> // for copy
4	#include <array>
5	#include <cmath> // for sqrt, floor, max
6	#include <iterator> // for back_inserter
7	#include <numeric> // for accumulate
8	#include <stdexcept> // for runtime_error
9	#include <string> // for string, stod
10
11	#include "openmc/error.h"
12	#include "openmc/math_functions.h"
13	#include "openmc/random_dist.h"
14	#include "openmc/random_lcg.h"
15	#include "openmc/xml_interface.h"
16
17	namespace openmc {
18
19	//==============================================================================
20	// DiscreteIndex implementation
21	//==============================================================================
22
23	DiscreteIndex::DiscreteIndex(pugi::xml_node node)	11,843✔
24	{
25	auto params = get_node_array<double>(node, "parameters");	11,843✔
26	std::size_t n = params.size() / 2;	11,843✔
27
28	assign({params.data() + n, n});	11,843✔
29	}	11,843✔
30
31	DiscreteIndex::DiscreteIndex(span<const double> p)	45,716✔
32	{
33	assign(p);	45,716✔
34	}	45,716✔
35
36	void DiscreteIndex::assign(span<const double> p)	64,529✔
37	{
38	prob_.assign(p.begin(), p.end());	64,529✔
39
40	this->init_alias();	64,529✔
41	}	64,529✔
42
43	void DiscreteIndex::init_alias()	64,529✔
44	{
45	normalize();	64,529✔
46
47	// The initialization and sampling method is based on Vose
48	// (DOI: 10.1109/32.92917)
49	// Vectors for large and small probabilities based on 1/n
50	vector<size_t> large;	64,529✔
51	vector<size_t> small;	64,529✔
52
53	size_t n = prob_.size();	64,529✔
54
55	// Set and allocate memory
56	alias_.assign(n, 0);	64,529✔
57
58	// Fill large and small vectors based on 1/n
59	for (size_t i = 0; i < n; i++) {	690,889✔
60	prob_[i] *= n;	626,360✔
61	if (prob_[i] > 1.0) {	626,360✔
62	large.push_back(i);	89,585✔
63	} else {
64	small.push_back(i);	536,775✔
65	}
66	}
67
68	while (!large.empty() && !small.empty()) {	552,504✔
69	int j = small.back();	487,975✔
70	int k = large.back();	487,975✔
71
72	// Remove last element of small
73	small.pop_back();	487,975✔
74
75	// Update probability and alias based on Vose's algorithm
76	prob_[k] += prob_[j] - 1.0;	487,975✔
77	alias_[j] = k;	487,975✔
78
79	// Move large index to small vector, if it is no longer large
80	if (prob_[k] < 1.0) {	487,975✔
81	small.push_back(k);	86,457✔
82	large.pop_back();	86,457✔
83	}
84	}
85	}	64,529✔
86
87	size_t DiscreteIndex::sample(uint64_t* seed) const	55,743,552✔
88	{
89	// Alias sampling of discrete distribution
90	size_t n = prob_.size();	55,743,552✔
91	if (n > 1) {	55,743,552✔
92	size_t u = prn(seed) * n;	12,645,394✔
93	if (prn(seed) < prob_[u]) {	12,645,394✔
94	return u;	7,246,961✔
95	} else {
96	return alias_[u];	5,398,433✔
97	}
98	} else {
99	return 0;	43,098,158✔
100	}
101	}
102
103	void DiscreteIndex::normalize()	64,529✔
104	{
105	// Renormalize density function so that it sums to unity. Note that we save
106	// the integral of the distribution so that if it is used as part of another
107	// distribution (e.g., Mixture), we know its relative strength.
108	integral_ = std::accumulate(prob_.begin(), prob_.end(), 0.0);	64,529✔
109	for (auto& p_i : prob_) {	690,889✔
110	p_i /= integral_;	626,360✔
111	}
112	}	64,529✔
113
114	//==============================================================================
115	// Discrete implementation
116	//==============================================================================
117
118	Discrete::Discrete(pugi::xml_node node) : di_(node)	11,843✔
119	{
120	auto params = get_node_array<double>(node, "parameters");	11,843✔
121
122	std::size_t n = params.size() / 2;	11,843✔
123
124	x_.assign(params.begin(), params.begin() + n);	11,843✔
125	}	11,843✔
126
127	Discrete::Discrete(const double* x, const double* p, size_t n) : di_({p, n})	45,716✔
128	{
129
130	x_.assign(x, x + n);	45,716✔
131	}	45,716✔
132
133	double Discrete::sample(uint64_t* seed) const	54,260,782✔
134	{
135	return x_[di_.sample(seed)];	54,260,782✔
136	}
137
138	//==============================================================================
139	// Uniform implementation
140	//==============================================================================
141
142	Uniform::Uniform(pugi::xml_node node)	245✔
143	{
144	auto params = get_node_array<double>(node, "parameters");	245✔
145	if (params.size() != 2) {	245✔
146	fatal_error("Uniform distribution must have two "	×
147	"parameters specified.");
148	}
149
150	a_ = params.at(0);	245✔
151	b_ = params.at(1);	245✔
152	}	245✔
153
154	double Uniform::sample(uint64_t* seed) const	285,006✔
155	{
156	return a_ + prn(seed) * (b_ - a_);	285,006✔
157	}
158
159	//==============================================================================
160	// PowerLaw implementation
161	//==============================================================================
162
163	PowerLaw::PowerLaw(pugi::xml_node node)	32✔
164	{
165	auto params = get_node_array<double>(node, "parameters");	32✔
166	if (params.size() != 3) {	32✔
167	fatal_error("PowerLaw distribution must have three "	×
168	"parameters specified.");
169	}
170
171	const double a = params.at(0);	32✔
172	const double b = params.at(1);	32✔
173	const double n = params.at(2);	32✔
174
175	offset_ = std::pow(a, n + 1);	32✔
176	span_ = std::pow(b, n + 1) - offset_;	32✔
177	ninv_ = 1 / (n + 1);	32✔
178	}	32✔
179
180	double PowerLaw::sample(uint64_t* seed) const	2,222✔
181	{
182	return std::pow(offset_ + prn(seed) * span_, ninv_);	2,222✔
183	}
184
185	//==============================================================================
186	// Maxwell implementation
187	//==============================================================================
188
189	Maxwell::Maxwell(pugi::xml_node node)	64✔
190	{
191	theta_ = std::stod(get_node_value(node, "parameters"));	64✔
192	}	64✔
193
194	double Maxwell::sample(uint64_t* seed) const	3,883✔
195	{
196	return maxwell_spectrum(theta_, seed);	3,883✔
197	}
198
199	//==============================================================================
200	// Watt implementation
201	//==============================================================================
202
203	Watt::Watt(pugi::xml_node node)	96✔
204	{
205	auto params = get_node_array<double>(node, "parameters");	96✔
206	if (params.size() != 2)	96✔
207	openmc::fatal_error("Watt energy distribution must have two "	×
208	"parameters specified.");
209
210	a_ = params.at(0);	96✔
211	b_ = params.at(1);	96✔
212	}	96✔
213
214	double Watt::sample(uint64_t* seed) const	9,256,124✔
215	{
216	return watt_spectrum(a_, b_, seed);	9,256,124✔
217	}
218
219	//==============================================================================
220	// Normal implementation
221	//==============================================================================
222	Normal::Normal(pugi::xml_node node)	×
223	{
224	auto params = get_node_array<double>(node, "parameters");	×
225	if (params.size() != 2) {	×
226	openmc::fatal_error("Normal energy distribution must have two "	×
227	"parameters specified.");
228	}
229
230	mean_value_ = params.at(0);	×
231	std_dev_ = params.at(1);	×
232	}
233
234	double Normal::sample(uint64_t* seed) const	×
235	{
236	return normal_variate(mean_value_, std_dev_, seed);	×
237	}
238
239	//==============================================================================
240	// Tabular implementation
241	//==============================================================================
242
243	Tabular::Tabular(pugi::xml_node node)	2,674✔
244	{
245	if (check_for_node(node, "interpolation")) {	2,674✔
246	std::string temp = get_node_value(node, "interpolation");	2,674✔
247	if (temp == "histogram") {	2,674✔
248	interp_ = Interpolation::histogram;	2,610✔
249	} else if (temp == "linear-linear") {	64✔
250	interp_ = Interpolation::lin_lin;	64✔
NEW 251	} else if (temp == "log-linear") {	×
NEW 252	interp_ = Interpolation::log_lin;	×
NEW 253	} else if (temp == "linear-log") {	×
NEW 254	interp_ = Interpolation::lin_log;	×
NEW 255	} else if (temp == "log-log") {	×
NEW 256	interp_ = Interpolation::log_log;	×
257	} else {
258	openmc::fatal_error(	×
259	"Unsupported interpolation type for distribution: " + temp);	×
260	}
261	} else {	2,674✔
262	interp_ = Interpolation::histogram;	×
263	}
264
265	// Read and initialize tabular distribution. If number of parameters is odd,
266	// add an extra zero for the 'p' array.
267	auto params = get_node_array<double>(node, "parameters");	2,674✔
268	if (params.size() % 2 != 0) {	2,674✔
269	params.push_back(0.0);	×
270	}
271	std::size_t n = params.size() / 2;	2,674✔
272	const double* x = params.data();	2,674✔
273	const double* p = x + n;	2,674✔
274	init(x, p, n);	2,674✔
275	}	2,674✔
276
277	Tabular::Tabular(const double* x, const double* p, int n, Interpolation interp,	42,679,390✔
278	const double* c)	42,679,390✔
279	: interp_ {interp}	42,679,390✔
280	{
281	init(x, p, n, c);	42,679,390✔
282	}	42,679,390✔
283
284	void Tabular::init(	42,682,064✔
285	const double* x, const double* p, std::size_t n, const double* c)
286	{
287	// Copy x/p arrays into vectors
288	std::copy(x, x + n, std::back_inserter(x_));	42,682,064✔
289	std::copy(p, p + n, std::back_inserter(p_));	42,682,064✔
290
291	// Calculate cumulative distribution function
292	if (c) {	42,682,064✔
293	std::copy(c, c + n, std::back_inserter(c_));	42,679,390✔
294	} else {
295	c_.resize(n);	2,674✔
296	c_[0] = 0.0;	2,674✔
297	for (int i = 1; i < n; ++i) {	35,160✔
298	if (interp_ == Interpolation::histogram) {	32,486✔
299	c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]);	32,358✔
300	} else if (interp_ == Interpolation::lin_lin) {	128✔
301	c_[i] = c_[i - 1] + 0.5 * (p_[i - 1] + p_[i]) * (x_[i] - x_[i - 1]);	128✔
NEW 302	} else if (interp_ == Interpolation::lin_log) {	×
NEW 303	double m = (p_[i] - p_[i - 1]) / std::log(x_[i] / x_[i - 1]);	×
NEW 304	c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]) +	×
NEW 305	m * (x_[i] * (std::log(x_[i] / x_[i - 1]) - 1.0) + x_[i - 1]);	×
NEW 306	} else if (interp_ == Interpolation::log_lin) {	×
NEW 307	double m = std::log(p_[i] / p_[i - 1]) / (x_[i] - x_[i - 1]);	×
NEW 308	c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]) *	×
NEW 309	exprel(m * (x_[i] - x_[i - 1]));	×
NEW 310	} else if (interp_ == Interpolation::log_log) {	×
NEW 311	double m = std::log((x_[i] * p_[i]) / (x_[i - 1] * p_[i - 1])) /	×
NEW 312	std::log(x_[i] / x_[i - 1]);	×
NEW 313	c_[i] = p_[i - 1] * std::log(x_[i] / x_[i - 1]) *	×
NEW 314	exprel(m * std::log(x_[i] / x_[i - 1]));	×
315	} else {
NEW 316	UNREACHABLE();	×
317	}
318	}
319	}
320
321	// Normalize density and distribution functions. Note that we save the
322	// integral of the distribution so that if it is used as part of another
323	// distribution (e.g., Mixture), we know its relative strength.
324	integral_ = c_[n - 1];	42,682,064✔
325	for (int i = 0; i < n; ++i) {	603,771,365✔
326	p_[i] = p_[i] / integral_;	561,089,301✔
327	c_[i] = c_[i] / integral_;	561,089,301✔
328	}
329	}	42,682,064✔
330
331	double Tabular::sample(uint64_t* seed) const	564,726,097✔
332	{
333	// Sample value of CDF
334	double c = prn(seed);	564,726,097✔
335
336	// Find first CDF bin which is above the sampled value
337	double c_i = c_[0];	564,726,097✔
338	int i;
339	std::size_t n = c_.size();	564,726,097✔
340	for (i = 0; i < n - 1; ++i) {	1,985,405,244✔
341	if (c <= c_[i + 1])	1,985,405,244✔
342	break;	564,726,097✔
343	c_i = c_[i + 1];	1,420,679,147✔
344	}
345
346	// Determine bounding PDF values
347	double x_i = x_[i];	564,726,097✔
348	double p_i = p_[i];	564,726,097✔
349
350	if (interp_ == Interpolation::histogram) {	564,726,097✔
351	// Histogram interpolation
352	if (p_i > 0.0) {	3,409,835✔
353	return x_i + (c - c_i) / p_i;	3,409,835✔
354	} else {
UNCOV 355	return x_i;	×
356	}
357	} else if (interp_ == Interpolation::lin_lin) {	561,316,262✔
358	// Linear-linear interpolation
359	double x_i1 = x_[i + 1];	561,316,262✔
360	double p_i1 = p_[i + 1];	561,316,262✔
361
362	double m = (p_i1 - p_i) / (x_i1 - x_i);	561,316,262✔
363	if (m == 0.0) {	561,316,262✔
364	return x_i + (c - c_i) / p_i;	304,439,514✔
365	} else {
366	return x_i +
367	(std::sqrt(std::max(0.0, p_i * p_i + 2 * m * (c - c_i))) - p_i) /	256,876,748✔
368	m;	256,876,748✔
369	}
UNCOV 370	} else if (interp_ == Interpolation::lin_log) {	×
371	// Linear-Log interpolation
NEW 372	double x_i1 = x_[i + 1];	×
NEW 373	double p_i1 = p_[i + 1];	×
374
NEW 375	double m = (p_i1 - p_i) / std::log(x_i1 / x_i);	×
NEW 376	double a = p_i / m - 1;	×
NEW 377	if (m > 0.0) {	×
NEW 378	return x_i * ((c - c_i) / (m * x_i) + a) /	×
NEW 379	lambert_w0(((c - c_i) / (m * x_i) + a) * std::exp(a));	×
NEW 380	} else if (m < 0.0) {	×
NEW 381	return x_i * ((c - c_i) / (m * x_i) + a) /	×
NEW 382	lambert_wm1(((c - c_i) / (m * x_i) + a) * std::exp(a));	×
383	} else {
NEW 384	return x_i + (c - c_i) / p_i;	×
385	}
NEW 386	} else if (interp_ == Interpolation::log_lin) {	×
387	// Log-linear interpolation
NEW 388	double x_i1 = x_[i + 1];	×
NEW 389	double p_i1 = p_[i + 1];	×
390
NEW 391	double m = std::log(p_i1 / p_i) / (x_i1 - x_i);	×
NEW 392	if (m == 0.0) {	×
NEW 393	return x_i + (c - c_i) / p_i;	×
394	} else {
NEW 395	return x_i + 1.0 / m * std::log1p(m * (c - c_i) / p_i);	×
396	}
NEW 397	} else if (interp_ == Interpolation::log_log) {	×
398	// Log-Log interpolation
NEW 399	double x_i1 = x_[i + 1];	×
NEW 400	double p_i1 = p_[i + 1];	×
401
NEW 402	double m = std::log((x_i1 * p_i1) / (x_i * p_i)) / std::log(x_i1 / x_i);	×
NEW 403	if (m == 0.0) {	×
NEW 404	return x_i * std::exp((c - c_i) / (p_i * x_i));	×
405	} else {
UNCOV 406	return x_i * std::pow(1.0 + m * (c - c_i) / (p_i * x_i), 1.0 / m);	×
407	}
408	} else {
UNCOV 409	UNREACHABLE();	×
410	}
411	}
412
413	//==============================================================================
414	// Equiprobable implementation
415	//==============================================================================
416
417	double Equiprobable::sample(uint64_t* seed) const	×
418	{
419	std::size_t n = x_.size();	×
420
421	double r = prn(seed);	×
UNCOV 422	int i = std::floor((n - 1) * r);	×
423
UNCOV 424	double xl = x_[i];	×
UNCOV 425	double xr = x_[i + i];	×
UNCOV 426	return xl + ((n - 1) * r - i) * (xr - xl);	×
427	}
428
429	//==============================================================================
430	// Mixture implementation
431	//==============================================================================
432
433	Mixture::Mixture(pugi::xml_node node)	48✔
434	{
435	double cumsum = 0.0;	48✔
436	for (pugi::xml_node pair : node.children("pair")) {	192✔
437	// Check that required data exists
438	if (!pair.attribute("probability"))	144✔
UNCOV 439	fatal_error("Mixture pair element does not have probability.");	×
440	if (!pair.child("dist"))	144✔
UNCOV 441	fatal_error("Mixture pair element does not have a distribution.");	×
442
443	// cummulative sum of probabilities
444	double p = std::stod(pair.attribute("probability").value());	144✔
445
446	// Save cummulative probability and distribution
447	auto dist = distribution_from_xml(pair.child("dist"));	144✔
448	cumsum += p * dist->integral();	144✔
449
450	distribution_.push_back(std::make_pair(cumsum, std::move(dist)));	144✔
451	}	144✔
452
453	// Save integral of distribution
454	integral_ = cumsum;	48✔
455
456	// Normalize cummulative probabilities to 1
457	for (auto& pair : distribution_) {	192✔
458	pair.first /= cumsum;	144✔
459	}
460	}	48✔
461
462	double Mixture::sample(uint64_t* seed) const	3,564✔
463	{
464	// Sample value of CDF
465	const double p = prn(seed);	3,564✔
466
467	// find matching distribution
468	const auto it = std::lower_bound(distribution_.cbegin(), distribution_.cend(),	3,564✔
469	p, [](const DistPair& pair, double p) { return pair.first < p; });	7,128✔
470
471	// This should not happen. Catch it
472	assert(it != distribution_.cend());	2,916✔
473
474	// Sample the chosen distribution
475	return it->second->sample(seed);	7,128✔
476	}
477
478	//==============================================================================
479	// Helper function
480	//==============================================================================
481
482	UPtrDist distribution_from_xml(pugi::xml_node node)	14,991✔
483	{
484	if (!check_for_node(node, "type"))	14,991✔
UNCOV 485	openmc::fatal_error("Distribution type must be specified.");	×
486
487	// Determine type of distribution
488	std::string type = get_node_value(node, "type", true, true);	14,991✔
489
490	// Allocate extension of Distribution
491	UPtrDist dist;	14,991✔
492	if (type == "uniform") {	14,991✔
493	dist = UPtrDist {new Uniform(node)};	245✔
494	} else if (type == "powerlaw") {	14,746✔
495	dist = UPtrDist {new PowerLaw(node)};	32✔
496	} else if (type == "maxwell") {	14,714✔
497	dist = UPtrDist {new Maxwell(node)};	64✔
498	} else if (type == "watt") {	14,650✔
499	dist = UPtrDist {new Watt(node)};	96✔
500	} else if (type == "normal") {	14,554✔
UNCOV 501	dist = UPtrDist {new Normal(node)};	×
502	} else if (type == "discrete") {	14,554✔
503	dist = UPtrDist {new Discrete(node)};	11,832✔
504	} else if (type == "tabular") {	2,722✔
505	dist = UPtrDist {new Tabular(node)};	2,674✔
506	} else if (type == "mixture") {	48✔
507	dist = UPtrDist {new Mixture(node)};	48✔
508	} else if (type == "muir") {	×
UNCOV 509	openmc::fatal_error(	×
510	"'muir' distributions are now specified using the openmc.stats.muir() "
511	"function in Python. Please regenerate your XML files.");
512	} else {
UNCOV 513	openmc::fatal_error("Invalid distribution type: " + type);	×
514	}
515	return dist;	29,982✔
516	}	14,991✔
517
518	} // namespace openmc

openmc-dev / openmc / 15568603301

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous