• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 20508376487

25 Dec 2025 05:10PM UTC coverage: 82.009% (+0.07%) from 81.94%
20508376487

Pull #3413

github

web-flow
Merge f7a0ff1de into 3f06a42ab
Pull Request #3413: More interpolation types in Tabular.

17072 of 23667 branches covered (72.13%)

Branch coverage included in aggregate %.

68 of 107 new or added lines in 3 files covered. (63.55%)

119 existing lines in 1 file now uncovered.

55231 of 64498 relevant lines covered (85.63%)

43348727.91 hits per line

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

74.94
/src/distribution.cpp
1
#include "openmc/distribution.h"
2

3
#include <algorithm> // for copy
4
#include <array>
5
#include <cmath>     // for sqrt, floor, max
6
#include <iterator>  // for back_inserter
7
#include <numeric>   // for accumulate
8
#include <stdexcept> // for runtime_error
9
#include <string>    // for string, stod
10

11
#include "openmc/error.h"
12
#include "openmc/math_functions.h"
13
#include "openmc/random_dist.h"
14
#include "openmc/random_lcg.h"
15
#include "openmc/xml_interface.h"
16

17
namespace openmc {
18

19
//==============================================================================
20
// DiscreteIndex implementation
21
//==============================================================================
22

23
DiscreteIndex::DiscreteIndex(pugi::xml_node node)
26,038✔
24
{
25
  auto params = get_node_array<double>(node, "parameters");
26,038✔
26
  std::size_t n = params.size() / 2;
26,038✔
27

28
  assign({params.data() + n, n});
26,038✔
29
}
26,038✔
30

31
DiscreteIndex::DiscreteIndex(span<const double> p)
46,867✔
32
{
33
  assign(p);
46,867✔
34
}
46,867✔
35

36
void DiscreteIndex::assign(span<const double> p)
80,806✔
37
{
38
  prob_.assign(p.begin(), p.end());
80,806✔
39

40
  this->init_alias();
80,806✔
41
}
80,806✔
42

43
void DiscreteIndex::init_alias()
80,806✔
44
{
45
  normalize();
80,806✔
46

47
  // The initialization and sampling method is based on Vose
48
  // (DOI: 10.1109/32.92917)
49
  // Vectors for large and small probabilities based on 1/n
50
  vector<size_t> large;
80,806✔
51
  vector<size_t> small;
80,806✔
52

53
  size_t n = prob_.size();
80,806✔
54

55
  // Set and allocate memory
56
  alias_.assign(n, 0);
80,806✔
57

58
  // Fill large and small vectors based on 1/n
59
  for (size_t i = 0; i < n; i++) {
1,551,214✔
60
    prob_[i] *= n;
1,470,408✔
61
    if (prob_[i] > 1.0) {
1,470,408✔
62
      large.push_back(i);
219,796✔
63
    } else {
64
      small.push_back(i);
1,250,612✔
65
    }
66
  }
67

68
  while (!large.empty() && !small.empty()) {
1,395,913✔
69
    int j = small.back();
1,315,107✔
70
    int k = large.back();
1,315,107✔
71

72
    // Remove last element of small
73
    small.pop_back();
1,315,107✔
74

75
    // Update probability and alias based on Vose's algorithm
76
    prob_[k] += prob_[j] - 1.0;
1,315,107✔
77
    alias_[j] = k;
1,315,107✔
78

79
    // Move large index to small vector, if it is no longer large
80
    if (prob_[k] < 1.0) {
1,315,107✔
81
      small.push_back(k);
211,879✔
82
      large.pop_back();
211,879✔
83
    }
84
  }
85
}
80,806✔
86

87
size_t DiscreteIndex::sample(uint64_t* seed) const
64,143,116✔
88
{
89
  // Alias sampling of discrete distribution
90
  size_t n = prob_.size();
64,143,116✔
91
  if (n > 1) {
64,143,116✔
92
    size_t u = prn(seed) * n;
15,730,684✔
93
    if (prn(seed) < prob_[u]) {
15,730,684✔
94
      return u;
9,295,909✔
95
    } else {
96
      return alias_[u];
6,434,775✔
97
    }
98
  } else {
99
    return 0;
48,412,432✔
100
  }
101
}
102

103
void DiscreteIndex::normalize()
80,806✔
104
{
105
  // Renormalize density function so that it sums to unity. Note that we save
106
  // the integral of the distribution so that if it is used as part of another
107
  // distribution (e.g., Mixture), we know its relative strength.
108
  integral_ = std::accumulate(prob_.begin(), prob_.end(), 0.0);
80,806✔
109
  for (auto& p_i : prob_) {
1,551,214✔
110
    p_i /= integral_;
1,470,408✔
111
  }
112
}
80,806✔
113

114
//==============================================================================
115
// Discrete implementation
116
//==============================================================================
117

118
Discrete::Discrete(pugi::xml_node node) : di_(node)
26,038✔
119
{
120
  auto params = get_node_array<double>(node, "parameters");
26,038✔
121

122
  std::size_t n = params.size() / 2;
26,038✔
123

124
  x_.assign(params.begin(), params.begin() + n);
26,038✔
125
}
26,038✔
126

127
Discrete::Discrete(const double* x, const double* p, size_t n) : di_({p, n})
46,867✔
128
{
129

130
  x_.assign(x, x + n);
46,867✔
131
}
46,867✔
132

133
double Discrete::sample(uint64_t* seed) const
59,379,915✔
134
{
135
  return x_[di_.sample(seed)];
59,379,915✔
136
}
137

138
//==============================================================================
139
// Uniform implementation
140
//==============================================================================
141

142
Uniform::Uniform(pugi::xml_node node)
256✔
143
{
144
  auto params = get_node_array<double>(node, "parameters");
256✔
145
  if (params.size() != 2) {
256!
146
    fatal_error("Uniform distribution must have two "
×
147
                "parameters specified.");
148
  }
149

150
  a_ = params.at(0);
256✔
151
  b_ = params.at(1);
256✔
152
}
256✔
153

154
double Uniform::sample(uint64_t* seed) const
283,488✔
155
{
156
  return a_ + prn(seed) * (b_ - a_);
283,488✔
157
}
158

159
//==============================================================================
160
// PowerLaw implementation
161
//==============================================================================
162

163
PowerLaw::PowerLaw(pugi::xml_node node)
32✔
164
{
165
  auto params = get_node_array<double>(node, "parameters");
32✔
166
  if (params.size() != 3) {
32!
167
    fatal_error("PowerLaw distribution must have three "
×
168
                "parameters specified.");
169
  }
170

171
  const double a = params.at(0);
32✔
172
  const double b = params.at(1);
32✔
173
  const double n = params.at(2);
32✔
174

175
  offset_ = std::pow(a, n + 1);
32✔
176
  span_ = std::pow(b, n + 1) - offset_;
32✔
177
  ninv_ = 1 / (n + 1);
32✔
178
}
32✔
179

180
double PowerLaw::sample(uint64_t* seed) const
2,222✔
181
{
182
  return std::pow(offset_ + prn(seed) * span_, ninv_);
2,222✔
183
}
184

185
//==============================================================================
186
// Maxwell implementation
187
//==============================================================================
188

189
Maxwell::Maxwell(pugi::xml_node node)
64✔
190
{
191
  theta_ = std::stod(get_node_value(node, "parameters"));
64✔
192
}
64✔
193

194
double Maxwell::sample(uint64_t* seed) const
3,883✔
195
{
196
  return maxwell_spectrum(theta_, seed);
3,883✔
197
}
198

199
//==============================================================================
200
// Watt implementation
201
//==============================================================================
202

203
Watt::Watt(pugi::xml_node node)
96✔
204
{
205
  auto params = get_node_array<double>(node, "parameters");
96✔
206
  if (params.size() != 2)
96!
207
    openmc::fatal_error("Watt energy distribution must have two "
×
208
                        "parameters specified.");
209

210
  a_ = params.at(0);
96✔
211
  b_ = params.at(1);
96✔
212
}
96✔
213

214
double Watt::sample(uint64_t* seed) const
12,138,075✔
215
{
216
  return watt_spectrum(a_, b_, seed);
12,138,075✔
217
}
218

219
//==============================================================================
220
// Normal implementation
221
//==============================================================================
222
Normal::Normal(pugi::xml_node node)
×
223
{
224
  auto params = get_node_array<double>(node, "parameters");
×
225
  if (params.size() != 2) {
×
226
    openmc::fatal_error("Normal energy distribution must have two "
×
227
                        "parameters specified.");
228
  }
229

230
  mean_value_ = params.at(0);
×
231
  std_dev_ = params.at(1);
×
232
}
×
233

234
double Normal::sample(uint64_t* seed) const
×
235
{
236
  return normal_variate(mean_value_, std_dev_, seed);
×
237
}
238

239
//==============================================================================
240
// Tabular implementation
241
//==============================================================================
242

243
Tabular::Tabular(pugi::xml_node node)
7,320✔
244
{
245
  if (check_for_node(node, "interpolation")) {
7,320!
246
    std::string temp = get_node_value(node, "interpolation");
7,320✔
247
    if (temp == "histogram") {
7,320✔
248
      interp_ = Interpolation::histogram;
7,256✔
249
    } else if (temp == "linear-linear") {
64!
250
      interp_ = Interpolation::lin_lin;
64✔
NEW
251
    } else if (temp == "log-linear") {
×
NEW
252
      interp_ = Interpolation::log_lin;
×
NEW
253
    } else if (temp == "log-log") {
×
NEW
254
      interp_ = Interpolation::log_log;
×
255
    } else {
256
      openmc::fatal_error(
×
257
        "Unsupported interpolation type for distribution: " + temp);
×
258
    }
259
  } else {
7,320✔
260
    interp_ = Interpolation::histogram;
×
261
  }
262

263
  // Read and initialize tabular distribution. If number of parameters is odd,
264
  // add an extra zero for the 'p' array.
265
  auto params = get_node_array<double>(node, "parameters");
7,320✔
266
  if (params.size() % 2 != 0) {
7,320!
267
    params.push_back(0.0);
×
268
  }
269
  std::size_t n = params.size() / 2;
7,320✔
270
  const double* x = params.data();
7,320✔
271
  const double* p = x + n;
7,320✔
272
  init(x, p, n);
7,320✔
273
}
7,320✔
274

275
Tabular::Tabular(const double* x, const double* p, int n, Interpolation interp,
46,274,153✔
276
  const double* c)
46,274,153✔
277
  : interp_ {interp}
46,274,153✔
278
{
279
  init(x, p, n, c);
46,274,153✔
280
}
46,274,153✔
281

282
void Tabular::init(
46,281,473✔
283
  const double* x, const double* p, std::size_t n, const double* c)
284
{
285
  // Copy x/p arrays into vectors
286
  std::copy(x, x + n, std::back_inserter(x_));
46,281,473✔
287
  std::copy(p, p + n, std::back_inserter(p_));
46,281,473✔
288

289
  // Calculate cumulative distribution function
290
  if (c) {
46,281,473✔
291
    std::copy(c, c + n, std::back_inserter(c_));
46,274,153✔
292
  } else {
293
    c_.resize(n);
7,320✔
294
    c_[0] = 0.0;
7,320✔
295
    for (int i = 1; i < n; ++i) {
90,912✔
296
      if (interp_ == Interpolation::histogram) {
83,592✔
297
        c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]);
83,464✔
298
      } else if (interp_ == Interpolation::lin_lin) {
128!
299
        c_[i] = c_[i - 1] + 0.5 * (p_[i - 1] + p_[i]) * (x_[i] - x_[i - 1]);
128✔
NEW
300
      } else if (interp_ == Interpolation::log_lin) {
×
NEW
301
        double m = std::log(p_[i] / p_[i - 1]) / (x_[i] - x_[i - 1]);
×
NEW
302
        c_[i] = c_[i - 1] + p_[i - 1] * (x_[i] - x_[i - 1]) *
×
NEW
303
                              exprel(m * (x_[i] - x_[i - 1]));
×
NEW
304
      } else if (interp_ == Interpolation::log_log) {
×
NEW
305
        double m = std::log((x_[i] * p_[i]) / (x_[i - 1] * p_[i - 1])) /
×
NEW
306
                   std::log(x_[i] / x_[i - 1]);
×
NEW
307
        c_[i] = c_[i - 1] + x_[i - 1] * p_[i - 1] *
×
NEW
308
                              std::log(x_[i] / x_[i - 1]) *
×
NEW
309
                              exprel(m * std::log(x_[i] / x_[i - 1]));
×
310
      } else {
NEW
311
        UNREACHABLE();
×
312
      }
313
    }
314
  }
315

316
  // Normalize density and distribution functions. Note that we save the
317
  // integral of the distribution so that if it is used as part of another
318
  // distribution (e.g., Mixture), we know its relative strength.
319
  integral_ = c_[n - 1];
46,281,473✔
320
  for (int i = 0; i < n; ++i) {
665,571,341✔
321
    p_[i] = p_[i] / integral_;
619,289,868✔
322
    c_[i] = c_[i] / integral_;
619,289,868✔
323
  }
324
}
46,281,473✔
325

326
double Tabular::sample(uint64_t* seed) const
691,517,353✔
327
{
328
  // Sample value of CDF
329
  double c = prn(seed);
691,517,353✔
330

331
  // Find first CDF bin which is above the sampled value
332
  double c_i = c_[0];
691,517,353✔
333
  int i;
334
  std::size_t n = c_.size();
691,517,353✔
335
  for (i = 0; i < n - 1; ++i) {
2,147,483,647!
336
    if (c <= c_[i + 1])
2,147,483,647✔
337
      break;
691,517,353✔
338
    c_i = c_[i + 1];
1,910,077,901✔
339
  }
340

341
  // Determine bounding PDF values
342
  double x_i = x_[i];
691,517,353✔
343
  double p_i = p_[i];
691,517,353✔
344

345
  if (interp_ == Interpolation::histogram) {
691,517,353✔
346
    // Histogram interpolation
347
    if (p_i > 0.0) {
3,421,748!
348
      return x_i + (c - c_i) / p_i;
3,421,748✔
349
    } else {
350
      return x_i;
×
351
    }
352
  } else if (interp_ == Interpolation::lin_lin) {
688,095,605!
353
    // Linear-linear interpolation
354
    double x_i1 = x_[i + 1];
688,095,605✔
355
    double p_i1 = p_[i + 1];
688,095,605✔
356

357
    double m = (p_i1 - p_i) / (x_i1 - x_i);
688,095,605✔
358
    if (m == 0.0) {
688,095,605✔
359
      return x_i + (c - c_i) / p_i;
321,979,701✔
360
    } else {
361
      return x_i +
362
             (std::sqrt(std::max(0.0, p_i * p_i + 2 * m * (c - c_i))) - p_i) /
366,115,904✔
363
               m;
366,115,904✔
364
    }
NEW
365
  } else if (interp_ == Interpolation::log_lin) {
×
366
    // Log-linear interpolation
NEW
367
    double x_i1 = x_[i + 1];
×
NEW
368
    double p_i1 = p_[i + 1];
×
369

NEW
370
    double m = std::log(p_i1 / p_i) / (x_i1 - x_i);
×
NEW
371
    double f = (c - c_i) / p_i;
×
NEW
372
    return x_i + f * log1prel(m * f);
×
NEW
373
  } else if (interp_ == Interpolation::log_log) {
×
374
    // Log-Log interpolation
NEW
375
    double x_i1 = x_[i + 1];
×
NEW
376
    double p_i1 = p_[i + 1];
×
377

NEW
378
    double m = std::log((x_i1 * p_i1) / (x_i * p_i)) / std::log(x_i1 / x_i);
×
NEW
379
    double f = (c - c_i) / (p_i * x_i);
×
NEW
380
    return x_i * std::exp(f * log1prel(m * f));
×
381
  } else {
NEW
382
    UNREACHABLE();
×
383
  }
384
}
385

386
//==============================================================================
387
// Equiprobable implementation
388
//==============================================================================
389

390
double Equiprobable::sample(uint64_t* seed) const
×
391
{
392
  std::size_t n = x_.size();
×
393

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

397
  double xl = x_[i];
×
398
  double xr = x_[i + i];
×
399
  return xl + ((n - 1) * r - i) * (xr - xl);
×
400
}
401

402
//==============================================================================
403
// Mixture implementation
404
//==============================================================================
405

406
Mixture::Mixture(pugi::xml_node node)
48✔
407
{
408
  double cumsum = 0.0;
48✔
409
  for (pugi::xml_node pair : node.children("pair")) {
192✔
410
    // Check that required data exists
411
    if (!pair.attribute("probability"))
144!
412
      fatal_error("Mixture pair element does not have probability.");
×
413
    if (!pair.child("dist"))
144!
414
      fatal_error("Mixture pair element does not have a distribution.");
×
415

416
    // cummulative sum of probabilities
417
    double p = std::stod(pair.attribute("probability").value());
144✔
418

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

423
    distribution_.push_back(std::make_pair(cumsum, std::move(dist)));
144✔
424
  }
144✔
425

426
  // Save integral of distribution
427
  integral_ = cumsum;
48✔
428

429
  // Normalize cummulative probabilities to 1
430
  for (auto& pair : distribution_) {
192✔
431
    pair.first /= cumsum;
144✔
432
  }
433
}
48✔
434

435
double Mixture::sample(uint64_t* seed) const
3,564✔
436
{
437
  // Sample value of CDF
438
  const double p = prn(seed);
3,564✔
439

440
  // find matching distribution
441
  const auto it = std::lower_bound(distribution_.cbegin(), distribution_.cend(),
3,564✔
442
    p, [](const DistPair& pair, double p) { return pair.first < p; });
7,128✔
443

444
  // This should not happen. Catch it
445
  assert(it != distribution_.cend());
2,916!
446

447
  // Sample the chosen distribution
448
  return it->second->sample(seed);
7,128✔
449
}
450

451
//==============================================================================
452
// Helper function
453
//==============================================================================
454

455
UPtrDist distribution_from_xml(pugi::xml_node node)
33,843✔
456
{
457
  if (!check_for_node(node, "type"))
33,843!
458
    openmc::fatal_error("Distribution type must be specified.");
×
459

460
  // Determine type of distribution
461
  std::string type = get_node_value(node, "type", true, true);
33,843✔
462

463
  // Allocate extension of Distribution
464
  UPtrDist dist;
33,843✔
465
  if (type == "uniform") {
33,843✔
466
    dist = UPtrDist {new Uniform(node)};
256✔
467
  } else if (type == "powerlaw") {
33,587✔
468
    dist = UPtrDist {new PowerLaw(node)};
32✔
469
  } else if (type == "maxwell") {
33,555✔
470
    dist = UPtrDist {new Maxwell(node)};
64✔
471
  } else if (type == "watt") {
33,491✔
472
    dist = UPtrDist {new Watt(node)};
96✔
473
  } else if (type == "normal") {
33,395!
474
    dist = UPtrDist {new Normal(node)};
×
475
  } else if (type == "discrete") {
33,395✔
476
    dist = UPtrDist {new Discrete(node)};
26,027✔
477
  } else if (type == "tabular") {
7,368✔
478
    dist = UPtrDist {new Tabular(node)};
7,320✔
479
  } else if (type == "mixture") {
48!
480
    dist = UPtrDist {new Mixture(node)};
48✔
481
  } else if (type == "muir") {
×
482
    openmc::fatal_error(
×
483
      "'muir' distributions are now specified using the openmc.stats.muir() "
484
      "function in Python. Please regenerate your XML files.");
485
  } else {
486
    openmc::fatal_error("Invalid distribution type: " + type);
×
487
  }
488
  return dist;
67,686✔
489
}
33,843✔
490

491
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc