v1.0.0-rc.128 (#35)

This is an UNSTABLE release. The release has significant modifications and bug fixes, if things go wrong, it is better to revert to 1.0.0-rc.124.

* jfjoch_broker: Handle properly reuse of image buffer locations
* jfjoch_broker: Fix bug in counting idle slots
* jfjoch_broker: Force obtuse angle for monoclinic cells
* jfjoch_process: Change scaling refinement tolerance

Reviewed-on: #35
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>
Co-committed-by: Filip Leonarski <filip.leonarski@psi.ch>

tools/jfjoch_process.cpp

@@ -44,8 +44,68 @@ void print_usage(Logger &logger) {
     logger.Info("   -M         Scale and merge (refine mosaicity) and write scaled.hkl + image.dat");
     logger.Info("   -P<txt>    Partiality refinement fixed|rot|unity (default: fixed)");
     logger.Info("   -A         Anomalous mode (don't merge Friedel pairs)");
+    logger.Info("   -C<cell>   Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)");
 }
 
+void trim_in_place(std::string& t) {
+    size_t b = 0;
+    while (b < t.size() && std::isspace(static_cast<unsigned char>(t[b]))) b++;
+    size_t e = t.size();
+    while (e > b && std::isspace(static_cast<unsigned char>(t[e - 1]))) e--;
+    t = t.substr(b, e - b);
+};
+
+std::optional<UnitCell> parse_unit_cell_arg(const char* arg)  {
+    if (!arg)
+        return std::nullopt;
+
+    std::string s(arg);
+
+
+    trim_in_place(s);
+
+    if (s.size() >= 2 && ((s.front() == '"' && s.back() == '"') || (s.front() == '\'' && s.back() == '\''))) {
+        s = s.substr(1, s.size() - 2);
+        trim_in_place(s);
+    }
+
+    std::vector<std::string> parts;
+    parts.reserve(6);
+    size_t start = 0;
+    while (true) {
+        size_t pos = s.find(',', start);
+        if (pos == std::string::npos) {
+            parts.push_back(s.substr(start));
+            break;
+        }
+        parts.push_back(s.substr(start, pos - start));
+        start = pos + 1;
+    }
+
+    if (parts.size() != 6)
+        return std::nullopt;
+
+    auto parse_float_strict = [](const std::string& t, float& out) -> bool {
+        try {
+            size_t idx = 0;
+            out = std::stof(t, &idx);
+            return idx == t.size();
+        } catch (...) {
+            return false;
+        }
+    };
+
+    UnitCell uc{};
+    if (!parse_float_strict(parts[0], uc.a)) return std::nullopt;
+    if (!parse_float_strict(parts[1], uc.b)) return std::nullopt;
+    if (!parse_float_strict(parts[2], uc.c)) return std::nullopt;
+    if (!parse_float_strict(parts[3], uc.alpha)) return std::nullopt;
+    if (!parse_float_strict(parts[4], uc.beta)) return std::nullopt;
+    if (!parse_float_strict(parts[5], uc.gamma)) return std::nullopt;
+
+    return uc;
+};
+
 int main(int argc, char **argv) {
     RegisterHDF5Filter();
 
@@ -66,6 +126,7 @@ int main(int argc, char **argv) {
     bool run_scaling = false;
     bool anomalous_mode = false;
     std::optional<int> space_group_number;
+    std::optional<UnitCell> fixed_reference_unit_cell;
 
     ScaleMergeOptions::PartialityModel partiality_model = ScaleMergeOptions::PartialityModel::Fixed;
 
@@ -78,7 +139,7 @@ int main(int argc, char **argv) {
     }
 
     int opt;
-    while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MP:AD:")) != -1) {
+    while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MP:AD:C:")) != -1) {
         switch (opt) {
             case 'o':
                 output_prefix = optarg;
@@ -121,6 +182,18 @@ int main(int argc, char **argv) {
             case 'A':
                 anomalous_mode = true;
                 break;
+            case 'C': {
+                auto uc = parse_unit_cell_arg(optarg);
+                if (!uc.has_value()) {
+                    logger.Error("Invalid -C unit cell. Expected: -C\"a,b,c,alpha,beta,gamma\" (6 floats, comma-separated, no spaces). Got: {}", optarg ? optarg : "<null>");
+                    print_usage(logger);
+                    exit(EXIT_FAILURE);
+                }
+                fixed_reference_unit_cell = uc;
+                logger.Info("Fixed reference unit cell set: a={:.3f} b={:.3f} c={:.3f} alpha={:.3f} beta={:.3f} gamma={:.3f}",
+                            uc->a, uc->b, uc->c, uc->alpha, uc->beta, uc->gamma);
+                break;
+            }
             case 'P':
                 if (strcmp(optarg, "unity") == 0)
                     partiality_model = ScaleMergeOptions::PartialityModel::Unity;
@@ -188,6 +261,10 @@ int main(int argc, char **argv) {
     experiment.OverwriteExistingFiles(true);
     experiment.PolarizationFactor(0.99);
 
+    if (fixed_reference_unit_cell.has_value()) {
+        experiment.SetUnitCell(*fixed_reference_unit_cell);
+    }
+
     // Configure Indexing
     IndexingSettings indexing_settings;
     if (use_fft)
@@ -368,8 +445,22 @@ int main(int argc, char **argv) {
             }
 
             // Progress log
-            if (current_idx_offset > 0 && current_idx_offset % 100 == 0)
-                logger.Info("Processed {} / {} images", current_idx_offset, images_to_process);
+            if (current_idx_offset > 0 && current_idx_offset % 100 == 0) {
+                std::optional<float> indexing_rate;
+                {
+                    std::lock_guard<std::mutex> lock(plots_mutex);
+                    indexing_rate = plots.GetIndexingRate();
+                }
+
+                if (indexing_rate.has_value()) {
+                    logger.Info("Processed {} / {} images (indexing rate {:.1f}%)",
+                                current_idx_offset, images_to_process,
+                                indexing_rate.value() * 100.0f);
+                } else {
+                    logger.Info("Processed {} / {} images (indexing rate N/A)",
+                                current_idx_offset, images_to_process);
+                }
+            }
         }
 
         // Finalize per-thread indexing (if any per-thread aggregation is needed, though pool handles most)
@@ -583,6 +674,13 @@ int main(int argc, char **argv) {
         logger.Info("Indexing rate:   {:.2f}%", end_msg.indexing_rate.value() * 100.0);
     }
     if (rotation_indexer_ret.has_value()) {
+        auto latt = rotation_indexer_ret->lattice;
+        if (auto axis_ = experiment.GetGoniometer()) {
+            const float angle_deg = axis_->GetAngle_deg(0);
+            const auto rot = axis_->GetTransformationAngle(angle_deg);
+            latt = latt.Multiply(rot);
+        }
+
         auto vec0 = rotation_indexer_ret->lattice.Vec0();
         auto vec1 = rotation_indexer_ret->lattice.Vec1();
         auto vec2 = rotation_indexer_ret->lattice.Vec2();