Indexing: let the FFT search reach past 500 A, and stop clipping the peak background
fft_max_unit_cell_A was both the default and an enforced check_max, so 500 A was the longest basis vector the FFT could ever return: FFTIndexer sizes its projected histogram from that value and the transform's last usable bin IS that length. Of the PDB's 206950 X-ray entries, 1091 (0.53%) have an axis longer than that and were unindexable by construction. The accepted range now goes to 1200 A, which leaves 8. The DEFAULT is unchanged at 500 - the histogram is sized from the value in use, so nothing pays for the wider range unless a caller asks for it. The peak picker's running-mean background was truncated at the ends of the spectrum rather than slid inward, so a peak within bg_half (~15 A) of either end - which is exactly where the longest cells sit - was judged on a one-sided background, biasing its prominence by however much the spectrum sloped there. Keep the window a constant width and slide it. Both bounds stay monotonically non-decreasing in j, so the GPU kernel's running sum is still valid. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
@@ -44,10 +44,16 @@ IndexingSettings &IndexingSettings::Algorithm(IndexingAlgorithmEnum input) {
|
||||
return *this;
|
||||
}
|
||||
|
||||
// The accepted range, not the default (which stays 500 A - see IndexingSettings.h). The FFT can
|
||||
// only recover a basis vector up to this length, so the ceiling is exactly the longest cell the
|
||||
// indexer can ever find; at 500 it excluded 1091 of the PDB's 206950 X-ray entries (0.53%) outright.
|
||||
// 1200 A leaves 8. Nothing pays for the wider range: the histogram is sized from the VALUE in use,
|
||||
// and only a caller that asks for more - a given cell that needs it, or the long-axis rescue - gets
|
||||
// a longer transform.
|
||||
IndexingSettings &IndexingSettings::FFT_MaxUnitCell_A(float input) {
|
||||
check_finite("FFT indexing max unit cell (A)", input);
|
||||
check_min("FFT indexing max unit cell (A)", input, 50);
|
||||
check_max("FFT indexing max unit cell (A)", input, 500);
|
||||
check_max("FFT indexing max unit cell (A)", input, fft_max_unit_cell_limit_A);
|
||||
fft_max_unit_cell_A = input;
|
||||
return *this;
|
||||
}
|
||||
|
||||
@@ -26,6 +26,14 @@ class IndexingSettings {
|
||||
float max_angle_from_ewald_deg = 2.0;
|
||||
float unit_cell_dist_tolerance_vs_reference = 0.05; // relative
|
||||
static constexpr float unit_cell_angle_tolerance_deg = 5.0; // degree
|
||||
public:
|
||||
// The longest cell the FFT search can be asked to reach, and so the longest it can ever find:
|
||||
// FFTIndexer sizes its histogram from fft_max_unit_cell_A and the transform's last usable bin IS
|
||||
// that length. Callers that widen the bound (a given cell, the long-axis rescue) must clamp to
|
||||
// this rather than let the setter throw - a rescue that recovers an implausible axis must not
|
||||
// take the whole run down with it.
|
||||
static constexpr float fft_max_unit_cell_limit_A = 1200.0;
|
||||
private:
|
||||
int64_t indexing_threads = 4;
|
||||
// Threads splitting the candidate-cell refinement WITHIN one indexer call. 1 (the default) is the
|
||||
// right answer whenever indexers already run one per image across all workers; it is raised only
|
||||
|
||||
@@ -110,8 +110,13 @@ void FFTIndexerCPU::ExecuteFFT(const std::vector<Coord> &coord, size_t nspots) {
|
||||
double len = len_coeff * static_cast<double>(j);
|
||||
if (len <= static_cast<double>(min_length_A)) continue;
|
||||
|
||||
const int lo = std::max(0, j - bg_half);
|
||||
const int hi = std::min(out_len, j + bg_half + 1);
|
||||
// Slide the background window inward at the ends rather than truncating it. A peak within
|
||||
// bg_half bins of either end - which is where the LONGEST cells sit, the last usable bin
|
||||
// being max_length_A itself - otherwise gets its background from a one-sided window, and
|
||||
// the prominence it is judged on is biased by however much the spectrum slopes there.
|
||||
int lo = j - bg_half, hi = j + bg_half + 1;
|
||||
if (lo < 0) { hi = std::min(out_len, hi - lo); lo = 0; }
|
||||
if (hi > out_len) { lo = std::max(0, lo - (hi - out_len)); hi = out_len; }
|
||||
const double background = (pref[hi] - pref[lo]) / static_cast<double>(hi - lo);
|
||||
const double prominence = mag[j] - background;
|
||||
|
||||
|
||||
@@ -40,9 +40,18 @@ __global__ void calculate_fft_result(
|
||||
FFTResult result{.magnitude = 0.0f, .direction = i, .length = -1};
|
||||
|
||||
for (int j = 0; j < out_len; ++j) {
|
||||
const int want_hi = min(j + bg_half, out_len - 1);
|
||||
// Constant-width window, slid inward at the ends instead of truncated (see
|
||||
// FFTIndexerCPU): a peak within bg_half of either end - where the LONGEST cells sit -
|
||||
// otherwise gets a one-sided background and a biased prominence. Both bounds stay
|
||||
// monotonically non-decreasing in j, so the running sum below is still valid.
|
||||
int want_lo = j - bg_half;
|
||||
int want_hi = j + bg_half;
|
||||
if (want_lo < 0) { want_hi = min(out_len - 1, want_hi - want_lo); want_lo = 0; }
|
||||
if (want_hi > out_len - 1) {
|
||||
want_lo = max(0, want_lo - (want_hi - (out_len - 1)));
|
||||
want_hi = out_len - 1;
|
||||
}
|
||||
while (whi < want_hi) { ++whi; winsum += complex_abs(d_output[offset + whi]); }
|
||||
const int want_lo = max(0, j - bg_half);
|
||||
while (wlo < want_lo) { winsum -= complex_abs(d_output[offset + wlo]); ++wlo; }
|
||||
|
||||
const float len = len_coeff * static_cast<float>(j);
|
||||
|
||||
Reference in New Issue
Block a user