Skip to content

Average equivalent q-vectors in explicit structure factor #479

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We鈥檒l occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mlund merged 1 commit into from
Mar 10, 2026
Merged

Average equivalent q-vectors in explicit structure factor #479

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 51 additions & 10 deletions src/scatter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -330,6 +330,36 @@ template <std::floating_point T> class SamplingPolicy
}
};

/**
* @brief Average values with equivalent keys (same rounded magnitude).
*
* Groups (key, value) pairs by rounded key and computes the mean value for each group.
* Used to average intensities from equivalent q-vectors with the same |q| magnitude.
*
* @see Equivalent to Rust's `average_duplicates` in pripps/src/explicit.rs
*/
template <std::floating_point T>
std::map<T, T> averageByMagnitude(const std::vector<std::pair<T, T>>& pairs,
T precision = T{10000})
{
Comment on lines +341 to +344

Copilot AI Mar 7, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

averageByMagnitude() duplicates the binning/rounding logic and hard-codes a default precision (10000) that is meant to match SamplingPolicy::precision. This creates a maintainability risk: if a different TSamplingPolicy is supplied (or SamplingPolicy's binning changes), duplicates may be grouped differently than addSampling() bins them, producing inconsistent results. Consider factoring the rounding/binning into a shared helper/constant, or plumb the precision/binner from the sampling policy into averageByMagnitude() so the grouping and accumulation always use the same binning rule.

Copilot uses AI. Check for mistakes.
struct Accumulator
{
T sum = T{0};
int count = 0;
};
Comment on lines +345 to +349

Copilot AI Mar 7, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Accumulator::count is an int, but it is counting elements from a std::vector and later cast to T. Using std::size_t (or at least an unsigned integral type) would better match container sizes and avoid potential overflow if pairs ever grows large.

Copilot uses AI. Check for mistakes.
std::map<T, Accumulator> bins;
for (const auto& [key, value] : pairs) {
const T rounded = std::round(key * precision) / precision;
bins[rounded].sum += value;
bins[rounded].count++;
}
std::map<T, T> result;
for (const auto& [key, acc] : bins) {
result[key] = acc.sum / static_cast<T>(acc.count);
}
return result;
}

/**
* @brief Calculate scattering intensity using explicit q averaging.
*
Expand Down Expand Up @@ -372,16 +402,22 @@ class StructureFactorPBC : private TSamplingPolicy
template <typename Tscatterers>
void sample(const Tscatterers& scatterers, const Point& boxlength)
{

Copilot AI Mar 7, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

p_max is an int coming from config (see analysis.cpp), and here it is cast to size_t to compute n. If p_max is negative, this will underflow to a huge size_t and attempt an enormous allocation even though the for (int p=1; p<=p_max; ++p) loop would do zero iterations. Consider validating q_multiplier/p_max > 0 in the constructor (throw) or early-return in sample() before computing n/allocating.

Suggested change
{
{
if (p_max <= 0) {
return;
}

Copilot uses AI. Check for mistakes.
const auto n = directions.size() * static_cast<size_t>(p_max);
std::vector<std::pair<T, T>> q_intensity(n);

#pragma omp parallel for collapse(2) default(shared)
for (size_t i = 0; i < directions.size(); ++i) { // openmp req. tradional loop
for (int p = 1; p <= p_max; ++p) { // loop over multiples of q
for (size_t i = 0; i < directions.size(); ++i) {
for (int p = 1; p <= p_max; ++p) {
const Point q =
2.0 * pc::pi * p * directions[i].cwiseQuotient(boxlength); // scattering vector
const auto intensity = calculateIntensity(scatterers, q);
#pragma omp critical // avoid race conditions when updating the map
addSampling(q.norm(), intensity, 1.0);
2.0 * pc::pi * p * directions[i].cwiseQuotient(boxlength);
q_intensity[i * static_cast<size_t>(p_max) + static_cast<size_t>(p - 1)] =
{static_cast<T>(q.norm()), calculateIntensity(scatterers, q)};
}
}

for (const auto& [q, intensity] : averageByMagnitude(q_intensity)) {
addSampling(q, intensity);
}
}

template <typename Tscatterers>
Expand Down Expand Up @@ -438,6 +474,9 @@ class StructureFactorIPBC : private TSamplingPolicy
template <typename Tscatterers>
void sample(const Tscatterers& scatterers, const Point& boxlength)
{
const auto n = directions.size() * static_cast<size_t>(p_max);
std::vector<std::pair<T, T>> q_intensity(n);

Comment on lines +477 to +479

Copilot AI Mar 7, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same issue as in StructureFactorPBC::sample(): p_max is an int and may be negative from configuration, but it is cast to size_t to compute n and size q_intensity. A negative p_max will underflow and can trigger an enormous allocation. Add validation (e.g., reject p_max <= 0) before casting/allocating.

Copilot uses AI. Check for mistakes.
// https://gcc.gnu.org/gcc-9/porting_to.html#ompdatasharing
// #pragma omp parallel for collapse(2) default(none) shared(directions, p_max, scatterers,
// boxlength)
Expand All @@ -462,18 +501,20 @@ class StructureFactorIPBC : private TSamplingPolicy
sum_f_cos += f * product;
sum_f_squared += f * f;
}
// collect average, `norm()` gives the scattering vector length
const T ipbc_factor =
std::pow(2, directions[i].count()); // 2 ^ number of non-zero elements
T intensity = T{0};
if (sum_f_squared != T{0}) {
intensity = (sum_f_cos * sum_f_cos) / sum_f_squared * ipbc_factor;
}
#pragma omp critical
// avoid race conditions when updating the map
addSampling(q.norm(), intensity, 1.0);
q_intensity[i * static_cast<size_t>(p_max) + static_cast<size_t>(p - 1)] =
{static_cast<T>(q_norm), intensity};
}
}

for (const auto& [q, intensity] : averageByMagnitude(q_intensity)) {
addSampling(q, intensity);
}
}

int getQMultiplier() { return p_max; }
Expand Down