15900894277

Committed 26 Jun 2025 11:42AM UTC coverage: 85.214% (-0.04%) from 85.252%

Build # 15900894277

Build Type

Pull #3452

github

Committed by

web-flow

Commit Message

Merge a81fd95f9 into 5c1021446

Pull Request Pull Request #3452: relaxed random ray constraints

Run Details

237 of 306 new or added lines in 8 files covered. (77.45%)

2 existing lines in 1 file now uncovered.

52791 of 61951 relevant lines covered (85.21%)

36510182.52 hits per line

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

79.43

/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/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_dist.h"
#include "openmc/random_lcg.h"
#include "openmc/vector.h"
#include "openmc/xml_interface.h"

namespace openmc {

//==============================================================================
// Distribution implementation
//==============================================================================

double Distribution::integral(double x0, double x1) const
{
  // Exit early if integral is over support
  auto sup = support();
  if ((x0 <= sup.first) && (sup.second <= x1))
    return 1.0;

  int32_t integral = 0;
  uint64_t seed = init_seed(0, 0);
  int32_t n = 1000000;
  for (auto i = 0; i < n; ++i) {
    double s = sample(&seed);
    if ((x0 <= s) && (s <= x1))
      ++integral;
  }
  return static_cast<float>(integral) / n;
}

//==============================================================================
// 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(vector<double>::iterator x) const
{
  // Alias sampling of discrete distribution
  size_t n = prob_.size();
  if (n > 1) {
    size_t u = *(x++) * n;
    if (*(x++) < 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);
  if (n > 1) {
    dims_ = 2;
  } else {
    dims_ = 0;
  }
}

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

  x_.assign(x, x + n);
  if (n > 1) {
    dims_ = 2;
  } else {
    dims_ = 0;
  }
}

double Discrete::sample(vector<double>::iterator x) const
{
  return x_[di_.sample(x)];
}

double Discrete::integral(double x0, double x1) const
{
  double integral = 0.0;
  size_t n = x_.size();
  for (int i = 0; i < n; ++i) {
    int j = di_.alias()[i];
    double p = di_.prob()[i];
    if ((x0 <= x_[i]) && (x_[i] <= x1))
      integral += p;
    if ((x0 <= x_[j]) && (x_[j] <= x1))
      integral += 1.0 - p;
  }
  return integral / n * di_.integral();
}

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

double Uniform::sample(vector<double>::iterator x) const
{
  return a_ + *(x++) * (b_ - a_);
}

double Uniform::integral(double x0, double x1) const
{
  return (std::min(x1, b_) - std::max(a_, x0)) / (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);
  dims_ = 1;
}

double PowerLaw::integral(double x0, double x1) const
{
  return monodiff(std::min(x1, b()), std::max(a(), x0), n() + 1.0) /
         monodiff(b(), a(), n() + 1.0);
}

double PowerLaw::sample(vector<double>::iterator x) const
{
  return std::pow(offset_ + *(x++) * 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);
}

double Maxwell::integral(double x0, double x1) const
{
  double y0 = std::sqrt(std::max(x0, 0.0) / theta_);
  double y1 = std::sqrt(x1 / theta_);
  const double ispi = 1.0 / SQRT_PI;
  return ((std::erf(y1) - std::erf(y0)) -
          2.0 * ispi * (y1 * std::exp(-y1 * y1) - y0 * std::exp(-y0 * y0)));
}

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

double Watt::integral(double x0, double x1) const
{
  double c = std::sqrt(a_ * b_) / 2.0;
  double y0 = std::sqrt(std::max(x0, 0.0) / a_);
  double y1 = std::sqrt(x1 / a_);
  const double ispi = 1.0 / SQRT_PI;
  return 0.5 * ((std::erf(y1 + c) - std::erf(y0 + c)) +
                 (std::erf(y1 - c) - std::erf(y0 - c)) +
                 ispi / c *
                   ((std::exp(-(c + y1) * (c + y1)) -
                      std::exp(-(c + y0) * (c + y0))) -
                     (std::exp(-(c - y1) * (c - y1)) -
                       std::exp(-(c - y0) * (c - y0)))));
}

//==============================================================================
// 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::integral(double x0, double x1) const
{
  double y0 = (x0 - mean_value_) / (std::sqrt(2.0) * std_dev_);
  double y1 = (x1 - mean_value_) / (std::sqrt(2.0) * std_dev_);
  return 0.5 * (std::erf(y1) - std::erf(y0));
}

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

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

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(vector<double>::iterator x) const
{
  // Sample value of CDF
  double c = *(x++);

  // 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::integral(double x0, double x1) const
{

  double integral = 0.0;

  std::size_t n = c_.size();

  double c_i = c_[0];
  int i;
  for (i = 0; i < n - 1; ++i) {
    if (x0 < x_[i + 1])
      break;
    c_i = c_[i + 1];
  }

  double c_j = c_[0];
  int j;
  for (j = 0; j < n - 1; ++j) {
    if (x1 < x_[j + 1])
      break;
    c_j = c_[j + 1];
  }

  integral += c_j - c_i;

  if (interp_ == Interpolation::histogram) {
    // Histogram interpolation
    integral -= (c_[i + 1] - c_[i]) * (x0 - x_[i]) / (x_[i + 1] - x_[i]);
    integral += (c_[j + 1] - c_[j]) * (x1 - x_[j]) / (x_[j + 1] - x_[j]);
  } else {
    // Linear-linear interpolation
    double m0 = (p_[i + 1] - p_[i]) / (x_[i + 1] - x_[i]);
    double m1 = (p_[j + 1] - p_[j]) / (x_[j + 1] - x_[j]);

    integral -= (p_[i] - m0 * x_[i]) * (x0 - x_[i]) +
                m0 / 2 * (x0 - x_[i]) * (x0 + x_[i]);
    integral += (p_[j] - m1 * x_[j]) * (x1 - x_[j]) +
                m1 / 2 * (x1 - x_[j]) * (x1 + x_[j]);
  }
  return integral;
}

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

double Equiprobable::sample(vector<double>::iterator x) const
{
  std::size_t n = x_.size();

  double r = *(x++);
  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::integral(double x0, double x1) const
{
  int32_t integral = 0;
  size_t n = x_.size();
  for (int i = 0; i < n; ++i) {
    if ((x0 <= x_[i]) && (x_[i] <= x1))
      ++integral;
  }
  return static_cast<double>(integral) / n;
}

//==============================================================================
// 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::integral(double x0, double x1) const
{
  double integral = 0.0;
  for (auto& pair : distribution_) {
    integral += pair.first * pair.second->integral(x0, x1);
  }
  return integral;
}

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

openmc-dev / openmc / 15900894277

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