leonarski_fandClaude Opus 4.8
9f23976e02
RfreeFlagsTest: determinism (pure function of the reflection), ~5% free fraction, Friedel/Bijvoet pairs never split, symmetry equivalents share a flag, resolution stratification (free set spans shells), and fraction=0 flags nothing. FrenchWilsonTest: strong reflections reduce to sqrt(I), weak/negative intensities get a positive amplitude, amplitudes are always finite and non-negative (incl. centric / epsilon>1 reflections in P4(3)2(1)2), and unusable sigma falls back to sqrt(max(I,0)). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
75 lines
3.0 KiB
C++
75 lines
3.0 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
|
|
#include <cmath>
|
|
#include <vector>
|
|
|
|
#include "../image_analysis/scale_merge/FrenchWilson.h"
|
|
|
|
namespace {
|
|
MergedReflection Refl(int h, int k, int l, float d, float I, float sigma) {
|
|
MergedReflection r;
|
|
r.h = h; r.k = k; r.l = l; r.d = d; r.I = I; r.sigma = sigma;
|
|
return r;
|
|
}
|
|
// A resolution-spread of ordinary reflections so a Wilson mean can be formed per shell.
|
|
std::vector<MergedReflection> Background() {
|
|
std::vector<MergedReflection> v;
|
|
for (int h = 1; h <= 12; ++h)
|
|
for (int k = 0; k <= 12; ++k)
|
|
for (int l = 0; l <= 12; ++l)
|
|
v.push_back(Refl(h, k, l, 40.0f / (1 + h * h + k * k + l * l), 800.0f, 20.0f));
|
|
return v;
|
|
}
|
|
}
|
|
|
|
TEST_CASE("French-Wilson: strong reflections reduce to sqrt(I)", "[french_wilson]") {
|
|
auto v = Background();
|
|
v.push_back(Refl(1, 0, 0, 25.0f, 40000.0f, 50.0f)); // I/sigma = 800, clearly strong
|
|
ApplyFrenchWilson(v, 1);
|
|
CHECK(v.back().F == Catch::Approx(std::sqrt(40000.0)).epsilon(0.02)); // ~200
|
|
CHECK(v.back().sigmaF >= 0.0f);
|
|
CHECK(std::isfinite(v.back().sigmaF));
|
|
}
|
|
|
|
TEST_CASE("French-Wilson: weak and negative intensities get a positive amplitude", "[french_wilson]") {
|
|
auto v = Background();
|
|
v.push_back(Refl(2, 0, 0, 20.0f, -40.0f, 50.0f)); // negative measured intensity
|
|
v.push_back(Refl(3, 0, 0, 15.0f, 10.0f, 50.0f)); // weak, I < sigma
|
|
ApplyFrenchWilson(v, 1);
|
|
const auto& neg = v[v.size() - 2];
|
|
const auto& weak = v.back();
|
|
CHECK(std::isfinite(neg.F));
|
|
CHECK(neg.F > 0.0f); // Bayesian estimate is positive (naive sqrt would give 0)
|
|
CHECK(std::isfinite(weak.F));
|
|
CHECK(weak.F > 0.0f);
|
|
}
|
|
|
|
TEST_CASE("French-Wilson: amplitudes are always finite and non-negative", "[french_wilson]") {
|
|
std::vector<MergedReflection> v;
|
|
// A deliberate mix: strong, weak, negative, tiny sigma, across resolution.
|
|
for (int i = 0; i < 300; ++i) {
|
|
const float d = 20.0f / (1 + 0.05f * i);
|
|
const float I = (i % 7 == 0) ? -30.0f : static_cast<float>((i % 50) * 40);
|
|
v.push_back(Refl(1 + i, 2, 3, d, I, 25.0f));
|
|
}
|
|
ApplyFrenchWilson(v, 96); // P4(3)2(1)2 (has centric reflections + epsilon>1 axes)
|
|
for (const auto& r : v) {
|
|
CHECK(std::isfinite(r.F));
|
|
CHECK(r.F >= 0.0f);
|
|
CHECK(std::isfinite(r.sigmaF));
|
|
CHECK(r.sigmaF >= 0.0f);
|
|
}
|
|
}
|
|
|
|
TEST_CASE("French-Wilson: unusable sigma falls back to sqrt(max(I,0))", "[french_wilson]") {
|
|
auto v = Background();
|
|
v.push_back(Refl(4, 0, 0, 12.0f, 144.0f, NAN)); // no sigma
|
|
v.push_back(Refl(5, 0, 0, 11.0f, -5.0f, NAN)); // no sigma, negative I
|
|
ApplyFrenchWilson(v, 1);
|
|
CHECK(v[v.size() - 2].F == Catch::Approx(12.0f)); // sqrt(144)
|
|
CHECK(v.back().F == Catch::Approx(0.0f)); // sqrt(max(-5,0))
|
|
}
|