Average equivalent q-vectors in explicit structure factor - #479
Conversation
StructureFactorPBC and StructureFactorIPBC were adding each direction's intensity individually to SamplingPolicy, bypassing the averaging of equivalent lattice vectors (e.g. [1,0,0], [0,1,0], [0,0,1]) within each configuration. This adds averageByMagnitude() to group and average intensities by |q| magnitude before accumulating across configurations, matching the approach used in pripps (Rust's average_duplicates). This also removes `#pragma omp critical` from the hot loop by writing to pre-allocated vector indices in the parallel section. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
There was a problem hiding this comment.
Pull request overview
This PR updates the explicit structure factor sampling to average intensities for equivalent q-vectors (same |q|) within each configuration before accumulating into the sampler, and it removes an OpenMP critical section by writing per-iteration results into a preallocated array.
Changes:
- Added
averageByMagnitude()helper to group(q_norm, intensity)pairs by rounded |q| and compute a per-bin mean. - Refactored
StructureFactorPBC::sample()andStructureFactorIPBC::sample()to compute intensities in parallel into a pre-sized vector, then average-by-|q| and calladdSampling(). - Removed
#pragma omp criticalfrom the parallel hot loop (map updates now happen after the parallel region).
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| @@ -372,16 +402,22 @@ class StructureFactorPBC : private TSamplingPolicy | |||
| template <typename Tscatterers> | |||
| void sample(const Tscatterers& scatterers, const Point& boxlength) | |||
| { | |||
There was a problem hiding this comment.
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.
| { | |
| { | |
| if (p_max <= 0) { | |
| return; | |
| } |
| const auto n = directions.size() * static_cast<size_t>(p_max); | ||
| std::vector<std::pair<T, T>> q_intensity(n); | ||
|
|
There was a problem hiding this comment.
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.
| template <std::floating_point T> | ||
| std::map<T, T> averageByMagnitude(const std::vector<std::pair<T, T>>& pairs, | ||
| T precision = T{10000}) | ||
| { |
There was a problem hiding this comment.
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.
| struct Accumulator | ||
| { | ||
| T sum = T{0}; | ||
| int count = 0; | ||
| }; |
There was a problem hiding this comment.
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.
Summary
StructureFactorPBCandStructureFactorIPBCwere adding each direction's intensity individually toSamplingPolicy, bypassing averaging of equivalent lattice vectors (e.g. [1,0,0], [0,1,0], [0,0,1]) within each configurationaverageByMagnitude()free function to group and average intensities by |q| magnitude before accumulating across configurations, matching the approach used in pripps (average_duplicatesin Rust)#pragma omp criticalfrom the parallel hot loop by writing to pre-allocated vector indicesTest plan
🤖 Generated with Claude Code