From f96b8639073a983e59d8f4ebb290e197c95017a8 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 29 Feb 2024 16:27:43 +0000 Subject: [PATCH 01/12] Add function to produce plots for IntegratePeaksMD results Could not make use of sliceviewer code as can import from mantidqt --- scripts/Diffraction/single_crystal/base_sx.py | 146 +++++++++++++++++- 1 file changed, 144 insertions(+), 2 deletions(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index 7c530cc6c602..3e2e87e162a6 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -9,11 +9,15 @@ from enum import Enum import mantid.simpleapi as mantid from mantid.api import FunctionFactory, AnalysisDataService as ADS -from mantid.kernel import logger +from mantid.kernel import logger, SpecialCoordinateSystem from FindGoniometerFromUB import getSignMaxAbsValInCol from mantid.geometry import CrystalStructure, SpaceGroupFactory, ReflectionGenerator, ReflectionConditionFilter from os import path - +from json import loads as json_loads +from matplotlib.pyplot import subplots, close +from matplotlib.patches import Circle +from matplotlib.colors import LogNorm +from matplotlib.backends.backend_pdf import PdfPages from abc import ABC, abstractmethod @@ -614,6 +618,144 @@ def remove_non_integrated_peaks(peaks): EnableLogging=False, ) + @staticmethod + def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, log_norm=True): + """ + :param wsMD: MD workspace to plot + :param peaks: integrated peaks using IntegratePeaksMD + :param filename: filename to store pdf output + :param nbins: number of bins along major radius of ellipsoid + :param extent: extent in units of largest outer background radius + :param log_norm: use log normalisation in colorscale + """ + wsMD = BaseSX.retrieve(wsMD) + peaks = BaseSX.retrieve(peaks) + + # find appropriate getter for peak centre given MD frame + frame = wsMD.getSpecialCoordinateSystem() + frame_to_peak_centre_attr = "getQLabFrame" + if frame == SpecialCoordinateSystem.QSample: + frame_to_peak_centre_attr = "getQSampleFrame" + elif frame == SpecialCoordinateSystem.HKL: + frame_to_peak_centre_attr = "getHKL" + + # loop over peaks and plot + try: + with PdfPages(filename) as pdf: + ws_cut = None + for ipk, pk in enumerate(peaks): + peak_shape = pk.getPeakShape() + print(peak_shape.shapeName().lower()) + if peak_shape.shapeName().lower() == "none": + continue + ws_cut, radii, bg_inner_radii, bg_outer_radii = BaseSX._bin_MD_around_peak( + wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_centre_attr + ) + + fig, axes = subplots(1, 3, figsize=(12, 4), subplot_kw={"projection": "mantid"}) + for iax, ax in enumerate(axes): + dims = list(range(3)) + dims.pop(iax) + # plot slice + im = ax.imshow(ws_cut.getSignalArray().sum(axis=iax)[:, ::-1].T) + im.set_extent([-1, 1, -1, 1]) # so that ellipsoid is a circle + if log_norm: + im.set_norm(LogNorm()) + # plot peak representation (circle given extents) + patch = Circle((0, 0), 1 / scale, facecolor="none", edgecolor=3 * [0.7]) # outer radius + ax.add_patch(patch) + patch = Circle( + (0, 0), bg_inner_radii[imax] / (bg_outer_radii[imax] * scale), facecolor="none", edgecolor=3 * [0.7], ls="--" + ) # inner radius + ax.add_patch(patch) + patch = Circle( + (0, 0), radii[imax] / (bg_outer_radii[imax] * scale), facecolor="none", edgecolor="r", ls="--" + ) # radius + ax.add_patch(patch) + # format axes + ax.set_aspect("equal") + ax.set_xlabel(ws_cut.getDimension(dims[0]).name) + ax.set_ylabel(ws_cut.getDimension(dims[1]).name) + fig.suptitle( + f"{ipk} ({','.join(str(np.round(pk.getHKL(), 2))[1:-1].split())})" + f" $I/\\sigma$={np.round(pk.getIntensityOverSigma(), 2)}\n" + rf"$\lambda$={np.round(pk.getWavelength(), 2)} $\AA$; " + rf"$2\theta={np.round(np.degrees(pk.getScattering()), 1)}^\circ$; " + rf"d={np.round(pk.getDSpacing(), 2)} $\AA$" + ) + fig.tight_layout() + pdf.savefig(fig) + close(fig) + + if ws_cut is not None: + mantid.DeleteWorkspace(ws_cut) + except OSError: + raise RuntimeError( + f"OutputFile ({output_file}) could not be opened - please check it is not open by " + f"another programme and that the user has permission to write to that directory." + ) + + @staticmethod + def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_centre_attr): + shape_info = json_loads(peak_shape.toJSON()) + if peak_shape.shapeName().lower() == "spherical": + BaseSX.convert_spherical_representation_to_ellipsoid(shape_info) + # get radii + radii = np.array([shape_info[f"radius{iax}"] for iax in range(3)]) + bg_inner_radii = np.array([shape_info[f"background_inner_radius{iax}"] for iax in range(3)]) + bg_outer_radii = np.array([shape_info[f"background_outer_radius{iax}"] for iax in range(3)]) + isort = np.argsort(radii) + imax = isort[-1] + # eignevectors of ellipsoid + evecs = np.zeros((3, 3)) + for iax in range(len(radii)): + evec = np.array([float(elem) for elem in shape_info[f"direction{iax}"].split()]) + evecs[:, iax] = evec / np.linalg.norm(evec) + # center and translation + translation = np.array([shape_info[f"translation{iax}"] for iax in range(3)]) + cen = getattr(pk, frame_to_peak_centre_attr)() + cen = np.matmul(evecs.T, np.array(cen) + translation) + # get extents + extents = np.vstack((cen - extent * bg_outer_radii, cen + extent * bg_outer_radii)) + # get nbins along each axis + nbins = [int(nbins_max * radii[iax] / radii[imax]) for iax in range(len(radii))] + # call BinMD + ws_cut = BinMD( + InputWorkspace=ws, + OutputWorkspace="__ws_cut", + AxisAligned=False, + BasisVector0=r"Q$_0$,unit," + ",".join(np.array2string(evecs[:, 0], precision=6).strip("[]").split()), + BasisVector1=r"Q$_1$,unit," + ",".join(np.array2string(evecs[:, 1], precision=6).strip("[]").split()), + BasisVector2=r"Q$_2$,unit," + ",".join(np.array2string(evecs[:, 2], precision=6).strip("[]").split()), + OutputExtents=extents.flatten(order="F"), + OutputBins=nbins, + EnableLogging=False, + ) + return ws_cut, radii, bg_inner_radii, bg_outer_radii + + @staticmethod + def _convert_spherical_representation_to_ellipsoid(shape_info): + # copied from mantidqt.widgets.sliceviewer.peaksviewer.representation.ellipsoid - can't import here though + # convert shape_info dict from sphere to ellipsoid for plotting + for key in ["radius", "background_inner_radius", "background_outer_radius"]: + if key in shape_info: + shape_info[f"{key}{0}"] = shape_info.pop(key) + else: + shape_info[f"{key}{0}"] = 0.0 # null value + for idim in [1, 2]: + shape_info[f"{key}{idim}"] = shape_info[f"{key}{0}"] + # add axes along basis vecs of frame and set 0 translation + shape_info.update( + { + "direction0": "1 0 0", + "direction1": "0 1 0", + "direction2": "0 0 1", + "translation0": 0.0, + "translation1": 0.0, + "translation2": 0.0, + } + ) + @staticmethod def retrieve(ws): if isinstance(ws, str): From b59a820fa7f12579106eef96f61403d5fbe68563 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 29 Feb 2024 16:41:44 +0000 Subject: [PATCH 02/12] Fix typo in variable name in OSError branch --- scripts/Diffraction/single_crystal/base_sx.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index 3e2e87e162a6..0ad28be5dc7e 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -691,7 +691,7 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo mantid.DeleteWorkspace(ws_cut) except OSError: raise RuntimeError( - f"OutputFile ({output_file}) could not be opened - please check it is not open by " + f"OutputFile ({filename}) could not be opened - please check it is not open by " f"another programme and that the user has permission to write to that directory." ) From 07acad2ee23bc82d390a329c52f250abe2d247ea Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 29 Feb 2024 17:04:43 +0000 Subject: [PATCH 03/12] Fix bug for integrated peaks without background shells --- scripts/Diffraction/single_crystal/base_sx.py | 33 ++++++++++--------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index 0ad28be5dc7e..dc0c84bed704 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -645,10 +645,9 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo ws_cut = None for ipk, pk in enumerate(peaks): peak_shape = pk.getPeakShape() - print(peak_shape.shapeName().lower()) if peak_shape.shapeName().lower() == "none": continue - ws_cut, radii, bg_inner_radii, bg_outer_radii = BaseSX._bin_MD_around_peak( + ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths = BaseSX._bin_MD_around_peak( wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_centre_attr ) @@ -657,21 +656,23 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo dims = list(range(3)) dims.pop(iax) # plot slice - im = ax.imshow(ws_cut.getSignalArray().sum(axis=iax)[:, ::-1].T) + im = ax.imshow(ws_cut.getSignalArray().sum(axis=iax)[:, ::-1].T) # reshape to match sliceviewer im.set_extent([-1, 1, -1, 1]) # so that ellipsoid is a circle if log_norm: im.set_norm(LogNorm()) # plot peak representation (circle given extents) - patch = Circle((0, 0), 1 / scale, facecolor="none", edgecolor=3 * [0.7]) # outer radius - ax.add_patch(patch) - patch = Circle( - (0, 0), bg_inner_radii[imax] / (bg_outer_radii[imax] * scale), facecolor="none", edgecolor=3 * [0.7], ls="--" - ) # inner radius - ax.add_patch(patch) - patch = Circle( - (0, 0), radii[imax] / (bg_outer_radii[imax] * scale), facecolor="none", edgecolor="r", ls="--" - ) # radius + patch = Circle((0, 0), radii[imax] / box_lengths[imax], facecolor="none", edgecolor="r", ls="--") ax.add_patch(patch) + # add background shell + if np.all(bg_outer_radii > 1e-10): + patch = Circle( + (0, 0), bg_outer_radii[imax] / box_lengths[imax], facecolor="none", edgecolor=3 * [0.7] + ) # outer radius + ax.add_patch(patch) + patch = Circle( + (0, 0), bg_inner_radii[imax] / box_lengths[imax], facecolor="none", edgecolor=3 * [0.7], ls="--" + ) # inner radius + ax.add_patch(patch) # format axes ax.set_aspect("equal") ax.set_xlabel(ws_cut.getDimension(dims[0]).name) @@ -716,11 +717,13 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c cen = getattr(pk, frame_to_peak_centre_attr)() cen = np.matmul(evecs.T, np.array(cen) + translation) # get extents - extents = np.vstack((cen - extent * bg_outer_radii, cen + extent * bg_outer_radii)) + box_lengths = bg_outer_radii if np.all(bg_outer_radii > 1e-10) else radii + box_lengths = box_lengths * extent + extents = np.vstack((cen - box_lengths, cen + box_lengths)) # get nbins along each axis nbins = [int(nbins_max * radii[iax] / radii[imax]) for iax in range(len(radii))] # call BinMD - ws_cut = BinMD( + ws_cut = mantid.BinMD( InputWorkspace=ws, OutputWorkspace="__ws_cut", AxisAligned=False, @@ -731,7 +734,7 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c OutputBins=nbins, EnableLogging=False, ) - return ws_cut, radii, bg_inner_radii, bg_outer_radii + return ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths @staticmethod def _convert_spherical_representation_to_ellipsoid(shape_info): From a157290350e27f0087aa938066c07cd406688d05 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 29 Feb 2024 17:14:00 +0000 Subject: [PATCH 04/12] Ensure minimum of 3 bins inside peak radius along each dim --- scripts/Diffraction/single_crystal/base_sx.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index dc0c84bed704..0e7352cce6d7 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -722,6 +722,10 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c extents = np.vstack((cen - box_lengths, cen + box_lengths)) # get nbins along each axis nbins = [int(nbins_max * radii[iax] / radii[imax]) for iax in range(len(radii))] + # ensure minimum of 3 bins inside radius along each dim + min_nbins_in_radius = np.min(nbins * radii / box_lengths) + if min_nbins_in_radius < 3: + nbins = nbins * 3 / min_nbins_in_radius # call BinMD ws_cut = mantid.BinMD( InputWorkspace=ws, @@ -731,7 +735,7 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c BasisVector1=r"Q$_1$,unit," + ",".join(np.array2string(evecs[:, 1], precision=6).strip("[]").split()), BasisVector2=r"Q$_2$,unit," + ",".join(np.array2string(evecs[:, 2], precision=6).strip("[]").split()), OutputExtents=extents.flatten(order="F"), - OutputBins=nbins, + OutputBins=nbins.astype(int), EnableLogging=False, ) return ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths From 44d21208a6e3e59db1057838b322e857f59f044f Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 1 Mar 2024 09:15:58 +0000 Subject: [PATCH 05/12] Add unit test for new BaseSX method --- .../single_crystal/test/test_base_sx.py | 40 +++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/scripts/Diffraction/single_crystal/test/test_base_sx.py b/scripts/Diffraction/single_crystal/test/test_base_sx.py index 605728f90258..6b2b56b851a2 100644 --- a/scripts/Diffraction/single_crystal/test/test_base_sx.py +++ b/scripts/Diffraction/single_crystal/test/test_base_sx.py @@ -5,14 +5,19 @@ LoadEmptyInstrument, AnalysisDataService, CreatePeaksWorkspace, + CreateMDWorkspace, SetUB, IndexPeaks, CreateSampleWorkspace, SetSample, TransformHKL, + IntegratePeaksMD, + FakeMDEventData, ) from mantid.dataobjects import Workspace2D from Diffraction.single_crystal.base_sx import BaseSX +import tempfile +import shutil base_sx_path = "Diffraction.single_crystal.base_sx" @@ -23,10 +28,12 @@ def setUpClass(cls): cls.ws = LoadEmptyInstrument(InstrumentName="SXD", OutputWorkspace="empty") axis = cls.ws.getAxis(0) axis.setUnit("TOF") + cls._test_dir = tempfile.mkdtemp() @classmethod def tearDownClass(cls): AnalysisDataService.clear() + shutil.rmtree(cls._test_dir) def test_retrieve(self): self.assertTrue(isinstance(BaseSX.retrieve("empty"), Workspace2D)) # ADS name of self.ws is "empty" @@ -151,6 +158,39 @@ def test_make_UB_consistent(self): allclose(peaks_ref.sample().getOrientedLattice().getUB().tolist(), peaks.sample().getOrientedLattice().getUB().tolist()) ) + def test_plot_integrated_peaks_MD(self): + # make fake dataset + ws = CreateMDWorkspace( + Dimensions="3", + Extents="-5,5,-5,5,-5,5", + Names="H,K,L", + Units="r.l.u.,r.l.u.,r.l.u.", + Frames="HKL,HKL,HKL", + SplitInto="2", + SplitThreshold="50", + ) + # generate fake data at (0,1,0) + FakeMDEventData(ws, EllipsoidParams="1e+03,0,1,0,1,0,0,0,1,0,0,0,1,0.01,0.04,0.02,1", RandomSeed="3873875") + # create peak at (0,1,0) and integrate + peaks = self._make_peaks_HKL(hs=[0], ks=[1], ls=[0], wsname="peaks_md") + peaks_int = IntegratePeaksMD( + InputWorkspace=ws, + PeakRadius=0.6, + BackgroundInnerRadius=0.6, + BackgroundOuterRadius=0.8, + PeaksWorkspace=peaks, + OutputWorkspace="peaks_int_sphere", + Ellipsoid=False, + UseOnePercentBackgroundCorrection=False, + ) + + # plot + out_file = path.join(self._test_dir, "out_plot_MD.pdf") + BaseSX.plot_integrated_peaks_MD(ws, peaks_int, out_file, nbins_max=21, extent=1.5, log_norm=True) + + # check output file saved + self.assertTrue(path.exists(out_file)) + # --- helper funcs --- def _make_peaks_HKL(self, hs=None, ks=[0], ls=[1], wsname="peaks_hkl"): From e349b70a604337a7e0808c22545e63956af0b2f1 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 1 Mar 2024 09:45:44 +0000 Subject: [PATCH 06/12] Add release note --- .../v6.10.0/Diffraction/Single_Crystal/New_features/36947.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/source/release/v6.10.0/Diffraction/Single_Crystal/New_features/36947.rst diff --git a/docs/source/release/v6.10.0/Diffraction/Single_Crystal/New_features/36947.rst b/docs/source/release/v6.10.0/Diffraction/Single_Crystal/New_features/36947.rst new file mode 100644 index 000000000000..3eb064040ada --- /dev/null +++ b/docs/source/release/v6.10.0/Diffraction/Single_Crystal/New_features/36947.rst @@ -0,0 +1 @@ +- Add method ``plot_integrated_peaks_MD`` to ``BaseSX`` to plot result of IntegratePeaksMD and save in pdf \ No newline at end of file From f5b7fe89df8c181b8748859cc7b0d5abd787a495 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 1 Mar 2024 12:19:52 +0000 Subject: [PATCH 07/12] Add marker for peak centre And fix bug by returning imax --- scripts/Diffraction/single_crystal/base_sx.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index 0e7352cce6d7..67f9dfa31319 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -647,7 +647,7 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo peak_shape = pk.getPeakShape() if peak_shape.shapeName().lower() == "none": continue - ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths = BaseSX._bin_MD_around_peak( + ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths, imax = BaseSX._bin_MD_around_peak( wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_centre_attr ) @@ -660,6 +660,8 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo im.set_extent([-1, 1, -1, 1]) # so that ellipsoid is a circle if log_norm: im.set_norm(LogNorm()) + # plot peak position + ax.plot(0, 0, "xr") # plot peak representation (circle given extents) patch = Circle((0, 0), radii[imax] / box_lengths[imax], facecolor="none", edgecolor="r", ls="--") ax.add_patch(patch) @@ -674,6 +676,8 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo ) # inner radius ax.add_patch(patch) # format axes + ax.set_xlim(-1, 1) + ax.set_ylim(-1, 1) ax.set_aspect("equal") ax.set_xlabel(ws_cut.getDimension(dims[0]).name) ax.set_ylabel(ws_cut.getDimension(dims[1]).name) @@ -721,14 +725,14 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c box_lengths = box_lengths * extent extents = np.vstack((cen - box_lengths, cen + box_lengths)) # get nbins along each axis - nbins = [int(nbins_max * radii[iax] / radii[imax]) for iax in range(len(radii))] + nbins = np.array([int(nbins_max * radii[iax] / radii[imax]) for iax in range(len(radii))]) # ensure minimum of 3 bins inside radius along each dim min_nbins_in_radius = np.min(nbins * radii / box_lengths) if min_nbins_in_radius < 3: nbins = nbins * 3 / min_nbins_in_radius # call BinMD ws_cut = mantid.BinMD( - InputWorkspace=ws, + InputWorkspace=wsMD, OutputWorkspace="__ws_cut", AxisAligned=False, BasisVector0=r"Q$_0$,unit," + ",".join(np.array2string(evecs[:, 0], precision=6).strip("[]").split()), @@ -738,7 +742,7 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c OutputBins=nbins.astype(int), EnableLogging=False, ) - return ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths + return ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths, imax @staticmethod def _convert_spherical_representation_to_ellipsoid(shape_info): From bbce042097193a74d589062575434b94ee80a2b7 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 6 Mar 2024 16:01:22 +0000 Subject: [PATCH 08/12] Split up large function, add doc string and avoid hard-coding 3 --- scripts/Diffraction/single_crystal/base_sx.py | 62 +++++++++++++------ 1 file changed, 42 insertions(+), 20 deletions(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index 67f9dfa31319..9c040bd56231 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -702,22 +702,25 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo @staticmethod def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_centre_attr): - shape_info = json_loads(peak_shape.toJSON()) - if peak_shape.shapeName().lower() == "spherical": - BaseSX.convert_spherical_representation_to_ellipsoid(shape_info) - # get radii - radii = np.array([shape_info[f"radius{iax}"] for iax in range(3)]) - bg_inner_radii = np.array([shape_info[f"background_inner_radius{iax}"] for iax in range(3)]) - bg_outer_radii = np.array([shape_info[f"background_outer_radius{iax}"] for iax in range(3)]) - isort = np.argsort(radii) - imax = isort[-1] - # eignevectors of ellipsoid - evecs = np.zeros((3, 3)) - for iax in range(len(radii)): - evec = np.array([float(elem) for elem in shape_info[f"direction{iax}"].split()]) - evecs[:, iax] = evec / np.linalg.norm(evec) - # center and translation - translation = np.array([shape_info[f"translation{iax}"] for iax in range(3)]) + """ + Bin MD workspace in peak region with projection axes given by the axes of the ellipsoid shape + :param wsMD: MD workspace (with 3 dims) + :param pk: Peak object that has been integrated + :param peak_shape: shape of integrated peak (spherical, ellipsoid or none) + :param nbins_max: number of bins along the major axis of the ellipsoid + :param extent: extent of output MD workspace along each dimension in units of the outer background radius + :param frame_to_peak_centre_attr: name of attribute to return peak centre in appropriate frame + :return ws_cut: MD workspace in peak region + :return radii: radii of the ellipsoid peak region + :return bg_inner_radii: inner background radii of the ellipsoid peak region + :return bg_outer_radii: outer background radii of the ellipsoid peak region + :return box_lengths: length of each dimension in the MD workspace + :return: imax: index of dimension with largest radius + """ + ndims = wsMD.getNumDims() + radii, bg_inner_radii, bg_outer_radii, evecs, translation = BaseSX._get_ellipsoid_params_from_peak(peak_shape, ndims) + imax = np.argmax(radii) + # calc center in frame of ellipsoid axes cen = getattr(pk, frame_to_peak_centre_attr)() cen = np.matmul(evecs.T, np.array(cen) + translation) # get extents @@ -725,11 +728,12 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c box_lengths = box_lengths * extent extents = np.vstack((cen - box_lengths, cen + box_lengths)) # get nbins along each axis - nbins = np.array([int(nbins_max * radii[iax] / radii[imax]) for iax in range(len(radii))]) + nbins = np.array([int(nbins_max * radii[iax] / radii[imax]) for iax in range(ndims)]) # ensure minimum of 3 bins inside radius along each dim - min_nbins_in_radius = np.min(nbins * radii / box_lengths) - if min_nbins_in_radius < 3: - nbins = nbins * 3 / min_nbins_in_radius + min_nbins_in_radius = 3 + nbins_in_radius = np.min(nbins * radii / box_lengths) # along most coarsely binned dimension + if nbins_in_radius < min_nbins_in_radius: + nbins = nbins * min_nbins_in_radius / nbins_in_radius # call BinMD ws_cut = mantid.BinMD( InputWorkspace=wsMD, @@ -744,6 +748,24 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c ) return ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths, imax + @staticmethod + def _get_ellipsoid_params_from_peak(peak_shape, ndims): + shape_info = json_loads(peak_shape.toJSON()) + if peak_shape.shapeName().lower() == "spherical": + BaseSX.convert_spherical_representation_to_ellipsoid(shape_info) + # get radii + radii = np.array([shape_info[f"radius{iax}"] for iax in range(ndims)]) + bg_inner_radii = np.array([shape_info[f"background_inner_radius{iax}"] for iax in range(ndims)]) + bg_outer_radii = np.array([shape_info[f"background_outer_radius{iax}"] for iax in range(ndims)]) + # eignevectors of ellipsoid + evecs = np.zeros((ndims, ndims)) + for iax in range(ndims): + evec = np.array([float(elem) for elem in shape_info[f"direction{iax}"].split()]) + evecs[:, iax] = evec / np.linalg.norm(evec) + # get translation in frame of wsMD + translation = np.array([shape_info[f"translation{iax}"] for iax in range(ndims)]) + return radii, bg_inner_radii, bg_outer_radii, evecs, translation + @staticmethod def _convert_spherical_representation_to_ellipsoid(shape_info): # copied from mantidqt.widgets.sliceviewer.peaksviewer.representation.ellipsoid - can't import here though From 0be5e08a408cffa45b5544f4eef10585d6f7680a Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 6 Mar 2024 16:07:39 +0000 Subject: [PATCH 09/12] Do not store cut MD workspace in ADS --- scripts/Diffraction/single_crystal/base_sx.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index 9c040bd56231..9c2f17b57159 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -642,7 +642,6 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo # loop over peaks and plot try: with PdfPages(filename) as pdf: - ws_cut = None for ipk, pk in enumerate(peaks): peak_shape = pk.getPeakShape() if peak_shape.shapeName().lower() == "none": @@ -691,9 +690,6 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo fig.tight_layout() pdf.savefig(fig) close(fig) - - if ws_cut is not None: - mantid.DeleteWorkspace(ws_cut) except OSError: raise RuntimeError( f"OutputFile ({filename}) could not be opened - please check it is not open by " @@ -745,6 +741,7 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c OutputExtents=extents.flatten(order="F"), OutputBins=nbins.astype(int), EnableLogging=False, + StoreInADS=False, ) return ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths, imax From e31e7707ac3c8602895042e0201326c0d16a1ad8 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 6 Mar 2024 16:19:24 +0000 Subject: [PATCH 10/12] Avoid hard-coding float tolerance and use 1 line if expression --- scripts/Diffraction/single_crystal/base_sx.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index 9c2f17b57159..e3bdf9340622 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -665,7 +665,7 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo patch = Circle((0, 0), radii[imax] / box_lengths[imax], facecolor="none", edgecolor="r", ls="--") ax.add_patch(patch) # add background shell - if np.all(bg_outer_radii > 1e-10): + if not np.allclose(bg_outer_radii, 0.0): patch = Circle( (0, 0), bg_outer_radii[imax] / box_lengths[imax], facecolor="none", edgecolor=3 * [0.7] ) # outer radius @@ -720,7 +720,7 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c cen = getattr(pk, frame_to_peak_centre_attr)() cen = np.matmul(evecs.T, np.array(cen) + translation) # get extents - box_lengths = bg_outer_radii if np.all(bg_outer_radii > 1e-10) else radii + box_lengths = bg_outer_radii if not np.allclose(bg_outer_radii, 0.0) else radii box_lengths = box_lengths * extent extents = np.vstack((cen - box_lengths, cen + box_lengths)) # get nbins along each axis @@ -768,10 +768,7 @@ def _convert_spherical_representation_to_ellipsoid(shape_info): # copied from mantidqt.widgets.sliceviewer.peaksviewer.representation.ellipsoid - can't import here though # convert shape_info dict from sphere to ellipsoid for plotting for key in ["radius", "background_inner_radius", "background_outer_radius"]: - if key in shape_info: - shape_info[f"{key}{0}"] = shape_info.pop(key) - else: - shape_info[f"{key}{0}"] = 0.0 # null value + shape_info[f"{key}{0}"] = shape_info.pop(key) if key in shape_info else 0.0 for idim in [1, 2]: shape_info[f"{key}{idim}"] = shape_info[f"{key}{0}"] # add axes along basis vecs of frame and set 0 translation From 65a45fa4a4915a868230ddded12ae9b4066d6ec8 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 6 Mar 2024 17:01:51 +0000 Subject: [PATCH 11/12] Add type hinting --- scripts/Diffraction/single_crystal/base_sx.py | 23 +++++++++++++------ 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index e3bdf9340622..0b7481d23f12 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -4,14 +4,14 @@ # NScD Oak Ridge National Laboratory, European Spallation Source # & Institut Laue - Langevin # SPDX - License - Identifier: GPL - 3.0 + -from typing import Sequence +from typing import Sequence, Optional import numpy as np from enum import Enum import mantid.simpleapi as mantid -from mantid.api import FunctionFactory, AnalysisDataService as ADS +from mantid.api import FunctionFactory, AnalysisDataService as ADS, IMDEventWorkspace, IPeaksWorkspace, IPeak from mantid.kernel import logger, SpecialCoordinateSystem from FindGoniometerFromUB import getSignMaxAbsValInCol -from mantid.geometry import CrystalStructure, SpaceGroupFactory, ReflectionGenerator, ReflectionConditionFilter +from mantid.geometry import CrystalStructure, SpaceGroupFactory, ReflectionGenerator, ReflectionConditionFilter, PeakShape from os import path from json import loads as json_loads from matplotlib.pyplot import subplots, close @@ -619,7 +619,14 @@ def remove_non_integrated_peaks(peaks): ) @staticmethod - def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, log_norm=True): + def plot_integrated_peaks_MD( + wsMD: IMDEventWorkspace, + peaks: IPeaksWorkspace, + filename: str, + nbins_max: Optional[int] = 21, + extent: Optional[float] = 1.5, + log_norm: Optional[bool] = True, + ): """ :param wsMD: MD workspace to plot :param peaks: integrated peaks using IntegratePeaksMD @@ -697,7 +704,9 @@ def plot_integrated_peaks_MD(wsMD, peaks, filename, nbins_max=21, extent=1.5, lo ) @staticmethod - def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_centre_attr): + def _bin_MD_around_peak( + wsMD: IMDEventWorkspace, pk: IPeak, peak_shape: PeakShape, nbins_max: int, extent: float, frame_to_peak_centre_attr: str + ): """ Bin MD workspace in peak region with projection axes given by the axes of the ellipsoid shape :param wsMD: MD workspace (with 3 dims) @@ -746,7 +755,7 @@ def _bin_MD_around_peak(wsMD, pk, peak_shape, nbins_max, extent, frame_to_peak_c return ws_cut, radii, bg_inner_radii, bg_outer_radii, box_lengths, imax @staticmethod - def _get_ellipsoid_params_from_peak(peak_shape, ndims): + def _get_ellipsoid_params_from_peak(peak_shape: PeakShape, ndims: int): shape_info = json_loads(peak_shape.toJSON()) if peak_shape.shapeName().lower() == "spherical": BaseSX.convert_spherical_representation_to_ellipsoid(shape_info) @@ -764,7 +773,7 @@ def _get_ellipsoid_params_from_peak(peak_shape, ndims): return radii, bg_inner_radii, bg_outer_radii, evecs, translation @staticmethod - def _convert_spherical_representation_to_ellipsoid(shape_info): + def _convert_spherical_representation_to_ellipsoid(shape_info: dict): # copied from mantidqt.widgets.sliceviewer.peaksviewer.representation.ellipsoid - can't import here though # convert shape_info dict from sphere to ellipsoid for plotting for key in ["radius", "background_inner_radius", "background_outer_radius"]: From 2b8f9a13af9758531016b5b390f516780e4f76d0 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 7 Mar 2024 12:50:13 +0000 Subject: [PATCH 12/12] Add more doc strings and remove OutputWorkspace arg --- scripts/Diffraction/single_crystal/base_sx.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index 0b7481d23f12..809b2dfe2829 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -628,6 +628,7 @@ def plot_integrated_peaks_MD( log_norm: Optional[bool] = True, ): """ + Function to plot summary of IntegratePeaksMD integration comprising 3 colorfill plots per peak saved in a pdf. :param wsMD: MD workspace to plot :param peaks: integrated peaks using IntegratePeaksMD :param filename: filename to store pdf output @@ -742,7 +743,6 @@ def _bin_MD_around_peak( # call BinMD ws_cut = mantid.BinMD( InputWorkspace=wsMD, - OutputWorkspace="__ws_cut", AxisAligned=False, BasisVector0=r"Q$_0$,unit," + ",".join(np.array2string(evecs[:, 0], precision=6).strip("[]").split()), BasisVector1=r"Q$_1$,unit," + ",".join(np.array2string(evecs[:, 1], precision=6).strip("[]").split()), @@ -756,6 +756,16 @@ def _bin_MD_around_peak( @staticmethod def _get_ellipsoid_params_from_peak(peak_shape: PeakShape, ndims: int): + """ + Extract ellipsoid parameters (eigenvectors, radii etc.) from PeakShape object + :param peak_shape: PeakShape object of a integrated peak + :param ndims: number of dimensions in the MD workspace + :return radii: array of ellipsoid radii (3 sigma) defining peak region + :return bg_inner_radii: array of inner radii for ellipsoid background shell + :return bg_outer_radii: array of outer radii for ellipsoid background shell + :return evecs: ndims x ndims array of eignevectors - each column corresponds to an ellipsoid axis + :return translation: translation of the ellipsoid center in the coordinates/frame of the MD workspace integrated + """ shape_info = json_loads(peak_shape.toJSON()) if peak_shape.shapeName().lower() == "spherical": BaseSX.convert_spherical_representation_to_ellipsoid(shape_info) @@ -774,6 +784,10 @@ def _get_ellipsoid_params_from_peak(peak_shape: PeakShape, ndims: int): @staticmethod def _convert_spherical_representation_to_ellipsoid(shape_info: dict): + """ + Update shape_info dict of a spherical peak shape to have same keys/fields as an ellipsoid + :param shape_info: dictionary of JSON shape representation + """ # copied from mantidqt.widgets.sliceviewer.peaksviewer.representation.ellipsoid - can't import here though # convert shape_info dict from sphere to ellipsoid for plotting for key in ["radius", "background_inner_radius", "background_outer_radius"]: