From de6646b95e2e0b4df240ed002e72c64f962828f3 Mon Sep 17 00:00:00 2001 From: takaba_k Date: Mon, 30 Mar 2026 16:21:34 +0200 Subject: [PATCH] jfjoch_process: refined logger and variable definition --- VERSION | 2 +- image_analysis/scale_merge/ScaleAndMerge.cpp | 48 ++++++++++---------- image_analysis/scale_merge/ScaleAndMerge.h | 4 ++ tools/jfjoch_process.cpp | 23 +++++++--- 4 files changed, 46 insertions(+), 31 deletions(-) diff --git a/VERSION b/VERSION index 31b74fae..91f33f63 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.0.0-rc.134 \ No newline at end of file +1.0.0-rc.133 \ No newline at end of file diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 57782ec6..d89d731d 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -18,6 +18,7 @@ #include #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 &obs, const ScaleMergeOptions &opt, std::vector &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::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 shell_acc(n_resolution_bins); + std::vector 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 &obs, - bool selection) { + bool selection, + Logger &logger) { ceres::Problem problem; std::vector 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 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 &corr_obs) { @@ -815,6 +812,9 @@ namespace { ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector > &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 Max spot count (default: 250)"); logger.Info(" -W HDF5 file with analysis results is written. 'l' or 'light' deactivates image-output"); logger.Info(" -T 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 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 filter_min_per_bin; std::optional max_spot_count_override; float sigma_spot_finding = 3.0; std::optional 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(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;