rugnux: add opt-in stills partiality (--still-partiality)

Serial stills currently treat every reflection as a full (partiality hardcoded to 1).
Add an opt-in Gaussian excitation-error partiality set at prediction time (CPU + CUDA):
p = exp(-dist_ewald^2 / (2*sigma^2)), sigma^2 = profile_radius^2 + (bandwidth_sigma*|recip_z|)^2,
with sigma = the per-image profile radius (ewald_dist_cutoff/2), so an edge-of-acceptance
reflection keeps p ~ exp(-2). Off by default; the merge weight (~p^2) then down-weights
far-from-Ewald partials instead of trusting them as fulls.

Validated: helps medium/strong stills (LOV R-free 0.336->0.329, lyso8 0.433->0.410, lowers
the systematic error-model b in both) but HURTS weak OCP (dividing by a small, uncertain p
amplifies orientation error -> high-res noise, resolution collapse), so it is left opt-in.
A static forward p explains only ~6-10% of the partiality scatter; the full win needs per-image
post-refinement (future work), for which this is the prediction-side groundwork.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

common/BraggIntegrationSettings.cpp

@@ -96,3 +96,12 @@ float BraggIntegrationSettings::GetDMinLimit_A() const {
 float BraggIntegrationSettings::GetMinimumSigmaInRegardsToI() const {
     return minimum_sigma_in_regards_to_i;
 }
+
+BraggIntegrationSettings &BraggIntegrationSettings::StillPartiality(bool input) {
+    still_partiality = input;
+    return *this;
+}
+
+bool BraggIntegrationSettings::GetStillPartiality() const {
+    return still_partiality;
+}

common/BraggIntegrationSettings.h

@@ -20,6 +20,7 @@ class BraggIntegrationSettings {
     float d_min_limit_A = 1.0;
     std::optional<float> fixed_profile_radius;
     float minimum_sigma_in_regards_to_i = 0.02;
+    bool still_partiality = false;   // experimental stills excitation-error partiality (rugnux --still-partiality)
 
 public:
     BraggIntegrationSettings& R1(float input);
@@ -28,6 +29,7 @@ public:
     BraggIntegrationSettings& DMinLimit_A(float input);
     BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional<float> input);
     BraggIntegrationSettings& Integrator(IntegratorMode input);
+    BraggIntegrationSettings& StillPartiality(bool input);
 
 
     [[nodiscard]] IntegratorMode GetIntegrator() const;
@@ -38,4 +40,5 @@ public:
     [[nodiscard]] float GetDMinLimit_A() const;
 
     [[nodiscard]] float GetMinimumSigmaInRegardsToI() const;
+    [[nodiscard]] bool GetStillPartiality() const;
 };

image_analysis/IndexAndRefine.cpp

@@ -324,7 +324,11 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
         .wedge_deg = std::fabs(wedge_deg),
         .mosaicity_deg = std::fabs(mos_deg),
         // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged.
-        .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f
+        .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f,
+        // Experimental stills partiality (off by default): sigma = ewald_dist_cutoff/2 = the per-image
+        // profile radius, so a reflection at the acceptance edge (dist_ewald ~ 2*sigma) keeps p ~ exp(-2).
+        .still_partiality = experiment.GetBraggIntegrationSettings().GetStillPartiality(),
+        .profile_radius_recipA = ewald_dist_cutoff * 0.5f
     };
 
     // Predict, then integrate with the selected integrator (box-sum or profile-fit).

image_analysis/bragg_prediction/BraggPrediction.cpp

@@ -156,6 +156,13 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                         continue;
 
                     float d = 1.0f / sqrtf(recip_sq);
+                    float partiality = 1.0f;
+                    if (settings.still_partiality && settings.profile_radius_recipA > 0.0f) {
+                        const float sig_bw = settings.bandwidth_sigma * std::fabs(recip_z);
+                        const float sigma2 = settings.profile_radius_recipA * settings.profile_radius_recipA
+                                             + sig_bw * sig_bw;
+                        partiality = std::exp(-0.5f * dist_ewald_sphere * dist_ewald_sphere / sigma2);
+                    }
                     reflections[i] = Reflection{
                         .h = h,
                         .k = k,
@@ -168,7 +175,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                         .d = d,
                         .dist_ewald = dist_ewald_sphere,
                         .rlp = 1.0,
-                        .partiality = 1.0,
+                        .partiality = partiality,
                         .zeta = 1.0,
                         .image_scale_corr = 1.0
                     };

image_analysis/bragg_prediction/BraggPrediction.h

@@ -21,9 +21,15 @@ struct BraggPredictionSettings {
     // Relative X-ray bandwidth Δλ/λ expressed as a Gaussian sigma (0 = monochromatic).
     // When > 0 the Ewald-shell acceptance is thickened radially per reflection by
     // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d²), so the 1/d² pink-beam smear no
-    // longer clips high-resolution reflections. Must stay the last member so existing
-    // designated initializers remain valid.
+    // longer clips high-resolution reflections.
     float bandwidth_sigma = 0.0f;
+    // Experimental stills partiality (rugnux --still-partiality). When still_partiality is set and
+    // profile_radius_recipA > 0, each reflection gets a Gaussian excitation-error partiality
+    // p = exp(-dist_ewald^2 / (2*sigma^2)), sigma^2 = profile_radius_recipA^2 + (bandwidth_sigma*|recip_z|)^2,
+    // instead of the fixed 1.0 (full). Off by default. Keep these two as the trailing members so the
+    // existing designated initializers (which stop at bandwidth_sigma) remain valid.
+    bool still_partiality = false;
+    float profile_radius_recipA = 0.0f;
 };
 
 class BraggPrediction {

image_analysis/bragg_prediction/BraggPredictionGPU.cu

@@ -127,7 +127,13 @@ namespace {
         out.d = 1.0f / sqrtf(recip_sq);
         out.dist_ewald = dist_ewald;
         out.rlp = 1.0f;
-        out.partiality = 1.0f;
+        float partiality = 1.0f;
+        if (C.still_partiality && C.profile_radius_recipA > 0.0f) {
+            const float sig_bw = C.bandwidth_sigma * fabsf(recip_z);
+            const float sigma2 = C.profile_radius_recipA * C.profile_radius_recipA + sig_bw * sig_bw;
+            partiality = expf(-0.5f * dist_ewald * dist_ewald / sigma2);
+        }
+        out.partiality = partiality;
         out.zeta = 1.0f;
         out.image_scale_corr = 1.0f;
         return true;
@@ -158,7 +164,9 @@ namespace {
                                           float high_res_A,
                                           float ewald_dist_cutoff,
                                           char centering,
-                                          float bandwidth_sigma) {
+                                          float bandwidth_sigma,
+                                          bool still_partiality,
+                                          float profile_radius_recipA) {
         KernelConsts kc{};
         auto geom = experiment.GetDiffractionGeometry();
         kc.det_width_pxl = static_cast<float>(experiment.GetXPixelsNum());
@@ -171,6 +179,8 @@ namespace {
         kc.one_over_wavelength = 1.0f / geom.GetWavelength_A();
         kc.ewald_cutoff = ewald_dist_cutoff;
         kc.bandwidth_sigma = bandwidth_sigma;
+        kc.still_partiality = still_partiality;
+        kc.profile_radius_recipA = profile_radius_recipA;
         kc.Astar = lattice.Astar();
         kc.Bstar = lattice.Bstar();
         kc.Cstar = lattice.Cstar();
@@ -193,7 +203,8 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment,
                              const BraggPredictionSettings &settings) {
     // Build constants on host
     KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff,
-                                        settings.centering, settings.bandwidth_sigma);
+                                        settings.centering, settings.bandwidth_sigma,
+                                        settings.still_partiality, settings.profile_radius_recipA);
     cudaMemcpyAsync(dK, &hK, sizeof(KernelConsts), cudaMemcpyHostToDevice, stream);
     cudaMemsetAsync(d_count, 0, sizeof(int), stream);
 

image_analysis/bragg_prediction/BraggPredictionGPU.h

@@ -19,6 +19,8 @@ struct KernelConsts {
     float one_over_dmax_sq;
     float ewald_cutoff;
     float bandwidth_sigma;   // relative Δλ/λ (sigma); 0 = monochromatic
+    bool still_partiality;        // experimental stills excitation-error partiality (off => p = 1, fulls)
+    float profile_radius_recipA;  // Gaussian sigma [1/A] for the stills partiality
     Coord Astar, Bstar, Cstar, S0;
     float rot[9];
     char centering;

rugnux/rugnux_cli.cpp

@@ -109,6 +109,7 @@ void print_usage() {
     std::cout << "   --bandwidth <num>               Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl;
     std::cout << "   --integration-radius <r>        Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl;
     std::cout << "   --integrator <txt>              Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl;
+    std::cout << "   --still-partiality              Experimental: weight stills reflections by a Gaussian excitation-error partiality exp(-dist_ewald^2/2sigma^2) instead of treating each as a full" << std::endl;
     std::cout << "   -q, --azim-q-spacing <num>      Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl;
     std::cout << "   --azim-min-q <num>              Azimuthal-integration minimum Q (1/A)" << std::endl;
     std::cout << "   --azim-max-q <num>              Azimuthal-integration maximum Q (1/A)" << std::endl;
@@ -153,6 +154,7 @@ enum {
     OPT_MODEL,
     OPT_DUMP_OBSERVATIONS,
     OPT_INTEGRATOR,
+    OPT_STILL_PARTIALITY,
     OPT_SCALE_FULLS,
     OPT_CAPTURE_UNCERTAINTY,
     OPT_MIN_CAPTURED_FRACTION,
@@ -246,6 +248,7 @@ static option long_options[] = {
     {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
     {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
     {"integrator", required_argument, nullptr, OPT_INTEGRATOR},
+    {"still-partiality", no_argument, nullptr, OPT_STILL_PARTIALITY},
     {"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS},
     {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS},
     {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF},
@@ -489,6 +492,7 @@ int main(int argc, char **argv) {
     std::optional<int> report_shell_count;                           // --resolution-shells
     std::optional<std::string> integration_radius_arg; // "r1" or "r1,r2,r3"
     std::optional<IntegratorMode> integrator_mode;     // --integrator boxsum|gaussian|empirical
+    bool still_partiality_flag = false;                // --still-partiality (experimental stills partiality)
     std::optional<double> outlier_reject_nsigma;       // merge per-observation outlier rejection
     std::optional<double> delta_cchalf_nsigma;         // per-crystal CC1/2-delta rejection
 
@@ -726,6 +730,9 @@ int main(int argc, char **argv) {
             else if (strcmp(optarg, "empirical") == 0) integrator_mode = IntegratorMode::ProfileEmpirical;
             else { logger.Error("--integrator expects boxsum|gaussian|empirical"); return 1; }
             break;
+        case OPT_STILL_PARTIALITY:
+            still_partiality_flag = true;
+            break;
         case OPT_REJECT_OUTLIERS:
             outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", logger);
             break;
@@ -1331,6 +1338,13 @@ int main(int argc, char **argv) {
                     : "profile (empirical)");
     }
 
+    if (still_partiality_flag) {
+        BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
+        bis.StillPartiality(true);
+        experiment.ImportBraggIntegrationSettings(bis);
+        logger.Info("Stills partiality enabled (experimental Gaussian excitation-error weighting)");
+    }
+
     SpotFindingSettings spot_settings;
     spot_settings.enable = true;
     spot_settings.indexing = true;