13183515203

Committed 06 Feb 2025 04:36PM UTC coverage: 82.601% (-2.3%) from 84.867%

Build # 13183515203

Build Type

Pull #3087

github

Committed by

web-flow

Commit Message

Merge d68c72d5e into 6e0f156d3

Pull Request Pull Request #3087: wheel building with scikit build core

Run Details

107123 of 129687 relevant lines covered (82.6%)

12608333.34 hits per line

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

85.91

/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 <gsl/gsl-lite.hpp>

#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(gsl::span<const double> p)
{
  assign(p);
}

void DiscreteIndex::assign(gsl::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 {
      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;
    }
  }
}

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

openmc-dev / openmc / 13183515203

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