mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 63 Wiki Activity
Files
6df657fe62960e1092ffde188e1ef3f4010723be
Jungfraujoch/tools/jfjoch_scale.cpp
T
leonarski_fandClaude Opus 4.8 7d5fb15413 Unify still processing on one --force-still flag; scaling follows rotation indexing 
Fold the two overlapping knobs into one. Rotation-vs-stills is now a single decision
(IsRotationIndexing() = goniometer && rotation-indexing) that drives BOTH indexing and
scaling/merging - there is no independent scaling-model switch:

- ScalingSettings: drop combine_3d and force_still_processing (and Combine3D/GetCombine3D,
  ForceStillProcessing/GetForceStillProcessing). The scaling stage reads
  experiment.IsRotationIndexing() instead - rotation indexing implies RotationScaleMerge,
  no rotation indexing implies per-image ScaleOnTheFly.
- jfjoch_process: merge --process-as-stills and --force-still-processing into one
  --force-still (turns rotation processing off entirely: still indexing + still scaling).
  The rotation-specific scaling defaults (scale-fulls, smooth-G, capture uncertainty,
  outlier rejection) now key on rotation_indexing.
- jfjoch_scale: replace -P fixed|rot3d with --force-still; rotation is auto-detected from the
  goniometer (RotationIndexing set accordingly) so it matches jfjoch_process.
- jfjoch_viewer: the "process as stills" toggle already sets RotationIndexing; drop the now
  redundant Combine3D/ForceStillProcessing calls.

Rotation path unchanged (lyso 16.4/99.6%/87.3%); --force-still routes to ScaleOnTheFly; the
-R + --force-still conflict still errors. CUDA + non-CUDA + viewer all build.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
2026-07-03 12:20:45 +02:00

391 lines
18 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <iostream>
#include <vector>
#include <string>
#include <future>
#include <chrono>
#include <fstream>
#include <getopt.h>
#include "../reader/JFJochHDF5Reader.h"
#include "../common/Logger.h"
#include "../common/DiffractionExperiment.h"
#include "../common/time_utc.h"
#include "../common/print_license.h"
#include "../image_analysis/MXAnalysisWithoutFPGA.h"
#include "../writer/FileWriter.h"
#include "../image_analysis/IndexAndRefine.h"
#include "../common/JFJochReceiverPlots.h"
#include "../compression/JFJochCompressor.h"
#include "../image_analysis/LoadFCalcFromMtz.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
#include "../image_analysis/WriteReflections.h"
#include "../image_analysis/UpdateReflectionResolution.h"
void print_usage() {
std::cout << "Usage ./jfjoch_scale {<options>} <input.h5>" << std::endl;
std::cout << "Options:" << std::endl;
std::cout << " -o, --output-prefix <txt> Output file prefix (default: output)" << std::endl;
std::cout << " -N, --threads <num> Number of threads (default: 1)" << std::endl;
std::cout << " -s, --start-image <num> Start image number (default: 0)" << std::endl;
std::cout << " -e, --end-image <num> End image number (default: all)" << std::endl;
std::cout << " -v, --verbose Verbose output" << std::endl;
std::cout << "" << std::endl;
std::cout << " Scaling and merging" << std::endl;
std::cout << " -S, --space-group <num> Space group number" << std::endl;
std::cout << " --scaling-high-resolution <num> High resolution limit for scaling/merging (default: 0.0; no limit)"
<< std::endl;
std::cout << " --force-still Scale as independent stills (ScaleOnTheFly) even for a rotation dataset "
"(default: rotation data uses RotationScaleMerge)" << std::endl;
std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl;
std::cout << " -B, --refine-bfactor Refine per image B-factor" << std::endl;
std::cout << " --min-partiality <num> Minimum partiality to accept reflection (default: 0.02)" <<
std::endl;
std::cout << " --min-image-cc <num> Per-image CC limit in percent (default: no limit)" << std::endl;
std::cout << " --scaling-iterations <num> Number of scaling iterations with no reference data (default: 3)"
<< std::endl;
std::cout << " --scaling-output <txt> Output format for scaling results mtz|cif|txt (default: cif)" << std::endl;
std::cout << " -z, --reference-mtz <file> Reference MTZ file" << std::endl;
std::cout << " --reference-column <label> Reference MTZ column to use (default: auto - F-model, else IMEAN/I/...)" << std::endl;
}
enum {
OPT_MIN_PARTIALITY = 1000,
OPT_MIN_IMAGE_CC,
OPT_SCALING_ITERATIONS,
OPT_SCALING_OUTPUT,
OPT_SCALING_HIGH_RESOLUTION,
OPT_REFERENCE_COLUMN,
OPT_FORCE_STILL
};
static option long_options[] = {
{"verbose", no_argument, nullptr, 'v'},
{"output-prefix", required_argument, nullptr, 'o'},
{"threads", required_argument, nullptr, 'N'},
{"start-image", required_argument, nullptr, 's'},
{"end-image", required_argument, nullptr, 'e'},
{"reference-mtz", required_argument, nullptr, 'z'},
{"reference-column", required_argument, nullptr, OPT_REFERENCE_COLUMN},
{"space-group", required_argument, nullptr, 'S'},
{"force-still", no_argument, nullptr, OPT_FORCE_STILL},
{"anomalous", no_argument, nullptr, 'A'},
{"refine-bfactor", no_argument, nullptr, 'B'},
{"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION},
{"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY},
{"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC},
{"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS},
{"scaling-output", required_argument, nullptr, OPT_SCALING_OUTPUT},
{nullptr, 0, nullptr, 0}
};
int main(int argc, char **argv) {
RegisterHDF5Filter();
print_license("jfjoch_scale");
Logger logger("jfjoch_scale");
std::string input_file;
std::string output_prefix = "output";
int nthreads = 0;
int start_image = 0;
int end_image = -1; // -1 indicates process until end
bool verbose = false;
bool anomalous_mode = false;
std::optional<int64_t> space_group_number;
bool refine_bfactor = false;
std::string ref_mtz;
std::string ref_column;
double min_partiality = 0.02;
double min_image_cc = 0.0;
int64_t scaling_iter = 3;
std::optional<IntensityFormat> intensity_format; // --scaling-output override; default lives in ScalingSettings (mmCIF)
bool force_still = false; // --force-still: scale a rotation dataset as stills
std::optional<float> d_min_scale_merge;
if (argc == 1) {
print_usage();
exit(EXIT_FAILURE);
}
int opt;
int option_index = 0;
const char *short_opts = "vo:N:s:e:z:S:AB";
while ((opt = getopt_long(argc, argv, short_opts, long_options, &option_index)) != -1) {
switch (opt) {
case 'o':
output_prefix = optarg;
break;
case 'v':
verbose = true;
break;
case 'N':
nthreads = atoi(optarg);
break;
case 's':
start_image = atoi(optarg);
break;
case 'e':
end_image = atoi(optarg);
break;
case 'z':
ref_mtz = optarg;
break;
case OPT_REFERENCE_COLUMN:
ref_column = optarg;
break;
case 'S':
space_group_number = atoi(optarg);
break;
case OPT_FORCE_STILL:
force_still = true;
break;
case 'A':
anomalous_mode = true;
break;
case 'B':
refine_bfactor = true;
break;
case OPT_SCALING_HIGH_RESOLUTION:
d_min_scale_merge = atof(optarg);
logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_scale_merge.value());
break;
case OPT_MIN_PARTIALITY:
min_partiality = std::stod(optarg);
break;
case OPT_MIN_IMAGE_CC:
min_image_cc = std::stod(optarg);
break;
case OPT_SCALING_OUTPUT:
if (strcmp(optarg, "mtz") == 0) {
intensity_format = IntensityFormat::MTZ;
} else if (strcmp(optarg, "cif") == 0) {
intensity_format = IntensityFormat::mmCIF;
} else if (strcmp(optarg, "txt") == 0) {
intensity_format = IntensityFormat::Text;
} else {
logger.Error("Invalid intensity format: {}", optarg);
exit(EXIT_FAILURE);
}
break;
case OPT_SCALING_ITERATIONS:
scaling_iter = atoi(optarg);
if (scaling_iter <= 0) {
logger.Error("Invalid scaling iteration count: {}", scaling_iter);
exit(EXIT_FAILURE);
}
break;
}
}
// Validate space group number early
const gemmi::SpaceGroup *space_group = nullptr;
if (space_group_number.has_value()) {
space_group = gemmi::find_spacegroup_by_number(space_group_number.value());
if (!space_group) {
logger.Error("Unknown space group number {}", space_group_number.value());
exit(EXIT_FAILURE);
}
logger.Info("Using space group {} (number {})", space_group->hm, space_group_number.value());
}
std::vector<IntegrationOutcome> reflections;
std::vector<float> mosaicity;
if (optind >= argc) {
logger.Error("Input file not specified");
print_usage();
exit(EXIT_FAILURE);
}
logger.Verbose(verbose);
logger.Info("Loading reflections");
// 1. Read Input files
JFJochHDF5Reader reader;
std::shared_ptr<const JFJochReaderDataset> dataset;
for (int i = optind; i < argc; i++) {
input_file = argv[optind];
try {
reader.ReadFile(input_file);
auto tmp_dataset = reader.GetDataset();
if (i == optind)
dataset = tmp_dataset;
uint64_t total_images_in_file = reader.GetNumberOfImages();
if (end_image < 0 || end_image >= total_images_in_file)
end_image = total_images_in_file - 1;
size_t mosaicity_size = mosaicity.size();
mosaicity.resize(mosaicity_size + end_image - start_image + 1);
for (int i = 0; i < end_image - start_image + 1; i++)
mosaicity[mosaicity_size + i] = dataset->mosaicity_deg[start_image + i];
reflections = reader.ReadReflections(start_image, end_image);
logger.Info("Loaded dataset from {}", input_file);
} catch (const std::exception &e) {
logger.Error("Error reading input file: {}", e.what());
// Hard exit is only for the first file
if (i == optind)
exit(EXIT_FAILURE);
}
}
std::vector<MergedReflection> reference_data;
if (!ref_mtz.empty()) {
try {
const auto reference = LoadReferenceMtz(
ref_mtz, ref_column.empty() ? std::nullopt : std::optional<std::string>(ref_column));
reference_data = reference.reflections;
logger.Info("Loaded {} reference reflections from {} (column {}{}{})",
reference_data.size(), ref_mtz, reference.used_column,
reference.squared ? ", squared to intensity" : "",
reference.default_column ? ", auto-selected" : ", user-specified");
if (reference.d_max > 0.0)
logger.Info("Reference resolution range {:.2f} - {:.2f} A", reference.d_max, reference.d_min);
if (reference.cell.has_value())
logger.Info("Reference unit cell: a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
reference.cell->a, reference.cell->b, reference.cell->c,
reference.cell->alpha, reference.cell->beta, reference.cell->gamma);
if (!reference.space_group_name.empty())
logger.Info("Reference space group: {} (number {})",
reference.space_group_name, reference.space_group_number.value_or(0));
const auto data_sg = space_group_number.has_value()
? std::optional<int>(static_cast<int>(*space_group_number))
: (dataset->experiment.GetSpaceGroupNumber().has_value()
? std::optional<int>(static_cast<int>(*dataset->experiment.GetSpaceGroupNumber()))
: std::nullopt);
const auto warning = ReferenceConsistencyWarning(reference, dataset->experiment.GetUnitCell(), data_sg);
if (!warning.empty())
logger.Warning("{}", warning);
} catch (const std::exception &e) {
logger.Error("Error reading reference MTZ {}: {}", ref_mtz, e.what());
exit(EXIT_FAILURE);
}
}
logger.Info("Starting analysis of {} images", reflections.size());
// 2. Setup Experiment & Components
DiffractionExperiment experiment(dataset->experiment);
experiment.BitDepthImage(32).Compression(CompressionAlgorithm::BSHUF_LZ4);
experiment.FilePrefix(output_prefix);
experiment.Mode(DetectorMode::Standard); // Ensure full image analysis
experiment.PixelSigned(true);
experiment.OverwriteExistingFiles(true);
experiment.PolarizationFactor(0.99);
experiment.SetFileWriterFormat(FileWriterFormat::NXmxLegacy);
// Keep the space group read from the input file (written by the full pipeline) unless -S
// overrides it.
if (space_group_number.has_value())
experiment.SpaceGroupNumber(space_group_number);
experiment.NumTriggers(1);
// Rotation indexing drives rotation scaling/merge downstream (IsRotationIndexing()); a rotation
// (goniometer) dataset uses RotationScaleMerge unless --force-still asks for per-image stills scaling.
IndexingSettings indexing_settings;
indexing_settings.RotationIndexing(experiment.GetGoniometer().has_value() && !force_still);
experiment.ImportIndexingSettings(indexing_settings);
ScalingSettings scaling_settings;
if (d_min_scale_merge)
scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
scaling_settings.MergeFriedel(!anomalous_mode);
scaling_settings.RefineB(refine_bfactor);
scaling_settings.MinPartiality(min_partiality);
scaling_settings.MinCCForImage(min_image_cc);
if (intensity_format)
scaling_settings.FileFormat(intensity_format.value());
experiment.ImportScalingSettings(scaling_settings);
if (!experiment.GetUnitCell()) {
logger.Error("Experiment unit cell not found, cannot update reflection resolution");
exit(EXIT_FAILURE);
}
auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections);
logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images);
experiment.ImagesPerTrigger(refl_stats.n_images);
auto scale_start = std::chrono::steady_clock::now();
std::vector<MergedReflection> merged_reflections;
MergeStatistics merged_statistics;
double error_model_isa = 0.0;
// Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine -> merge
// (on the GPU when present). It does not support external-reference scaling, B-factor refinement or
// wedge refinement. Everything else (stills, reference scaling) uses ScaleOnTheFly + MergeOnTheFly.
const bool is_rotation = experiment.IsRotationIndexing();
if (is_rotation) {
if (!reference_data.empty() || experiment.GetScalingSettings().GetRefineB()
|| experiment.GetRefineRotationWedgeInScaling()
|| experiment.GetScalingSettings().GetRotationWedgeForScaling().has_value())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Rotation scaling/merging (RotationScaleMerge) does not support reference "
"scaling, B-factor refinement or wedge refinement");
RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(),
scaling_iter, 0.0f, nthreads, logger);
rsm.Ingest();
auto r = rsm.Run(false);
merged_reflections = std::move(r.merged);
merged_statistics = std::move(r.statistics);
error_model_isa = r.isa;
} else {
for (int i = 0; i < scaling_iter; i++) {
if (reference_data.empty())
ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads);
else
ScaleOnTheFly(experiment, reference_data).Scale(reflections, nthreads);
}
MergeOnTheFly merge_engine(experiment);
merge_engine.ReferenceCell(experiment.GetUnitCell());
for (size_t i = 0; i < reflections.size(); ++i)
merge_engine.AddImage(reflections[i], static_cast<int64_t>(i));
merged_reflections = merge_engine.ExportReflections();
merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data);
error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0;
}
logger.Info("Scale + merge completed in {:.2f} s ({} unique reflections)",
std::chrono::duration<double>(std::chrono::steady_clock::now() - scale_start).count(),
merged_reflections.size());
// Print resolution-shell statistics table
std::cout << merged_statistics;
// Space-group determination lives in the full jfjoch_process pipeline, where the lattice search
// and a consistent integration are available; this re-scaling tool only consumes a space group
// (from the file's /entry/sample/space_group_number or -S) and merges in it.
const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();
if (!fixed_space_group)
logger.Warning("No space group in the input file or on the command line - merged in P1. "
"Re-run jfjoch_process (which determines and stores the space group) or pass "
"-S to scale and merge in the correct symmetry.");
const auto twin_sg_number = experiment.GetSpaceGroupNumber();
const gemmi::SpaceGroup *twin_sg = twin_sg_number
? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr;
std::cout << std::endl << TwinningAnalysisToText(AnalyzeTwinning(merged_reflections, twin_sg)) << std::endl;
if (!output_prefix.empty())
WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics,
error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?",
output_prefix);
}
Reference in New Issue View Git Blame Copy Permalink