ScaleAndMerge: Increase function tolerance for faster stop

common/CrystalLattice.cpp

@@ -139,35 +139,6 @@ void CrystalLattice::ReorderABEqual() {
 }
 
 void CrystalLattice::ReorderMonoclinic() {
-    float alpha = angle_deg(vec[1], vec[2]);
-    float beta = angle_deg(vec[0], vec[2]);
-    float gamma = angle_deg(vec[0], vec[1]);
-
-    float da = std::abs(alpha - 90.0f);
-    float db = std::abs(beta - 90.0f);
-    float dg = std::abs(gamma - 90.0f);
-
-    Coord a = vec[0];
-    Coord b = vec[1];
-    Coord c = vec[2];
-
-    if (da > db && da > dg) {
-        // alpha most different -> [b, c, a]
-        vec[0] = b;
-        vec[1] = a;
-        vec[2] = c;
-    } else if (dg > db && dg > da) {
-        // gamma most different -> [c, a, b]
-        vec[0] = a;
-        vec[1] = c;
-        vec[2] = b;
-    } // else beta most different -> keep [a, b, c]
-
-    if (vec[0].Length() > vec[2].Length())
-        std::swap(vec[0], vec[2]);
-
-    FixHandeness();
-
     // Enforce obtuse beta (>= 90°). Beta is the angle between a and c.
     // Flip signs of a and b simultaneously to keep handedness and lengths unchanged,
     // which maps beta -> 180° - beta.

image_analysis/IndexAndRefine.cpp

@@ -82,6 +82,7 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data
             .crystal_system = sym_result.system
         };
         outcome.lattice_candidate = sym_result.conventional;
+
     }
 
     return outcome;
@@ -125,6 +126,9 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In
 
     outcome.lattice_candidate = data.latt;
 
+    if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Monoclinic)
+        outcome.lattice_candidate->ReorderMonoclinic();
+
     if (outcome.beam_center_updated) {
         msg.beam_corr_x = data.beam_corr_x;
         msg.beam_corr_y = data.beam_corr_y;

image_analysis/RotationIndexer.cpp

@@ -76,8 +76,8 @@ void RotationIndexer::TryIndex() {
     if (!indexer_result.lattice.empty()) {
         // Find lattice type
         search_result_ = LatticeSearch(indexer_result.lattice[0]);
-        // Run refinement
 
+        // Run refinement
         DiffractionExperiment experiment_copy(experiment);
         XtalOptimizerData data{
             .geom = experiment_copy.GetDiffractionGeometry(),
@@ -95,6 +95,8 @@ void RotationIndexer::TryIndex() {
         if (data.crystal_system == gemmi::CrystalSystem::Trigonal)
             data.crystal_system = gemmi::CrystalSystem::Hexagonal;
 
+        if (data.crystal_system == gemmi::CrystalSystem::Monoclinic)
+            data.latt.ReorderMonoclinic();
         if (XtalOptimizer(data, v_sel)) {
             indexed_lattice = data.latt;
             updated_geom_ = data.geom;
@@ -132,6 +134,7 @@ std::optional<RotationIndexerResult> RotationIndexer::ProcessImage(int64_t image
     }
     if (!indexed_lattice)
         return {};
+
     return RotationIndexerResult{
         .lattice = indexed_lattice.value(),
         .search_result = search_result_,

image_analysis/scale_merge/ScaleAndMerge.cpp

@@ -388,6 +388,7 @@ namespace {
         options.max_num_iterations = opt.max_num_iterations;
         options.max_solver_time_in_seconds = opt.max_solver_time_s;
         options.num_threads = static_cast<int>(hw);
+        options.function_tolerance = 1e-4;
 
         ceres::Solver::Summary summary;
         ceres::Solve(options, &problem, &summary);

tests/CrystalLatticeTest.cpp

@@ -67,10 +67,8 @@ TEST_CASE("CrystalLattice_Sort") {
 
 TEST_CASE("CrystalLattice_ReorderMonoclinic") {
     std::vector<CrystalLattice> latt = {
-        {85,70,60, 85, 90, 90},
-        {60,70,85, 90, 90, 85},
-        {60,85,70, 90, 85, 90},
-        {70,60,85, 90, 90, 85}
+        {60, 85, 70, 90, 70, 90},
+        {60, 85, 70, 90, 110, 90},
     };
     for (const auto &l_in :latt) {
         CrystalLattice l = l_in;
@@ -79,7 +77,7 @@ TEST_CASE("CrystalLattice_ReorderMonoclinic") {
         CHECK(l.Vec0().Length() == Catch::Approx(60));
         CHECK(l.Vec1().Length() == Catch::Approx(85));
         CHECK(l.Vec2().Length() == Catch::Approx(70));
-        CHECK(l.GetUnitCell().beta == Catch::Approx(95.0));
+        CHECK(l.GetUnitCell().beta == Catch::Approx(110.0));
     }
 }
 

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();