17361657497

Committed 31 Aug 2025 07:54PM UTC coverage: 84.729%. First build

Build # 17361657497

Build Type

Pull #3550

github

Committed by

web-flow

Commit Message

Merge 98ca32e81 into 00edc7769

Pull Request Pull Request #3550: [Point Detector] Add distribution get_pdf functionality

Run Details

7 of 356 new or added lines in 10 files covered. (1.97%)

52947 of 62490 relevant lines covered (84.73%)

38803853.85 hits per line

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

72.14

/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/search.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)];
}

// Not implemented
double Discrete::get_pdf(double x) const
{
  return -1;
}

//==============================================================================
// 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_);
}

double Uniform::get_pdf(double x) const
{
  if (x <= b_ && x >= a_)
    return 1 / (b_ - a_);
  else
    return 0;
}

//==============================================================================
// 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_);
}

double PowerLaw::get_pdf(double x) const
{
  // Recover the lower/upper bounds from offset_ and span_
  double a = std::pow(offset_, ninv_);         // since offset_ = a^(n+1)
  double b = std::pow(offset_ + span_, ninv_); // since offset_+span_ = b^(n+1)

  if (x < a || x > b) {
    return 0.0; // outside support
  }

  double n = 1.0 / ninv_ - 1.0; // since ninv_ = 1/(n+1)

  return (n + 1.0) * std::pow(x, n) / (std::pow(b, n + 1) - std::pow(a, n + 1));
}

//==============================================================================
// 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);
}

// Not implemented
double Maxwell::get_pdf(double x) const
{
  return -1;
}

//==============================================================================
// 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);
}

// Not implemented
double Watt::get_pdf(double x) const
{
  return -1;
}

//==============================================================================
// 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);
}

double Normal::get_pdf(double x) const
{
  double exponent = -0.5 * std::pow((x - mean_value_) / std_dev_, 2);
  double coefficient = 1 / (std_dev_ * std::sqrt(2 * PI));
  return coefficient * std::exp(exponent);
}

//==============================================================================
// 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 {
      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_));

  // 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_));
  } 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]);
      }
    }
  }

  // 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 {
    // 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;
    }
  }
}

double Tabular::get_pdf(double x) const
{
  // get PDF value at x

  int i;
  std::size_t n = x_.size();
  if (x < x_[0]) {
    return 0;
  } else if (x > x_[n - 1]) {
    return 0;
  } else {
    i = lower_bound_index(x_.begin(), x_.end(), x);
  }

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

  if (interp_ == Interpolation::histogram) {
    // Histogram interpolation
    return p_i;
  } else {
    // 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 p_i;
    } else {
      return p_i + (x - x_i) * m;
    }
  }
}

//==============================================================================
// 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);
}

double Equiprobable::get_pdf(double x) const
{
  return -1;
}

//==============================================================================
// 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);
}

double Mixture::get_pdf(double x) const
{
  return -1;
}

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

openmc-dev / openmc / 17361657497

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