diff --git a/build-deb.sh b/build-deb.sh index d771e06..e32e7ab 100755 --- a/build-deb.sh +++ b/build-deb.sh @@ -68,12 +68,12 @@ then bookworm) debian_version=12 ;; - trixie) - debian_version=13 - ;; - sid) - debian_version=13 - ;; + trixie) + debian_version=13 + ;; + sid) + debian_version=13 + ;; esac fi @@ -109,6 +109,7 @@ optional arguments: --debian10 Simulate a debian 10 Buster system --debian11 Simulate a debian 11 Bullseye system --debian12 Simulate a debian 12 Bookworm system + --debian13 Simulate a debian 13 Trixie system " install=0 diff --git a/doc/source/dahu.rst b/doc/source/dahu.rst index d94ed72..c9d5684 100644 --- a/doc/source/dahu.rst +++ b/doc/source/dahu.rst @@ -11,15 +11,15 @@ The *dahu* server executes **jobs**: * The job (de-) serializes JSON strings coming from/returning to Tango * Jobs are executed asynchronously, the request for calculation is answered instantaneously with a *jobid* (an integer, unique for the process). * The *jobid* can be used to poll the server for the status of the job or for manual synchronization (mind that Tango can time-out!). -* When jobs are finished, the client is notified via Tango events about the status +* When jobs are finished, the client is notified via **Tango events** about the status change * Results can be retrieved after the job has finished. Jobs execute **plugin**: ------------------------ -* Plugins are written in Python (extension in Cython or OpenCL are common) +* Plugins are written in Python (extensions in Cython or OpenCL are common) * Plugins can be classes or simple functions -* The input and output MUST be JSON-seriablisable as simple dictionnaries +* The input and output MUST be JSON-serializable as simple dictionaries * Plugins are dynamically loaded from Python modules * Plugins can be profiled for performance analysis diff --git a/plugins/bm29/__init__.py b/plugins/bm29/__init__.py index 20fc1a4..23e75b0 100644 --- a/plugins/bm29/__init__.py +++ b/plugins/bm29/__init__.py @@ -11,7 +11,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "03/12/2024" +__date__ = "05/05/2025" __status__ = "development" __version__ = "0.2.0" @@ -19,6 +19,8 @@ from .integrate import IntegrateMultiframe from .subtracte import SubtractBuffer from .hplc import HPLC +from .mesh import Mesh register(IntegrateMultiframe, fqn="bm29.integratemultiframe") register(SubtractBuffer, fqn="bm29.subtractbuffer") register(HPLC, fqn="bm29.hplc") +register(Mesh, fqn="bm29.mesh") \ No newline at end of file diff --git a/plugins/bm29/hplc.py b/plugins/bm29/hplc.py index 5e3a8da..1efa7f1 100644 --- a/plugins/bm29/hplc.py +++ b/plugins/bm29/hplc.py @@ -10,7 +10,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "27/05/2025" +__date__ = "23/02/2026" __status__ = "development" __version__ = "0.3.0" @@ -21,6 +21,7 @@ from math import log, pi import posixpath import copy +import zipfile from collections import namedtuple from urllib3.util import parse_url from dahu.plugin import Plugin @@ -29,11 +30,12 @@ logger = logging.getLogger("bm29.hplc") import numpy import h5py -import pyFAI, pyFAI.azimuthalIntegrator, pyFAI.units +import pyFAI, pyFAI.integrator.azimuthal, pyFAI.units from pyFAI.method_registry import IntegrationMethod import freesas, freesas.cormap, freesas.invariants from freesas.autorg import auto_gpa, autoRg, auto_guinier from freesas.bift import BIFT +from freesas.app.extract_ascii import write_ascii from scipy.optimize import minimize import scipy.signal import scipy.ndimage @@ -132,6 +134,41 @@ def build_background(I, std=None, keep=0.3): return bg_avg, bg_std, to_keep +def save_zip(filename, config, I, sigma): + """Save a stack of I into a zipfile with each frames in a dat-file. + + :param filename: name of the zip-file + :param confif: this is some NexusJuice namedtuple. we use only q and the sample description. + :param I: 2D array with the intensity of the stack of curves + :param sigma: 2D array with the uncertainties of the stack of frames + :return: nothing + """ + basename = os.path.basename(filename) + base = os.path.splitext(basename)[0] + destz = base + "_%04i.dat" + common = {"q": config.q} + if config.sample: + sample = config.sample + if sample.name: + common["sample"]: sample.name + if sample.buffer: + common["buffer"] = sample.buffer + if sample.temperature_env: + common["storage temperature"] = sample.temperature_env + if sample.temperature: + common["exposure temperature"] = sample.temperature + if sample.concentration: + common["concentration"] = sample.concentration + res = [] + for i, s in zip(I, sigma): + r = copy.copy(common) + r["I"] = i + r["std"] = s + res.append(r) + with zipfile.ZipFile(filename, "w") as z: + for idx, frame in enumerate(res): + z.writestr(destz % idx, write_ascii(frame)) + class HPLC(Plugin): """ Rebuild the complete chromatogram and perform basic analysis on it. @@ -291,13 +328,15 @@ def create_nexus(self): time_ds.attrs["interpretation"] = "spectrum" time_ds.attrs["long_name"] = "Time stamps (s)" - integration_data = nxs.new_class(chroma_grp, "results", "NXdata") + integration_data = nxs.new_class(chroma_grp, "result", "NXdata") chroma_grp.attrs["title"] = str(self.juices[0].sample) int_ds = integration_data.create_dataset("I", data=numpy.ascontiguousarray(I, dtype=numpy.float32)) std_ds = integration_data.create_dataset("errors", data=numpy.ascontiguousarray(sigma, dtype=numpy.float32)) q_ds = integration_data.create_dataset("q", data=self.juices[0].q) q_ds.attrs["interpretation"] = "spectrum" + q_ds.attrs["unit"] = unit_name + q_ds.attrs["long_name"] = "Scattering vector q (nm⁻¹)" integration_data.attrs["signal"] = "I" integration_data.attrs["axes"] = [".", "q"] integration_data.attrs["SILX_style"] = SAXS_STYLE @@ -309,6 +348,9 @@ def create_nexus(self): int_ds.attrs["scale"] = "log" std_ds.attrs["interpretation"] = "spectrum" + save_zip(os.path.splitext(self.output_file)[0]+".zip", + self.juices[0], I, sigma) + # Process 2: SVD decomposition svd_grp = nxs.new_class(entry_grp, "2_SVD", "NXprocess") svd_grp["sequence_index"] = self.sequence_index() @@ -385,7 +427,7 @@ def create_nexus(self): self.to_pyarch["buffer_frames"] = to_keep self.to_pyarch["buffer_I"] = bg_avg self.to_pyarch["buffer_Stdev"] = bg_std - bg_data = nxs.new_class(bg_grp, "results", "NXdata") + bg_data = nxs.new_class(bg_grp, "result", "NXdata") bg_data.attrs["signal"] = "I" bg_data.attrs["SILX_style"] = SAXS_STYLE bg_data.attrs["axes"] = radial_unit @@ -493,16 +535,15 @@ def one_fraction(self, fraction, index, nxs, top_grp): guinier_autorg = nxs.new_class(guinier_grp, "autorg", "NXcollection") guinier_gpa = nxs.new_class(guinier_grp, "gpa", "NXcollection") guinier_guinier = nxs.new_class(guinier_grp, "guinier", "NXcollection") - guinier_data = nxs.new_class(guinier_grp, "results", "NXdata") + guinier_data = nxs.new_class(guinier_grp, "result", "NXdata") guinier_data.attrs["SILX_style"] = NORMAL_STYLE guinier_data.attrs["title"] = "Guinier analysis" # Stage4 processing: autorg and auto_gpa sasm = numpy.vstack((q, I_frc, sigma_frc)).T - try: gpa = auto_gpa(sasm) except Exception as error: - guinier_gpa["Failed"] = "%s: %s" % (error.__class__.__name__, error) + guinier_gpa["Failed"] = f"{error.__class__.__name__}: {error}" gpa = None else: # "Rg sigma_Rg I0 sigma_I0 start_point end_point quality aggregated" @@ -520,7 +561,7 @@ def one_fraction(self, fraction, index, nxs, top_grp): try: guinier = auto_guinier(sasm) except Exception as error: - guinier_guinier["Failed"] = "%s: %s" % (error.__class__.__name__, error) + guinier_guinier["Failed"] = f"{error.__class__.__name__}: {error}" guinier = None else: # "Rg sigma_Rg I0 sigma_I0 start_point end_point quality aggregated" @@ -540,7 +581,7 @@ def one_fraction(self, fraction, index, nxs, top_grp): try: autorg = autoRg(sasm) except Exception as err: - guinier_autorg["Failed"] = "%s: %s" % (err.__class__.__name__, err) + guinier_autorg["Failed"] = f"{err.__class__.__name__}: {err}" autorg = None else: if autorg.Rg < 0: @@ -594,7 +635,7 @@ def one_fraction(self, fraction, index, nxs, top_grp): dlogI = err[mask] / logI q2_ds = guinier_data.create_dataset("q2", data=q2.astype(numpy.float32)) q2_ds.attrs["unit"] = radius_unit + "⁻²" - q2_ds.attrs["long_name"] = "q² (%s⁻²)" % radius_unit + q2_ds.attrs["long_name"] = f"q² ({radius_unit}⁻²)" q2_ds.attrs["interpretation"] = "spectrum" lnI_ds = guinier_data.create_dataset("logI", data=logI.astype(numpy.float32)) lnI_ds.attrs["long_name"] = "log(I)" @@ -623,7 +664,7 @@ def one_fraction(self, fraction, index, nxs, top_grp): kratky_grp["program"] = "freesas.autorg" kratky_grp["version"] = freesas.version kratky_grp["date"] = get_isotime() - kratky_data = nxs.new_class(kratky_grp, "results", "NXdata") + kratky_data = nxs.new_class(kratky_grp, "result", "NXdata") kratky_data.attrs["SILX_style"] = NORMAL_STYLE kratky_data.attrs["title"] = "Dimensionless Kratky plots" kratky_grp.attrs["default"] = posixpath.relpath(kratky_data.name, kratky_grp.name) @@ -637,21 +678,21 @@ def one_fraction(self, fraction, index, nxs, top_grp): qRg_ds = kratky_data.create_dataset("qRg", data=xdata.astype(numpy.float32)) qRg_ds.attrs["interpretation"] = "spectrum" qRg_ds.attrs["long_name"] = "q·Rg (unit-less)" - k_ds = kratky_data.create_dataset("q2Rg2I/I0", data=ydata.astype(numpy.float32)) + k_ds = kratky_data.create_dataset("q2Rg2I÷I0", data=ydata.astype(numpy.float32)) k_ds.attrs["interpretation"] = "spectrum" k_ds.attrs["long_name"] = "q²Rg²I(q)/I₀" ke_ds = kratky_data.create_dataset("errors", data=dy.astype(numpy.float32)) ke_ds.attrs["interpretation"] = "spectrum" kratky_data_attrs = kratky_data.attrs - kratky_data_attrs["signal"] = "q2Rg2I/I0" - kratky_data_attrs["axes"] = "qRg" + kratky_data_attrs["signal"] = k_ds.name + kratky_data_attrs["axes"] = qRg_ds.name # stage 6: Rambo-Tainer invariant rti_grp = nxs.new_class(f_grp, "4_invariants", "NXprocess") rti_grp["sequence_index"] = self.sequence_index() rti_grp["program"] = "freesas.invariants" rti_grp["version"] = freesas.version - rti_data = nxs.new_class(rti_grp, "results", "NXdata") + rti_data = nxs.new_class(rti_grp, "result", "NXdata") # average_data.attrs["SILX_style"] = SAXS_STYLE # average_data.attrs["signal"] = "intensity_normed" # Rambo_Tainer @@ -688,7 +729,7 @@ def one_fraction(self, fraction, index, nxs, top_grp): bift_grp["program"] = "freesas.bift" bift_grp["version"] = freesas.version bift_grp["date"] = get_isotime() - bift_data = nxs.new_class(bift_grp, "results", "NXdata") + bift_data = nxs.new_class(bift_grp, "result", "NXdata") bift_data.attrs["SILX_style"] = NORMAL_STYLE bift_data.attrs["title"] = "Pair distance distribution function p(r)" @@ -726,11 +767,11 @@ def one_fraction(self, fraction, index, nxs, top_grp): res = minimize(bo.opti_evidence, (Dmax, log(alpha)), args=(npt, use_wisdom), method="powell") cfg_grp["Powell_steps"] = res.nfev cfg_grp["Monte-Carlo_steps"] = 0 + stats = bo.calc_stats() except Exception as error: - bift_grp["Failed"] = "%s: %s" % (error.__class__.__name__, error) + bift_grp["Failed"] = f"{error.__class__.__name__}: {error}" bo = None else: - stats = bo.calc_stats() bift_grp["alpha"] = stats.alpha_avg bift_grp["alpha_error"] = stats.alpha_std self.Dmax = bift_grp["Dmax"] = stats.Dmax_avg @@ -750,7 +791,7 @@ def one_fraction(self, fraction, index, nxs, top_grp): r_ds.attrs["interpretation"] = "spectrum" r_ds.attrs["unit"] = radius_unit - r_ds.attrs["long_name"] = "radius r(%s)" % radius_unit + r_ds.attrs["long_name"] = f"radius r({radius_unit})" p_ds = bift_data.create_dataset("p(r)", data=stats.density_avg.astype(numpy.float32)) p_ds.attrs["interpretation"] = "spectrum" bift_data["errors"] = stats.density_std @@ -882,13 +923,13 @@ def read_nexus(filename): entry_name = nxsr.h5.attrs["default"] entry_grp = nxsr.h5[entry_name] h5path = entry_grp.name - nxdata_grp = nxsr.h5[entry_grp.attrs["default"]] + nxdata_grp = entry_grp[entry_grp.attrs["default"]] assert nxdata_grp.name.endswith("hplc") # we are reading HPLC data signal = nxdata_grp.attrs["signal"] axis = nxdata_grp.attrs["axes"] Isum = nxdata_grp[signal][()] idx = nxdata_grp[axis][()] - integrated = nxdata_grp.parent["results"] + integrated = nxdata_grp.parent["result"] signal = integrated.attrs["signal"] I = integrated[signal][()] axes = integrated.attrs["axes"][-1] @@ -956,6 +997,7 @@ def send_to_icat(self): raw=os.path.dirname(os.path.abspath(self.input_files[0])), path=os.path.dirname(os.path.abspath(self.output_file)), data=to_icat, + dataset="HPLC", gallery=gallery, metadata=metadata) @@ -974,3 +1016,4 @@ def save_csv(self, filename, sum_I, Rg): + diff --git a/plugins/bm29/icat.py b/plugins/bm29/icat.py index e6eac7a..2957345 100644 --- a/plugins/bm29/icat.py +++ b/plugins/bm29/icat.py @@ -11,7 +11,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "21/02/2025" +__date__ = "05/05/2025" __status__ = "development" version = "0.3.0" @@ -46,7 +46,7 @@ def send_icat(proposal=None, beamline=None, sample=None, dataset=None, path=None :param proposal: mx1324 :param beamline: name of the beamline :param sample: sample name as registered in icat - :param dataset: name given by BLISS + :param dataset: name of the dataset: integration, subtraction, HPLC, ... :param path: directory name where processed data are staying :param raw: list of directory name of the raw data (not the processed ones) :param data: dict with all data sent to iCat @@ -55,37 +55,44 @@ def send_icat(proposal=None, beamline=None, sample=None, dataset=None, path=None :return: data sent to icat as a dict """ gallery = _ensure_gallery(gallery) + print(gallery) tmp = gallery.strip("/").split("/") - idx_process = [i for i,j in enumerate(tmp) if j.lower().startswith("process")][-1] - if tmp[idx_process] == "processed": - assert idx_process>=6 - if proposal is None: - proposal = tmp[idx_process-6] - if beamline is None: - beamline = tmp[idx_process-5] - if sample is None: - sample = tmp[idx_process-2] - if dataset is None: - dataset = tmp[idx_process+1] - if path is None: - path = os.path.dirname(gallery) - if raw is None: - raw = os.path.abspath(gallery[:gallery.lower().index("process")]) - elif tmp[idx_process] == "PROCESSED_DATA": - if proposal is None: - proposal = tmp[idx_process-3] - if beamline is None: - beamline = tmp[idx_process-2] - if sample is None: - sample = tmp[idx_process+1] - if dataset is None: - dataset = tmp[idx_process+2] - if path is None: - path = os.path.dirname(gallery) - if raw is None: - raw = os.path.dirname(os.path.dirname(os.path.abspath(gallery.replace("PROCESSED_DATA", "RAW_DATA")))) + idx_process = [i for i,j in enumerate(tmp) if j.lower().startswith("process")] + if idx_process: + idx_process=idx_process[-1] + if tmp[idx_process] == "processed": + assert idx_process>=6 + if proposal is None: + proposal = tmp[idx_process-6] + if beamline is None: + beamline = tmp[idx_process-5] + if sample is None: + sample = tmp[idx_process-2] + if dataset is None: + dataset = tmp[idx_process+1] + if path is None: + path = os.path.dirname(gallery) + if raw is None: + raw = os.path.abspath(gallery[:gallery.lower().index("process")]) + elif tmp[idx_process] == "PROCESSED_DATA": + if proposal is None: + proposal = tmp[idx_process-3] + if beamline is None: + beamline = tmp[idx_process-2] + if sample is None: + sample = tmp[idx_process+1] + if dataset is None: + dataset = tmp[idx_process+2] + if path is None: + path = os.path.dirname(gallery) + if raw is None: + raw = os.path.dirname(os.path.dirname(os.path.abspath(gallery.replace("PROCESSED_DATA", "RAW_DATA")))) + else: + logger.error("Unrecognized path layout") + return else: - logger.error("Unrecognized path layout") + logger.error("No gallery provided") + return if metadata is None: metadata = {} diff --git a/plugins/bm29/integrate.py b/plugins/bm29/integrate.py index c447015..b631757 100644 --- a/plugins/bm29/integrate.py +++ b/plugins/bm29/integrate.py @@ -11,9 +11,9 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "04/06/2025" +__date__ = "05/02/2026" __status__ = "development" -__version__ = "0.3.0" +__version__ = "0.3.1" import os import time @@ -30,7 +30,7 @@ import numpy import h5py import pyFAI -import pyFAI.azimuthalIntegrator +from pyFAI import integrator import freesas import freesas.cormap @@ -189,6 +189,7 @@ def setup(self, kwargs=None): self.monitor_values = numpy.array(self.input.get("monitor_values", 1), dtype=numpy.float64) if self.input.get("average_out_monitor_values"): self.monitor_values = numpy.zeros_like(self.monitor_values) + self.monitor_values.mean() + self.log_warning("Averaging-out the monitor values !") self.normalization_factor = float(self.input.get("normalization_factor", 1)) self.scale_factor = float(self.input.get("exposure_time", 1)) / self.normalization_factor @@ -197,6 +198,7 @@ def teardown(self): logger.debug("IntegrateMultiframe.teardown") # export the output file location self.output["output_file"] = self.output_file + self.output["monitor_noise_%"] = 100 * self.monitor_values.std() / self.monitor_values.mean() if self.nxs is not None: self.nxs.close() if self.ai is not None: @@ -249,7 +251,7 @@ def wait_file(self, filename, timeout=None): """ timeout = self.timeout if timeout is None else timeout end_time = time.perf_counter() + timeout - dirname = os.path.dirname(filename) + dirname = os.path.dirname(filename) or "." while not os.path.isdir(dirname): if time.perf_counter() > end_time: self.log_error(f"Filename {filename} did not appear in {timeout} seconds") @@ -273,8 +275,10 @@ def wait_file(self, filename, timeout=None): def create_nexus(self): "create the nexus result file with basic structure" - if not os.path.isdir(os.path.dirname(self.output_file)): - os.makedirs(os.path.dirname(self.output_file)) + dirname = os.path.dirname(self.output_file) + if dirname: + if not os.path.isdir(dirname): + os.makedirs(dirname) creation_time = os.stat(self.input_file).st_ctime nxs = self.nxs = Nexus(self.output_file, mode="w", creator="dahu") @@ -395,15 +399,15 @@ def create_nexus(self): pol_ds = cfg_grp.create_dataset("polarization_factor", data=polarization_factor) pol_ds.attrs["comment"] = "Between -1 and +1, 0 for circular" cfg_grp.create_dataset("integration_method", data=json.dumps(method.method._asdict())) - integration_data = nxs.new_class(integration_grp, "results", "NXdata") + integration_data = nxs.new_class(integration_grp, "result", "NXdata") integration_grp.attrs["title"] = str(self.sample) # Stage 1 processing: Integration frame per frame - integrate1_results = self.process1_integration(self.input_frames) + integrate1_result = self.process1_integration(self.input_frames) radial_unit, unit_name = str(self.unit).split("_", 1) - q = numpy.ascontiguousarray(integrate1_results.radial, numpy.float32) - I = numpy.ascontiguousarray(integrate1_results.intensity, dtype=numpy.float32) - sigma = numpy.ascontiguousarray(integrate1_results.sigma, dtype=numpy.float32) + q = numpy.ascontiguousarray(integrate1_result.radial, numpy.float32) + I = numpy.ascontiguousarray(integrate1_result.intensity, dtype=numpy.float32) + sigma = numpy.ascontiguousarray(integrate1_result.sigma, dtype=numpy.float32) self.to_memcached[radial_unit] = q self.to_memcached["I"] = I @@ -428,7 +432,7 @@ def create_nexus(self): hplc_data = nxs.new_class(integration_grp, "hplc", "NXdata") hplc_data.attrs["title"] = "Chromatogram" - sum_ds = hplc_data.create_dataset("sum", data=numpy.ascontiguousarray(integrate1_results.intensity.sum(axis=-1), dtype=numpy.float32)) + sum_ds = hplc_data.create_dataset("sum", data=numpy.ascontiguousarray(integrate1_result.intensity.sum(axis=-1), dtype=numpy.float32)) sum_ds.attrs["interpretation"] = "spectrum" sum_ds.attrs["long_name"] = "Summed Intensity" hplc_data["frame_ids"] = frame_ds @@ -449,7 +453,7 @@ def create_nexus(self): cormap_grp["program"] = "freesas.cormap" cormap_grp["version"] = freesas.version cormap_grp["date"] = get_isotime() - cormap_data = nxs.new_class(cormap_grp, "results", "NXdata") + cormap_data = nxs.new_class(cormap_grp, "result", "NXdata") cormap_data.attrs["SILX_style"] = NORMAL_STYLE cfg_grp = nxs.new_class(cormap_grp, "configuration", "NXcollection") @@ -459,33 +463,33 @@ def create_nexus(self): cfg_grp["fidelity_rel"] = fidelity_rel # Stage 2 processing - cormap_results = self.process2_cormap(integrate1_results.intensity, fidelity_abs, fidelity_rel) + cormap_result = self.process2_cormap(integrate1_result.intensity, fidelity_abs, fidelity_rel) cormap_data.attrs["signal"] = "probability" - cormap_ds = cormap_data.create_dataset("probability", data=cormap_results.probability) + cormap_ds = cormap_data.create_dataset("probability", data=cormap_result.probability) cormap_ds.attrs["interpretation"] = "image" cormap_ds.attrs["long_name"] = "Probability to be the same" - count_ds = cormap_data.create_dataset("count", data=cormap_results.count) + count_ds = cormap_data.create_dataset("count", data=cormap_result.count) count_ds.attrs["interpretation"] = "image" count_ds.attrs["long_name"] = "Longest sequence where curves do not cross each other" - to_merge_ds = cormap_data.create_dataset("to_merge", data=numpy.arange(*cormap_results.tomerge, dtype=numpy.uint16)) + to_merge_ds = cormap_data.create_dataset("to_merge", data=numpy.arange(*cormap_result.tomerge, dtype=numpy.uint16)) to_merge_ds.attrs["long_name"] = "Index of equivalent frames" cormap_grp.attrs["default"] = posixpath.relpath(cormap_data.name, cormap_grp.name) if self.ispyb.url: - self.to_pyarch["merged"] = cormap_results.tomerge + self.to_pyarch["merged"] = cormap_result.tomerge # Process 3: time average and standard deviation average_grp = nxs.new_class(entry_grp, "3_time_average", "NXprocess") average_grp["sequence_index"] = 3 average_grp["program"] = fully_qualified_name(self.__class__) average_grp["version"] = __version__ - average_data = nxs.new_class(average_grp, "results", "NXdata") + average_data = nxs.new_class(average_grp, "result", "NXdata") average_data.attrs["SILX_style"] = SAXS_STYLE average_data.attrs["signal"] = "intensity_normed" # Stage 3 processing - res3 = self.process3_average(cormap_results.tomerge) + res3 = self.process3_average(cormap_result.tomerge) Iavg = numpy.ascontiguousarray(res3.average, dtype=numpy.float32) sigma_avg = numpy.ascontiguousarray(res3.deviation, dtype=numpy.float32) @@ -514,7 +518,7 @@ def create_nexus(self): ai2_grp["program"] = "pyFAI" ai2_grp["version"] = pyFAI.version ai2_grp["date"] = get_isotime() - ai2_data = nxs.new_class(ai2_grp, "results", "NXdata") + ai2_data = nxs.new_class(ai2_grp, "result", "NXdata") ai2_data.attrs["signal"] = "I" ai2_data.attrs["axes"] = radial_unit ai2_data.attrs["SILX_style"] = SAXS_STYLE @@ -656,6 +660,7 @@ def send_to_icat(self): raw=os.path.dirname(os.path.dirname(os.path.abspath(self.input_file))), path=os.path.dirname(os.path.abspath(self.output_file)), data=to_icat, + dataset = "integrate", gallery=self.ispyb.gallery or os.path.join(os.path.dirname(os.path.abspath(self.output_file)), "gallery"), metadata=metadata) diff --git a/plugins/bm29/ispyb.py b/plugins/bm29/ispyb.py index d27b7ae..ef589ef 100644 --- a/plugins/bm29/ispyb.py +++ b/plugins/bm29/ispyb.py @@ -228,7 +228,17 @@ def _mk_filename(self, index, path, basename="frame", ext=".dat"): os.makedirs(dest) except Exception as err: logger.error("Unable to create directory %s: %s: %s", dest, type(err), err) - os.stat(dest) #this is to enforce the mounting of the directory + else: + # Once the directory is pretendily created, write something in is an delete it. + delete_me = os.path.join(dest, "delete.me") + res = os.system(f"touch {delete_me}") + if res: + logger.error(f"`touch {delete_me}` return error {res}, directory creation did probably not work as expected !") + try: + os.remove(delete_me) + except FileNotFoundError: + logger.error(f"`rm {delete_me}` raised FileNotFoundError, directory creation did probably not work as expected !") + os.stat(dest) # this is to enforce the mounting of the directory if isinstance(index, int): filename = os.path.join(dest, "%s_%04d%s" % (basename, index, ext)) else: diff --git a/plugins/bm29/memcached.py b/plugins/bm29/memcached.py index 4bd8432..a7a95dd 100644 --- a/plugins/bm29/memcached.py +++ b/plugins/bm29/memcached.py @@ -10,7 +10,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "21/02/2025" +__date__ = "22/04/2025" __status__ = "development" __version__ = "0.3.0" @@ -29,5 +29,7 @@ def to_memcached(dico): mc = memcache.Client([(SERVER, 11211)]) rc["server"] = socket.getfqdn()+":11211" for k, v in dico.items(): + if len(k)>250: + k = k[-250:] rc[k] = mc.set(k, v) return rc diff --git a/plugins/bm29/mesh.py b/plugins/bm29/mesh.py new file mode 100644 index 0000000..cfbd26d --- /dev/null +++ b/plugins/bm29/mesh.py @@ -0,0 +1,547 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +"""Data Analysis plugin for BM29: BioSaxs + +* Mesh mode: Rebuild the complete map and performs basic analysis on it. +""" + +__authors__ = ["Jérôme Kieffer"] +__contact__ = "Jerome.Kieffer@ESRF.eu" +__license__ = "MIT" +__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" +__date__ = "04/12/2025" +__status__ = "development" +__version__ = "0.1.0" + +import os +import posixpath +import json +import glob +import collections +from dataclasses import dataclass, fields, asdict +import numpy +from dahu.plugin import Plugin +import h5py +import matplotlib +from matplotlib.pyplot import subplots +import pyFAI +from pyFAI.method_registry import IntegrationMethod +from pyFAI.io.ponifile import PoniFile +from pyFAI.io.diffmap_config import DiffmapConfig, WorkerConfig, MotorRange, ListDataSet, DataSet +from .common import Sample, Ispyb, get_equivalent_frames, cmp_float, get_integrator, KeyCache, \ + polarization_factor, method, Nexus, get_isotime, SAXS_STYLE, NORMAL_STYLE, \ + create_nexus_sample +matplotlib.use("Agg") +NexusJuice = collections.namedtuple("NexusJuice", "filename h5path npt unit idx Isum q I sigma poni mask energy polarization method sample timestamps") +Position = collections.namedtuple('Position', 'index slow fast') + +@dataclass(slots=True) +class Scan: + """Class describing 2D mesh-scan""" + fast_motor_name: str = "fast" + fast_motor_start: float = 0.0 + fast_motor_stop: float = 1.0 + fast_motor_step: int = 0 # bliss-like steps, actually step+1 points + slow_motor_name: str = "slow" + slow_motor_start: float = 0.0 + slow_motor_stop: float = 1.0 + slow_motor_step: int = 0 # bliss-like steps, actually step+1 points + backnforth: bool = False + + def __repr__(self): + return json.dumps(self.as_dict(), indent=4) + + def as_dict(self): + """Like asdict, without extra features: + + :return: dict which can be JSON-serialized + """ + dico = {} + for key, value in asdict(self).items(): + dico[key] = value + return dico + + @classmethod + def from_dict(cls, dico): + """Alternative constructor, + :param dico: dict with the config + :return: instance of the dataclass + """ + to_init = dico + self = cls(**to_init) + return self + + def save(self, filename): + """Dump the content of the dataclass as JSON file""" + with open(filename, "w") as w: + w.write(json.dumps(self.as_dict(), indent=2)) + + def get_pos(self, idx=None): + """ + Calculate the position in the mesh scan fram according to its index + + :param idx: index of current frame + :return: namedtuple Position=(frame_index, slow_motor_position, fast_motor_position) + if valid else None + """ + width = self.fast_motor_step + 1 + line = idx // width + if self.backnforth and line % 2 == 1: + row = self.fast_motor_step - (idx % width) + else: + row = idx % width + if line<=self.slow_motor_step and row<=self.fast_motor_step: + return Position(idx, line, row) + else: + return None + + @property + def shape(self): + return (self.slow_motor_step + 1, self.fast_motor_step + 1) + + @classmethod + def parse(cls, text): + """Alternative constructor, + :param text: string containing the bliss command (starting with `amesh`) + :return: instance of the dataclass + """ + res = None + if text.startswith("amesh"): + words = text.split() + if len(words) >= 9: + res = cls(words[1], float(words[2]),float(words[3]), int(words[4]), + words[5], float(words[6]),float(words[7]), int(words[8]), + True) + return res + +def input_from_master(master_file): + """Convert a bliss masterfile containing one or multiple NXentry + with a 2D scan into a set of plugins to be launched + + :param master_file: path to a bliss masterfile + :return: list of job-input-dicts. + """ + result = [] + + with Nexus(master_file, mode="r", pure=True) as master: + for entry in master.get_entries(): + job = {"plugin_name":"bm29.mesh"} + name = posixpath.split(entry.name)[-1] + filename = os.path.abspath(entry.file.filename) + dirtree = filename.split(os.sep)[:-1] + raw_idx = dirtree.index("RAW_DATA") + dirtree[raw_idx] = "PROCESSED_DATA" + job["output_file"] = os.sep.join(dirtree + ["mesh", "mesh.h5"]) + wildcard = os.sep.join(dirtree + ["integrate", "*.h5"]) + input_files = glob.glob(wildcard) + input_files.sort() + job["integrated_files"] = input_files + + title = entry.get("title", "") + if isinstance(title, h5py.Dataset): + title = title[()] + if isinstance(title, bytes): + title = title.decode() + if title: + scan = Scan.parse(title) + if scan is None: + continue + else: + job["scan"] = scan.as_dict() + result.append(job) + return result + + +def mesh_plot(mesh, + title="mesh-scan", + x_label="fast motor", + y_label="slow motor", + x_range=None, + y_range=None, + filename=None, + img_format="png", + ax=None, + labelsize=None, + fontsize=None,): + """ + Generate an image of the mesh-scan + + :param mesh: sum of integrated data + :param filename: name of the file where the cuve should be saved + :param img_format: image image format + :param ax: subplotib where to plot in + :return: the matplotlib figure + """ + if ax: + fig = ax.figure + else: + fig, ax = subplots(figsize=(12, 10)) + ax.imshow(mesh, cmap="binary") + ax.set_xlabel(x_label, fontsize=fontsize) + ax.set_ylabel(y_label, fontsize=fontsize) + ax.set_title(title) + if x_range is not None: + ax.xaxis.set_major_formatter(lambda x, pos: f"{numpy.interp(x, numpy.arange(x_range.size), x_range):.2f}") + if y_range is not None: + ax.yaxis.set_major_formatter(lambda x, pos: f"{numpy.interp(x, numpy.arange(y_range.size), y_range):.2f}") + + ax.tick_params(axis="x", labelsize=labelsize) + ax.tick_params(axis="y", labelsize=labelsize) + + if filename: + if img_format: + fig.savefig(filename, format=img_format) + else: + fig.savefig(filename) + return fig + + + +class Mesh(Plugin): + """ Rebuild the complete map and perform basic analysis on it. + + Typical JSON file: + { + "integrated_files": ["img_001.h5", "img_002.h5"], + "output_file": "mesh.h5" + "ispyb": { + "url": "http://ispyb.esrf.fr:1234", + "pyarch": "/data/pyarch/mx1234/sample", + "measurement_id": -1, + "collection_id": -1 + }, + "scan": { + "fast_motor_name": "chipz", + "fast_motor_start": -2.2, + "fast_motor_stop": -3.2, + "fast_motor_step": 3, + "slow_motor_name": "chipy", + "slow_motor_start": -5.2, + "slow_motor_stop": -12.2, + "slow_motor_step": 7, + "backnforth": False + } + "transpose": False, + "wait_for": [jobid_img001, jobid_img002], + "plugin_name": "bm29.mesh" + } + """ + + def __init__(self): + Plugin.__init__(self) + self.input_files = [] + self.nxs = None + self.output_file = None + self.scan = Scan() + self.transpose = None + self.juices = [] + self.to_pyarch = {} + self.ispyb = None + self._pid = 0 + + def sequence_index(self): + value = self._pid + self._pid += 1 + return value + + def setup(self): + Plugin.setup(self) + + for job_id in self.input.get("wait_for", []): + self.wait_for(job_id) + + self.input_files = [os.path.abspath(i) for i in self.input.get("integrated_files", "")] + self.output_file = self.input.get("output_file") + if not self.output_file: + dirname, basename = os.path.split(os.path.commonprefix(self.input_files) + "_mesh.h5") + dirname = os.path.dirname(dirname) + dirname = os.path.join(dirname, "mesh") + self.output_file = os.path.join(dirname, basename) + if not os.path.isdir(dirname): + try: + os.makedirs(dirname) + except Exception as err: + self.log_warning(f"Unable to create dir {dirname}. {type(err)}: {err}") + + self.log_warning(f"No output file provided, using {self.output_file}") + + #Manage gallery here + dirname = os.path.dirname(self.output_file) + gallery = os.path.join(dirname, "gallery") + if not os.path.isdir(gallery): + try: + os.makedirs(gallery) + except Exception as err: + self.log_warning(f"Unable to create dir {gallery}. {type(err)}: {err}") + ispydict = self.input.get("ispyb", {}) + ispydict["gallery"] = gallery + self.ispyb = Ispyb._fromdict(ispydict) + self.scan = Scan.from_dict(self.input.get("scan", {})) + + def process(self): + self.create_nexus() + self.to_pyarch["hdf5_filename"] = self.output_file + # self.to_pyarch["chunk_size"] = self.juices[0].Isum.size + self.to_pyarch["id"] = os.path.commonprefix(self.input_files) + self.to_pyarch["sample_name"] = self.juices[0].sample.name + if not self.input.get("no_ispyb"): + self.send_to_ispyb() + # self.output["icat"] = + self.send_to_icat() + + def teardown(self): + Plugin.teardown(self) + # logger.debug("HPLC.teardown") + # export the output file location + self.output["output_file"] = self.output_file + if self.nxs is not None: + self.nxs.close() + self.to_pyarch = None + self.ispyb = None + + def create_nexus(self): + nxs = Nexus(self.output_file, mode="w") + entry_grp = nxs.new_entry("entry", self.input.get("plugin_name", "dahu"), + title='BioSaxs Mesh experiment', + force_time=get_isotime()) + entry_grp["version"] = __version__ + nxs.h5.attrs["default"] = entry_grp.name.strip("/") + + # Configuration + cfg_grp = nxs.new_class(entry_grp, "configuration", "NXnote") + cfg_grp.create_dataset("data", data=json.dumps(self.input, indent=2, separators=(",\r\n", ":\t"))) + cfg_grp.create_dataset("format", data="text/json") + + # Process 0: Measurement group + input_grp = nxs.new_class(entry_grp, "0_measurement", "NXcollection") + input_grp["sequence_index"] = self.sequence_index() + + for idx, filename in enumerate(self.input_files): + juice = self.read_nexus(filename) + if juice is not None: + rel_path = os.path.relpath(os.path.abspath(filename), os.path.dirname(os.path.abspath(self.output_file))) + input_grp["LImA_%04i" % idx] = h5py.ExternalLink(rel_path, juice.h5path) + self.juices.append(juice) + + assert self.juices + q = self.juices[0].q + unit = self.juices[0].unit + radial_unit, unit_name = str(unit).split("_", 1) + + # Sample: outsourced ! + create_nexus_sample(nxs, entry_grp, self.juices[0].sample) + + # Process 1: Mesh + mesh_grp = nxs.new_class(entry_grp, "1_mesh", "NXprocess") + mesh_grp["sequence_index"] = self.sequence_index() + mesh_src = mesh_grp.create_group("sources") + input_dataset = ListDataSet() + if self.juices: + for idx, juice in enumerate(self.juices): + with h5py.File(juice.filename) as h: + grp = h[juice.h5path] + meas = grp["0_measurement"] + ds = meas["images"] + src = os.path.abspath(ds.file.filename) + name = ds.name + nframes, *img_shape = ds.shape + rel_path = os.path.relpath(src, os.path.dirname(os.path.abspath(self.output_file))) + mesh_src[f"images_{idx:04d}"] = h5py.ExternalLink(rel_path, name) + input_dataset.append(DataSet(src, name, nframes, img_shape)) + poni = PoniFile(juice.poni) + # mask = juice.mask + polarization = juice.polarization + method = juice.method + else: + poni = mask = energy = polarization = method = None + + nbin = q.size + # Creates a configuration NXnote in the NXProcess like diffmap would do""" + diffmap_grp = nxs.new_class(mesh_grp, "configuration", "NXnote") + diffmap_grp["type"] = "text/json" + worker = WorkerConfig(poni=poni, + nbpt_rad=nbin, + nbpt_azim=1) + worker.unit = unit + worker.method = (method.split, method.algorithm, method.implementation) + worker.opencl_device = method.target + worker.polarization_factor = polarization + + diffmap = DiffmapConfig(experiment_title="bm29.mesh", + slow_motor=MotorRange(start=self.scan.slow_motor_start, + stop=self.scan.slow_motor_stop, + points=self.scan.slow_motor_step+1, + name=self.scan.slow_motor_name), + fast_motor=MotorRange(start=self.scan.fast_motor_start, + stop=self.scan.fast_motor_stop, + points=self.scan.fast_motor_step+1, + name=self.scan.fast_motor_name), + offset=0, + zigzag_scan=self.scan.backnforth, + ai=worker, + input_data=input_dataset, + output_file=self.output_file) + # print(diffmap.as_dict()) + diffmap_grp.create_dataset("data", + data=json.dumps(diffmap.as_dict(), + indent=2, + separators=(",\r\n", ":\t"))) + + shape = self.scan.shape + (nbin,) + I = numpy.zeros(shape, dtype=numpy.float32) + slow = numpy.linspace(self.scan.slow_motor_start, + self.scan.slow_motor_stop, + self.scan.slow_motor_step+1, + endpoint=True).astype("float32") + fast = numpy.linspace(self.scan.fast_motor_start, + self.scan.fast_motor_stop, + self.scan.fast_motor_step+1, + endpoint=True).astype("float32") + + sigma = numpy.zeros(shape, dtype=numpy.float32) + Isum = numpy.zeros(self.scan.shape, dtype=numpy.float32) + indices = numpy.zeros(self.scan.shape, dtype=numpy.uint32) + + timestamps = [] + + for juice in self.juices: + timestamps.append(juice.timestamps) + # print(juice.idx) + for j, i in enumerate(juice.idx): + p = self.scan.get_pos(i) + if p is None: + continue + # print(p) + indices[p.slow, p.fast] = p.index + I[p.slow, p.fast] = juice.I[j] + Isum[p.slow, p.fast] = juice.Isum[j] + sigma[p.slow, p.fast] = juice.sigma[j] + + timestamps = self.to_pyarch["time"] = numpy.concatenate(timestamps) + + mesh_data = nxs.new_class(mesh_grp, "mesh", "NXdata") + mesh_data.attrs["title"] = "Mesh scan" + sum_ds = mesh_data.create_dataset("sum", data=Isum, dtype=numpy.float32) + sum_ds.attrs["interpretation"] = "image" + sum_ds.attrs["long_name"] = "Summed Intensity" + frame_ds = mesh_data.create_dataset("frame_ids", data=indices, dtype=numpy.uint32) + frame_ds.attrs["interpretation"] = "image" + frame_ds.attrs["long_name"] = "frame index" + slow_motor_ds = mesh_data.create_dataset("slow_motor", data =slow) + slow_motor_ds.attrs["interpretation"] = "spectrum" + slow_motor_ds.attrs["long_name"] = self.scan.slow_motor_name + fast_motor_ds = mesh_data.create_dataset("fast_motor", data =fast) + fast_motor_ds.attrs["interpretation"] = "spectrum" + fast_motor_ds.attrs["long_name"] = self.scan.fast_motor_name + mesh_data.attrs["signal"] = "sum" + mesh_data.attrs["axes"] = ["slow_motor", "fast_motor"] + + mesh_grp.attrs["default"] = posixpath.relpath(mesh_data.name, mesh_grp.name) + entry_grp.attrs["default"] = posixpath.relpath(mesh_data.name, entry_grp.name) + time_ds = mesh_data.create_dataset("timestamps", data=timestamps, dtype=numpy.uint32) + time_ds.attrs["interpretation"] = "spectrum" + time_ds.attrs["long_name"] = "Time stamps (s)" + + integration_data = nxs.new_class(mesh_grp, "result", "NXdata") + mesh_grp.attrs["title"] = str(self.juices[0].sample) + integration_data["map_ptr"] = frame_ds + int_ds = integration_data.create_dataset("I", data=numpy.ascontiguousarray(I, dtype=numpy.float32)) + std_ds = integration_data.create_dataset("errors", data=numpy.ascontiguousarray(sigma, dtype=numpy.float32)) + q_ds = integration_data.create_dataset("q", data=self.juices[0].q) + slow_motor_ds = integration_data.create_dataset("slow_motor", data =slow) + slow_motor_ds.attrs["interpretation"] = "spectrum" + slow_motor_ds.attrs["long_name"] = self.scan.slow_motor_name + fast_motor_ds = integration_data.create_dataset("fast_motor", data =fast) + fast_motor_ds.attrs["interpretation"] = "spectrum" + fast_motor_ds.attrs["long_name"] = self.scan.fast_motor_name + q_ds.attrs["interpretation"] = "spectrum" + q_ds.attrs["unit"] = unit_name + q_ds.attrs["long_name"] = "Scattering vector q (nm⁻¹)" + integration_data.attrs["signal"] = "I" + integration_data.attrs["axes"] = ["slow_motor", "fast_motor", "q"] + integration_data.attrs["SILX_style"] = SAXS_STYLE + + int_ds.attrs["interpretation"] = "spectrum" + int_ds.attrs["units"] = "arbitrary" + int_ds.attrs["long_name"] = "Intensity (absolute, normalized on water)" + # int_ds.attrs["uncertainties"] = "errors" This does not work + int_ds.attrs["scale"] = "log" + std_ds.attrs["interpretation"] = "spectrum" + + mesh_plot(Isum, + title="Mesh scan", + x_label=self.scan.fast_motor_name, + y_label=self.scan.slow_motor_name, + x_range=fast, + y_range=slow, + filename=os.path.join(self.ispyb.gallery, "mesh.png")) + # save_zip(os.path.splitext(self.output_file)[0]+".zip", + # self.juices[0], I, sigma) + + def read_nexus(self, filename): + "return some NexusJuice from a HDF5 file " + with Nexus(filename, "r") as nxsr: + entry_name = nxsr.h5.attrs["default"] + entry_grp = nxsr.h5[entry_name] + h5path = entry_grp.name + nxdata_grp = entry_grp[entry_grp.attrs["default"]] + # assert nxdata_grp.name.endswith("hplc") # we are reading HPLC data + signal = nxdata_grp.attrs["signal"] + axis = nxdata_grp.attrs["axes"] + Isum = nxdata_grp[signal][()] + idx = nxdata_grp[axis][()] + try: + integrated = nxdata_grp.parent["result"] + except KeyError: + integrated = nxdata_grp.parent["results"] + self.log_warning(f"Parsing old file {filename} !") + signal = integrated.attrs["signal"] + I = integrated[signal][()] + axes = integrated.attrs["axes"][-1] + q = integrated[axes][()] + sigma = integrated["errors"][()] + + npt = len(q) + unit = pyFAI.units.to_unit(axes + "_" + integrated[axes].attrs["units"]) + integration_grp = nxdata_grp.parent + poni = str(integration_grp["configuration/file_name"][()]).strip() + if not os.path.exists(poni): + poni = integration_grp["configuration/data"][()] + if isinstance(poni, bytes): + poni = poni.decode() + poni = json.loads(poni) + polarization = integration_grp["configuration/polarization_factor"][()] + method = IntegrationMethod.select_method(**json.loads(integration_grp["configuration/integration_method"][()]))[0] + instrument_grp = nxsr.get_class(entry_grp, class_type="NXinstrument")[0] + detector_grp = nxsr.get_class(instrument_grp, class_type="NXdetector")[0] + mask = detector_grp["pixel_mask"].attrs["filename"] + mono_grp = nxsr.get_class(instrument_grp, class_type="NXmonochromator")[0] + energy = mono_grp["energy"][()] + + # Read the sample description: + sample_grp = nxsr.get_class(entry_grp, class_type="NXsample")[0] + sample_name = posixpath.split(sample_grp.name)[-1] + + buffer = sample_grp["buffer"][()] if "buffer" in sample_grp else "" + concentration = sample_grp["concentration"][()] if "concentration" in sample_grp else "" + description = sample_grp["description"][()] if "description" in sample_grp else "" + hplc = sample_grp["hplc"][()] if "hplc" in sample_grp else "" + temperature = sample_grp["temperature"][()] if "temperature" in sample_grp else "" + temperature_env = sample_grp["temperature_env"][()] if "temperature_env" in sample_grp else "" + sample = Sample(sample_name, description, buffer, concentration, hplc, temperature_env, temperature) + meas_grp = nxsr.get_class(entry_grp, class_type="NXdata")[0] + timestamps = [] + for ts_name in ("timestamps", "time-stamps"): + if ts_name in meas_grp: + timestamps = meas_grp[ts_name][()] + break + return NexusJuice(filename, h5path, npt, unit, idx, Isum, q, I, sigma, poni, mask, energy, polarization, method, sample, timestamps) + "filename h5path npt unit idx Isum q I sigma poni mask energy polarization method sample timestamps" + + def send_to_ispyb(self): + self.log_warning("send_to_ispyb: unimplemented") + + def send_to_icat(self): + self.log_warning("send_to_icat: unimplemented") + diff --git a/plugins/bm29/meson.build b/plugins/bm29/meson.build index 3c08e41..8a7730b 100644 --- a/plugins/bm29/meson.build +++ b/plugins/bm29/meson.build @@ -7,7 +7,8 @@ py.install_sources( [ 'nexus.py', 'subtracte.py', 'memcached.py', - 'icat.py' + 'icat.py', + 'mesh.py' ], pure: false, # Will be installed next to binaries subdir: 'dahu/plugins/bm29' # Folder relative to site-packages to install to diff --git a/plugins/bm29/nexus.py b/plugins/bm29/nexus.py index 411f949..b19e677 100644 --- a/plugins/bm29/nexus.py +++ b/plugins/bm29/nexus.py @@ -70,6 +70,15 @@ def is_hdf5(filename): return sig == signature +def fully_qualified_name(o): + """Return the fully qualified name of the class""" + klass = o.__class__ + module = klass.__module__ + if module == 'builtins': + return klass.__qualname__ # avoid outputs like 'builtins.str' + return module + '.' + klass.__qualname__ + + class Nexus: """ Writer class to handle Nexus/HDF5 data @@ -88,7 +97,8 @@ class Nexus: def __init__(self, filename, mode=None, creator=None, timeout=None, - start_time=None): + start_time=None, + pure=False): """ Constructor @@ -97,6 +107,8 @@ def __init__(self, filename, mode=None, :param creator: set as attr of the NXroot :param timeout: retry for that amount of time (in seconds) :param start_time: set as attr of the NXroot + :param pure: use pure h5py mode. Unless, try to be clever when + accessing write-opened files (can breaks external links) """ self.filename = os.path.abspath(filename) self.mode = mode @@ -104,18 +116,18 @@ def __init__(self, filename, mode=None, logger.error("h5py module missing: NeXus not supported") raise RuntimeError("H5py module is missing") - pre_existing = os.path.exists(self.filename) or "w" in mode - if self.mode is None: - if pre_existing: - self.mode = "r" - else: - self.mode = "a" if timeout: end = time.perf_counter() + timeout while time.perf_counter() < end : + pre_existing = os.path.exists(self.filename) + if self.mode is None: + if pre_existing: + mode = "r" + else: + mode = "a" try: - if self.mode == "r": + if mode == "r" and not pure: self.file_handle = open(self.filename, mode="rb") self.h5 = h5py.File(self.file_handle, mode="r") else: @@ -125,24 +137,36 @@ def __init__(self, filename, mode=None, os.stat(os.path.dirname(self.filename)) time.sleep(1) else: + self.mode = mode break else: raise OSError(f"Unable to open HDF5 file {self.filename}") else: - if self.mode == "r": + pre_existing = os.path.exists(self.filename) + if self.mode is None: + if pre_existing: + self.mode = "r" + else: + self.mode = "a" + + if not pure and self.mode == "r" and h5py.version.version_tuple >= (2, 9): self.file_handle = open(self.filename, mode=self.mode + "b") self.h5 = h5py.File(self.file_handle, mode=self.mode) else: self.file_handle = None self.h5 = h5py.File(self.filename, mode=self.mode) self.to_close = [] - if not pre_existing: + + if not pre_existing or "w" in self.mode: self.h5.attrs["NX_class"] = "NXroot" self.h5.attrs["file_time"] = get_isotime(start_time) self.h5.attrs["file_name"] = self.filename self.h5.attrs["HDF5_Version"] = h5py.version.hdf5_version self.h5.attrs["creator"] = creator or self.__class__.__name__ + def __repr__(self): + return f"<{fully_qualified_name(self)} file on {self.h5}>" + def __del__(self): self.close() @@ -202,13 +226,14 @@ def get_entries(self): :return: list of HDF5 groups """ - entries = [(grp, from_isotime(self.h5[grp + "/start_time"][()])) - for grp in self.h5 - if isinstance(self.h5[grp], h5py.Group) and - ("start_time" in self.h5[grp]) and - self.get_attr(self.h5[grp], "NX_class") == "NXentry"] - entries.sort(key=lambda a: a[1], reverse=True) # sort entries in decreasing time - return [self.h5[i[0]] for i in entries] + entries = [(name, grp, from_isotime(self.h5[name + "/start_time"][()])) + for name, grp in self.h5.items() + if isinstance(grp, h5py.Group) and + ("start_time" in grp) and + self.get_attr(grp, "NX_class") == "NXentry"] + # print(entries) + entries.sort(key=lambda a: a[-1], reverse=True) # sort entries in decreasing time + return [i[1] for i in entries] def find_detector(self, all=False): """ @@ -243,11 +268,14 @@ def new_entry(self, entry="entry", program_name="pyFAI", if not force_name: nb_entries = len(self.get_entries()) entry = "%s_%04i" % (entry, nb_entries) - entry_grp = self.h5.require_group(entry) + entry_grp = self.h5 + for i in entry.split("/"): + if i: + entry_grp = entry_grp.require_group(i) self.h5.attrs["default"] = entry_grp.name.strip("/") entry_grp.attrs["NX_class"] = "NXentry" entry_grp["title"] = str(title) - entry_grp["program_name"] = program_name + entry_grp["program_name"] = str(program_name) if isinstance(force_time, str): entry_grp["start_time"] = force_time else: diff --git a/plugins/bm29/subtracte.py b/plugins/bm29/subtracte.py index 4ea0cbe..5a19885 100644 --- a/plugins/bm29/subtracte.py +++ b/plugins/bm29/subtracte.py @@ -11,14 +11,15 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "27/05/2025" +__date__ = "25/06/2025" __status__ = "development" -__version__ = "0.3.0" +__version__ = "0.4.0" import os import posixpath import json import copy +import zipfile from math import log, pi from collections import namedtuple from urllib3.util import parse_url @@ -33,12 +34,13 @@ logger.error("Numexpr is not installed, falling back on numpy's implementations") numexpr = None import h5py -import pyFAI, pyFAI.azimuthalIntegrator +import pyFAI, pyFAI.integrator.azimuthal from pyFAI.containers import Integrate1dResult from pyFAI.method_registry import IntegrationMethod import freesas, freesas.cormap, freesas.invariants from freesas.autorg import auto_gpa, autoRg, auto_guinier from freesas.bift import BIFT +from freesas.app.extract_ascii import write_ascii from scipy.optimize import minimize from .common import Sample, Ispyb, get_equivalent_frames, cmp_float, get_integrator, KeyCache, \ polarization_factor, method, Nexus, get_isotime, SAXS_STYLE, NORMAL_STYLE, \ @@ -47,7 +49,72 @@ from .memcached import to_memcached from .icat import send_icat -NexusJuice = namedtuple("NexusJuice", "filename h5path npt unit q I sigma poni mask energy polarization method signal2d error2d normalization sample") + +NexusJuice = namedtuple("NexusJuice", "filename h5path npt unit " + "q I sigma poni mask energy polarization method signal2d " + "error2d normalization sample " + "I_all, sigma_all") + + +def save_zip(filename, sample_juice, buffer_juices): + """Save a stack of I into a zipfile with each frames in a dat-file. + + :param filename: name of the zip-file + :param sample_juice: + :param buffer_juices: list of buffer juice + :return: nothing + """ + basename = os.path.basename(filename) + base = os.path.splitext(basename)[0] + destz_sample = "sample/" + destz_buffer = "buffer_%1i/" + common = {"q": sample_juice.q} + if sample_juice.sample: + sample = sample_juice.sample + if sample.name: + common["sample"]: sample.name + destz_sample += sample.name + else: + destz_sample += "sample" + + if sample.buffer: + common["buffer"] = sample.buffer + destz_buffer += sample.buffer if isinstance(sample.buffer, str) else sample.buffer.decode() + else: + destz_buffer += "buffer" + + if sample.temperature_env: + common["storage temperature"] = sample.temperature_env + + if sample.temperature: + common["exposure temperature"] = sample.temperature + + if sample.concentration: + common["concentration"] = sample.concentration + destz_sample += "_%04i.dat" + destz_buffer += "_%04i.dat" + res = {} + # sample + idx = 0 + for i, s in zip(sample_juice.I_all, sample_juice.sigma_all): + r = copy.copy(common) + r["I"] = i + r["std"] = s + res[destz_sample % idx] = r + idx+=1 + # buffers + for buffer_idx, buffer in enumerate(buffer_juices): + idx = 0 + for i, s in zip(buffer.I_all, buffer.sigma_all): + r = copy.copy(common) + r["I"] = i + r["std"] = s + res[destz_buffer % (buffer_idx, idx)] = r + idx+=1 + + with zipfile.ZipFile(filename, "w") as z: + for name, frame in res.items(): + z.writestr(name, write_ascii(frame)) class SubtractBuffer(Plugin): @@ -209,7 +276,7 @@ def create_nexus(self): entry_grp = nxs.new_entry("entry", self.input.get("plugin_name", "dahu"), title='BioSaxs buffer subtraction', force_time=get_isotime()) - nxs.h5.attrs["default"] = entry_grp.name.strip["/"] + nxs.h5.attrs["default"] = entry_grp.name.strip("/") # Configuration cfg_grp = nxs.new_class(entry_grp, "configuration", "NXnote") @@ -232,13 +299,18 @@ def create_nexus(self): # Sample: outsourced ! create_nexus_sample(nxs, entry_grp, self.sample_juice.sample) + #save input curves as zipfile: TODO Check that this is is working: + save_zip(os.path.splitext(self.output_file)[0]+".zip", + self.sample_juice, + self.buffer_juices) + # Process 1: CorMap cormap_grp = nxs.new_class(entry_grp, "1_correlation_mapping", "NXprocess") cormap_grp["sequence_index"] = 1 cormap_grp["program"] = "freesas.cormap" cormap_grp["version"] = freesas.version cormap_grp["date"] = get_isotime() - cormap_data = nxs.new_class(cormap_grp, "results", "NXdata") + cormap_data = nxs.new_class(cormap_grp, "result", "NXdata") cormap_data.attrs["SILX_style"] = NORMAL_STYLE cfg_grp = nxs.new_class(cormap_grp, "configuration", "NXcollection") @@ -278,7 +350,7 @@ def create_nexus(self): average_grp["sequence_index"] = 2 average_grp["program"] = fully_qualified_name(self.__class__) average_grp["version"] = __version__ - average_data = nxs.new_class(average_grp, "results", "NXdata") + average_data = nxs.new_class(average_grp, "result", "NXdata") average_data.attrs["SILX_style"] = SAXS_STYLE average_data.attrs["signal"] = "intensity_normed" # Stage 2 processing @@ -341,7 +413,7 @@ def create_nexus(self): ai2_grp["version"] = pyFAI.version ai2_grp["date"] = get_isotime() radial_unit, unit_name = str(key_cache.unit).split("_", 1) - ai2_data = nxs.new_class(ai2_grp, "results", "NXdata") + ai2_data = nxs.new_class(ai2_grp, "result", "NXdata") ai2_data.attrs["SILX_style"] = SAXS_STYLE ai2_data.attrs["title"] = "%s, subtracted" % self.sample_juice.sample.name ai2_data.attrs["signal"] = "I" @@ -404,7 +476,7 @@ def create_nexus(self): guinier_autorg = nxs.new_class(guinier_grp, "autorg", "NXcollection") guinier_gpa = nxs.new_class(guinier_grp, "gpa", "NXcollection") guinier_guinier = nxs.new_class(guinier_grp, "guinier", "NXcollection") - guinier_data = nxs.new_class(guinier_grp, "results", "NXdata") + guinier_data = nxs.new_class(guinier_grp, "result", "NXdata") guinier_data.attrs["SILX_style"] = NORMAL_STYLE guinier_data.attrs["title"] = "Guinier analysis" # Stage4 processing: autorg and auto_gpa @@ -538,7 +610,7 @@ def create_nexus(self): kratky_grp["program"] = "freesas.autorg" kratky_grp["version"] = freesas.version kratky_grp["date"] = get_isotime() - kratky_data = nxs.new_class(kratky_grp, "results", "NXdata") + kratky_data = nxs.new_class(kratky_grp, "result", "NXdata") kratky_data.attrs["SILX_style"] = NORMAL_STYLE kratky_data.attrs["title"] = "Dimensionless Kratky plots" kratky_grp.attrs["default"] = posixpath.relpath(kratky_data.name, kratky_grp.name) @@ -554,21 +626,21 @@ def create_nexus(self): qRg_ds.attrs["long_name"] = "q·Rg (unit-less)" #Nota the "/" hereafter is chr(8725), the division sign and not the usual slash - k_ds = kratky_data.create_dataset("q2Rg2I∕I0", data=ydata.astype(numpy.float32)) + k_ds = kratky_data.create_dataset("q2Rg2I÷I0", data=ydata.astype(numpy.float32)) k_ds.attrs["interpretation"] = "spectrum" k_ds.attrs["long_name"] = "q²Rg²I(q)/I₀" ke_ds = kratky_data.create_dataset("errors", data=dy.astype(numpy.float32)) ke_ds.attrs["interpretation"] = "spectrum" kratky_data_attrs = kratky_data.attrs - kratky_data_attrs["signal"] = "q2Rg2I∕I0" - kratky_data_attrs["axes"] = "qRg" + kratky_data_attrs["signal"] = k_ds.name + kratky_data_attrs["axes"] = qRg_ds.name # stage 6: Rambo-Tainer invariant rti_grp = nxs.new_class(entry_grp, "6_invariants", "NXprocess") rti_grp["sequence_index"] = 6 rti_grp["program"] = "freesas.invariants" rti_grp["version"] = freesas.version - rti_data = nxs.new_class(rti_grp, "results", "NXdata") + rti_data = nxs.new_class(rti_grp, "result", "NXdata") # average_data.attrs["SILX_style"] = SAXS_STYLE # average_data.attrs["signal"] = "intensity_normed" # Rambo_Tainer @@ -606,7 +678,7 @@ def create_nexus(self): bift_grp["program"] = "freesas.bift" bift_grp["version"] = freesas.version bift_grp["date"] = get_isotime() - bift_data = nxs.new_class(bift_grp, "results", "NXdata") + bift_data = nxs.new_class(bift_grp, "result", "NXdata") bift_data.attrs["SILX_style"] = NORMAL_STYLE bift_data.attrs["title"] = "Pair distance distribution function p(r)" @@ -690,7 +762,7 @@ def read_nexus(filename): with Nexus(filename, "r") as nxsr: entry_grp = nxsr.get_entries()[0] h5path = entry_grp.name - nxdata_grp = nxsr.h5[entry_grp.attrs["default"]] + nxdata_grp = entry_grp[entry_grp.attrs["default"]] signal = nxdata_grp.attrs["signal"] axis = nxdata_grp.attrs["axes"] I = nxdata_grp[signal][()] @@ -700,13 +772,9 @@ def read_nexus(filename): unit = pyFAI.units.to_unit(axis + "_" + nxdata_grp[axis].attrs["units"]) integration_grp = nxdata_grp.parent poni = integration_grp["configuration/file_name"][()] - if isinstance(poni, bytes): - poni = poni.decode() - else: - poni = str(poni) - poni = poni.strip() + poni = str_(poni).strip() if not os.path.exists(poni): - poni = str(integration_grp["configuration/data"][()]).strip() + poni = str_(integration_grp["configuration/data"][()]).strip() polarization = integration_grp["configuration/polarization_factor"][()] method = IntegrationMethod.select_method(**json.loads(integration_grp["configuration/integration_method"][()]))[0] instrument_grp = nxsr.get_class(entry_grp, class_type="NXinstrument")[0] @@ -722,7 +790,7 @@ def read_nexus(filename): sample_grp = nxsr.get_class(entry_grp, class_type="NXsample")[0] sample_name = posixpath.basename(sample_grp.name) - buffer = sample_grp["buffer"][()] if "buffer" in sample_grp else "" + buffer = str_(sample_grp["buffer"][()] if "buffer" in sample_grp else "") concentration = sample_grp["concentration"][()] if "concentration" in sample_grp else "" description = sample_grp["description"][()] if "description" in sample_grp else "" hplc = sample_grp["hplc"][()] if "hplc" in sample_grp else "" @@ -730,7 +798,15 @@ def read_nexus(filename): temperature_env = sample_grp["temperature_env"][()] if "temperature_env" in sample_grp else "" sample = Sample(sample_name, description, buffer, concentration, hplc, temperature_env, temperature) - return NexusJuice(filename, h5path, npt, unit, q, I, sigma, poni, mask, energy, polarization, method, image2d, error2d, norm, sample) + if "1_integration" in entry_grp: + I_all = entry_grp["1_integration/result/I"][()] + sigma_all = entry_grp["1_integration/result/errors"][()] + else: + I_all = [] + sigma_all = [] + + return NexusJuice(filename, h5path, npt, unit, q, I, sigma, poni, mask, energy, polarization, + method, image2d, error2d, norm, sample, I_all, sigma_all) def send_to_ispyb(self): if self.ispyb.url and parse_url(self.ispyb.url).host: @@ -753,6 +829,7 @@ def send_to_icat(self): raw=raw, path=os.path.dirname(os.path.abspath(self.output_file)), data=to_icat, + dataset="subtraction", gallery=self.ispyb.gallery or os.path.join(os.path.dirname(os.path.abspath(self.output_file)), "gallery"), metadata=metadata) @@ -766,3 +843,8 @@ def send_to_memcached(self): return to_memcached(dico) +def str_(smth): + if isinstance(smth, bytes): + return smth.decode() + else: + return str(smth) diff --git a/src/dahu/server.py b/src/dahu/server.py index 3d1ca6f..885e528 100644 --- a/src/dahu/server.py +++ b/src/dahu/server.py @@ -3,7 +3,7 @@ from __future__ import with_statement, print_function, absolute_import, division """ -Data Analysis RPC server over Tango: +Data Analysis RPC server over Tango: Tango device server """ @@ -11,7 +11,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "09/07/2021" +__date__ = "04/06/2025" __status__ = "production" __docformat__ = 'restructuredtext' @@ -120,7 +120,7 @@ def listPlugins(self): res = ["List of all plugin currently loaded (use initPlugin to loaded additional plugins):"] plugins = list(plugin_factory.registry.keys()) plugins.sort() - return os.linesep.join(res + [" %s : %s" % (i, plugin_factory.registry[i].__doc__.split("\n")[0]) for i in plugins]) + return os.linesep.join(res + [f' {i} : {plugin_factory.registry[i].__doc__.split("\n")[0]}' for i in plugins]) def initPlugin(self, name): """ @@ -134,9 +134,9 @@ def initPlugin(self, name): err = "plugin %s failed to be instanciated: %s" % (name, error) logger.error(err) if plugin is None or err: - return "Plugin not found: %s, err" % (name, err) + return f"Plugin not found: {name}, {err}" else: - return "Plugin loaded: %s%s%s" % (name, os.linesep, plugin.__doc__) + return f"Plugin loaded: {name}{os.linesep}{plugin.__doc__}" def abort(self, jobId): """ @@ -283,7 +283,7 @@ def waitJob(self, jobId): Wait for a job to be finished and returns the status. May cause Tango timeout if too slow to finish .... May do polling to wait the job actually started - + @param jobId: identifier of the job (int) @return: status of the job """ diff --git a/version.py b/version.py index 13364f1..582e92e 100755 --- a/version.py +++ b/version.py @@ -64,7 +64,7 @@ "final": 15} MAJOR = 2026 -MINOR = 1 +MINOR = 2 MICRO = 0 RELEV = "dev" # <16 SERIAL = 0 # <16