jfjoch_process: refined logger and variable definition

VERSION

@@ -1 +1 @@
-1.0.0-rc.134
+1.0.0-rc.133

image_analysis/scale_merge/ScaleAndMerge.cpp

@@ -18,6 +18,7 @@
 #include <vector>
 
 #include "../../common/ResolutionShells.h"
+#include "../../common/Logger.h"
 
 namespace {
     struct HKLKey {
@@ -248,10 +249,7 @@ namespace {
     void select_reflections_by_quasi_random(const std::vector<ObsRef> &obs,
         const ScaleMergeOptions &opt,
         std::vector<bool> &hkl_selected,
-        double Isigma_cutoff = 1.0,
-        int min_per_bin = 2500, // needs optimization
-        int max_per_bin = 80000,
-        int n_resolution_bins = 20) {
+        Logger &logger) {
         float stat_d_min = std::numeric_limits<float>::max();
         float stat_d_max = 0.0f;
 
@@ -269,7 +267,7 @@ namespace {
         int reflection_above_cutoff = 0;
 
         for (const auto &o: obs) {
-            if (o.r->I / o.r->sigma < Isigma_cutoff) {
+            if (o.r->I / o.r->sigma < opt.filter_sigma_cutoff) {
                 o.selected = false;
                 continue;
             }
@@ -287,16 +285,16 @@ namespace {
             }
         }
 
-        std::cout << "Reflections of I/sigma > " << Isigma_cutoff << " : " << reflection_above_cutoff << std::endl;
-        if (reflection_above_cutoff < min_per_bin * n_resolution_bins) {
-            std::cout << "No additional selection applied before scaling" << std::endl;
+        logger.Info("Reflections of I/sigma > {} : {}", opt.filter_sigma_cutoff, reflection_above_cutoff);
+        if (reflection_above_cutoff < opt.filter_min_per_bin * opt.filter_n_resolution_bins) {
+            logger.Info("No additional selection applied before scaling");
             return;
         }
 
         if (stat_d_min < stat_d_max && stat_d_min > 0.0f) {
             const float d_min_pad = stat_d_min * 0.999f;
             const float d_max_pad = stat_d_max * 1.001f;
-            ResolutionShells scaling_shells(d_min_pad, d_max_pad, n_resolution_bins);
+            ResolutionShells scaling_shells(d_min_pad, d_max_pad, opt.filter_n_resolution_bins);
 
             for (int h = 0; h < nhkl; ++h) {
                 const auto d = per_hkl[h].d;
@@ -315,11 +313,11 @@ namespace {
                 int obs_total = 0;
                 bool selected = true;
             };
-            std::vector<ShellAccum> shell_acc(n_resolution_bins);
+            std::vector<ShellAccum> shell_acc(opt.filter_n_resolution_bins);
 
             for (int h = 0; h < nhkl; ++h) {
                 auto &sa = shell_acc[per_hkl[h].shell_id];
-                if (sa.obs_unique  * n_operator  > min_per_bin * 1.2 || sa.obs_total > max_per_bin)
+                if (sa.obs_unique  * n_operator  > opt.filter_min_per_bin * 1.2 || sa.obs_total > opt.filter_max_per_bin)
                     hkl_selected[h] = false;
                 else
                     sa.obs_unique += 1;
@@ -328,15 +326,13 @@ namespace {
 
             const auto shell_min_res = scaling_shells.GetShellMinRes();
 
-            std::cout << "| d-mean | n_refl_uni | n_refl_tot | selection |" << std::endl;
-            for (int n=0; n < n_resolution_bins; ++n) {
-                if (shell_acc[n].obs_unique * n_operator < min_per_bin)
+            logger.Info("| d-mean | n_refl_uni | n_refl_tot | selection |");
+            for (int n=0; n < opt.filter_n_resolution_bins; ++n) {
+                if (shell_acc[n].obs_unique * n_operator < opt.filter_min_per_bin)
                     shell_acc[n].selected = false;
                 if (shell_min_res[n] <= 0.0f) continue;
-                std::cout << std::setw(8) << std::fixed << std::setprecision(3) << shell_min_res[n]
-                          << std::setw(12) << std::fixed << std::setprecision(0) << shell_acc[n].obs_unique
-                          << std::setw(12) << shell_acc[n].obs_total
-                          << "      " << shell_acc[n].selected << std::endl;
+                logger.Info("| {:6.3f} | {:10d} | {:10d} | {:9d} |",
+                    shell_min_res[n], shell_acc[n].obs_unique, shell_acc[n].obs_total, shell_acc[n].selected);
             }
 
             for (int h = 0; h < nhkl; ++h) {
@@ -355,7 +351,8 @@ namespace {
                bool rotation_crystallography,
                size_t nhkl,
                const std::vector<ObsRef> &obs,
-               bool selection) {
+               bool selection,
+               Logger &logger) {
         ceres::Problem problem;
 
         std::vector<double> Itrue(nhkl, 0.0);
@@ -383,7 +380,7 @@ namespace {
 
         double wedge = opt.wedge_deg.value_or(0.0);
 
-        select_reflections_by_quasi_random(obs, opt, hkl_selected);
+        if (selection) select_reflections_by_quasi_random(obs, opt, hkl_selected, logger);
 
         std::vector<bool> is_valid_hkl_slot(nhkl, false);
         for (const auto &o: obs) {
@@ -508,9 +505,9 @@ namespace {
         options.function_tolerance = 1e-4;
 
         ceres::Solver::Summary summary;
-        std::cout << "Now start the ceres-solver with residual blocks: " << problem.NumResidualBlocks() << std::endl;
+        logger.Info("Now start the ceres-solver with residual blocks: {}", problem.NumResidualBlocks());
         ceres::Solve(options, &problem, &summary);
-        std::cout << summary.FullReport() << std::endl;
+        logger.Info(summary.FullReport());
     }
 
     void merge(size_t nhkl, ScaleMergeResult &out, const std::vector<CorrectedObs> &corr_obs) {
@@ -815,6 +812,9 @@ namespace {
 
 ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection> > &observations,
                                                const ScaleMergeOptions &opt) {
+
+    Logger logger("ScaleAndMergeReflections");
+
     if (opt.image_cluster <= 0)
         throw std::invalid_argument("image_cluster must be positive");
 
@@ -857,8 +857,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
         }
     }
 
-    std::cout << "Now scale the reflections: " << nrefl << std::endl;
-    scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs, selection);
+    logger.Info("Now scale the reflections: {}", nrefl);
+    scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs, selection, logger);
 
     ScaleMergeResult out;
 

image_analysis/scale_merge/ScaleAndMerge.h

@@ -51,6 +51,10 @@ struct ScaleMergeOptions {
     bool refine_wedge = false;
 
     bool selection = true;
+    double filter_sigma_cutoff = 1.0;
+    int filter_min_per_bin = 10000; // needs optimization
+    int filter_max_per_bin = 80000;
+    int filter_n_resolution_bins = 20;
 
     enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed;
 };

tools/jfjoch_process.cpp

@@ -52,8 +52,8 @@ void print_usage(Logger &logger) {
     logger.Info("   -c<num>    Max spot count (default: 250)");
     logger.Info("   -W<txt>    HDF5 file with analysis results is written. 'l' or 'light' deactivates image-output");
     logger.Info("   -T<num>    Noise sigma level for spot finding (default: 3.0)");
-    logger.Info("   -a         Use all reflection for scaling without filtering (default: false)");
-    logger.Info("   -w         Reflection wedge in scaling (default: false)");
+    logger.Info("   -m<num>    Min unique reflections per bin used for scaling. Use all when no value is specified (default: 10000)");
+    logger.Info("   -w         Refine wedge in scaling (default: false)");
 }
 
 void trim_in_place(std::string& t) {
@@ -168,6 +168,7 @@ int main(int argc, char **argv) {
     bool write_output = false;
     bool write_output_noimage = false;
     bool filtering = true;
+    std::optional<int16_t> filter_min_per_bin;
     std::optional<int64_t> max_spot_count_override;
     float sigma_spot_finding = 3.0;
     std::optional<float> merging_threshold;
@@ -185,7 +186,7 @@ int main(int argc, char **argv) {
     }
 
     int opt;
-    while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:M::P:AD:C:T:W:aw")) != -1) {
+    while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:M::P:AD:C:T:W:m::w")) != -1) {
         switch (opt) {
             case 'o':
                 output_prefix = optarg;
@@ -247,8 +248,10 @@ int main(int argc, char **argv) {
             case 'x':
                 refine_beam_center = false;
                 break;
-            case 'a':
+            case 'm':
                 filtering = false;
+                if (optarg)
+                    filter_min_per_bin = atoi(optarg);
                 break;
             case 'w':
                 refine_wedge = true;
@@ -648,8 +651,15 @@ int main(int argc, char **argv) {
         scale_opts.max_solver_time_s = 240.0; // generous cutoff for now
         scale_opts.merge_friedel = !anomalous_mode;
         scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0);
-        scale_opts.selection = filtering;
         scale_opts.refine_wedge = refine_wedge;
+        if (filter_min_per_bin.has_value()) {
+            scale_opts.selection = true;
+            scale_opts.filter_min_per_bin = filter_min_per_bin.value();
+        } else {
+            scale_opts.selection = filtering;
+        }
+        if (rotation_indexing)
+            scale_opts.filter_min_per_bin = scale_opts.filter_min_per_bin * 0.25;
 
         const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();
 
@@ -666,6 +676,7 @@ int main(int argc, char **argv) {
         double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
 
         if (scale_result) print_statistics(logger, scale_result->statistics);
+//        if (scale_opts.wedge_deg.has_value()) logger.Info("Refined wedge: {:.3f}", scale_opts.wedge_deg.value());
 
         if (scale_result && !fixed_space_group) {
             logger.Info("Searching for space group from P1-merged reflections ...");
@@ -698,10 +709,10 @@ int main(int argc, char **argv) {
             logger.Info("");
 
             if (sg_search.best_space_group.has_value()) {
+                logger.Info("SG-search wall-clock time: {:.2f} s", sg_search_time);
                 if (sg_search.best_space_group->number != 0) {
                     if (sg_search.best_space_group->number != scale_opts.space_group->number) {
                         logger.Info("Re-running scaling in detected space group {}", sg_search.best_space_group->short_name());
-                        logger.Info("SG-search wall-clock time: {:.2f} s", sg_search_time);
 
                         scale_opts.space_group = *sg_search.best_space_group;