Merge mtz-amp into rc181: FWT/DELFWT non-negative, phase turned by 180 degrees where negative

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

image_analysis/structure_refinement/ModelValidation.cpp

@@ -1835,15 +1835,24 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
             anom_by_hkl[hkl_key(hv.hkl)] = hv.value;
         for (size_t i = 0; i < terms.size(); i++) {
             // [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg
-            // returns (-180, 180] and a P column is not supposed to.
-            const double phi_wrapped = terms[i].phi * 180.0 / PI;
-            const auto phi_deg = static_cast<float>(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped);
+            // returns (-180, 180] and a P column is not supposed to - nor 360 itself, which the rounding to float
+            // can give.
+            const double deg = terms[i].phi * 180.0 / PI;
+            const auto wrap = [](double d) {
+                const auto f = static_cast<float>(d - 360.0 * std::floor(d / 360.0));
+                return f < 360.0f ? f : 0.0f;
+            };
+            const float phi_deg = wrap(deg);
+            // A negative coefficient is written as its magnitude on the phase turned by 180 degrees: an MTZ
+            // amplitude column is non-negative by convention, and not every reader accepts otherwise.
+            const float phwt = wrap(fwt[i] < 0 ? deg + 180.0 : deg);
+            const float phdelwt = wrap(delfwt[i] < 0 ? deg + 180.0 : deg);
             const float row[] = {static_cast<float>(terms[i].hkl[0]), static_cast<float>(terms[i].hkl[1]),
                                  static_cast<float>(terms[i].hkl[2]),
                                  static_cast<float>(terms[i].fo), static_cast<float>(terms[i].fc),
                                  phi_deg,
-                                 fwt[i], phi_deg,
-                                 delfwt[i], phi_deg,
+                                 std::abs(fwt[i]), phwt,
+                                 std::abs(delfwt[i]), phdelwt,
                                  static_cast<float>(sigma_a.weight[i].m),
                                  terms[i].free ? 0.0f : 1.0f};
             for (size_t c = 0; c < mtz.columns.size(); c++)

image_analysis/structure_refinement/Refinement.cpp

@@ -1150,6 +1150,12 @@ double BalancedADPWeight(Context &ctx, const RestraintSet &b_restraints, const s
     return nx > 0 ? std::sqrt(ng / nx) : 1.0;
 }
 
+// A phase in degrees as an MTZ 'P' column holds it: in [0, 360), also after the rounding to float.
+float MtzPhase(double deg) {
+    const auto f = static_cast<float>(deg - 360.0 * std::floor(deg / 360.0));
+    return f < 360.0f ? f : 0.0f;
+}
+
 // The model's map coefficients and its difference map, weighted by the refinement's own likelihood: m its figure
 // of merit and D its D (2mFo - DFc, mFo - DFc). With a twin law, Fo is detwinned by the model in proportion to its
 // own share of the twinned intensity, Fo |F_model(h)| / |F_twin(h)|, and Fc is the model's own |F_model(h)|.
@@ -1200,14 +1206,14 @@ void MakeMaps(Context &ctx, const std::vector<MergedReflection> &merged, double
         const double fwt = ctx.refl.centric[m] ? mw * fo : 2 * mw * fo - d * fc;
         const double delfwt = mw * fo - d * fc;
         const double phi = std::arg(fmodel[m]);
-        double phi_deg = phi * 180.0 / PI;
-        if (phi_deg < 0)
-            phi_deg += 360.0;
+        const double phi_deg = phi * 180.0 / PI;
+        // A negative coefficient is written as its magnitude on the phase turned by 180 degrees: an MTZ amplitude
+        // column is non-negative by convention, and not every reader accepts otherwise.
         const float row[] = {static_cast<float>(ctx.refl.hkl[m][0]), static_cast<float>(ctx.refl.hkl[m][1]),
                              static_cast<float>(ctx.refl.hkl[m][2]), static_cast<float>(fo), static_cast<float>(fc),
-                             static_cast<float>(phi_deg), static_cast<float>(fwt), static_cast<float>(phi_deg),
-                             static_cast<float>(delfwt), static_cast<float>(phi_deg), static_cast<float>(mw),
-                             ctx.refl.free[m] ? 0.0f : 1.0f};
+                             MtzPhase(phi_deg), static_cast<float>(std::abs(fwt)), MtzPhase(fwt < 0 ? phi_deg + 180.0 : phi_deg),
+                             static_cast<float>(std::abs(delfwt)), MtzPhase(delfwt < 0 ? phi_deg + 180.0 : phi_deg),
+                             static_cast<float>(mw), ctx.refl.free[m] ? 0.0f : 1.0f};
         for (size_t c = 0; c < mtz.columns.size(); c++)
             mtz.columns[c].values.push_back(row[c]);
         diff.v.push_back({ctx.refl.hkl[m], std::polar(static_cast<float>(delfwt), static_cast<float>(phi))});

tests/ModelValidationTest.cpp

@@ -613,6 +613,16 @@ TEST_CASE("ModelValidation_MapCoefficientsCarryTheAnomalousMap", "[ModelValidati
     REQUIRE(plain.ok);
     REQUIRE(plain.maps);
     CHECK(plain.maps->coefficients.columns.size() == 12);
+    // FWT and DELFWT as amplitude and phase: non-negative, their phases in [0, 360)
+    for (size_t c : {6, 8}) {
+        const auto &f = plain.maps->coefficients.columns[c], &phi = plain.maps->coefficients.columns[c + 1];
+        CHECK(f.label == (c == 6 ? "FWT" : "DELFWT"));
+        for (size_t i = 0; i < f.values.size(); i++) {
+            CHECK(f.values[i] >= 0.0f);
+            CHECK(phi.values[i] >= 0.0f);
+            CHECK(phi.values[i] < 360.0f);
+        }
+    }
 
     // With it: ANOM / PHANOM beside the others, missing on the reflections without a pair.
     for (size_t i = 0; i < obs.size(); i++)

tests/RefinementTest.cpp

@@ -834,6 +834,17 @@ TEST_CASE("Refinement_restrained_recovers_model", "[Refinement]") {
             rms += refined[i]->pos.dist_sq(ref[i]->pos);
         rms = std::sqrt(rms / ref.size());
         CHECK(rms < 0.15);
+        // the map coefficients as amplitude and phase: FWT and DELFWT non-negative, their phases in [0, 360)
+        REQUIRE(r.maps);
+        const auto &columns = r.maps->coefficients.columns;
+        REQUIRE(columns[6].label == "FWT");
+        REQUIRE(columns[8].label == "DELFWT");
+        for (size_t c : {6, 8})
+            for (size_t i = 0; i < columns[c].values.size(); i++) {
+                CHECK(columns[c].values[i] >= 0.0f);
+                CHECK(columns[c + 1].values[i] >= 0.0f);
+                CHECK(columns[c + 1].values[i] < 360.0f);
+            }
 
         // one trial at a time instead of three side by side (sharing a card's tables, mask and scale fit): the same
         // model, bit for bit