from coseda.logging_utils import log_print
import os
import glob
import subprocess
import multiprocessing
import configparser
import datetime
from coseda.initialize import get_parent_ini_file, shoutout
from coseda.export import hkl2mtz
from coseda.nexus.process import write_nxprocess_merging
import re
import csv
def _extract_h5_paths_from_stream(stream_path):
if not stream_path or not os.path.isfile(stream_path):
return []
base_dir = os.path.dirname(stream_path)
paths = []
seen = set()
try:
with open(stream_path, "r", encoding="utf-8", errors="replace") as handle:
for raw_line in handle:
line = raw_line.strip()
if not line:
continue
if line.startswith("Image filename:"):
token = line.split(":", 1)[1].strip()
elif line.startswith("Image filename ="):
token = line.split("=", 1)[1].strip()
else:
continue
if not token:
continue
token = token.strip("\"'")
token = token.split()[0].strip("\"'")
if token.startswith("file://"):
token = token[7:]
if "//" in token:
token = token.split("//", 1)[0]
if not os.path.isabs(token):
token = os.path.abspath(os.path.join(base_dir, token))
if token in seen:
continue
seen.add(token)
paths.append(token)
except Exception:
return []
return paths
[docs]
def run_merging(input_path):
"""
Read merging parameters from the INI file and execute partialator merging.
If 'merging_threads' is not set, defaults to the number of available CPU cores.
"""
# Determine the INI file path
if isinstance(input_path, str) and input_path.lower().endswith('.ini'):
configfile_path = input_path
else:
configfile_path = get_parent_ini_file(input_path)
shoutout(configfile_path)
# Load configuration
config = configparser.ConfigParser()
try:
with open(configfile_path, 'r', encoding='utf-8') as f:
config.read_file(f)
except UnicodeDecodeError:
with open(configfile_path, 'r', encoding='latin-1') as f:
config.read_file(f)
# Determine base directory and output directory from INI
config_dir = os.path.dirname(configfile_path)
output_folder = config.get("Paths", "outputfolder", fallback="")
if output_folder:
base_output_dir = os.path.join(config_dir, output_folder)
else:
base_output_dir = config_dir
# Create a new timestamped merge directory inside the output folder
timestamp = datetime.datetime.now().strftime("%Y%m%d_%H%M%S")
merge_dir = os.path.join(base_output_dir, f"merge_{timestamp}")
os.makedirs(merge_dir, exist_ok=True)
# --- Read Paths ---
if config.has_option('Paths', 'hkl_file'):
hkl_filename = config.get('Paths', 'hkl_file')
else:
raise KeyError("Missing 'hkl_file' in [Paths] section of INI")
# Output directory is based on output_folder setting
output_dir = base_output_dir
# --- Locate the .stream file to merge ---
# Always read 'streamfile' from INI; if missing, fallback to glob
if config.has_option('Paths', 'streamfile'):
stream_file = config.get('Paths', 'streamfile')
else:
stream_glob = os.path.join(output_dir, '*.stream')
stream_files = glob.glob(stream_glob)
if not stream_files:
raise FileNotFoundError(f"No .stream file found in directory {output_dir}")
stream_file = stream_files[0]
# --- Read Parameters ---
if not config.has_section('Parameters'):
raise KeyError("Missing [Parameters] section in INI")
# Required parameters
try:
pointgroup = config.get('Parameters', 'merging_pointgroup')
except (configparser.NoOptionError, configparser.NoSectionError):
raise KeyError("Missing 'merging_pointgroup' in [Parameters]")
try:
min_res = config.get('Parameters', 'merging_min_res')
except (configparser.NoOptionError, configparser.NoSectionError):
raise KeyError("Missing 'merging_min_res' in [Parameters]")
try:
iterations = config.get('Parameters', 'merging_iterations')
except (configparser.NoOptionError, configparser.NoSectionError):
raise KeyError("Missing 'merging_iterations' in [Parameters]")
# Optional threads, default to CPU count
if config.has_option('Parameters', 'merging_threads'):
threads = config.getint('Parameters', 'merging_threads')
else:
threads = multiprocessing.cpu_count()
# New optional merging parameters with defaults
model = config.get('Parameters', 'merging_model') if config.has_option('Parameters', 'merging_model') else "offset"
polarisation = config.get('Parameters', 'merging_polarisation') if config.has_option('Parameters', 'merging_polarisation') else "none"
min_measurements = config.get('Parameters', 'merging_min_measurements') if config.has_option('Parameters', 'merging_min_measurements') else "2"
max_adu = config.get('Parameters', 'merging_max_adu') if config.has_option('Parameters', 'merging_max_adu') else "inf"
push_res = config.get('Parameters', 'merging_push_res') if config.has_option('Parameters', 'merging_push_res') else "inf"
no_Bscale = config.getboolean('Parameters', 'merging_no_Bscale') if config.has_option('Parameters', 'merging_no_Bscale') else True
output_every_cycle = True
save_partialator_logs = config.getboolean('Parameters', 'merging_save_partialator_logs') if config.has_option('Parameters', 'merging_save_partialator_logs') else False
# Read optional unique axis
if config.has_option('Parameters', 'merging_unique_axis'):
unique_axis = config.get('Parameters', 'merging_unique_axis')
else:
unique_axis = None
# Construct merging command
merging_cmd = [
'partialator',
'-i', str(stream_file),
f'--model={model}',
'-j', str(threads),
'-o', os.path.join(merge_dir, hkl_filename),
'-y', pointgroup,
f'--polarisation={polarisation}',
f'--min-measurements={min_measurements}',
f'--max-adu={max_adu}',
f'--min-res={min_res}',
f'--push-res={push_res}'
]
# Add or omit flags based on booleans
if no_Bscale:
merging_cmd.append('--no-Bscale')
merging_cmd.append('--output-every-cycle')
# Only append a standalone --unique-axis if pointgroup does NOT already contain "_u"
if unique_axis and '_u' not in pointgroup:
merging_cmd.append(f'--unique-axis={unique_axis}')
# Append iterations and harvest-file. Partialator diagnostic logs are very
# file-heavy, so keep them opt-in.
merging_cmd += [
f'--iterations={iterations}',
'--harvest-file=' + os.path.join(merge_dir, "parameters.json"),
]
if save_partialator_logs:
merging_cmd.append('--log-folder=' + os.path.join(merge_dir, "pr-logs"))
log_print(f"Merging command: {' '.join(merging_cmd)}")
if save_partialator_logs:
os.makedirs(os.path.join(merge_dir, "pr-logs"), exist_ok=True)
# Run merging, capture stdout and stderr
stdout_path = os.path.join(merge_dir, "stdout.log")
stderr_path = os.path.join(merge_dir, "stderr.log")
try:
params = {
"pointgroup": pointgroup,
"min_res": min_res,
"iterations": iterations,
"threads": threads,
"model": model,
"polarisation": polarisation,
"min_measurements": min_measurements,
"max_adu": max_adu,
"push_res": push_res,
"no_Bscale": no_Bscale,
"output_every_cycle": output_every_cycle,
"save_partialator_logs": save_partialator_logs,
"unique_axis": unique_axis,
"stream_file": stream_file,
"merge_dir": merge_dir,
"command": " ".join(merging_cmd),
"output_hkl": os.path.join(merge_dir, hkl_filename),
"harvest_file": os.path.join(merge_dir, "parameters.json"),
}
if save_partialator_logs:
params["log_folder"] = os.path.join(merge_dir, "pr-logs")
for h5_path in _extract_h5_paths_from_stream(stream_file):
if not h5_path or not os.path.exists(h5_path):
continue
try:
import h5py
with h5py.File(h5_path, "r+") as h5file:
write_nxprocess_merging(
h5file,
program="coseda.merge",
method="crystfel.partialator",
parameters=params,
input_path=stream_file,
output_path=os.path.join(merge_dir, hkl_filename),
)
except Exception as exc:
log_print(f"Warning: failed to write NXprocess merging for {h5_path}: {exc}")
except Exception as exc:
log_print(f"Warning: failed to prepare NXprocess merging: {exc}")
with open(stdout_path, "w") as stdout, open(stderr_path, "w") as stderr:
log_print(f"Running merging for point group {pointgroup}, min_res={min_res}, iterations={iterations}")
subprocess.run(merging_cmd, stdout=stdout, stderr=stderr, check=True)
log_print("Merging completed successfully.")
# --- Convert merged HKL to MTZ ---
# Path to the merged HKL file
merged_hkl = os.path.join(merge_dir, hkl_filename)
# Determine cell file path from INI: it should reside relative to the INI's directory
if config.has_option('Paths', 'cellfile'):
cellfile_rel = config.get('Paths', 'cellfile')
cellfile_path = os.path.normpath(os.path.join(config_dir, cellfile_rel))
else:
cellfile_path = None
try:
# Call hkl2mtz with the merge directory as the output_dir
hkl2mtz(merged_hkl, merge_dir, cellfile_path)
except Exception as e:
log_print(f"HKL to MTZ conversion failed: {e}")
[docs]
def parse_scaling_file(input_filepath, output_filepath):
"""
Parse the scaling/refinement log output from partialator and save a CSV summary.
"""
cycle_pattern = re.compile(r'Scaling and refinement cycle (\d+) of \d+')
cchalf_pattern = re.compile(r'Overall CChalf = ([\d.]+) % \((\d+) reflections\)')
delta_cchalf_pattern = re.compile(r'deltaCChalf = ([\d.]+) ± ([\d.]+) %')
bad_crystals_pattern = re.compile(r'(\d+) bad crystals:')
ok_crystals_pattern = re.compile(r'(\d+) OK')
negative_delta_pattern = re.compile(r'(\d+) negative delta CC½')
reflections_input_pattern = re.compile(r'(\d+) reflections in input.')
data = []
current_cycle = None
current_cchalf = None
current_reflections = None
current_delta_cchalf = None
current_delta_cchalf_error = None
current_bad_crystals = None
current_ok_crystals = None
current_negative_delta = None
reflections_input = None
with open(input_filepath, 'r') as file:
for line in file:
line = line.strip()
cycle_match = cycle_pattern.match(line)
if cycle_match:
if current_cycle is not None:
data.append([
current_cycle, current_cchalf, current_reflections, current_delta_cchalf,
current_delta_cchalf_error, current_bad_crystals, current_ok_crystals, current_negative_delta
])
current_cycle = cycle_match.group(1)
current_cchalf = None
current_reflections = None
current_delta_cchalf = None
current_delta_cchalf_error = None
current_bad_crystals = None
current_ok_crystals = None
current_negative_delta = None
cchalf_match = cchalf_pattern.match(line)
if cchalf_match:
current_cchalf = cchalf_match.group(1)
current_reflections = cchalf_match.group(2)
delta_cchalf_match = delta_cchalf_pattern.match(line)
if delta_cchalf_match:
current_delta_cchalf = delta_cchalf_match.group(1)
current_delta_cchalf_error = delta_cchalf_match.group(2)
bad_crystals_match = bad_crystals_pattern.match(line)
if bad_crystals_match:
current_bad_crystals = bad_crystals_match.group(1)
ok_crystals_match = ok_crystals_pattern.match(line)
if ok_crystals_match:
current_ok_crystals = ok_crystals_match.group(1)
negative_delta_match = negative_delta_pattern.match(line)
if negative_delta_match:
current_negative_delta = negative_delta_match.group(1)
reflections_input_match = reflections_input_pattern.match(line)
if reflections_input_match:
reflections_input = reflections_input_match.group(1)
if current_cycle is not None:
data.append([
current_cycle, current_cchalf, current_reflections, current_delta_cchalf,
current_delta_cchalf_error, current_bad_crystals, current_ok_crystals, current_negative_delta
])
header = [
'Cycle', 'Overall CChalf (%)', 'Reflections', 'deltaCChalf', 'deltaCChalf Error (%)',
'Bad Crystals', 'OK Crystals', 'Negative delta CC½'
]
log_print("\t".join(header))
for row in data:
log_print("\t".join(str(value) if value is not None else '' for value in row))
with open(output_filepath, 'w', newline='') as csvfile:
csv_writer = csv.writer(csvfile)
csv_writer.writerow(header)
for row in data:
csv_writer.writerow(row)
if reflections_input is not None:
log_print(f"Reflections in input: {reflections_input}")