diff --git a/sasdata/data.py b/sasdata/data.py index 788f3d7f5..73e64e612 100644 --- a/sasdata/data.py +++ b/sasdata/data.py @@ -39,7 +39,9 @@ def __init__( @property def ordinate(self) -> Quantity: match self.dataset_type: - case dataset_types.one_dim | dataset_types.two_dim: + case (dataset_types.one_dim | + dataset_types.two_dim | + dataset_types.angle_dim): return self._data_contents["I"] case dataset_types.sesans: return self._data_contents["Depolarisation"] @@ -70,6 +72,8 @@ def abscissae(self) -> Quantity: # probably want to avoid creating a new Quantity but at the moment I # can't see a way around it. return Quantity(data_contents, reference_data_content.units, name=self._data_contents["Qx"].name, id_header=self._data_contents["Qx"]._id_header) + case dataset_types.angle_dim: + return self._data_contents["Phi"] case dataset_types.sesans: return self._data_contents["SpinEchoLength"] case _: @@ -148,3 +152,38 @@ def default(self, obj): } case _: return super().default(obj) + + +def sasdata_reader2D_converter(data2d: SasData | None = None) -> SasData: + """ + convert old 2d format opened by IhorReader or danse_reader + to new 2D SasData format + This is mainly used by the Readers + + :param data2d: SasData object with 2D arrays + :return: SasData object with 1D arrays + + """ + if data2d._data_contents["I"] is None or data2d.x_bins is None or data2d.y_bins is None: + raise ValueError("Can't convert this data: data=None...") + new_x = np.tile(data2d.x_bins, (len(data2d.y_bins), 1)) + new_y = np.tile(data2d.y_bins, (len(data2d.x_bins), 1)) + new_y = new_y.swapaxes(0, 1) + + new_data = data2d._data_contents["I"].value.flatten() + qx_data = new_x.flatten() + qy_data = new_y.flatten() + err_data = np.sqrt(data2d._data_contents["I"].variance.value) + if not data2d._data_contents["I"].has_variance or np.any(err_data <= 0): + new_err_data = np.sqrt(np.abs(new_data)) + else: + new_err_data = err_data.flatten() + mask = np.ones(len(new_data), dtype=bool) + + data2d._data_contents["I"].value = new_data + data2d._data_contents["I"].variance.value = new_err_data ** 2 + data2d._data_contents["Qx"].value = qx_data + data2d._data_contents["Qy"].value = qy_data + data2d.mask = mask + + return data2d diff --git a/sasdata/data_util/averaging.py b/sasdata/data_util/averaging.py index 2e9f0f1b0..386543bc1 100644 --- a/sasdata/data_util/averaging.py +++ b/sasdata/data_util/averaging.py @@ -1,20 +1,21 @@ -""" -This module contains various data processors used by Sasview's slicers. -""" +"""This module contains various data processors used by Sasview's slicers.""" import math import numpy as np import numpy.typing as npt +from sasdata.data import SasData from sasdata.data_util.binning import DirectionalAverage from sasdata.data_util.interval import IntervalType from sasdata.data_util.roi import CartesianROI, PolarROI -from sasdata.dataloader.data_info import Data1D, Data2D +from sasdata.dataset_types import angle_dim, one_dim from sasdata.quantities.constants import Pi, TwoPi +from sasdata.quantities.quantity import Quantity +from sasdata.quantities.units import radians -def get_dq_data(data2d: Data2D) -> npt.NDArray[np.floating]: +def get_dq_data(data2d: SasData) -> npt.NDArray[np.floating]: """ Get the dq for resolution averaging The pinholes and det. pix contribution present @@ -23,23 +24,31 @@ def get_dq_data(data2d: Data2D) -> npt.NDArray[np.floating]: q = 0. Note This method works on only pinhole geometry. Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. """ - z_max = max(data2d.q_data) - z_min = min(data2d.q_data) - dqx_at_z_max = data2d.dqx_data[np.argmax(data2d.q_data)] - dqx_at_z_min = data2d.dqx_data[np.argmin(data2d.q_data)] - dqy_at_z_max = data2d.dqy_data[np.argmax(data2d.q_data)] - dqy_at_z_min = data2d.dqy_data[np.argmin(data2d.q_data)] + + q_data = np.sqrt(data2d._data_contents["Qx"].value**2 + data2d._data_contents["Qy"].value**2) + + z_max = max(q_data) + z_min = min(q_data) + + dqx_data = np.sqrt(data2d._data_contents["Qx"].variance.value) + dqy_data = np.sqrt(data2d._data_contents["Qy"].variance.value) + + dqx_at_z_max = dqx_data[np.argmax(q_data)] + dqx_at_z_min = dqx_data[np.argmin(q_data)] + dqy_at_z_max = dqy_data[np.argmax(q_data)] + dqy_at_z_min = dqy_data[np.argmin(q_data)] + # Find qdx at q = 0 dq_overlap_x = (dqx_at_z_min * z_max - dqx_at_z_max * z_min) / (z_max - z_min) # when extrapolation goes wrong - if dq_overlap_x > min(data2d.dqx_data): - dq_overlap_x = min(data2d.dqx_data) + if dq_overlap_x > min(dqx_data): + dq_overlap_x = min(dqx_data) dq_overlap_x *= dq_overlap_x # Find qdx at q = 0 dq_overlap_y = (dqy_at_z_min * z_max - dqy_at_z_max * z_min) / (z_max - z_min) # when extrapolation goes wrong - if dq_overlap_y > min(data2d.dqy_data): - dq_overlap_y = min(data2d.dqy_data) + if dq_overlap_y > min(dqy_data): + dq_overlap_y = min(dqy_data) # get dq at q=0. dq_overlap_y *= dq_overlap_y @@ -47,8 +56,8 @@ def get_dq_data(data2d: Data2D) -> npt.NDArray[np.floating]: # Final protection of dq if dq_overlap < 0: dq_overlap = dqy_at_z_min - dqx_data = data2d.dqx_data[np.isfinite(data2d.data)] - dqy_data = data2d.dqy_data[np.isfinite(data2d.data)] - dq_overlap + dqx_data = dqx_data[np.isfinite(data2d.ordinate.value)] + dqy_data = dqy_data[np.isfinite(data2d.ordinate.value)] - dq_overlap # def; dqx_data = dq_r dqy_data = dq_phi # Convert dq 2D to 1D here dq_data = np.sqrt(dqx_data**2 + dqy_data**2) @@ -56,9 +65,7 @@ def get_dq_data(data2d: Data2D) -> npt.NDArray[np.floating]: class Boxsum(CartesianROI): - """ - Compute the sum of the intensity within a rectangular Region Of Interest. - """ + """Compute the sum of the intensity within a rectangular Region Of Interest.""" def __init__(self, qx_range: tuple[float, float] = (0.0, 0.0), qy_range: tuple[float, float] = (0.0, 0.0)) -> None: """ @@ -70,11 +77,11 @@ def __init__(self, qx_range: tuple[float, float] = (0.0, 0.0), qy_range: tuple[f """ super().__init__(qx_range=qx_range, qy_range=qy_range) - def __call__(self, data2d: Data2D) -> tuple[float, float, float]: + def __call__(self, data2d: SasData = None) -> tuple[float, float, float]: """ Coordinate data processing operations and return the results. - :param data2d: The Data2D object for which the sum is calculated. + :param data2d: The SasData object for which the sum is calculated. """ self.validate_and_assign_data(data2d) total_sum, error, count = self._sum() @@ -97,7 +104,7 @@ def _sum(self) -> tuple[float, float, float]: data = weights * self.data # Not certain that the weights should be squared here, I'm just copying # how it was done in the old manipulations.py - err_squared = weights * weights * self.err_data * self.err_data + err_squared = (weights * self.err_data) ** 2 total_sum = np.sum(data) total_errors_squared = np.sum(err_squared) @@ -107,25 +114,22 @@ def _sum(self) -> tuple[float, float, float]: class Boxavg(Boxsum): - """ - Compute the average intensity within a rectangular Region Of Interest. - """ + """Compute the average intensity within a rectangular Region Of Interest.""" def __init__(self, qx_range: tuple[float, float] = (0.0, 0.0), qy_range: tuple[float, float] = (0.0, 0.0)) -> None: """ Set up the Region of Interest and its boundaries. - The units of these parameters are A^-1 :param qx_range: Bounds of the ROI along the Q_x direction. :param qy_range: Bounds of the ROI along the Q_y direction. """ super().__init__(qx_range=qx_range, qy_range=qy_range) - def __call__(self, data2d: Data2D) -> tuple[float, float]: + def __call__(self, data2d: SasData) -> tuple[float, float]: """ Coordinate data processing operations and return the results. - :param data2d: The Data2D object for which the average is calculated. + :param data2d: The SasData object for which the average is calculated. """ self.validate_and_assign_data(data2d) total_sum, error, count = super()._sum() @@ -138,7 +142,7 @@ class SlabX(CartesianROI): Average I(Q_x, Q_y) along the y direction (within a ROI), giving I(Q_x). This class is initialised by specifying the boundaries of the ROI and is - called by supplying a Data2D object. It returns a Data1D object. + called by supplying a SasData object. It returns a SasData object. The averaging process can also be thought of as projecting 2D -> 1D. There also exists the option to "fold" the ROI, where Q data on opposite @@ -157,7 +161,6 @@ def __init__( """ Set up the ROI boundaries, the binning of the output 1D data, and fold. - The units of these parameters are A^-1 :param qx_range: Bounds of the ROI along the Q_x direction. :param qy_range: Bounds of the ROI along the Q_y direction. :param nbins: The number of bins data is sorted into along Q_x. @@ -169,12 +172,12 @@ def __init__( self.fold: bool = fold self.base: float | None = base - def __call__(self, data2d: Data2D) -> Data1D: + def __call__(self, data2d: SasData = None) -> SasData: """ Compute the 1D average of 2D data, projecting along the Q_x axis. - :param data2d: The Data2D object for which the average is computed. - :return: Data1D object for plotting. + :param data2d: The SasData object for which the average is computed. + :return: SasData object for plotting. """ self.validate_and_assign_data(data2d) @@ -197,9 +200,14 @@ def __call__(self, data2d: Data2D) -> Data1D: nbins=self.nbins, base=self.base, ) + qx_data, intensity, error = directional_average(data=self.data, err_data=self.err_data) - return Data1D(x=qx_data, y=intensity, dy=error) + data_contents = { + "Q": Quantity(qx_data, data2d._data_contents["Qx"].units, None), + "I": Quantity(intensity, data2d.ordinate.units, error), + } + return SasData(f"{data2d.name}: Slab X Average", data_contents, one_dim, data2d.metadata) class SlabY(CartesianROI): @@ -207,7 +215,7 @@ class SlabY(CartesianROI): Average I(Q_x, Q_y) along the x direction (within a ROI), giving I(Q_y). This class is initialised by specifying the boundaries of the ROI and is - called by supplying a Data2D object. It returns a Data1D object. + called by supplying a SasData object. It returns a SasData object. The averaging process can also be thought of as projecting 2D -> 1D. There also exists the option to "fold" the ROI, where Q data on opposite @@ -226,7 +234,6 @@ def __init__( """ Set up the ROI boundaries, the binning of the output 1D data, and fold. - The units of these parameters are A^-1 :param qx_range: Bounds of the ROI along the Q_x direction. :param qy_range: Bounds of the ROI along the Q_y direction. :param qy_max: Upper bound of the ROI along the Q_y direction. @@ -239,12 +246,12 @@ def __init__( self.fold: bool = fold self.base: float | None = base - def __call__(self, data2d: Data2D) -> Data1D: + def __call__(self, data2d: SasData = None) -> SasData: """ Compute the 1D average of 2D data, projecting along the Q_y axis. - :param data2d: The Data2D object for which the average is computed. - :return: Data1D object for plotting. + :param data2d: The SasData object for which the average is computed. + :return: SasData object for plotting. """ self.validate_and_assign_data(data2d) @@ -269,7 +276,11 @@ def __call__(self, data2d: Data2D) -> Data1D: ) qy_data, intensity, error = directional_average(data=self.data, err_data=self.err_data) - return Data1D(x=qy_data, y=intensity, dy=error) + data_contents = { + "Q": Quantity(qy_data, data2d._data_contents["Qy"].units, None), + "I": Quantity(intensity, data2d.ordinate.units, error), + } + return SasData(f"{data2d.name}: Slab Y Average", data_contents, one_dim, data2d.metadata) class CircularAverage(PolarROI): @@ -279,7 +290,7 @@ class CircularAverage(PolarROI): This class is initialised by specifying lower and upper limits on the magnitude of Q values to consider during the averaging, though currently SasView always calls this class using the full range of data. When called, - this class is supplied with a Data2D object. It returns a Data1D object + this class is supplied with a SasData object. It returns a SasData object where intensity is given as a function of Q only. """ @@ -293,7 +304,6 @@ def __init__( """ Set up the lower and upper radial limits as well as the number of bins. - The units are A^-1 for the radial parameters. :param r_min: Lower limit for |Q| values to use during averaging. :param r_max: Upper limit for |Q| values to use during averaging. :param nbins: The number of bins data is sorted into along |Q| the axis @@ -302,27 +312,27 @@ def __init__( self.nbins: int = nbins self.base: float | None = base - def __call__(self, data2D: Data2D, ismask: bool = False) -> Data1D: + def __call__(self, data2D: SasData, ismask: bool = False) -> SasData: """ Perform circular averaging on the data. Uses DirectionalAverage for bin construction and weights, and computes dx (d_q) using get_dq_data averaged with those weights so behavior matches the legacy implementation. - :param data2D: Data2D object + :param data2D: SasData object :param ismask: If True, respect data2D.mask (skip masked points). If False, ignore mask. - :return: Data1D object with x (bin centers), y (intensity), dy and dx (if available) + :return: SasData object with x (bin centers), y (intensity), dy and dx (if available) """ # Work on unmasked finite arrays first (matches legacy filtering) - finite_mask = np.isfinite(data2D.data) + finite_mask = np.isfinite(data2D.ordinate.value) if not np.any(finite_mask): raise RuntimeError(f"Circular averaging: invalid q_data: {data2D.q_data}") - data_all = data2D.data[finite_mask] - q_all = data2D.q_data[finite_mask] - qx_all = data2D.qx_data[finite_mask] - qy_all = data2D.qy_data[finite_mask] - err_all = data2D.err_data[finite_mask] if data2D.err_data is not None else None - mask_all = data2D.mask[finite_mask] + data_all = data2D.ordinate.value[finite_mask] + qx_all = data2D._data_contents["Qx"].value[finite_mask] + qy_all = data2D._data_contents["Qy"].value[finite_mask] + q_all = np.sqrt(data2D._data_contents["Qx"].value**2 + data2D._data_contents["Qy"].value**2)[finite_mask] + err_all = np.sqrt(data2D.ordinate.variance.value)[finite_mask] + mask_all = (data2D.mask if data2D.mask is not None else np.ones_like(data2D.ordinate.value, dtype=bool))[finite_mask] # Optional mask handling: legacy used an ismask flag to optionally skip masked points if ismask: @@ -334,24 +344,22 @@ def __call__(self, data2D: Data2D, ismask: bool = False) -> Data1D: major_axis = q_all[sel] phi_axis = np.arctan2(qy_all[sel], qx_all[sel]) data_vals = data_all[sel] - err_vals = err_all[sel] if err_all is not None else None + err_vals = err_all[sel] # Prepare dq_data if available, aligned to the finite mask and selection dq_vals = None - if getattr(data2D, "dqx_data", None) is not None and getattr(data2D, "dqy_data", None) is not None: + if data2D._data_contents["Qx"].has_variance and data2D._data_contents["Qy"].has_variance: dq_full = get_dq_data(data2D) # already uses np.isfinite(data2D.data) dq_vals = dq_full[sel] # Set up DirectionalAverage with full-circle phi range major_lims = (self.r_min, self.r_max) minor_lims = (0.0, TwoPi) - directional_average = DirectionalAverage( - major_axis=major_axis, - minor_axis=phi_axis, - lims=(major_lims, minor_lims), - nbins=self.nbins, - base=self.base, - ) + directional_average = DirectionalAverage(major_axis=major_axis, + minor_axis=phi_axis, + lims=(major_lims, minor_lims), + nbins=self.nbins, + base=self.base) # Compute weights, then produce averaged intensity/error via DirectionalAverage weights = directional_average.compute_weights() @@ -371,11 +379,16 @@ def __call__(self, data2D: Data2D, ismask: bool = False) -> Data1D: counts = np.sum(weights, axis=1) with np.errstate(divide="ignore", invalid="ignore"): dx_full = dq_weighted / counts - dx = dx_full[populated] + dQ = dx_full[populated] else: - dx = None + dQ = None + + data_contents = { + "Q": Quantity(x, data2D._data_contents["Qx"].units, dQ), + "I": Quantity(intensity, data2D.ordinate.units, error), + } + return SasData(f"{data2D.name}: Circular Average", data_contents, one_dim, data2D.metadata) - return Data1D(x=x, y=intensity, dy=error, dx=dx) class Ring(PolarROI): @@ -384,9 +397,8 @@ class Ring(PolarROI): This class is initialised by specifying lower and upper limits on the magnitude of Q values to consider during the averaging. When called, - this class is supplied with a Data2D object. It returns a Data1D object. - This Data1D object gives intensity as a function of the angle from the - positive x-axis, φ, only. + this class is supplied with a SasData object. It returns a SasData object + which gives intensity as a function of the angle from the positive x-axis, φ, only. """ def __init__( @@ -399,7 +411,6 @@ def __init__( """ Set up the lower and upper radial limits as well as the number of bins. - The units are A^-1 for the radial parameters. :param r_min: Lower limit for |Q| values to use during averaging. :param r_max: Upper limit for |Q| values to use during averaging. :param nbins: The number of bins data is sorted into along Phi the axis @@ -411,27 +422,32 @@ def __init__( self.nbins: int = nbins self.base: float | None = base - def __call__(self, data2D: Data2D) -> Data1D: + def __call__(self, data2D: SasData) -> SasData: """ Apply the ring to the data set. Returns the angular distribution for a given q range - :param data2D: Data2D object + :param data2D: SasData object - :return: Data1D object + :return: SasData object """ - if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: - raise RuntimeError("Ring averaging only take plottable_2D objects") + if not isinstance(data2D, SasData): + msg = "Data supplied for ring averaging must be of type SasData." + raise RuntimeError(msg) + if not ("Qx" in data2D._data_contents and + "Qy" in data2D._data_contents): + msg = "SasData object for ring averaging must contain 'Qx' and 'Qy' data." + raise RuntimeError(msg) # Get data - data = data2D.data[np.isfinite(data2D.data)] - q_data = data2D.q_data[np.isfinite(data2D.data)] - err_data = None - if data2D.err_data is not None: - err_data = data2D.err_data[np.isfinite(data2D.data)] - qx_data = data2D.qx_data[np.isfinite(data2D.data)] - qy_data = data2D.qy_data[np.isfinite(data2D.data)] - mask_data = data2D.mask[np.isfinite(data2D.data)] + valid_data = np.isfinite(data2D.ordinate.value) + + data = data2D.ordinate.value[valid_data] + err_data = np.sqrt(data2D.ordinate.variance.value)[valid_data] + qx_data = data2D._data_contents["Qx"].value[valid_data] + qy_data = data2D._data_contents["Qy"].value[valid_data] + q_data = np.sqrt(data2D._data_contents["Qx"].value ** 2 + data2D._data_contents["Qy"].value ** 2)[valid_data] + mask_data = (data2D.mask if data2D.mask is not None else np.ones_like(data2D.ordinate.value, dtype=bool))[valid_data] # Set space for 1d outputs phi_bins = np.zeros(self.nbins) @@ -486,66 +502,11 @@ def __call__(self, data2D: Data2D) -> Data1D: msg = "Average Error: No points inside ROI to average..." raise ValueError(msg) - return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx]) - - ''' - def __call__(self, data2d: Data2D = None) -> Data1D: - """ - Compute the 1D average of 2D data, projecting along the Phi axis. - - :param data2d: The Data2D object for which the average is computed. - :return: Data1D object for plotting. - """ - self.validate_and_assign_data(data2d) - - # half-bin shift so the first bin is centered at zero - phi_shift = np.pi / self.nbins - shifted_phi = (self.phi_data +Pi+ phi_shift) % (TwoPi) - - self.phi_data=self.phi_data+Pi - self.phi_min = 0.0 - self.phi_max = TwoPi - - minor_lims = (self.r_min, self.r_max) - major_lims = (self.phi_min, self.phi_max) - - - # major_lims is None because a full-circle angular range is used - directional_average = DirectionalAverage(major_axis=shifted_phi, - minor_axis=self.q_data, - lims=(major_lims,minor_lims), - nbins=self.nbins, base=self.base) - # Reuse DirectionalAverage's weights, then compute the same sums/divisions - weights = directional_average.compute_weights() - - # Projected x values (averaged shifted phi per bin) -- not used as final x, - # but computed here to mirror DirectionalAverage internal behaviour. - x_axis_values = np.sum(weights * shifted_phi, axis=1) - - intensity = np.sum(weights * self.data, axis=1) - errs_squared = np.sum((weights * self.err_data) ** 2, axis=1) - - bin_counts = np.sum(weights, axis=1) - if not np.any(bin_counts > 0): - raise ValueError("Average Error: No bins inside ROI to average...") - - errors = np.sqrt(errs_squared) - - # Only compute divisions where bin_counts > 0 (others will become NaN/inf) - with np.errstate(divide='ignore', invalid='ignore'): - x_axis_values = x_axis_values / bin_counts - intensity = intensity / bin_counts - errors = errors / bin_counts - - # Legacy reported x values are the unshifted bin starts (i.e. 2*pi*i/nbins) - phi_values = np.linspace(0.0, TwoPi, self.nbins, endpoint=False) - - finite = np.isfinite(intensity) - if not finite.any(): - raise ValueError("Average Error: No points inside ROI to average...") - - return Data1D(x=phi_values[finite], y=intensity[finite], dy=errors[finite]) - ''' + data_contents = { + "Phi": Quantity(phi_values[idx], radians, None), + "I": Quantity(phi_bins[idx], data2D.ordinate.units, phi_err[idx]), + } + return SasData(f"{data2D.name}: Ring Average", data_contents, angle_dim, data2D.metadata) class SectorQ(PolarROI): @@ -566,8 +527,8 @@ class SectorQ(PolarROI): the data from the two regions are graphed separeately, with the secondary ROI data labelled using negative Q values. - When called, this class is supplied with a Data2D object. It returns a - Data1D object where intensity is given as a function of Q only. + When called, this class is supplied with a SasData object. It returns a + SasData object where intensity is given as a function of Q only. """ def __init__( @@ -582,9 +543,8 @@ def __init__( """ Set up the ROI boundaries, the binning of the output 1D data, and fold. - The units are A^-1 for radial parameters, and radians for anglar ones. :param r_range: Tuple (r_min, r_max) defining limits for |Q| values to use during averaging. - :param phi_range: Tuple (phi_min, phi_max) defining limits for φ in radians (in the primary ROI). + :param phi_range: Tuple (phi_min, phi_max) defining limits for φ in the primary ROI. :Defaults to full circle (0, 2*pi). :param nbins: The number of bins data is sorted into along the |Q| axis :param fold: Whether the primary and secondary ROIs should be folded @@ -596,12 +556,12 @@ def __init__( self.fold: bool = fold self.base: float | None = base - def __call__(self, data2d: Data2D) -> Data1D: + def __call__(self, data2d: SasData = None) -> SasData: """ Compute the 1D average of 2D data, projecting along the Q_y axis. - :param data2d: The Data2D object for which the average is computed. - :return: Data1D object for plotting. + :param data2d: The SasData object for which the average is computed. + :return: SasData object for plotting. """ self.validate_and_assign_data(data2d) @@ -643,8 +603,10 @@ def __call__(self, data2d: Data2D) -> Data1D: base=self.base, ) - primary_q, primary_I, primary_err = primary_region(data=self.data, err_data=self.err_data) - secondary_q, secondary_I, secondary_err = secondary_region(data=self.data, err_data=self.err_data) + primary_q, primary_I, primary_err = primary_region(data=self.data, + err_data=self.err_data) + secondary_q, secondary_I, secondary_err = secondary_region(data=self.data, + err_data=self.err_data) if self.fold: # Combining the two regions requires re-binning; the q value @@ -676,15 +638,22 @@ def __call__(self, data2d: Data2D) -> Data1D: finite = np.isfinite(average_intensity) - data1d = Data1D(x=combined_q[finite], y=average_intensity[finite], dy=combined_err[finite]) + data_contents = { + "Q": Quantity(combined_q[finite], data2d._data_contents["Qx"].units, None), + "I": Quantity(average_intensity[finite], data2d.ordinate.units, combined_err[finite]), + } else: # The secondary ROI is labelled with negative Q values. combined_q = np.append(np.flip(-1 * secondary_q), primary_q) combined_intensity = np.append(np.flip(secondary_I), primary_I) combined_error = np.append(np.flip(secondary_err), primary_err) - data1d = Data1D(x=combined_q, y=combined_intensity, dy=combined_error) - return data1d + data_contents = { + "Q": Quantity(combined_q, data2d._data_contents["Qx"].units, None), + "I": Quantity(combined_intensity, data2d.ordinate.units, combined_error), + } + + return SasData(f"{data2d.name}:SectorQ Average", data_contents, one_dim, data2d.metadata) class WedgeQ(PolarROI): @@ -696,9 +665,9 @@ class WedgeQ(PolarROI): the positive x-axis. This class is initialised by specifying lower and upper limits on both the - magnitude of Q and the angle φ. - When called, this class is supplied with a Data2D object. It returns a - Data1D object where intensity is given as a function of Q only. + magnitude of Q and the angle φ. When called, this class is supplied with a + SasData object. It returns a sasData object where intensity is given as a + function of Q only. """ def __init__( @@ -712,9 +681,8 @@ def __init__( """ Set up the ROI boundaries, and the binning of the output 1D data. - The units are A^-1 for radial parameters, and radians for anglar ones. :param r_range: Tuple (r_min, r_max) defining limits for |Q| values to use during averaging. - :param phi_range: Tuple (phi_min, phi_max) defining limits for φ in radians (in the primary ROI). + :param phi_range: Tuple (phi_min, phi_max) defining limits for φ in the primary ROI. :Defaults to full circle (0, 2*pi). :param nbins: The number of bins data is sorted into along the |Q| axis """ @@ -722,12 +690,12 @@ def __init__( self.nbins: int = nbins self.base: float | None = base - def __call__(self, data2d: Data2D) -> Data1D: + def __call__(self, data2d: SasData = None) -> SasData: """ Compute the 1D average of 2D data, projecting along the Q_y axis. - :param data2d: The Data2D object for which the average is computed. - :return: Data1D object for plotting. + :param data2d: The SasData object for which the average is computed. + :return: SasData object for plotting. """ self.validate_and_assign_data(data2d) @@ -764,9 +732,14 @@ def __call__(self, data2d: Data2D) -> Data1D: nbins=self.nbins, base=self.base, ) + q_data, intensity, error = directional_average(data=self.data, err_data=self.err_data) - return Data1D(x=q_data, y=intensity, dy=error) + data_contents = { + "Q": Quantity(q_data, data2d._data_contents["Qx"].units, None), + "I": Quantity(intensity, data2d.ordinate.units, error), + } + return SasData(f"{data2d.name}: Wedge Q Average", data_contents, one_dim, data2d.metadata) class WedgePhi(PolarROI): @@ -778,9 +751,8 @@ class WedgePhi(PolarROI): This class is initialised by specifying lower and upper limits on both the magnitude of Q and the angle φ, measured anticlockwise from the positive - x-axis. - When called, this class is supplied with a Data2D object. It returns a - Data1D object where intensity is given as a function of Q only. + x-axis. When called, this class is supplied with a SasData object. It returns + a SasData object where intensity is given as a function of φ only. """ def __init__( @@ -794,7 +766,6 @@ def __init__( """ Set up the ROI boundaries, and the binning of the output 1D data. - The units are A^-1 for radial parameters, and radians for anglar ones. :param r_range: Tuple (r_min, r_max) defining limits for |Q| values to use during averaging. :param phi_range: Tuple (phi_min, phi_max) defining angular bounds in radians. Defaults to full circle (0, 2*pi). @@ -805,12 +776,12 @@ def __init__( self.nbins: int = nbins self.base: float | None = base - def __call__(self, data2d: Data2D) -> Data1D: + def __call__(self, data2d: SasData = None) -> SasData: """ Compute the 1D average of 2D data, projecting along the Q_y axis. - :param data2d: The Data2D object for which the average is computed. - :return: Data1D object for plotting. + :param data2d: The SasData object for which the average is computed. + :return: SasData object for plotting. """ self.validate_and_assign_data(data2d) @@ -829,8 +800,8 @@ def __call__(self, data2d: Data2D) -> Data1D: # Remember to transform back afterward as we're plotting against phi. phi_offset = self.phi_min self.phi_min = 0.0 - self.phi_max = (self.phi_max - phi_offset) % (TwoPi) - self.phi_data = (self.phi_data - phi_offset) % (TwoPi) + self.phi_max = (self.phi_max - phi_offset) % TwoPi + self.phi_data = (self.phi_data - phi_offset) % TwoPi # Averaging takes place between angular and radial limits # When phi_max and phi_min have the same angle, ROI is a full circle. @@ -847,7 +818,7 @@ def __call__(self, data2d: Data2D) -> Data1D: nbins=self.nbins, base=self.base, ) - phi_data, intensity, error = directional_average(data=self.data, err_data=self.err_data) + _, intensity, error = directional_average(data=self.data, err_data=self.err_data) # Compute phi bin starts to match legacy behaviour (Ring / old SectorPhi) # phi_min has been normalized to 0 earlier; phi_offset stores original start. @@ -858,7 +829,7 @@ def __call__(self, data2d: Data2D) -> Data1D: full_phi = np.linspace(self.phi_min, self.phi_max, self.nbins, endpoint=False) # Shift back to original phi range - full_phi = (full_phi + phi_offset) % (TwoPi) + full_phi = (full_phi + phi_offset) % TwoPi # Determine which bins were populated using the weights (preserves full bin index space) weights = directional_average.compute_weights() @@ -872,21 +843,11 @@ def __call__(self, data2d: Data2D) -> Data1D: phi_centers = full_phi[populated] + directional_average.bin_widths[populated] / 2.0 # intensity and error returned by DirectionalAverage are already filtered to the populated/finite bins - return Data1D(x=phi_centers, y=intensity, dy=error) - - """ - # Convert angular data back to the original phi range - phi_data += phi_offset - # In the old manipulations.py, we also had this shift to plot the data - # at the centre of the bins. I'm not sure why it's only angular binning - # which gets this treatment. - # TODO: Update this once non-linear binning options are implemented - weights = directional_average.compute_weights() - populated = np.sum(weights, axis=1) > 0 - phi_data += directional_average.bin_widths[populated] / 2 - - return Data1D(x=phi_data, y=intensity, dy=error) - """ + data_contents = { + "Phi": Quantity(phi_centers, radians, None), + "I": Quantity(intensity, data2d.ordinate.units, error), + } + return SasData(f"{data2d.name}: Wedge Phi Average", data_contents, angle_dim, data2d.metadata) class SectorPhi(WedgePhi): @@ -915,23 +876,20 @@ def __init__( nbins: int = 100, ) -> None: # Forward to WedgePhi using the tuple-based it expects. - super().__init__(r_range=(r_min, r_max), phi_range=(phi_min, phi_max), center=center, nbins=nbins) - ################################################################################ class Ringcut(PolarROI): """ - Defines a ring on a 2D data set. - The ring is defined by r_min, r_max, and + Defines a ring on a 2D data set. The ring is defined by r_min, r_max, and the position of the center of the ring. - The data returned is the region inside the ring + The data returned is the region inside the ring. Phi_min and phi_max should be defined between 0 and 2*pi - in anti-clockwise starting from the x- axis on the left-hand side + in anti-clockwise starting from the x-axis on the left-hand side """ def __init__( @@ -942,19 +900,19 @@ def __init__( ) -> None: super().__init__(r_range, phi_range, center) - def __call__(self, data2D: Data2D) -> npt.NDArray[np.bool_]: + def __call__(self, data2D: SasData) -> npt.NDArray[np.bool_]: """ Apply the ring to the data set. - Returns the angular distribution for a given q range + Returns the angular distribution for a given q range. - :param data2D: Data2D object + :param data2D: SasData object :return: index array in the range """ - super().validate_and_assign_data(data2D) + self.validate_and_assign_data(data2D) # Calculate q_data using unmasked qx_data and qy_data - q_data = np.sqrt(data2D.qx_data * data2D.qx_data + data2D.qy_data * data2D.qy_data) + q_data = np.sqrt(self.qx_data**2 + self.qy_data**2) # check whether each data point is inside ROI out = (self.r_min <= q_data) & (self.r_max >= q_data) @@ -962,25 +920,23 @@ def __call__(self, data2D: Data2D) -> npt.NDArray[np.bool_]: class Boxcut(CartesianROI): - """ - Find a rectangular 2D region of interest. - """ + """Find a rectangular 2D region of interest.""" def __init__(self, qx_range: tuple[float, float] = (0.0, 0.0), qy_range: tuple[float, float] = (0.0, 0.0)) -> None: super().__init__(qx_range=qx_range, qy_range=qy_range) - def __call__(self, data2D: Data2D) -> npt.NDArray[np.bool_]: + def __call__(self, data2D: SasData) -> npt.NDArray[np.bool_]: """ - Find a rectangular 2D region of interest where data points inside the ROI are True, and False otherwise + Find a rectangular 2D region of interest where data points inside the ROI are True, and False otherwise. - :param data2D: Data2D object + :param data2D: SasData object :return: mask, 1d array (len = len(data)) """ - super().validate_and_assign_data(data2D) + self.validate_and_assign_data(data2D) # check whether each data point is inside ROI - outx = (self.qx_min <= data2D.qx_data) & (self.qx_max > data2D.qx_data) - outy = (self.qy_min <= data2D.qy_data) & (self.qy_max > data2D.qy_data) + outx = (self.qx_min <= self.qx_data) & (self.qx_max > self.qx_data) + outy = (self.qy_min <= self.qy_data) & (self.qy_max > self.qy_data) return outx & outy @@ -999,24 +955,22 @@ class Sectorcut(PolarROI): def __init__(self, phi_range: tuple[float, float] = (0.0, Pi), center: tuple[float, float] = (0.0, 0.0)) -> None: super().__init__(r_range=(0, np.inf), phi_range=phi_range, center=center) - def __call__(self, data2D: Data2D) -> npt.NDArray[np.bool_]: + def __call__(self, data2D: SasData) -> npt.NDArray[np.bool_]: """ - Find a rectangular 2D region of interest where data points inside the ROI are True, and False otherwise + Find a rectangular 2D region of interest where data points inside the ROI are True, and False otherwise - :param data2D: Data2D object + :param data2D: SasData object :return: mask, 1d array (len = len(data)) """ - super().validate_and_assign_data(data2D) + self.validate_and_assign_data(data2D) # Ensure unmasked data is used for the phi_data calculation to ensure data sizes match self.phi_data = np.arctan2(data2D.qy_data, data2D.qx_data) - # Calculate q_data using unmasked qx_data and qy_data to ensure data sizes match - q_data = np.sqrt(data2D.qx_data * data2D.qx_data + data2D.qy_data * data2D.qy_data) phi_offset = self.phi_min self.phi_min = 0.0 - self.phi_max = (self.phi_max - phi_offset) % (TwoPi) - self.phi_data = (self.phi_data - phi_offset) % (TwoPi) + self.phi_max = (self.phi_max - phi_offset) % TwoPi + self.phi_data = (self.phi_data - phi_offset) % TwoPi phi_shifted = self.phi_data - Pi # Determine angular bounds for both upper and lower half of image diff --git a/sasdata/data_util/binning.py b/sasdata/data_util/binning.py index fdacda103..737e8f044 100644 --- a/sasdata/data_util/binning.py +++ b/sasdata/data_util/binning.py @@ -11,7 +11,7 @@ class DirectionalAverage: Average along one coordinate axis of 2D data and return data for a 1D plot. This can also be thought of as a projection onto the major axis: 2D -> 1D. - This class operates on a decomposed Data2D object, and returns data needed + This class operates on a decomposed SasData object, and returns data needed to construct a Data1D object. The class is instantiated with two arrays of orthogonal coordinate data (depending on the coordinate system, these may have undergone some pre-processing) and two corresponding two-element @@ -20,7 +20,7 @@ class DirectionalAverage: the other is divided into bins and becomes the dependent variable of the eventual 1D plot. These are called the minor and major axes respectively. When a class instance is called, it is passed the intensity and error data - from the original Data2D object. These should not have undergone any + from the original SasData object. These should not have undergone any coordinate system dependent pre-processing. Note that the old version of manipulations.py had an option for logarithmic @@ -203,7 +203,7 @@ def __call__(self, data, err_data): """ Compute the directional average of the supplied intensity & error data. - :param data: intensity data from the origninal Data2D object. + :param data: intensity data from the original SasData object. :param err_data: the corresponding errors for the intensity data. """ weights = self.compute_weights() diff --git a/sasdata/data_util/manipulations.py b/sasdata/data_util/manipulations.py index 315824612..5b0d9285f 100644 --- a/sasdata/data_util/manipulations.py +++ b/sasdata/data_util/manipulations.py @@ -1,11 +1,11 @@ """ Data manipulations for 2D data sets. -Using the meta data information, various types of averaging are performed in Q-space +Using the meta data information, various types of averaging are performed in Q-space. To test this module use: ``` cd test -PYTHONPATH=../src/ python2 -m sasmanipulations.test.utest_averaging DataInfoTests.test_sectorphi_quarter +PYTHONPATH=../src/ python2 -m sasmanipulations.test.utest_averaging DataInfoTests.test_sectorphi_quarter ``` """ ##################################################################### @@ -23,6 +23,8 @@ import numpy as np +from sasdata.data import SasData + ################################################################################ # Backwards-compatible wrappers that delegate to the new implementations # in averaging.py. @@ -42,10 +44,13 @@ SlabY, WedgePhi, WedgeQ, + get_dq_data, ) -from sasdata.dataloader.data_info import Data1D, Data2D +from sasdata.dataloader.data_info import Data2D from sasdata.dataloader.data_info import reader2D_converter as _di_reader2D_converter +from sasdata.dataset_types import one_dim from sasdata.quantities.constants import Pi, TwoPi +from sasdata.quantities.quantity import Quantity warn("sasdata.data_util.manipulations is deprecated. Unless otherwise noted, update your import to " "sasdata.data_util.averaging.", DeprecationWarning, stacklevel=2) @@ -70,10 +75,10 @@ def deduce_qz(qx: float, qy: float, wavelength: float) -> float: def position_and_wavelength_to_q(dx: float, dy: float, detector_distance: float, wavelength: float) -> float: """ - :param dx: x-distance from beam center [mm] - :param dy: y-distance from beam center [mm] - :param detector_distance: sample to detector distance [mm] - :param wavelength: neutron wavelength [nm] + :param dx: x-distance from beam center + :param dy: y-distance from beam center + :param detector_distance: sample to detector distance + :param wavelength: neutron wavelength :return: q-value at the given position """ # Distance from beam center in the plane of detector @@ -233,55 +238,13 @@ def get_pixel_fraction(q_max: float, q_00: float, q_01: float, q_10: float, q_11 return frac_max -def get_dq_data(data2d: Data2D) -> np.array: - ''' - Get the dq for resolution averaging - The pinholes and det. pix contribution present - in both direction of the 2D which must be subtracted when - converting to 1D: dq_overlap should be calculated ideally at - q = 0. Note This method works on only pinhole geometry. - Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. - ''' - z_max = max(data2d.q_data) - z_min = min(data2d.q_data) - dqx_at_z_max = data2d.dqx_data[np.argmax(data2d.q_data)] - dqx_at_z_min = data2d.dqx_data[np.argmin(data2d.q_data)] - dqy_at_z_max = data2d.dqy_data[np.argmax(data2d.q_data)] - dqy_at_z_min = data2d.dqy_data[np.argmin(data2d.q_data)] - # Find qdx at q = 0 - dq_overlap_x = (dqx_at_z_min * z_max - dqx_at_z_max * z_min) / (z_max - z_min) - # when extrapolation goes wrong - if dq_overlap_x > min(data2d.dqx_data): - dq_overlap_x = min(data2d.dqx_data) - dq_overlap_x *= dq_overlap_x - # Find qdx at q = 0 - dq_overlap_y = (dqy_at_z_min * z_max - dqy_at_z_max * z_min) / (z_max - z_min) - # when extrapolation goes wrong - if dq_overlap_y > min(data2d.dqy_data): - dq_overlap_y = min(data2d.dqy_data) - # get dq at q=0. - dq_overlap_y *= dq_overlap_y - - dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) - # Final protection of dq - if dq_overlap < 0: - dq_overlap = dqy_at_z_min - dqx_data = data2d.dqx_data[np.isfinite(data2d.data)] - dqy_data = data2d.dqy_data[np.isfinite( - data2d.data)] - dq_overlap - # def; dqx_data = dq_r dqy_data = dq_phi - # Convert dq 2D to 1D here - dq_data = np.sqrt(dqx_data**2 + dqy_data**2) - return dq_data - ################################################################################ def reader2D_converter(data2d: Data2D | None = None) -> Data2D: """ - convert old 2d format opened by IhorReader or danse_reader - to new Data2D format - This is mainly used by the Readers + Convert old 2d format opened by IhorReader or danse_reader to new Data2D format. + This is mainly used by the Readers. :param data2d: 2d array of Data2D object :return: 1d arrays of Data2D object @@ -358,8 +321,10 @@ def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, width = max(abs(x_max - x_min), 1e-12) nbins = int(math.ceil(width / abs(bin_width))) if bin_width != 0 else 1 self.nbins=nbins - super().__init__(qx_range=(x_min, x_max), qy_range=(y_min, y_max), - nbins=nbins, fold=fold) + super().__init__(qx_range=(x_min, x_max), + qy_range=(y_min, y_max), + nbins=nbins, + fold=fold) class SlabY(SlabY): @@ -382,8 +347,10 @@ def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, height = max(abs(y_max - y_min), 1e-12) nbins = int(math.ceil(height / abs(bin_width))) if bin_width != 0 else 1 self.nbins=nbins - super().__init__(qx_range=(x_min, x_max), qy_range=(y_min, y_max), - nbins=nbins, fold=fold) + super().__init__(qx_range=(x_min, x_max), + qy_range=(y_min, y_max), + nbins=nbins, + fold=fold) ################################################################################ @@ -391,13 +358,11 @@ def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, class Boxsum(Boxsum): def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): - super().__init__(qx_range=(x_min, x_max), qy_range=(y_min, y_max)) class Boxavg(Boxavg): def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): - super().__init__(qx_range=(x_min, x_max), qy_range=(y_min, y_max)) ################################################################################ @@ -422,23 +387,22 @@ def __call__(self, data2D, ismask=False): """ Perform circular averaging on the data - :param data2D: Data2D object - :return: Data1D object + :param data2D: SasData object + :return: SasData object """ # Get data W/ finite values - data = data2D.data[np.isfinite(data2D.data)] - q_data = data2D.q_data[np.isfinite(data2D.data)] - err_data = None - if data2D.err_data is not None: - err_data = data2D.err_data[np.isfinite(data2D.data)] - mask_data = data2D.mask[np.isfinite(data2D.data)] + finite_mask = np.isfinite(data2D.ordinate.value) + data = data2D.ordinate.value[finite_mask] + q_data = np.sqrt(data2D._data_contents["Qx"].value**2 + data2D._data_contents["Qy"].value**2)[finite_mask] + err_data = np.sqrt(data2D.ordinate.variance.value)[finite_mask] + mask_data = (data2D.mask if data2D.mask is not None else np.ones_like(data2D.ordinate.value, dtype=bool))[finite_mask] dq_data = None - if data2D.dqx_data is not None and data2D.dqy_data is not None: + if data2D._data_contents["Qx"].has_variance and data2D._data_contents["Qy"].has_variance: dq_data = get_dq_data(data2D) if len(q_data) == 0: - msg = "Circular averaging: invalid q_data: %g" % data2D.q_data + msg = "Circular averaging: invalid q_data: %g" % q_data raise RuntimeError(msg) # Build array of Q intervals @@ -508,25 +472,26 @@ def __call__(self, data2D, ismask=False): idx = (np.isfinite(y)) & (np.isfinite(x)) if err_x is not None: - d_x = err_x[idx] / y_counts[idx] + dQ = err_x[idx] / y_counts[idx] else: - d_x = None + dQ = None if not idx.any(): msg = "Average Error: No points inside ROI to average..." raise ValueError(msg) - return Data1D(x=x[idx], y=y[idx], dy=err_y[idx], dx=d_x) + data_contents = { + "Q": Quantity(x[idx], data2D._data_contents["Qx"].units, dQ), + "I": Quantity(y[idx], data2D.ordinate.units, err_y[idx]), + } + return SasData("Circular Average", data_contents, one_dim, data2D.metadata) ################################################################################ class Ring(Ring): - """ - Wrapper for new Ring. - """ - + """Wrapper for new Ring.""" @property def nbins_phi(self): @@ -537,12 +502,9 @@ def nbins_phi(self, value): self.nbins = value def __init__(self, r_min=0.0, r_max=0.0, center_x=0.0, center_y=0.0, nbins=36): - - super().__init__( - r_range=(r_min, r_max), - center=(center_x, center_y), - nbins=nbins - ) + super().__init__(r_range=(r_min, r_max), + center=(center_x, center_y), + nbins=nbins) class _Sector: """ @@ -556,12 +518,10 @@ class _Sector: starting from the negative x-axis. """ - def __init__(self, r_min, r_max, phi_min=0, phi_max=2 * math.pi, nbins=20, - base=None): - ''' - :param base: must be a valid base for an algorithm, i.e., - a positive number - ''' + def __init__(self, r_min, r_max, phi_min=0, phi_max=2*math.pi, nbins=20, base=None): + """ + :param base: must be a valid base for an algorithm, i.e., a positive number + """ self.r_min = r_min self.r_max = r_max self.phi_min = phi_min @@ -576,26 +536,28 @@ def _agv(self, data2D, run='phi'): """ Perform sector averaging. - :param data2D: Data2D object - :param run: define the varying parameter ('phi' , or 'sector') + :param data2D: SasData object + :param run: define the varying parameter ('phi' or 'sector') - :return: Data1D object + :return: SasData object """ - if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: - raise RuntimeError("Ring averaging only take plottable_2D objects") + if not ("Qx" in data2D._data_contents and + "Qy" in data2D._data_contents): + raise RuntimeError("For averaging the SasData object must contain 'Qx' and 'Qy' data.") # Get all the data & info - data = data2D.data[np.isfinite(data2D.data)] - q_data = data2D.q_data[np.isfinite(data2D.data)] - err_data=None - if data2D.err_data is not None: - err_data = data2D.err_data[np.isfinite(data2D.data)] - qx_data = data2D.qx_data[np.isfinite(data2D.data)] - qy_data = data2D.qy_data[np.isfinite(data2D.data)] - mask_data = data2D.mask[np.isfinite(data2D.data)] + finite_mask = np.isfinite(data2D.ordinate.value) + data = data2D.ordinate.value[finite_mask] + err_data = np.sqrt(data2D.ordinate.variance.value)[finite_mask] + qx_data = data2D._data_contents["Qx"].value[finite_mask] + qy_data = data2D._data_contents["Qy"].value[finite_mask] + q_data = np.sqrt(data2D._data_contents["Qx"].value ** 2 + + data2D._data_contents["Qy"].value ** 2 + )[finite_mask] + mask_data = (data2D.mask if data2D.mask is not None else np.ones_like(data2D.ordinate.value, dtype=bool))[finite_mask] dq_data = None - if data2D.dqx_data is not None and data2D.dqy_data is not None: + if data2D._data_contents["Qx"].has_variance and data2D._data_contents["Qy"].has_variance: dq_data = get_dq_data(data2D) # set space for 1d outputs @@ -753,13 +715,19 @@ def _agv(self, data2D, run='phi'): idx = (np.isfinite(y) & np.isfinite(y_err)) if x_err is not None: - d_x = x_err[idx] / y_counts[idx] + dQ = x_err[idx] / y_counts[idx] else: - d_x = None + dQ = None if not idx.any(): msg = "Average Error: No points inside sector of ROI to average..." raise ValueError(msg) - return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x) + + data_contents = { + "Q": Quantity(x[idx], data2D._data_contents["Qx"].units, dQ), + "I": Quantity(y[idx], data2D.ordinate.units, y_err[idx]), + } + return SasData("agv", data_contents, one_dim, data2D.metadata) + class SectorPhi(_Sector): @@ -775,8 +743,8 @@ def __call__(self, data2D): """ Perform sector average and return I(phi). - :param data2D: Data2D object - :return: Data1D object + :param data2D: SasData object + :return: SasData object """ return self._agv(data2D, 'phi') @@ -798,38 +766,28 @@ def __call__(self, data2D): """ Perform sector average and return I(Q). - :param data2D: Data2D object + :param data2D: SasData object - :return: Data1D object + :return: SasData object """ return self._agv(data2D, 'sector') class WedgePhi(WedgePhi): - """ - Wrapper for new WedgePhi (behaviour matches legacy WedgePhi expectations). - """ + """Wrapper for new WedgePhi (behaviour matches legacy WedgePhi expectations).""" def __init__(self, r_min, r_max, phi_min=0, phi_max=TwoPi, center_x=0.0, center_y=0.0, nbins=10): - - super().__init__( - r_range=(r_min, r_max), - phi_range=(phi_min, phi_max), - center=(center_x, center_y), - nbins=nbins - ) + super().__init__(r_range=(r_min, r_max), + phi_range=(phi_min, phi_max), + center=(center_x, center_y), + nbins=nbins) class WedgeQ(WedgeQ): - """ - Wrapper for new WedgeQ (behaviour matches legacy WedgeQ expectations). - """ + """Wrapper for new WedgeQ (behaviour matches legacy WedgeQ expectations).""" def __init__(self, r_min, r_max, phi_min=0, phi_max=TwoPi, center_x=0.0, center_y=0.0, nbins=10): - - super().__init__( - r_range=(r_min, r_max), - phi_range=(phi_min, phi_max), - center=(center_x, center_y), - nbins=nbins - ) + super().__init__(r_range=(r_min, r_max), + phi_range=(phi_min, phi_max), + center=(center_x, center_y), + nbins=nbins) ################################################################################ @@ -837,12 +795,9 @@ def __init__(self, r_min, r_max, phi_min=0, phi_max=TwoPi, center_x=0.0, center_ class Ringcut(Ringcut): def __init__(self, r_min=0.0, r_max=0.0, center_x=0.0, center_y=0.0): # center_x, center_y ignored for compatibility - - super().__init__( - r_range=(r_min, r_max), - phi_range=(0.0, TwoPi), - center=(center_x, center_y) - ) + super().__init__(r_range=(r_min, r_max), + phi_range=(0.0, TwoPi), + center=(center_x, center_y)) ################################################################################ diff --git a/sasdata/data_util/roi.py b/sasdata/data_util/roi.py index 73cee5c2d..eb2752a13 100644 --- a/sasdata/data_util/roi.py +++ b/sasdata/data_util/roi.py @@ -1,18 +1,18 @@ import numpy as np -from sasdata.dataloader.data_info import Data2D +from sasdata.data import SasData from sasdata.quantities.constants import TwoPi class GenericROI: """ - Base class used to set up the data from a Data2D object for processing. + Base class used to set up the data from a SasData object for processing. This class performs any coordinate system independent setup and validation. """ def __init__(self, center: tuple[float, float] = (0.0, 0.0)): """ - Assign the variables used to label the properties of the Data2D object. + Assign the variables used to label the properties of the SasData object. In classes inheriting from GenericROI, the variables used to define the boundaries of the Region Of Interest are also set up during __init__. @@ -26,106 +26,100 @@ def __init__(self, center: tuple[float, float] = (0.0, 0.0)): self.qx_data = None self.qy_data = None - def validate_and_assign_data(self, data2d: Data2D = None) -> None: + def validate_and_assign_data(self, data2d: SasData = None) -> None: """ Check that the data supplied is valid and assign data to variables. - This method must be executed before any further data processing happens + This method must be executed before any further data processing happens. - :param data2d: A Data2D object which is the target of a child class' + :param data2d: A SasData object which is the target of a child class' data manipulations. """ - # Check that the supplied data2d is valid and usable. - if not isinstance(data2d, Data2D): - msg = "Data supplied must be of type Data2D." + # Check that the supplied data is valid and usable. + if not isinstance(data2d, SasData): + msg = "Data supplied must be of type SasData." raise TypeError(msg) - if len(data2d.detector) > 1: - msg = f"Invalid number of detectors: {len(data2d.detector)}" + if not ("Qx" in data2d._data_contents and + "Qy" in data2d._data_contents): + msg = "SasData object must contain 'Qx' and 'Qy' data." + raise TypeError(msg) + if len(data2d.metadata.instrument.detector) > 1: + msg = (f"Invalid number of detectors: {len(data2d.metadata.instrument.detector)}." + "Cannot have more than 1 detector.") raise ValueError(msg) # Only use data which is finite and not masked off - valid_data = np.isfinite(data2d.data) & data2d.mask + if data2d.mask is not None: + valid_data = np.isfinite(data2d.ordinate.value) & data2d.mask + else: + valid_data = np.isfinite(data2d.ordinate.value) - # Assign properties of the Data2D object to variables for reference + # Assign properties of the SasData object to variables for reference # during data processing. - self.data = data2d.data[valid_data] - self.err_data = data2d.err_data[valid_data] + self.data = data2d.ordinate.value[valid_data] + self.err_data = np.sqrt(data2d.ordinate.variance.value)[valid_data] - self.qx_data = data2d.qx_data[valid_data]-self.center_x - self.qy_data = data2d.qy_data[valid_data]-self.center_y - self.q_data = np.sqrt(self.qx_data * self.qx_data + self.qy_data * self.qy_data) + self.qx_data = data2d._data_contents["Qx"].value[valid_data] - self.center_x + self.qy_data = data2d._data_contents["Qy"].value[valid_data] - self.center_y + self.q_data = np.sqrt(self.qx_data ** 2 + self.qy_data ** 2) # Compute phi in the legacy convention: atan2(qy,qx) + pi # (legacy code used this origin; keeping it here makes all polar - # averaging implementations agree and restores the tests). + # averaging implementations agree and restores the tests). self.phi_data = np.arctan2(self.qy_data, self.qx_data) + np.pi # No points should have zero error, if they do then assume the error is # the square root of the data. This code was added to replicate # previous functionality. It's a bit dodgy, so feel free to remove. - self.err_data[self.err_data == 0] = \ - np.sqrt(np.abs(self.data[self.err_data == 0])) + self.err_data[self.err_data == 0] = np.sqrt(np.abs(self.data[self.err_data == 0])) class CartesianROI(GenericROI): - """ - Base class for data manipulators with a Cartesian (rectangular) ROI. - """ + """Base class for data manipulators with a Cartesian (rectangular) ROI.""" def __init__(self, qx_range: tuple[float, float] = (0.0, 0.0), qy_range: tuple[float, float] = (0.0, 0.0)) -> None: """ - Assign the variables used to label the properties of the Data2D object. + Assign the variables used to label the properties of the SasData object. Also establish the upper and lower bounds defining the ROI. - The units of these parameters are A^-1 :param qx_range: Bounds of the ROI along the Q_x direction. :param qy_range: Bounds of the ROI along the Q_y direction. """ + super().__init__() qx_min, qx_max = qx_range qy_min, qy_max = qy_range - super().__init__() self.qx_min = qx_min self.qx_max = qx_max self.qy_min = qy_min self.qy_max = qy_max class PolarROI(GenericROI): - """ - Base class for data manipulators with a polar ROI. - """ + """Base class for data manipulators with a polar ROI.""" - def __init__(self, r_range: tuple[float, float], phi_range: tuple[float, float] = (0.0, TwoPi), center: tuple[float, float] = (0.0, 0.0)) -> None: + def __init__(self, + r_range: tuple[float, float], + phi_range: tuple[float, float] = (0.0, TwoPi), + center: tuple[float, float] = (0.0, 0.0) + ) -> None: """ - Assign the variables used to label the properties of the Data2D object. + Assign the variables used to label the properties of the SasData object. Also establish the upper and lower bounds defining the ROI. - The units are A^-1 for radial parameters, and radians for anglar ones. :param r_range: Tuple (r_min, r_max) defining limits for |Q| values to use during averaging. - :param phi_range: Tuple (phi_min, phi_max) defining limits for φ in radians (in the ROI). + :param phi_range: Tuple (phi_min, phi_max) defining limits for φ in the ROI. Note that Phi is measured anti-clockwise from the positive x-axis. """ - r_min, r_max = r_range - phi_min, phi_max = phi_range super().__init__(center = center) self.phi_data = None + r_min, r_max = r_range + phi_min, phi_max = phi_range + if r_min >= r_max: msg = "Minimum radius cannot be greater than maximum radius." raise ValueError(msg) - # Units A^-1 for radii, radians for angles + self.r_min = r_min self.r_max = r_max self.phi_min = phi_min self.phi_max = phi_max - - def validate_and_assign_data(self, data2d: Data2D = None) -> None: - """ - Check that the data supplied valid and assign data variables. - This method must be executed before any further data processing happens - - :param data2d: A Data2D object which is the target of a child class' - data manipulations. - """ - - # Most validation and pre-processing is taken care of by GenericROI. - super().validate_and_assign_data(data2d) diff --git a/sasdata/dataset_types.py b/sasdata/dataset_types.py index cbe36825d..fc1a3169b 100644 --- a/sasdata/dataset_types.py +++ b/sasdata/dataset_types.py @@ -42,13 +42,21 @@ class DatasetType: ["Qx", "Qy", "Qz", "I", "dI"], ["Qx", "Qy", "Qz", "dQx", "dQy", "dQz", "I", "dI"]]) +angle_dim = DatasetType( + name="I vs Phi", + required=["Phi", "I"], + optional=["dI", "dPhi", "Shadowfactor"], + expected_orders=[ + ["Phi", "I", "dI"], + ["Phi", "dPhi", "I", "dI"]]) + sesans = DatasetType( name="SESANS", required=["SpinEchoLength", "Depolarisation", "Wavelength"], optional=["Transmission", "Polarisation"], expected_orders=[["z", "G"]]) -dataset_types = {dataset.name for dataset in [one_dim, two_dim, sesans]} +dataset_types = {dataset.name for dataset in [one_dim, two_dim, angle_dim, sesans]} # diff --git a/test/sasmanipulations/helper.py b/test/sasmanipulations/helper.py index 21d0a288b..cc1c8025b 100644 --- a/test/sasmanipulations/helper.py +++ b/test/sasmanipulations/helper.py @@ -4,35 +4,41 @@ import numpy as np from scipy import integrate -from sasdata.dataloader import data_info +from sasdata.data import SasData +from sasdata.dataset_types import two_dim +from sasdata.metadata import Instrument, Metadata, Source from sasdata.quantities.constants import TwoPi +from sasdata.quantities.quantity import Quantity +from sasdata.quantities.units import angstroms, per_angstrom, per_centimeter def make_dd_from_func(func, matrix_size=201): """ - Create a MatrixToData2D from a function of (x, y). Returns the MatrixToData2D - instance and matrix_size for convenience. + Create a MatrixToSasData from a function of (x, y). + Returns the MatrixToSasData instance and matrix_size for convenience. func should accept (x, y) meshgrid arrays and return a 2D array. """ x, y = np.meshgrid(np.linspace(-1, 1, matrix_size), np.linspace(-1, 1, matrix_size)) mat = func(x, y) - return MatrixToData2D(data2d=mat), matrix_size + return MatrixToSasData(data2d=mat), matrix_size def make_uniform_dd(shape=(100, 100), value=1.0): mat = np.full(shape, value, dtype=float) - return MatrixToData2D(data2d=mat) + return MatrixToSasData(data2d=mat) def integrate_1d_output(output, method="simpson"): """ Integrate output from an averager consistently. - - If output is a Data1D-like object with .x and .y -> integrate y(x) + - If output is a SasData-like object with "Q" and "I" -> integrate I(Q) - If output is a tuple (result, error[, npoints]) -> return numeric result """ - if hasattr(output, "x") and hasattr(output, "y"): + if (hasattr(output, "_data_contents") and + "Q" in output._data_contents and + "I" in output._data_contents): if method == "trapezoid": - return integrate.trapezoid(output.y, output.x) - return integrate.simpson(output.y, output.x) + return integrate.trapezoid(output._data_contents["I"].value, output._data_contents["Q"].value) + return integrate.simpson(output._data_contents["I"].value, output._data_contents["Q"].value) if isinstance(output, tuple) and len(output) >= 1: return output[0] raise TypeError("Unsupported averager output type: %r" % type(output)) @@ -53,13 +59,13 @@ def expected_slaby_area(qx_min, qx_max, qy_min, qy_max): return x_part_avg * y_part_integ def make_uniform_dd(shape=(100, 100), value=1.0): - """Convenience for tests that need a constant matrix Data2D.""" + """Convenience for tests that need a constant matrix SasData.""" mat = np.full(shape, value, dtype=float) - return MatrixToData2D(data2d=mat) + return MatrixToSasData(data2d=mat) def run_and_integrate(averager, dd, integrator="simpson"): """ - Run an averager (callable) with a Data2D container returned by MatrixToData2D + Run an averager (callable) with a SasData container returned by MatrixToSasData and return the integrated result (scalar area / sum) consistently. """ out = averager(dd.data) @@ -90,9 +96,9 @@ def expected_boxavg_and_err(matrix, slice_rows=None, slice_cols=None): return avg, err -class MatrixToData2D: +class MatrixToSasData: """ - Create Data2D objects from supplied 2D arrays of data. + Create 2D SasData objects from supplied 2D arrays of data. Error data can also be included. Adapted from sasdata.data_util.manipulations.reader_2D_converter @@ -105,18 +111,43 @@ def __init__(self, data2d=None, err_data=None): # qmax can be any number, 1 just makes things simple. self.qmax = 1 - # Creating a Data2D object to use for testing the averagers. - self.data = data_info.Data2D(data=data_flat, err_data=err_flat, - qx_data=qx_data, qy_data=qy_data, - q_data=q_data, mask=mask) + + # Create a SasData object to use for testing the averagers. + data_contents = { + "Qx": Quantity(qx_data, per_angstrom), + "Qy": Quantity(qy_data, per_angstrom), + "I": Quantity(data_flat, per_centimeter), + "dI": Quantity(err_flat, per_centimeter) + } + + wavelength = Quantity(1., angstroms) + source = Source(radiation=None, + beam_shape=None, + beam_size=None, + wavelength=wavelength, + wavelength_max=None, + wavelength_min=None, + wavelength_spread=None) + instrument = Instrument(collimations=[], + source=source, + detector=[]) + metadata=Metadata(title=None, + run=[], + definition=None, + process=[], + sample=None, + instrument=instrument, + raw=None) + + self.data = SasData("Matrix Data", data_contents, two_dim, metadata) def _validate_and_convert_inputs(self, data2d, err_data): """Validate inputs and coerce to numpy arrays. Returns (matrix, err_data_or_None).""" if data2d is None: - raise ValueError("Data must be supplied to convert to Data2D") + raise ValueError("Data must be supplied to convert to SasData") matrix = np.asarray(data2d) if matrix.ndim != 2: - raise ValueError("Supplied array must have 2 dimensions to convert to Data2D") + raise ValueError("Supplied array must have 2 dimensions to convert to SasData") if err_data is not None: err_arr = np.asarray(err_data) @@ -133,16 +164,6 @@ def _compute_bins(self, matrix_shape): qx_bins = np.linspace(start=-1 * 1, stop=1, num=cols, endpoint=True) qy_bins = np.linspace(start=-1 * 1, stop=1, num=rows, endpoint=True) return qx_bins, qy_bins - # qmax can be any number, 1 just makes things simple. - self.qmax = 1 - qx_bins = np.linspace(start=-1 * self.qmax, - stop=self.qmax, - num=matrix.shape[1], - endpoint=True) - qy_bins = np.linspace(start=-1 * self.qmax, - stop=self.qmax, - num=matrix.shape[0], - endpoint=True) def _build_flat_arrays(self, matrix, err_arr, qx_bins, qy_bins): """Flatten matrix and build qx, qy, q arrays plus mask and error handling.""" diff --git a/test/sasmanipulations/utest_averaging.py b/test/sasmanipulations/utest_averaging.py index c8281f6ee..5b6b02eaa 100644 --- a/test/sasmanipulations/utest_averaging.py +++ b/test/sasmanipulations/utest_averaging.py @@ -3,7 +3,7 @@ import numpy as np -import sasdata.dataloader.data_info as data_info +from sasdata.data import SasData, sasdata_reader2D_converter from sasdata.data_util.manipulations import ( Boxavg, Boxsum, @@ -14,10 +14,14 @@ SlabX, SlabY, position_and_wavelength_to_q, - reader2D_converter, ) -from sasdata.dataloader.loader import Loader +from sasdata.dataset_types import two_dim +from sasdata.metadata import Detector, Instrument, Metadata, Source, Vec3 from sasdata.quantities.constants import Pi, TwoPi +from sasdata.quantities.quantity import Quantity +from sasdata.quantities.units import angstroms, millimeters, none, per_angstrom, per_centimeter +from sasdata.temp_ascii_reader import load_data_default_params as ascii_load_data +from sasdata.temp_hdf5_reader import load_data as hdf_load_data def find(filename): @@ -25,9 +29,7 @@ def find(filename): class Averaging(unittest.TestCase): - """ - Test averaging manipulations on a flat distribution - """ + """Test averaging manipulations on a flat distribution.""" def setUp(self): """ @@ -37,25 +39,45 @@ def setUp(self): x_0 = np.ones([100, 100]) dx_0 = np.ones([100, 100]) - self.data = data_info.Data2D(data=x_0, err_data=dx_0) - detector = data_info.Detector() - detector.distance = 1000.0 # mm - detector.pixel_size.x = 1.0 # mm - detector.pixel_size.y = 1.0 # mm - - # center in pixel position = (len(x_0)-1)/2 - detector.beam_center.x = (len(x_0) - 1) / 2 # pixel number - detector.beam_center.y = (len(x_0) - 1) / 2 # pixel number - self.data.detector.append(detector) - - source = data_info.Source() - source.wavelength = 10.0 # A - self.data.source = source + data_contents = { + "Qx": Quantity(np.arange(100), per_angstrom), + "Qy": Quantity(np.arange(100), per_angstrom), + "I": Quantity(x_0, per_centimeter), + "dI": Quantity(dx_0, per_centimeter) + } + + wavelength = Quantity(10.0, angstroms) + source = Source(radiation=None, + beam_shape=None, + beam_size=None, + wavelength=wavelength, + wavelength_max=None, + wavelength_min=None, + wavelength_spread=None) + detector = Detector(name = None, + distance = Quantity(1000.0, millimeters), + offset = None, + orientation = None, + beam_center = Vec3(x=Quantity(0.5 * (len(x_0) - 1), none), y=Quantity(0.5 * (len(x_0) - 1), none), z=None), + pixel_size = Vec3(x=Quantity(1.0, millimeters), y=Quantity(1.0, millimeters), z=None), + slit_length = None) + instrument = Instrument(collimations=[], + source=source, + detector=[detector]) + metadata=Metadata(title=None, + run=[], + definition=None, + process=[], + sample=None, + instrument=instrument, + raw=None) + + self.data = SasData("Test Averaging", data_contents, two_dim, metadata) # get_q(dx, dy, det_dist, wavelength) where units are mm,mm,mm,and A # respectively. - self.qmin = position_and_wavelength_to_q(1.0, 1.0, detector.distance, source.wavelength) - self.qmax = position_and_wavelength_to_q(49.5, 49.5, detector.distance, source.wavelength) + self.qmin = position_and_wavelength_to_q(1.0, 1.0, detector.distance.value, source.wavelength.value) + self.qmax = position_and_wavelength_to_q(49.5, 49.5, detector.distance.value, source.wavelength.value) self.qstep = len(x_0) x = np.linspace(start=-1 * self.qmax, @@ -68,35 +90,31 @@ def setUp(self): endpoint=True) self.data.x_bins = x self.data.y_bins = y - self.data = reader2D_converter(self.data) + self.data = sasdata_reader2D_converter(self.data) def test_ring_flat_distribution(self): - """ - Test ring averaging - """ - r = Ring(r_min=2 * self.qmin, r_max=5 * self.qmin, - center_x=self.data.detector[0].beam_center.x, - center_y=self.data.detector[0].beam_center.y) + """Test ring averaging.""" + r = Ring(r_min=2*self.qmin, + r_max=5*self.qmin, + center_x=self.data.metadata.instrument.detector[0].beam_center.x.value, + center_y=self.data.metadata.instrument.detector[0].beam_center.y.value) r.nbins_phi = 20 o = r(self.data) for i in range(20): - self.assertEqual(o.y[i], 1.0) + self.assertEqual(o._data_contents["I"].value[i], 1.0) def test_sectorphi_full(self): - """ - Test sector averaging - """ + """Test sector averaging.""" r = SectorPhi(r_min=self.qmin, r_max=3 * self.qmin, phi_min=0, phi_max=TwoPi) r.nbins_phi = 20 o = r(self.data) for i in range(7): - self.assertEqual(o.y[i], 1.0) + self.assertEqual(o._data_contents["I"].value[i], 1.0) def test_sectorphi_partial(self): - """ - """ + """Test sector averaging.""" phi_max = Pi * 1.5 r = SectorPhi(r_min=self.qmin, r_max=3 * self.qmin, phi_min=0, phi_max=phi_max) @@ -105,36 +123,42 @@ def test_sectorphi_partial(self): o = r(self.data) self.assertEqual(r.phi_max, phi_max) for i in range(17): - self.assertEqual(o.y[i], 1.0) + self.assertEqual(o._data_contents["I"].value[i], 1.0) class DataInfoTests(unittest.TestCase): def setUp(self): filepath = find('MAR07232_rest.h5') - self.data_list = Loader().load(filepath) - self.data = self.data_list[0] + data_dict = hdf_load_data(filepath) + self.data = data_dict['sasentry01'] def test_ring(self): - """ - Test ring averaging - """ + """Test ring averaging.""" + if beam_center := self.data.metadata.instrument.detector[0].beam_center: + center_x = beam_center.x.value + center_y = beam_center.y.value + else: + center_x = 0.0 + center_y = 0.0 + r = Ring(r_min=.005, r_max=.01, - center_x=self.data.detector[0].beam_center.x, - center_y=self.data.detector[0].beam_center.y, + center_x=center_x, + center_y=center_y, nbins=20) r.nbins_phi = 20 o = r(self.data) filepath = find('ring_testdata.txt') - answer_list = Loader().load(filepath) + answer_list = ascii_load_data(filepath) answer = answer_list[0] self.assertEqual(len(answer_list), 1) for i in range(r.nbins_phi - 1): - self.assertAlmostEqual(o.x[i + 1], answer.x[i], 4) - self.assertAlmostEqual(o.y[i + 1], answer.y[i], 4) - self.assertAlmostEqual(o.dy[i + 1], answer.dy[i], 4) + # Current ascii reader implementation assumes file data is "one_dim" + self.assertAlmostEqual(o._data_contents["Phi"].value[i], answer._data_contents["Q"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4) def test_circularavg(self): """ @@ -148,11 +172,11 @@ def test_circularavg(self): o = r(self.data) filepath = find('avg_testdata.txt') - answer = Loader().load(filepath)[0] + answer = ascii_load_data(filepath)[0] for i in range(r.nbins_phi): - self.assertAlmostEqual(o.x[i], answer.x[i], delta=1e-4) - self.assertAlmostEqual(o.y[i], answer.y[i], delta=1e-4) - self.assertAlmostEqual(o.dy[i], answer.dy[i],delta=1e-4) + self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], delta=1e-4) + self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], delta=1e-4) + self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], delta=1e-4) def test_box(self): """ @@ -183,11 +207,11 @@ def test_slabX(self): o = r(self.data) filepath = find('slabx_testdata.txt') - answer = Loader().load(filepath)[0] - for i in range(len(o.x)): - self.assertAlmostEqual(o.x[i], answer.x[i], 4) - self.assertAlmostEqual(o.y[i], answer.y[i], 4) - self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) + answer = ascii_load_data(filepath)[0] + for i in range(len(o._data_contents["Q"].value)): + self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4) def test_slabY(self): """ @@ -201,11 +225,11 @@ def test_slabY(self): o = r(self.data) filepath = find('slaby_testdata.txt') - answer = Loader().load(filepath)[0] - for i in range(len(o.x)): - self.assertAlmostEqual(o.x[i], answer.x[i], 4) - self.assertAlmostEqual(o.y[i], answer.y[i], 4) - self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) + answer = ascii_load_data(filepath)[0] + for i in range(len(o._data_contents["Q"].value)): + self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4) def test_sectorphi_full(self): """ @@ -227,11 +251,11 @@ def test_sectorphi_full(self): o = r(self.data) filepath = find('ring_testdata.txt') - answer = Loader().load(filepath)[0] - for i in range(len(o.x)): - self.assertAlmostEqual(o.x[i], answer.x[i], 4) - self.assertAlmostEqual(o.y[i], answer.y[i], 4) - self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) + answer = ascii_load_data(filepath)[0] + for i in range(len(o._data_contents["Q"].value)-1): + self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i+1], 4) + self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i+1], 4) + self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i+1], 4) def test_sectorphi_quarter(self): """ @@ -244,11 +268,11 @@ def test_sectorphi_quarter(self): o = r(self.data) filepath = find('sectorphi_testdata.txt') - answer = Loader().load(filepath)[0] - for i in range(len(o.x)): - self.assertAlmostEqual(o.x[i], answer.x[i], 4) - self.assertAlmostEqual(o.y[i], answer.y[i], 4) - self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) + answer = ascii_load_data(filepath)[0] + for i in range(len(o._data_contents["Q"].value)): + self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4) def test_sectorq_full(self): """ @@ -261,11 +285,11 @@ def test_sectorq_full(self): o = r(self.data) filepath = find('sectorq_testdata.txt') - answer = Loader().load(filepath)[0] - for i in range(len(o.x)): - self.assertAlmostEqual(o.x[i], answer.x[i], 4) - self.assertAlmostEqual(o.y[i], answer.y[i], 4) - self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) + answer = ascii_load_data(filepath)[0] + for i in range(len(o._data_contents["Q"].value)): + self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4) + self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4) def test_sectorq_log(self): """ @@ -278,8 +302,8 @@ def test_sectorq_log(self): o = r(self.data) expected_binning = np.logspace(np.log10(0.005), np.log10(0.01), 20, base=10) - for i in range(len(o.x)): - self.assertAlmostEqual(o.x[i], expected_binning[i], 3) + for i in range(len(o._data_contents["Q"].value)): + self.assertAlmostEqual(o._data_contents["Q"].value[i], expected_binning[i], 3) # TODO: Test for Y values (o.y) # print len(self.data.x_bins) diff --git a/test/sasmanipulations/utest_averaging_box.py b/test/sasmanipulations/utest_averaging_box.py index f3b4ec4a4..806d5ccad 100644 --- a/test/sasmanipulations/utest_averaging_box.py +++ b/test/sasmanipulations/utest_averaging_box.py @@ -6,9 +6,9 @@ import numpy as np from sasdata.data_util.averaging import Boxavg, Boxsum -from sasdata.dataloader import data_info +from sasdata.metadata import Detector from test.sasmanipulations.helper import ( - MatrixToData2D, + MatrixToSasData, expected_boxavg_and_err, expected_boxsum_and_err, make_uniform_dd, @@ -43,10 +43,18 @@ def test_boxsum_multiple_detectors(self): Test Boxsum raises an error when there are multiple detectors. """ dd = make_uniform_dd((100, 100), value=1.0) - detector1 = data_info.Detector() - detector2 = data_info.Detector() - dd.data.detector.append(detector1) - dd.data.detector.append(detector2) + + detector = Detector( + name = None, + distance = None, + offset = None, + orientation = None, + beam_center = None, + pixel_size = None, + slit_length = None) + + dd.data.metadata.instrument.detector.append(detector) + dd.data.metadata.instrument.detector.append(detector) box_object = Boxsum() self.assertRaises(ValueError, box_object, dd.data) @@ -58,7 +66,7 @@ def test_boxsum_total(self): # Creating a 100x100 matrix for a distribution which is flat in y # and linear in x. test_data = np.tile(np.arange(100), (100, 1)) - dd = MatrixToData2D(data2d=test_data) + dd = MatrixToSasData(data2d=test_data) box_object = Boxsum(qx_range=(-1 * dd.qmax, dd.qmax), qy_range=(-1 * dd.qmax, dd.qmax)) result, error, npoints = box_object(dd.data) @@ -74,7 +82,7 @@ def test_boxsum_subset_total(self): # Creating a 100x100 matrix for a distribution which is flat in y # and linear in x. test_data = np.tile(np.arange(100), (100, 1)) - dd = MatrixToData2D(data2d=test_data) + dd = MatrixToSasData(data2d=test_data) # region corresponds to central 50x50 in original test box_object = Boxsum(qx_range=(-0.5 * dd.qmax, 0.5 * dd.qmax), qy_range=(-0.5 * dd.qmax, 0.5 * dd.qmax)) @@ -91,7 +99,7 @@ def test_boxsum_zero_sum(self): """ test_data = np.ones([100, 100]) test_data[25:75, 25:75] = 0 - dd = MatrixToData2D(data2d=test_data) + dd = MatrixToSasData(data2d=test_data) box_object = Boxsum(qx_range=(-0.5 * dd.qmax, 0.5 * dd.qmax), qy_range=(-0.5 * dd.qmax, 0.5 * dd.qmax)) result, error, npoints = box_object(dd.data) @@ -127,10 +135,18 @@ def test_boxavg_multiple_detectors(self): Test Boxavg raises an error when there are multiple detectors. """ dd = make_uniform_dd((100, 100), value=1.0) - detector1 = data_info.Detector() - detector2 = data_info.Detector() - dd.data.detector.append(detector1) - dd.data.detector.append(detector2) + + detector = Detector( + name = None, + distance = None, + offset = None, + orientation = None, + beam_center = None, + pixel_size = None, + slit_length = None) + + dd.data.metadata.instrument.detector.append(detector) + dd.data.metadata.instrument.detector.append(detector) box_object = Boxavg() self.assertRaises(ValueError, box_object, dd.data) @@ -142,7 +158,7 @@ def test_boxavg_total(self): # Creating a 100x100 matrix for a distribution which is flat in y # and linear in x. test_data = np.tile(np.arange(100), (100, 1)) - dd = MatrixToData2D(data2d=test_data) + dd = MatrixToSasData(data2d=test_data) box_object = Boxavg(qx_range=(-1 * dd.qmax, dd.qmax), qy_range=(-1 * dd.qmax, dd.qmax)) result, error = box_object(dd.data) @@ -158,7 +174,7 @@ def test_boxavg_subset_total(self): # Creating a 100x100 matrix for a distribution which is flat in y # and linear in x. test_data = np.tile(np.arange(100), (100, 1)) - dd = MatrixToData2D(data2d=test_data) + dd = MatrixToSasData(data2d=test_data) box_object = Boxavg(qx_range=(-0.5 * dd.qmax, 0.5 * dd.qmax), qy_range=(-0.5 * dd.qmax, 0.5 * dd.qmax)) result, error = box_object(dd.data) @@ -175,7 +191,7 @@ def test_boxavg_zero_average(self): test_data = np.ones([100, 100]) # Make a hole in the middle with zeros test_data[25:75, 25:75] = np.zeros([50, 50]) - dd = MatrixToData2D(data2d=test_data) + dd = MatrixToSasData(data2d=test_data) box_object = Boxavg(qx_range=(-0.5 * dd.qmax, 0.5 * dd.qmax), qy_range=(-0.5 * dd.qmax, 0.5 * dd.qmax)) result, error = box_object(dd.data) diff --git a/test/sasmanipulations/utest_averaging_circle.py b/test/sasmanipulations/utest_averaging_circle.py index 36723c203..adcecb2a2 100644 --- a/test/sasmanipulations/utest_averaging_circle.py +++ b/test/sasmanipulations/utest_averaging_circle.py @@ -8,7 +8,7 @@ from sasdata.data_util.averaging import CircularAverage, Ring, SectorQ, WedgePhi, WedgeQ from sasdata.quantities.constants import Pi -from test.sasmanipulations.helper import CircularTestingMatrix, MatrixToData2D +from test.sasmanipulations.helper import CircularTestingMatrix, MatrixToSasData # TODO - also check the errors are being calculated correctly @@ -56,7 +56,7 @@ def test_circularaverage_check_q_data(self): Check CircularAverage ensures the data supplied has `q_data` populated """ # test_data = np.ones([100, 100]) - # averager_data = DataMatrixToData2D(test_data) + # averager_data = DataMatrixToSasData(test_data) # # Overwrite q_data so it's empty # averager_data.data.q_data = np.array([]) # circ_object = CircularAverage() @@ -76,7 +76,7 @@ def test_circularaverage_no_points_to_average(self): Test CircularAverage raises ValueError when the ROI contains no data """ test_data = np.ones([100, 100]) - averager_data = MatrixToData2D(test_data) + averager_data = MatrixToSasData(test_data) # Region of interest well outside region with data circ_object = CircularAverage(r_range=(2 * averager_data.qmax,3 * averager_data.qmax)) @@ -88,7 +88,7 @@ def test_circularaverage_averages_circularly(self): """ test_data = CircularTestingMatrix(frequency=2, matrix_size=201, major_axis='Q') - averager_data = MatrixToData2D(test_data.matrix) + averager_data = MatrixToSasData(test_data.matrix) # Test the ability to average over a subsection of the data r_min = averager_data.qmax * 0.25 @@ -99,7 +99,7 @@ def test_circularaverage_averages_circularly(self): data1d = circ_object(averager_data.data) expected_area = test_data.area_under_region(r_min=r_min, r_max=r_max) - actual_area = integrate.trapezoid(data1d.y, data1d.x) + actual_area = integrate.trapezoid(data1d._data_contents["I"].value, data1d._data_contents["Q"].value) # This used to be able to pass with a precision of 3 d.p. with the old # manipulations.py - I'm not sure why it doesn't anymore. @@ -146,7 +146,7 @@ def test_ring_no_points_to_average(self): Test Ring raises ValueError when the ROI contains no data """ test_data = np.ones([100, 100]) - averager_data = MatrixToData2D(test_data) + averager_data = MatrixToSasData(test_data) # Region of interest well outside region with data ring_object = Ring(r_range=(2 * averager_data.qmax, 3 * averager_data.qmax)) @@ -158,7 +158,7 @@ def test_ring_averages_azimuthally(self): """ test_data = CircularTestingMatrix(frequency=1, matrix_size=201, major_axis='Phi') - averager_data = MatrixToData2D(test_data.matrix) + averager_data = MatrixToSasData(test_data.matrix) # Test the ability to average over a subsection of the data r_min = 0.25 * averager_data.qmax @@ -169,7 +169,7 @@ def test_ring_averages_azimuthally(self): data1d = ring_object(averager_data.data) expected_area = test_data.area_under_region(r_min=r_min, r_max=r_max) - actual_area = integrate.simpson(data1d.y, data1d.x) + actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Phi"].value) self.assertAlmostEqual(actual_area, expected_area, 1) @@ -219,13 +219,13 @@ def test_sectorq_averaging_without_fold(self): """ test_data = CircularTestingMatrix(frequency=1, matrix_size=201, major_axis='Q') - averager_data = MatrixToData2D(test_data.matrix) + averager_data = MatrixToSasData(test_data.matrix) r_min = 0 r_max = 0.9 * averager_data.qmax - phi_min = Pi/6 - phi_max = 5*Pi/6 - nbins = int(test_data.matrix_size * np.sqrt(2)/4) # usually reliable + phi_min = Pi / 6.0 + phi_max = 5.0 * Pi / 6.0 + nbins = int(0.25 * test_data.matrix_size * np.sqrt(2)) # usually reliable wedge_object = SectorQ(r_range=(r_min, r_max), phi_range=(phi_min,phi_max), nbins=nbins) # Explicitly set fold to False - results span full +/- range @@ -241,7 +241,7 @@ def test_sectorq_averaging_without_fold(self): expected_area += test_data.area_under_region(r_min=r_min, r_max=r_max, phi_min=phi_min+Pi, phi_max=phi_max+Pi) - actual_area = integrate.simpson(data1d.y, data1d.x) + actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Q"].value) self.assertAlmostEqual(actual_area, expected_area, 1) @@ -252,13 +252,13 @@ def test_sectorq_averaging_with_fold(self): """ test_data = CircularTestingMatrix(frequency=1, matrix_size=201, major_axis='Q') - averager_data = MatrixToData2D(test_data.matrix) + averager_data = MatrixToSasData(test_data.matrix) r_min = 0 r_max = 0.9 * averager_data.qmax - phi_min = Pi/6 - phi_max = 5*Pi/6 - nbins = int(test_data.matrix_size * np.sqrt(2)/4) # usually reliable + phi_min = Pi / 6.0 + phi_max = 5.0 * Pi / 6.0 + nbins = int(0.25 * test_data.matrix_size * np.sqrt(2)) # usually reliable wedge_object = SectorQ(r_range=(r_min, r_max), phi_range=(phi_min,phi_max), nbins=nbins) # Explicitly set fold to True - points either side of 0,0 are averaged @@ -275,7 +275,7 @@ def test_sectorq_averaging_with_fold(self): phi_min=phi_min+Pi, phi_max=phi_max+Pi) expected_area /= 2 - actual_area = integrate.simpson(data1d.y, data1d.x) + actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Q"].value) self.assertAlmostEqual(actual_area, expected_area, 1) @@ -313,21 +313,20 @@ def test_wedgeq_averaging(self): """ test_data = CircularTestingMatrix(frequency=3, matrix_size=201, major_axis='Q') - averager_data = MatrixToData2D(test_data.matrix) + averager_data = MatrixToSasData(test_data.matrix) r_min = 0.1 * averager_data.qmax r_max = 0.9 * averager_data.qmax - phi_min = Pi/6 - phi_max = 5*Pi/6 + phi_min = Pi / 6.0 + phi_max = 5.0 * Pi / 6.0 nbins = int(test_data.matrix_size * np.sqrt(2)/4) # usually reliable wedge_object = WedgeQ(r_range=(r_min, r_max), phi_range=(phi_min,phi_max), nbins=nbins) data1d = wedge_object(averager_data.data) expected_area = test_data.area_under_region(r_min=r_min, r_max=r_max, - phi_min=phi_min, - phi_max=phi_max) - actual_area = integrate.simpson(data1d.y, data1d.x) + phi_min=phi_min, phi_max=phi_max) + actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Q"].value) self.assertAlmostEqual(actual_area, expected_area, 1) @@ -351,8 +350,6 @@ def test_wedgephi_init(self): nbins = 100 # base = 10 - # wedge_object = WedgePhi(r_min=r_min, r_max=r_max, phi_min=phi_min, - # phi_max=phi_max, nbins=nbins, base=base) wedge_object = WedgePhi(r_range=(r_min, r_max), phi_range=(phi_min,phi_max), nbins=nbins) self.assertEqual(wedge_object.r_min, r_min) @@ -376,21 +373,21 @@ def test_wedgephi_averaging(self): """ test_data = CircularTestingMatrix(frequency=1, matrix_size=201, major_axis='Phi') - averager_data = MatrixToData2D(test_data.matrix) + averager_data = MatrixToSasData(test_data.matrix) r_min = 0.1 * averager_data.qmax r_max = 0.9 * averager_data.qmax - phi_min = Pi/6 - phi_max = 5*Pi/6 - nbins = int(test_data.matrix_size * np.sqrt(2)/4) # usually reliable + phi_min = Pi / 6.0 + phi_max = 5.0 * Pi / 6.0 + nbins = int(0.25 *test_data.matrix_size * np.sqrt(2)) # usually reliable wedge_object = WedgePhi(r_range=(r_min, r_max), phi_range=(phi_min,phi_max), nbins=nbins) data1d = wedge_object(averager_data.data) expected_area = test_data.area_under_region(r_min=r_min, r_max=r_max, - phi_min=phi_min, - phi_max=phi_max) - actual_area = integrate.simpson(data1d.y, data1d.x) + phi_min=phi_min, phi_max=phi_max) + + actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Phi"].value) self.assertAlmostEqual(actual_area, expected_area, 1) diff --git a/test/sasmanipulations/utest_averaging_directional.py b/test/sasmanipulations/utest_averaging_directional.py index 7d8dc5174..7fbb5607c 100644 --- a/test/sasmanipulations/utest_averaging_directional.py +++ b/test/sasmanipulations/utest_averaging_directional.py @@ -9,7 +9,7 @@ import numpy as np from sasdata.data_util.binning import DirectionalAverage -from test.sasmanipulations.helper import MatrixToData2D +from test.sasmanipulations.helper import MatrixToSasData # TODO - also check the errors are being calculated correctly @@ -73,7 +73,7 @@ def setUp(self): x, y = np.meshgrid(self.qx_data, self.qy_data) # quadratic in x, linear in y data = x * x * y - self.data2d = MatrixToData2D(data) + self.data2d = MatrixToSasData(data) # ROI is the first quadrant only. Same limits for both axes. self.lims = (0.0, 1.0) @@ -85,9 +85,10 @@ def setUp(self): self.bin_width = (self.lims[1] - self.lims[0]) / self.nbins self.directional_average = \ - DirectionalAverage(major_axis=self.data2d.data.qx_data, - minor_axis=self.data2d.data.qy_data, - lims=(self.lims,self.lims), nbins=self.nbins) + DirectionalAverage(major_axis=self.data2d.data._data_contents["Qx"].value, + minor_axis=self.data2d.data._data_contents["Qy"].value, + lims=(self.lims,self.lims), + nbins=self.nbins) def test_bin_width(self): """ @@ -140,8 +141,8 @@ def test_directional_averaging(self): the bins. """ x_axis_values, intensity, errors = \ - self.directional_average(data=self.data2d.data.data, - err_data=self.data2d.data.err_data) + self.directional_average(data=self.data2d.data._data_contents["I"].value, + err_data=self.data2d.data._data_contents["dI"].value) expected_x = self.qx_data[self.in_roi] expected_intensity = np.mean(self.qy_data[self.in_roi]) * expected_x**2 @@ -156,8 +157,10 @@ def test_no_points_in_roi(self): # move the region of interest to outside the range of the data self.directional_average.major_lims = (2, 3) self.directional_average.minor_lims = (2, 3) - self.assertRaises(ValueError, self.directional_average, - self.data2d.data.data, self.data2d.data.err_data) + self.assertRaises(ValueError, + self.directional_average, + self.data2d.data._data_contents["I"].value, + self.data2d.data._data_contents["dI"].value) if __name__ == '__main__': unittest.main() diff --git a/test/sasmanipulations/utest_averaging_slab.py b/test/sasmanipulations/utest_averaging_slab.py index da6485867..b1a701dc1 100644 --- a/test/sasmanipulations/utest_averaging_slab.py +++ b/test/sasmanipulations/utest_averaging_slab.py @@ -6,9 +6,9 @@ import numpy as np from sasdata.data_util.averaging import SlabX, SlabY -from sasdata.dataloader import data_info +from sasdata.metadata import Detector from test.sasmanipulations.helper import ( - MatrixToData2D, + MatrixToSasData, expected_slabx_area, expected_slaby_area, integrate_1d_output, @@ -47,11 +47,19 @@ def test_slabx_multiple_detectors(self): """ Test that SlabX raises an error when there are multiple detectors """ - averager_data = MatrixToData2D(np.ones([100, 100])) - detector1 = data_info.Detector() - detector2 = data_info.Detector() - averager_data.data.detector.append(detector1) - averager_data.data.detector.append(detector2) + averager_data = MatrixToSasData(np.ones([100, 100])) + + detector = Detector( + name = None, + distance = None, + offset = None, + orientation = None, + beam_center = None, + pixel_size = None, + slit_length = None) + + averager_data.data.metadata.instrument.detector.append(detector) + averager_data.data.metadata.instrument.detector.append(detector) slab_object = SlabX() self.assertRaises(ValueError, slab_object, averager_data.data) @@ -61,7 +69,7 @@ def test_slabx_no_points_to_average(self): Test SlabX raises ValueError when the ROI contains no data """ test_data = np.ones([100, 100]) - averager_data = MatrixToData2D(data2d=test_data) + averager_data = MatrixToSasData(data2d=test_data) # Region of interest well outside region with data qx_min = 2 * averager_data.qmax @@ -155,11 +163,20 @@ def test_slaby_multiple_detectors(self): """ Test that SlabY raises an error when there are multiple detectors """ - averager_data = MatrixToData2D(np.ones([100, 100])) - detector1 = data_info.Detector() - detector2 = data_info.Detector() - averager_data.data.detector.append(detector1) - averager_data.data.detector.append(detector2) + averager_data = MatrixToSasData(np.ones([100, 100])) + + detector = Detector( + name = None, + distance = None, + offset = None, + orientation = None, + beam_center = None, + pixel_size = None, + slit_length = None) + + averager_data.data.metadata.instrument.detector.append(detector) + averager_data.data.metadata.instrument.detector.append(detector) + slab_object = SlabY() self.assertRaises(ValueError, slab_object, averager_data.data) @@ -168,7 +185,7 @@ def test_slaby_no_points_to_average(self): Test SlabY raises ValueError when the ROI contains no data """ test_data = np.ones([100, 100]) - averager_data = MatrixToData2D(data2d=test_data) + averager_data = MatrixToSasData(data2d=test_data) # Region of interest well outside region with data qx_min = 2 * averager_data.qmax diff --git a/test/utest_new_sasdata.py b/test/utest_new_sasdata.py index cc497d9ae..5ee073b32 100644 --- a/test/utest_new_sasdata.py +++ b/test/utest_new_sasdata.py @@ -2,11 +2,12 @@ from sasdata.data import SasData from sasdata.data_backing import Group -from sasdata.dataset_types import one_dim, three_dim, two_dim +from sasdata.dataset_types import angle_dim, one_dim, three_dim, two_dim from sasdata.metadata import Instrument, Metadata, Source from sasdata.postprocess import deduce_qz +from sasdata.quantities.constants import Pi from sasdata.quantities.quantity import Quantity -from sasdata.quantities.units import angstroms, per_angstrom, per_centimeter +from sasdata.quantities.units import angstroms, per_angstrom, per_centimeter, radians def test_1d(): @@ -105,3 +106,20 @@ def test_3d(): deduce_qz(data) assert (data._data_contents['Qz'].value != (0*data._data_contents['Qx'].value)).all() + +def test_angle(): + phi = [0.4 * Pi, 0.8 * Pi, 1.2 * Pi, 1.6 * Pi, 2 * Pi] + i = [5, 4, 3, 2, 1] + + phi_quantity = Quantity(np.array(phi), radians) + i_quantity = Quantity(np.array(i), per_centimeter) + + data_contents = { + 'Phi': phi_quantity, + 'I': i_quantity + } + + data = SasData('TestData', data_contents, angle_dim, Group('root', {}), True) + + assert all(data.abscissae.value == np.array(phi)) + assert all(data.ordinate.value == np.array(i))