convolution

Convolution

Bases: NumericalConvolutionBase

Convolution class that combines analytical and numerical convolution methods to efficiently perform convolutions of ComponentCollections with ResolutionComponents. Supports analytical convolution for pairs of analytical model components (DeltaFunction, Gaussian, Lorentzian, Voigt), while using numerical convolution for other components. If temperature is provided, detailed balance correction is applied to the sample model. In this case, all convolutions are handled numerically. Includes a setting to normalize the detailed balance correction. Includes optional upsampling and extended range to improve accuracy of the numerical convolutions. Also warns about numerical instabilities if peaks are very wide or very narrow.

energy : np.ndarray or scipp.Variable 1D array of energy values where the convolution is evaluated. sample_components : ComponentCollection or ModelComponent The sample components to be convolved. resolution_components : ComponentCollection or ModelComponent The resolution components to convolve with. upsample_factor : int, optional The factor by which to upsample the input data before convolution. Default is 5. extension_factor : float, optional The factor by which to extend the input data range before convolution. Default is 0.2. temperature : Parameter, float, or None, optional The temperature to use for detailed balance correction. Default is None. temperature_unit : str or sc.Unit, optional The unit of the temperature parameter. Default is 'K'. energy_unit : str or sc.Unit, optional The unit of the energy. Default is 'meV'. normalize_detailed_balance : bool, optional Whether to normalize the detailed balance correction. Default is True.

Source code in src/easydynamics/convolution/convolution.py

class Convolution(NumericalConvolutionBase):
    """Convolution class that combines analytical and numerical
    convolution methods to efficiently perform convolutions of
    ComponentCollections with ResolutionComponents. Supports analytical
    convolution for pairs of analytical model components (DeltaFunction,
    Gaussian, Lorentzian, Voigt), while using numerical convolution for
    other components. If temperature is provided, detailed balance
    correction is applied to the sample model. In this case, all
    convolutions are handled numerically. Includes a setting to
    normalize the detailed balance correction. Includes optional
    upsampling and extended range to improve accuracy of the numerical
    convolutions. Also warns about numerical instabilities if peaks are
    very wide or very narrow.

    Args:
    energy : np.ndarray or scipp.Variable
        1D array of energy values where the convolution is evaluated.
    sample_components : ComponentCollection or ModelComponent
        The sample components to be convolved.
    resolution_components : ComponentCollection or ModelComponent
        The resolution components to convolve with.
    upsample_factor : int, optional
        The factor by which to upsample the input data before
        convolution. Default is 5.
    extension_factor : float, optional
        The factor by which to extend the input data range before
        convolution. Default is 0.2.
    temperature : Parameter, float, or None, optional
        The temperature to use for detailed balance correction.
        Default is None.
    temperature_unit : str or sc.Unit, optional
        The unit of the temperature parameter. Default is 'K'.
    energy_unit : str or sc.Unit, optional
        The unit of the energy. Default is 'meV'.
    normalize_detailed_balance : bool, optional
        Whether to normalize the detailed balance correction.
        Default is True.
    """

    # When these attributes are changed, the convolution plan
    # needs to be rebuilt
    _invalidate_plan_on_change = {
        'energy',
        '_energy',
        '_energy_grid',
        '_sample_components',
        '_resolution_components',
        '_temperature',
        '_upsample_factor',
        '_extension_factor',
        '_energy_unit',
        '_normalize_detailed_balance',
    }

    def __init__(
        self,
        energy: np.ndarray | sc.Variable,
        sample_components: ComponentCollection | ModelComponent,
        resolution_components: ComponentCollection | ModelComponent,
        upsample_factor: Numerical = 5,
        extension_factor: Numerical = 0.2,
        temperature: Parameter | Numerical | None = None,
        temperature_unit: str | sc.Unit = 'K',
        energy_unit: str | sc.Unit = 'meV',
        normalize_detailed_balance: bool = True,
    ):
        self._convolution_plan_is_valid = False
        self._reactions_enabled = False
        super().__init__(
            energy=energy,
            sample_components=sample_components,
            resolution_components=resolution_components,
            upsample_factor=upsample_factor,
            extension_factor=extension_factor,
            temperature=temperature,
            temperature_unit=temperature_unit,
            energy_unit=energy_unit,
            normalize_detailed_balance=normalize_detailed_balance,
        )

        self._reactions_enabled = True
        # Separate sample model components into pairs that can be
        # handled analytically, delta functions, and the rest
        # Also initialize analytical and numerical convolvers based on s
        # ample model component
        self._build_convolution_plan()

    def convolution(
        self,
    ) -> np.ndarray:
        """Perform convolution using analytical convolutions where
        possible, and numerical convolutions for the remaining
        components.

        Returns:
            np.ndarray
                The convolved values evaluated at energy.
        """
        if not self._convolution_plan_is_valid:
            self._build_convolution_plan()
        total = np.zeros_like(self.energy.values, dtype=float)

        # Analytical convolution
        if self._analytical_convolver is not None:
            total += self._analytical_convolver.convolution()

        # Numerical convolution
        if self._numerical_convolver is not None:
            total += self._numerical_convolver.convolution()

        # Delta function components
        if self._delta_sample_components.components:
            total += self._convolve_delta_functions()

        return total

    def _convolve_delta_functions(self) -> np.ndarray:
        "Convolve delta function components of the sample model with"
        'the resolution components.'
        'No detailed balance correction is applied to delta functions.'
        return sum(
            delta.area.value
            * self._resolution_components.evaluate(self.energy.values - delta.center.value)
            for delta in self._delta_sample_components.components
        )

    def _check_if_pair_is_analytic(
        self,
        sample_component: ModelComponent,
        resolution_component: ModelComponent,
    ) -> bool:
        """Check if the convolution of the given component pair can be
        handled analytically.

        Args:
            sample_component : ModelComponent
                The sample component to be convolved.
            resolution_component : ModelComponent
                The resolution component to convolve with.
        Returns:
            bool
                True if the component pair can be handled analytically,
                False otherwise.
        """

        if not isinstance(sample_component, ModelComponent):
            raise TypeError(
                f'`sample_component` is an instance of {type(sample_component).__name__}, \
                but must be a ModelComponent.'
            )

        if not isinstance(resolution_component, ModelComponent):
            raise TypeError(
                f'`resolution_component` is an instance of {type(resolution_component).__name__}, \
                    but must be a ModelComponent.'
            )

        if isinstance(resolution_component, DeltaFunction):
            raise TypeError(
                'resolution components contains delta functions. This is not supported.'
            )

        analytical_types = (Gaussian, Lorentzian, Voigt)
        if isinstance(sample_component, analytical_types) and isinstance(
            resolution_component, analytical_types
        ):
            return True

        return False

    def _build_convolution_plan(self) -> None:
        """Separate sample model components into analytical pairs, delta
        functions, and the rest.
        """

        analytical_sample_components = ComponentCollection()
        delta_sample_components = ComponentCollection()
        numerical_sample_components = ComponentCollection()

        for sample_component in self._sample_components.components:
            # If delta function, put in delta sample model and go to the
            # next component
            if isinstance(sample_component, DeltaFunction):
                delta_sample_components.append_component(sample_component)
                continue

            # If temperature is set, all other components go to
            # numerical sample model
            if self.temperature is not None:
                numerical_sample_components.append_component(sample_component)
                continue

            # If temperature is not set, check if all
            # resolution components can be convolved analytically with
            # this sample component
            pair_is_analytic = []
            for resolution_component in self._resolution_components.components:
                pair_is_analytic.append(
                    self._check_if_pair_is_analytic(sample_component, resolution_component)
                )
            # If all resolution components can be convolved analytically
            # with this sample component, add it to analytical
            # sample model. If not, it goes to numerical sample model.
            if all(pair_is_analytic):
                analytical_sample_components.append_component(sample_component)
            else:
                numerical_sample_components.append_component(sample_component)

        self._analytical_sample_components = analytical_sample_components
        self._delta_sample_components = delta_sample_components
        self._numerical_sample_components = numerical_sample_components
        self._convolution_plan_is_valid = True

        # Update convolvers
        self._set_convolvers()

    def _set_convolvers(self) -> None:
        """Initialize analytical and numerical convolvers based on
        sample model components.

        There is no delta function convolver, as delta functions are
        handled directly in the convolution method.
        """

        if self._analytical_sample_components.components:
            self._analytical_convolver = AnalyticalConvolution(
                energy=self.energy,
                sample_components=self._analytical_sample_components,
                resolution_components=self._resolution_components,
            )
        else:
            self._analytical_convolver = None

        if self._numerical_sample_components.components:
            self._numerical_convolver = NumericalConvolution(
                energy=self.energy,
                sample_components=self._numerical_sample_components,
                resolution_components=self._resolution_components,
                upsample_factor=self.upsample_factor,
                extension_factor=self.extension_factor,
                temperature=self.temperature,
                temperature_unit=self._temperature_unit,
                normalize_detailed_balance=self.normalize_detailed_balance,
            )
        else:
            self._numerical_convolver = None

    # Update some setters so the internal sample models are updated
    def __setattr__(self, name, value):
        """Custom setattr to invalidate convolution plan on relevant
        attribute changes, and build a new plan.

        The new plan is only built after initialization (when
        _reactions_enabled is True) to avoid issues during __init__.
        """
        super().__setattr__(name, value)

        if name in self._invalidate_plan_on_change:
            self._convolution_plan_is_valid = False

        if getattr(self, '_reactions_enabled', False) and name in self._invalidate_plan_on_change:
            self._build_convolution_plan()

__setattr__(name, value)

Custom setattr to invalidate convolution plan on relevant attribute changes, and build a new plan.

The new plan is only built after initialization (when _reactions_enabled is True) to avoid issues during init.

Source code in src/easydynamics/convolution/convolution.py

def __setattr__(self, name, value):
    """Custom setattr to invalidate convolution plan on relevant
    attribute changes, and build a new plan.

    The new plan is only built after initialization (when
    _reactions_enabled is True) to avoid issues during __init__.
    """
    super().__setattr__(name, value)

    if name in self._invalidate_plan_on_change:
        self._convolution_plan_is_valid = False

    if getattr(self, '_reactions_enabled', False) and name in self._invalidate_plan_on_change:
        self._build_convolution_plan()

convolution()

Perform convolution using analytical convolutions where possible, and numerical convolutions for the remaining components.

Returns:

Type	Description
ndarray	np.ndarray The convolved values evaluated at energy.

Source code in src/easydynamics/convolution/convolution.py

def convolution(
    self,
) -> np.ndarray:
    """Perform convolution using analytical convolutions where
    possible, and numerical convolutions for the remaining
    components.

    Returns:
        np.ndarray
            The convolved values evaluated at energy.
    """
    if not self._convolution_plan_is_valid:
        self._build_convolution_plan()
    total = np.zeros_like(self.energy.values, dtype=float)

    # Analytical convolution
    if self._analytical_convolver is not None:
        total += self._analytical_convolver.convolution()

    # Numerical convolution
    if self._numerical_convolver is not None:
        total += self._numerical_convolver.convolution()

    # Delta function components
    if self._delta_sample_components.components:
        total += self._convolve_delta_functions()

    return total

analytical_convolution

AnalyticalConvolution

Bases: ConvolutionBase

Analytical convolution of a ModelComponent or ComponentCollection with a ResolutionModel.

Possible analytical convolutions are any combination of delta functions, Gaussians, Lorentzians and Voigt profiles. Args: energy : np.ndarray or scipp.Variable 1D array of energy values where the convolution is evaluated. sample_components : ComponentCollection or ModelComponent The sample model to be convolved. resolution_components : ComponentCollection or ModelComponent The resolution model to convolve with.

Source code in src/easydynamics/convolution/analytical_convolution.py

class AnalyticalConvolution(ConvolutionBase):
    """Analytical convolution of a ModelComponent or ComponentCollection
    with a ResolutionModel.

    Possible analytical convolutions are any combination of
    delta functions, Gaussians, Lorentzians and Voigt profiles.
    Args:
        energy : np.ndarray or scipp.Variable
            1D array of energy values where the convolution is
            evaluated.
        sample_components : ComponentCollection or ModelComponent
            The sample model to be convolved.
        resolution_components : ComponentCollection or ModelComponent
            The resolution model to convolve with.
    """

    # Mapping of supported component type pairs to convolution methods.
    # Delta functions are handled separately.
    _CONVOLUTIONS = {
        ('Gaussian', 'Gaussian'): '_convolute_gaussian_gaussian',
        ('Gaussian', 'Lorentzian'): '_convolute_gaussian_lorentzian',
        ('Gaussian', 'Voigt'): '_convolute_gaussian_voigt',
        ('Lorentzian', 'Lorentzian'): '_convolute_lorentzian_lorentzian',
        ('Lorentzian', 'Voigt'): '_convolute_lorentzian_voigt',
        ('Voigt', 'Voigt'): '_convolute_voigt_voigt',
    }

    def __init__(
        self,
        energy: np.ndarray | sc.Variable,
        energy_unit: str | sc.Unit = 'meV',
        sample_components: ComponentCollection | ModelComponent | None = None,
        resolution_components: ComponentCollection | ModelComponent | None = None,
    ):
        super().__init__(
            energy=energy,
            energy_unit=energy_unit,
            sample_components=sample_components,
            resolution_components=resolution_components,
        )

    def convolution(
        self,
    ) -> np.ndarray:
        """Convolve sample with resolution analytically if possible.
        Accepts ComponentCollection or single ModelComponent for each.
        Possible analytical convolutions are any combination of delta
        functions, Gaussians, Lorentzians and Voigt profiles.

        Returns:
            np.ndarray
                The convolution of the sample_components and resolution_
                components values evaluated at energy.

        Raises:
            ValueError
                If resolution_components contains delta functions.
            ValueError
                If component pair cannot be handled analytically.
        """

        # prepare list of components
        if isinstance(self.sample_components, ComponentCollection):
            sample_components = self.sample_components.components
        else:
            sample_components = [self.sample_components]

        if isinstance(self.resolution_components, ComponentCollection):
            resolution_components = self.resolution_components.components
        else:
            resolution_components = [self.resolution_components]

        total = np.zeros_like(self.energy.values, dtype=float)

        for sample_component in sample_components:
            # Go through resolution components,
            # adding analytical contributions
            for resolution_component in resolution_components:
                contrib = self._convolute_analytic_pair(
                    sample_component=sample_component,
                    resolution_component=resolution_component,
                )
                total += contrib

        return total

    def _convolute_analytic_pair(
        self,
        sample_component: ModelComponent,
        resolution_component: ModelComponent,
    ) -> np.ndarray:
        """Analytic convolution for component pair (sample_component,
        resolution_component). The convolution of two gaussian
        components results in another gaussian component with width
        sqrt(w1^2 + w2^2). The convolution of two lorentzian components
        results in another lorentzian component with width w1 + w2. The
        convolution of a gaussian and a lorentzian results in a voigt
        profile. The convolution of a gaussian and a voigt profile
        results in another voigt profile, with the lorentzian width
        unchanged and the gaussian widths summed in quadrature. The
        convolution of a lorentzian and a voigt profile results in
        another voigt profile, with the gaussian width unchanged and the
        lorentzian widths summed. The convolution of two voigt profiles
        results in another voigt profile, with the gaussian widths
        summed in quadrature and the lorentzian widths summed. The
        convolution of a delta function with any component or
        ComponentCollection results in the same component or
        ComponentCollection shifted by the delta center. All areas are
        multiplied.

        Args:
            sample_component : ModelComponent
                The sample component to be convolved.
            resolution_component : ModelComponent
                The resolution component to convolve with.

        Returns:
            np.ndarray: The convolution result

        Raises:
            ValueError:
            If the component pair cannot be handled analytically.
        """

        if isinstance(resolution_component, DeltaFunction):
            raise ValueError(
                'Analytical convolution with a delta function \
                    in the resolution model is not supported.'
            )

        # Delta function + anything -->
        # anything, shifted by delta center with area A1 * A2
        if isinstance(sample_component, DeltaFunction):
            return self._convolute_delta_any(
                sample_component,
                resolution_component,
            )

        pair = (type(sample_component).__name__, type(resolution_component).__name__)
        swapped = False

        if pair not in self._CONVOLUTIONS:
            # Try reversing the pair
            pair = (
                type(resolution_component).__name__,
                type(sample_component).__name__,
            )
            swapped = True

        func_name = self._CONVOLUTIONS.get(pair)

        if func_name is None:
            raise ValueError(
                f'Analytical convolution not supported for component pair: '
                f'{type(sample_component).__name__}, {type(resolution_component).__name__}'
            )

        # Call the corresponding method
        if swapped:
            return getattr(self, func_name)(resolution_component, sample_component)
        else:
            return getattr(self, func_name)(sample_component, resolution_component)

    def _convolute_delta_any(
        self,
        sample_component: DeltaFunction,
        resolution_components: ComponentCollection | ModelComponent,
    ):
        """Convolution of delta function with any ModelComponent or
        ComponentCollection results in the same component or
        ComponentCollection shifted by the delta center. The areas are
        multiplied.

        Args:
            sample_component : DeltaFunction
                The sample component to be convolved.
            resolution_components : ComponentCollection | ModelComponent
                The resolution model to convolve with.
        Returns:
            np.ndarray
                The evaluated convolution values at self.energy.
        """
        return sample_component.area.value * resolution_components.evaluate(
            self.energy.values - sample_component.center.value
        )

    def _convolute_gaussian_gaussian(
        self,
        sample_component: Gaussian,
        resolution_component: Gaussian,
    ) -> np.ndarray:
        """Convolution of two Gaussian components results in another
        Gaussian component with width sqrt(w1^2 + w2^2). The areas are
        multiplied.

        Args:
            sample_component : Gaussian
                The sample Gaussian component to be convolved.
            resolution_component : Gaussian
                The resolution Gaussian component to convolve with.

        Returns:
            np.ndarray

                The evaluated convolution values at self.energy.
        """

        width = np.sqrt(sample_component.width.value**2 + resolution_component.width.value**2)

        area = sample_component.area.value * resolution_component.area.value

        center = sample_component.center.value + resolution_component.center.value

        return self._gaussian_eval(area=area, center=center, width=width)

    def _convolute_gaussian_lorentzian(
        self,
        sample_component: Gaussian,
        resolution_component: Lorentzian,
    ) -> np.ndarray:
        """Convolution of a Gaussian and a Lorentzian results in a Voigt
        profile. The areas are multiplied.

        Args:
            sample_component : Gaussian
                The sample Gaussian component to be convolved.
            resolution_component : Lorentzian
                The resolution Lorentzian component to convolve with.

        Returns:
            np.ndarray
                The evaluated convolution values at self.energy.
        """
        center = sample_component.center.value + resolution_component.center.value
        area = sample_component.area.value * resolution_component.area.value

        return self._voigt_eval(
            area=area,
            center=center,
            gaussian_width=sample_component.width.value,
            lorentzian_width=resolution_component.width.value,
        )

    def _convolute_gaussian_voigt(
        self,
        sample_component: Gaussian,
        resolution_component: Voigt,
    ) -> np.ndarray:
        """Convolution of a Gaussian and a Voigt profile results in
        another Voigt profile. The Lorentzian width remains unchanged,
        while the Gaussian widths are summed in quadrature. The areas
        are multiplied.

        Args:
            sample_component : Gaussian
                The sample Gaussian component to be convolved.
            resolution_component : Voigt
                The resolution Voigt component to convolve with.

        Returns:
            np.ndarray
                The evaluated convolution values at self.energy.
        """
        area = sample_component.area.value * resolution_component.area.value

        center = sample_component.center.value + resolution_component.center.value

        gaussian_width = np.sqrt(
            sample_component.width.value**2 + resolution_component.gaussian_width.value**2
        )

        lorentzian_width = resolution_component.lorentzian_width.value

        return self._voigt_eval(
            area=area,
            center=center,
            gaussian_width=gaussian_width,
            lorentzian_width=lorentzian_width,
        )

    def _convolute_lorentzian_lorentzian(
        self,
        sample_component: Lorentzian,
        resolution_component: Lorentzian,
    ) -> np.ndarray:
        """Convolution of two Lorentzian components results in another
        Lorentzian component with width w1 + w2. The areas are
        multiplied.

        Args:
            sample_component : Lorentzian
                The sample Lorentzian component to be convolved.
            resolution_component : Lorentzian
                The resolution Lorentzian component to convolve with.
        Returns:
            np.ndarray
                The evaluated convolution values at self.energy.
        """
        area = sample_component.area.value * resolution_component.area.value

        center = sample_component.center.value + resolution_component.center.value

        width = sample_component.width.value + resolution_component.width.value

        return self._lorentzian_eval(area=area, center=center, width=width)

    def _convolute_lorentzian_voigt(
        self,
        sample_component: Lorentzian,
        resolution_component: Voigt,
    ) -> np.ndarray:
        """Convolution of a Lorentzian and a Voigt profile results in
        another Voigt profile.

        The Gaussian width remains unchanged, while the Lorentzian
        widths are summed.
        The areas are multiplied.
        Args:
            sample_component : Lorentzian
                The sample Lorentzian component to be convolved.
            resolution_component : Voigt
                The resolution Voigt component to convolve with.
        Returns:
            np.ndarray
                The evaluated convolution values at self.energy.
        """
        area = sample_component.area.value * resolution_component.area.value

        center = sample_component.center.value + resolution_component.center.value

        gaussian_width = resolution_component.gaussian_width.value

        lorentzian_width = (
            sample_component.width.value + resolution_component.lorentzian_width.value
        )

        return self._voigt_eval(
            area=area,
            center=center,
            gaussian_width=gaussian_width,
            lorentzian_width=lorentzian_width,
        )

    def _convolute_voigt_voigt(
        self,
        sample_component: Voigt,
        resolution_component: Voigt,
    ) -> np.ndarray:
        """Convolution of two Voigt profiles results in another Voigt
        profile.

        The Gaussian widths are summed in quadrature,
        while the Lorentzian widths are summed.
        The areas are multiplied.
        Args:
            sample_component : Voigt
                The sample Voigt component to be convolved.
            resolution_component : Voigt
                The resolution Voigt component to convolve with.
        Returns:
            np.ndarray
                The evaluated convolution values at self.energy.
        """
        area = sample_component.area.value * resolution_component.area.value

        center = sample_component.center.value + resolution_component.center.value

        gaussian_width = np.sqrt(
            sample_component.gaussian_width.value**2 + resolution_component.gaussian_width.value**2
        )

        lorentzian_width = (
            sample_component.lorentzian_width.value + resolution_component.lorentzian_width.value
        )
        return self._voigt_eval(
            area=area,
            center=center,
            gaussian_width=gaussian_width,
            lorentzian_width=lorentzian_width,
        )

    def _gaussian_eval(
        self,
        area: float,
        center: float,
        width: float,
    ) -> np.ndarray:
        """Evaluate a Gaussian function. y = (area/(sqrt(2pi) *
        width))*exp(-0.5*((x-center) / width)^2) All checks are handled
        in the calling function.

        Args:
            area : float
                The area under the Gaussian curve.
            center : float
                The center of the Gaussian.
            width : float
                The width (sigma) of the Gaussian.
        Returns:
            np.ndarray
                The evaluated Gaussian values at self.energy.
        """

        normalization = 1 / (np.sqrt(2 * np.pi) * width)
        exponent = -0.5 * ((self.energy.values - center) / width) ** 2

        return area * normalization * np.exp(exponent)

    def _lorentzian_eval(self, area: float, center: float, width: float) -> np.ndarray:
        """
        Evaluate a Lorentzian function.
        y = (area * width / pi) / ((x - center)^2 + width^2).
        All checks are handled in the calling function.

        Args:
            area : float
                The area under the Lorentzian.
            center : float
                The center of the Lorentzian.
            width : float
                The width (HWHM) of the Lorentzian.
        Returns:
            np.ndarray
                The evaluated Lorentzian values at self.energy.
        """

        normalization = width / np.pi
        denominator = (self.energy.values - center) ** 2 + width**2

        return area * normalization / denominator

    def _voigt_eval(
        self,
        area: float,
        center: float,
        gaussian_width: float,
        lorentzian_width: float,
    ) -> np.ndarray:
        """Evaluate a Voigt profile function using scipy's
        voigt_profile.

        Args:
            area : float
                The area under the Voigt profile.
            center : float
                The center of the Voigt profile.
            gaussian_width : float
                The Gaussian width (sigma) of the Voigt profile.
            lorentzian_width : float
                The Lorentzian width (HWHM) of the Voigt profile.
        Returns:
            np.ndarray
                The evaluated Voigt profile values at self.energy.
        """

        return area * voigt_profile(self.energy.values - center, gaussian_width, lorentzian_width)

convolution()

Convolve sample with resolution analytically if possible. Accepts ComponentCollection or single ModelComponent for each. Possible analytical convolutions are any combination of delta functions, Gaussians, Lorentzians and Voigt profiles.

Returns:

Type	Description
ndarray	np.ndarray The convolution of the sample_components and resolution_ components values evaluated at energy.

Source code in src/easydynamics/convolution/analytical_convolution.py

def convolution(
    self,
) -> np.ndarray:
    """Convolve sample with resolution analytically if possible.
    Accepts ComponentCollection or single ModelComponent for each.
    Possible analytical convolutions are any combination of delta
    functions, Gaussians, Lorentzians and Voigt profiles.

    Returns:
        np.ndarray
            The convolution of the sample_components and resolution_
            components values evaluated at energy.

    Raises:
        ValueError
            If resolution_components contains delta functions.
        ValueError
            If component pair cannot be handled analytically.
    """

    # prepare list of components
    if isinstance(self.sample_components, ComponentCollection):
        sample_components = self.sample_components.components
    else:
        sample_components = [self.sample_components]

    if isinstance(self.resolution_components, ComponentCollection):
        resolution_components = self.resolution_components.components
    else:
        resolution_components = [self.resolution_components]

    total = np.zeros_like(self.energy.values, dtype=float)

    for sample_component in sample_components:
        # Go through resolution components,
        # adding analytical contributions
        for resolution_component in resolution_components:
            contrib = self._convolute_analytic_pair(
                sample_component=sample_component,
                resolution_component=resolution_component,
            )
            total += contrib

    return total

convolution

Convolution

Bases: NumericalConvolutionBase

Convolution class that combines analytical and numerical convolution methods to efficiently perform convolutions of ComponentCollections with ResolutionComponents. Supports analytical convolution for pairs of analytical model components (DeltaFunction, Gaussian, Lorentzian, Voigt), while using numerical convolution for other components. If temperature is provided, detailed balance correction is applied to the sample model. In this case, all convolutions are handled numerically. Includes a setting to normalize the detailed balance correction. Includes optional upsampling and extended range to improve accuracy of the numerical convolutions. Also warns about numerical instabilities if peaks are very wide or very narrow.

energy : np.ndarray or scipp.Variable 1D array of energy values where the convolution is evaluated. sample_components : ComponentCollection or ModelComponent The sample components to be convolved. resolution_components : ComponentCollection or ModelComponent The resolution components to convolve with. upsample_factor : int, optional The factor by which to upsample the input data before convolution. Default is 5. extension_factor : float, optional The factor by which to extend the input data range before convolution. Default is 0.2. temperature : Parameter, float, or None, optional The temperature to use for detailed balance correction. Default is None. temperature_unit : str or sc.Unit, optional The unit of the temperature parameter. Default is 'K'. energy_unit : str or sc.Unit, optional The unit of the energy. Default is 'meV'. normalize_detailed_balance : bool, optional Whether to normalize the detailed balance correction. Default is True.

Source code in src/easydynamics/convolution/convolution.py

class Convolution(NumericalConvolutionBase):
    """Convolution class that combines analytical and numerical
    convolution methods to efficiently perform convolutions of
    ComponentCollections with ResolutionComponents. Supports analytical
    convolution for pairs of analytical model components (DeltaFunction,
    Gaussian, Lorentzian, Voigt), while using numerical convolution for
    other components. If temperature is provided, detailed balance
    correction is applied to the sample model. In this case, all
    convolutions are handled numerically. Includes a setting to
    normalize the detailed balance correction. Includes optional
    upsampling and extended range to improve accuracy of the numerical
    convolutions. Also warns about numerical instabilities if peaks are
    very wide or very narrow.

    Args:
    energy : np.ndarray or scipp.Variable
        1D array of energy values where the convolution is evaluated.
    sample_components : ComponentCollection or ModelComponent
        The sample components to be convolved.
    resolution_components : ComponentCollection or ModelComponent
        The resolution components to convolve with.
    upsample_factor : int, optional
        The factor by which to upsample the input data before
        convolution. Default is 5.
    extension_factor : float, optional
        The factor by which to extend the input data range before
        convolution. Default is 0.2.
    temperature : Parameter, float, or None, optional
        The temperature to use for detailed balance correction.
        Default is None.
    temperature_unit : str or sc.Unit, optional
        The unit of the temperature parameter. Default is 'K'.
    energy_unit : str or sc.Unit, optional
        The unit of the energy. Default is 'meV'.
    normalize_detailed_balance : bool, optional
        Whether to normalize the detailed balance correction.
        Default is True.
    """

    # When these attributes are changed, the convolution plan
    # needs to be rebuilt
    _invalidate_plan_on_change = {
        'energy',
        '_energy',
        '_energy_grid',
        '_sample_components',
        '_resolution_components',
        '_temperature',
        '_upsample_factor',
        '_extension_factor',
        '_energy_unit',
        '_normalize_detailed_balance',
    }

    def __init__(
        self,
        energy: np.ndarray | sc.Variable,
        sample_components: ComponentCollection | ModelComponent,
        resolution_components: ComponentCollection | ModelComponent,
        upsample_factor: Numerical = 5,
        extension_factor: Numerical = 0.2,
        temperature: Parameter | Numerical | None = None,
        temperature_unit: str | sc.Unit = 'K',
        energy_unit: str | sc.Unit = 'meV',
        normalize_detailed_balance: bool = True,
    ):
        self._convolution_plan_is_valid = False
        self._reactions_enabled = False
        super().__init__(
            energy=energy,
            sample_components=sample_components,
            resolution_components=resolution_components,
            upsample_factor=upsample_factor,
            extension_factor=extension_factor,
            temperature=temperature,
            temperature_unit=temperature_unit,
            energy_unit=energy_unit,
            normalize_detailed_balance=normalize_detailed_balance,
        )

        self._reactions_enabled = True
        # Separate sample model components into pairs that can be
        # handled analytically, delta functions, and the rest
        # Also initialize analytical and numerical convolvers based on s
        # ample model component
        self._build_convolution_plan()

    def convolution(
        self,
    ) -> np.ndarray:
        """Perform convolution using analytical convolutions where
        possible, and numerical convolutions for the remaining
        components.

        Returns:
            np.ndarray
                The convolved values evaluated at energy.
        """
        if not self._convolution_plan_is_valid:
            self._build_convolution_plan()
        total = np.zeros_like(self.energy.values, dtype=float)

        # Analytical convolution
        if self._analytical_convolver is not None:
            total += self._analytical_convolver.convolution()

        # Numerical convolution
        if self._numerical_convolver is not None:
            total += self._numerical_convolver.convolution()

        # Delta function components
        if self._delta_sample_components.components:
            total += self._convolve_delta_functions()

        return total

    def _convolve_delta_functions(self) -> np.ndarray:
        "Convolve delta function components of the sample model with"
        'the resolution components.'
        'No detailed balance correction is applied to delta functions.'
        return sum(
            delta.area.value
            * self._resolution_components.evaluate(self.energy.values - delta.center.value)
            for delta in self._delta_sample_components.components
        )

    def _check_if_pair_is_analytic(
        self,
        sample_component: ModelComponent,
        resolution_component: ModelComponent,
    ) -> bool:
        """Check if the convolution of the given component pair can be
        handled analytically.

        Args:
            sample_component : ModelComponent
                The sample component to be convolved.
            resolution_component : ModelComponent
                The resolution component to convolve with.
        Returns:
            bool
                True if the component pair can be handled analytically,
                False otherwise.
        """

        if not isinstance(sample_component, ModelComponent):
            raise TypeError(
                f'`sample_component` is an instance of {type(sample_component).__name__}, \
                but must be a ModelComponent.'
            )

        if not isinstance(resolution_component, ModelComponent):
            raise TypeError(
                f'`resolution_component` is an instance of {type(resolution_component).__name__}, \
                    but must be a ModelComponent.'
            )

        if isinstance(resolution_component, DeltaFunction):
            raise TypeError(
                'resolution components contains delta functions. This is not supported.'
            )

        analytical_types = (Gaussian, Lorentzian, Voigt)
        if isinstance(sample_component, analytical_types) and isinstance(
            resolution_component, analytical_types
        ):
            return True

        return False

    def _build_convolution_plan(self) -> None:
        """Separate sample model components into analytical pairs, delta
        functions, and the rest.
        """

        analytical_sample_components = ComponentCollection()
        delta_sample_components = ComponentCollection()
        numerical_sample_components = ComponentCollection()

        for sample_component in self._sample_components.components:
            # If delta function, put in delta sample model and go to the
            # next component
            if isinstance(sample_component, DeltaFunction):
                delta_sample_components.append_component(sample_component)
                continue

            # If temperature is set, all other components go to
            # numerical sample model
            if self.temperature is not None:
                numerical_sample_components.append_component(sample_component)
                continue

            # If temperature is not set, check if all
            # resolution components can be convolved analytically with
            # this sample component
            pair_is_analytic = []
            for resolution_component in self._resolution_components.components:
                pair_is_analytic.append(
                    self._check_if_pair_is_analytic(sample_component, resolution_component)
                )
            # If all resolution components can be convolved analytically
            # with this sample component, add it to analytical
            # sample model. If not, it goes to numerical sample model.
            if all(pair_is_analytic):
                analytical_sample_components.append_component(sample_component)
            else:
                numerical_sample_components.append_component(sample_component)

        self._analytical_sample_components = analytical_sample_components
        self._delta_sample_components = delta_sample_components
        self._numerical_sample_components = numerical_sample_components
        self._convolution_plan_is_valid = True

        # Update convolvers
        self._set_convolvers()

    def _set_convolvers(self) -> None:
        """Initialize analytical and numerical convolvers based on
        sample model components.

        There is no delta function convolver, as delta functions are
        handled directly in the convolution method.
        """

        if self._analytical_sample_components.components:
            self._analytical_convolver = AnalyticalConvolution(
                energy=self.energy,
                sample_components=self._analytical_sample_components,
                resolution_components=self._resolution_components,
            )
        else:
            self._analytical_convolver = None

        if self._numerical_sample_components.components:
            self._numerical_convolver = NumericalConvolution(
                energy=self.energy,
                sample_components=self._numerical_sample_components,
                resolution_components=self._resolution_components,
                upsample_factor=self.upsample_factor,
                extension_factor=self.extension_factor,
                temperature=self.temperature,
                temperature_unit=self._temperature_unit,
                normalize_detailed_balance=self.normalize_detailed_balance,
            )
        else:
            self._numerical_convolver = None

    # Update some setters so the internal sample models are updated
    def __setattr__(self, name, value):
        """Custom setattr to invalidate convolution plan on relevant
        attribute changes, and build a new plan.

        The new plan is only built after initialization (when
        _reactions_enabled is True) to avoid issues during __init__.
        """
        super().__setattr__(name, value)

        if name in self._invalidate_plan_on_change:
            self._convolution_plan_is_valid = False

        if getattr(self, '_reactions_enabled', False) and name in self._invalidate_plan_on_change:
            self._build_convolution_plan()

__setattr__(name, value)

Custom setattr to invalidate convolution plan on relevant attribute changes, and build a new plan.

The new plan is only built after initialization (when _reactions_enabled is True) to avoid issues during init.

Source code in src/easydynamics/convolution/convolution.py

def __setattr__(self, name, value):
    """Custom setattr to invalidate convolution plan on relevant
    attribute changes, and build a new plan.

    The new plan is only built after initialization (when
    _reactions_enabled is True) to avoid issues during __init__.
    """
    super().__setattr__(name, value)

    if name in self._invalidate_plan_on_change:
        self._convolution_plan_is_valid = False

    if getattr(self, '_reactions_enabled', False) and name in self._invalidate_plan_on_change:
        self._build_convolution_plan()

convolution()

Perform convolution using analytical convolutions where possible, and numerical convolutions for the remaining components.

Returns:

Type	Description
ndarray	np.ndarray The convolved values evaluated at energy.

Source code in src/easydynamics/convolution/convolution.py

def convolution(
    self,
) -> np.ndarray:
    """Perform convolution using analytical convolutions where
    possible, and numerical convolutions for the remaining
    components.

    Returns:
        np.ndarray
            The convolved values evaluated at energy.
    """
    if not self._convolution_plan_is_valid:
        self._build_convolution_plan()
    total = np.zeros_like(self.energy.values, dtype=float)

    # Analytical convolution
    if self._analytical_convolver is not None:
        total += self._analytical_convolver.convolution()

    # Numerical convolution
    if self._numerical_convolver is not None:
        total += self._numerical_convolver.convolution()

    # Delta function components
    if self._delta_sample_components.components:
        total += self._convolve_delta_functions()

    return total

convolution_base

ConvolutionBase

Base class for convolutions of sample and resolution models. This base class has no convolution functionality.

energy : np.ndarray or scipp.Variable 1D array of energy values where the convolution is evaluated. sample_components : ComponentCollection or ModelComponent The sample model to be convolved. resolution_components : ComponentCollection or ModelComponent The resolution model to convolve with. energy_unit : str or sc.Unit, optional The unit of the energy. Default is 'meV'.

Source code in src/easydynamics/convolution/convolution_base.py

class ConvolutionBase:
    """Base class for convolutions of sample and resolution models. This
    base class has no convolution functionality.

    Args:
    energy : np.ndarray or scipp.Variable
        1D array of energy values where the convolution is evaluated.
    sample_components : ComponentCollection or ModelComponent
        The sample model to be convolved.
    resolution_components : ComponentCollection or ModelComponent
        The resolution model to convolve with.
    energy_unit : str or sc.Unit, optional
        The unit of the energy. Default is 'meV'.
    """

    def __init__(
        self,
        energy: np.ndarray | sc.Variable,
        sample_components: ComponentCollection | ModelComponent = None,
        resolution_components: ComponentCollection | ModelComponent = None,
        energy_unit: str | sc.Unit = 'meV',
    ):
        if isinstance(energy, Numerical):
            energy = np.array([float(energy)])

        if not isinstance(energy, (np.ndarray, sc.Variable)):
            raise TypeError('Energy must be a numpy ndarray or a scipp Variable.')

        if not isinstance(energy_unit, (str, sc.Unit)):
            raise TypeError('Energy_unit must be a string or sc.Unit.')

        if isinstance(energy, np.ndarray):
            energy = sc.array(dims=['energy'], values=energy, unit=energy_unit)

        self._energy = energy
        self._energy_unit = energy_unit

        if sample_components is not None and not (
            isinstance(sample_components, ComponentCollection)
            or isinstance(sample_components, ModelComponent)
        ):
            raise TypeError(
                f'`sample_components` is an instance of {type(sample_components).__name__}, but must be a ComponentCollection or ModelComponent.'  # noqa: E501
            )
        self._sample_components = sample_components

        if resolution_components is not None and not (
            isinstance(resolution_components, ComponentCollection)
            or isinstance(resolution_components, ModelComponent)
        ):
            raise TypeError(
                f'`resolution_components` is an instance of {type(resolution_components).__name__}, but must be a ComponentCollection or ModelComponent.'  # noqa: E501
            )
        self._resolution_components = resolution_components

    @property
    def energy(self) -> sc.Variable:
        """Get the energy."""

        return self._energy

    @energy.setter
    def energy(self, energy: np.ndarray) -> None:
        """Set the energy.
         Args:
            energy : np.ndarray or scipp.Variable
                1D array of energy values where the convolution is
                evaluated.

        Raises:
            TypeError: If energy is not a numpy ndarray or a
            scipp Variable.
        """

        if isinstance(energy, Numerical):
            energy = np.array([float(energy)])

        if not isinstance(energy, (np.ndarray, sc.Variable)):
            raise TypeError('Energy must be a Number, a numpy ndarray or a scipp Variable.')

        if isinstance(energy, np.ndarray):
            self._energy = sc.array(dims=['energy'], values=energy, unit=self._energy.unit)

        if isinstance(energy, sc.Variable):
            self._energy = energy
            self._energy_unit = energy.unit

    @property
    def energy_unit(self) -> str:
        """Get the energy unit."""
        return self._energy_unit

    @energy_unit.setter
    def energy_unit(self, unit_str: str) -> None:
        raise AttributeError(
            (
                f'Unit is read-only. Use convert_unit to change the unit between allowed types '
                f'or create a new {self.__class__.__name__} with the desired unit.'
            )
        )  # noqa: E501

    def convert_energy_unit(self, energy_unit: str | sc.Unit) -> None:
        """Convert the energy to the specified unit
        Args:
            energy_unit : str or sc.Unit
                The unit of the energy.

        Raises:
            TypeError: If energy_unit is not a string or scipp unit.
        """
        if not isinstance(energy_unit, (str, sc.Unit)):
            raise TypeError('Energy unit must be a string or scipp unit.')

        self.energy = sc.to_unit(self.energy, energy_unit)
        self._energy_unit = energy_unit

    @property
    def sample_components(self) -> ComponentCollection | ModelComponent:
        """Get the sample model."""
        return self._sample_components

    @sample_components.setter
    def sample_components(self, sample_components: ComponentCollection | ModelComponent) -> None:
        """Set the sample model.
        Args:
            sample_components : ComponentCollection or ModelComponent
                The sample model to be convolved.

        Raises:
            TypeError: If sample_components is not a ComponentCollection
            or ModelComponent.
        """
        if not isinstance(sample_components, (ComponentCollection, ModelComponent)):
            raise TypeError(
                f'`sample_components` is an instance of {type(sample_components).__name__}, but must be a ComponentCollection or ModelComponent.'  # noqa: E501
            )
        self._sample_components = sample_components

    @property
    def resolution_components(self) -> ComponentCollection | ModelComponent:
        """Get the resolution model."""
        return self._resolution_components

    @resolution_components.setter
    def resolution_components(
        self, resolution_components: ComponentCollection | ModelComponent
    ) -> None:
        """Set the resolution model.
        Args:
            resolution_components : ComponentCollection or
            ModelComponent
                The resolution model to convolve with.

        Raises:
            TypeError: If resolution_components is not a
            ComponentCollection or ModelComponent.
        """
        if not isinstance(resolution_components, (ComponentCollection, ModelComponent)):
            raise TypeError(
                f'`resolution_components` is an instance of {type(resolution_components).__name__}, but must be a ComponentCollection or ModelComponent.'  # noqa: E501
            )
        self._resolution_components = resolution_components

convert_energy_unit(energy_unit)

Convert the energy to the specified unit Args: energy_unit : str or sc.Unit The unit of the energy.

Raises:

Type	Description
TypeError	If energy_unit is not a string or scipp unit.

Source code in src/easydynamics/convolution/convolution_base.py

def convert_energy_unit(self, energy_unit: str | sc.Unit) -> None:
    """Convert the energy to the specified unit
    Args:
        energy_unit : str or sc.Unit
            The unit of the energy.

    Raises:
        TypeError: If energy_unit is not a string or scipp unit.
    """
    if not isinstance(energy_unit, (str, sc.Unit)):
        raise TypeError('Energy unit must be a string or scipp unit.')

    self.energy = sc.to_unit(self.energy, energy_unit)
    self._energy_unit = energy_unit

energy property writable

Get the energy.

energy_unit property writable

Get the energy unit.

resolution_components property writable

Get the resolution model.

sample_components property writable

Get the sample model.

energy_grid

EnergyGrid dataclass

Container for the dense energy grid and related metadata.

Attributes:

Name	Type	Description
energy_dense		np.ndarray The upsampled and extended energy array.
energy_dense_centered		np.ndarray The centered version of energy_dense (used for resolution evaluation).
energy_dense_step		float The spacing of energy_dense (used for width checks and normalization).
energy_span_dense		float The total span of energy_dense. (used for width checks).
energy_even_length_offset		float The offset to apply if energy_dense has even length (used for convolution alignment).

Source code in src/easydynamics/convolution/energy_grid.py

@dataclass(frozen=True)
class EnergyGrid:
    """Container for the dense energy grid and related metadata.

    Attributes:
        energy_dense : np.ndarray
            The upsampled and extended energy array.
        energy_dense_centered : np.ndarray
            The centered version of energy_dense
            (used for resolution evaluation).
        energy_dense_step : float
            The spacing of energy_dense
            (used for width checks and normalization).
        energy_span_dense : float
            The total span of energy_dense.
            (used for width checks).
        energy_even_length_offset : float
            The offset to apply if energy_dense has even length
            (used for convolution alignment).
    """

    energy_dense: np.ndarray
    energy_dense_centered: np.ndarray
    energy_dense_step: float
    energy_span_dense: float
    energy_even_length_offset: float

numerical_convolution

NumericalConvolution

Bases: NumericalConvolutionBase

Numerical convolution of a ComponentCollection with a ComponentCollection using FFT. Includes optional upsampling and extended range to improve accuracy. Warns about very wide or very narrow peaks in the models. If temperature is provided, detailed balance correction is applied to the sample model.

energy : np.ndarray or scipp.Variable 1D array of energy values where the convolution is evaluated. sample_components : ComponentCollection or ModelComponent The sample model to be convolved. resolution_components : ComponentCollection or ModelComponent The resolution model to convolve with. upsample_factor : int, optional The factor by which to upsample the input data before convolution. Default is 5. extension_factor : float, optional The factor by which to extend the input data range before convolution. Default is 0.2. temperature : Parameter, float, or None, optional The temperature to use for detailed balance correction. Default is None. temperature_unit : str or sc.Unit, optional The unit of the temperature parameter. Default is 'K'. energy_unit : str or sc.Unit, optional The unit of the energy. Default is 'meV'. normalize_detailed_balance : bool, optional Whether to normalize the detailed balance correction. Default is True.

Source code in src/easydynamics/convolution/numerical_convolution.py

class NumericalConvolution(NumericalConvolutionBase):
    """Numerical convolution of a ComponentCollection with a
    ComponentCollection using FFT. Includes optional upsampling and
    extended range to improve accuracy. Warns about very wide or very
    narrow peaks in the models. If temperature is provided, detailed
    balance correction is applied to the sample model.

    Args:
    energy : np.ndarray or scipp.Variable
        1D array of energy values where the convolution is evaluated.
    sample_components : ComponentCollection or ModelComponent
        The sample model to be convolved.
    resolution_components : ComponentCollection or ModelComponent
        The resolution model to convolve with.
    upsample_factor : int, optional
        The factor by which to upsample the input data
        before convolution.
        Default is 5.
    extension_factor : float, optional
        The factor by which to extend the input data range
        before convolution.
        Default is 0.2.
    temperature : Parameter, float, or None, optional
        The temperature to use for detailed balance correction.
          Default is None.
    temperature_unit : str or sc.Unit, optional
        The unit of the temperature parameter. Default is 'K'.
    energy_unit : str or sc.Unit, optional
        The unit of the energy. Default is 'meV'.
    normalize_detailed_balance : bool, optional
        Whether to normalize the detailed balance correction.
        Default is True.
    """

    def __init__(
        self,
        energy: np.ndarray | sc.Variable,
        sample_components: ComponentCollection | ModelComponent,
        resolution_components: ComponentCollection | ModelComponent,
        upsample_factor: Numerical = 5,
        extension_factor: float = 0.2,
        temperature: Parameter | float | None = None,
        temperature_unit: str | sc.Unit = 'K',
        energy_unit: str | sc.Unit = 'meV',
        normalize_detailed_balance: bool = True,
    ):
        super().__init__(
            energy=energy,
            sample_components=sample_components,
            resolution_components=resolution_components,
            upsample_factor=upsample_factor,
            extension_factor=extension_factor,
            temperature=temperature,
            temperature_unit=temperature_unit,
            energy_unit=energy_unit,
            normalize_detailed_balance=normalize_detailed_balance,
        )

    def convolution(
        self,
    ) -> np.ndarray:
        """Calculate the convolution of the sample and resolution models
        at the values given in energy. Includes detailed balance
        correction if temperature is provided.

        Returns:
            np.ndarray
                The convolved values evaluated at energy.
        """

        # Give warnings if peaks are very wide or very narrow
        self._check_width_thresholds(
            model=self.sample_components,
            model_name='sample model',
        )
        self._check_width_thresholds(
            model=self.resolution_components,
            model_name='resolution model',
        )

        # Evaluate sample model. If called via the Convolution class,
        # delta functions are already filtered out.
        sample_vals = self.sample_components.evaluate(
            self._energy_grid.energy_dense - self._energy_grid.energy_even_length_offset
        )

        # Detailed balance correction
        if self.temperature is not None:
            detailed_balance_factor_correction = detailed_balance_factor(
                energy=self._energy_grid.energy_dense,
                temperature=self.temperature,
                energy_unit=self.energy.unit,
                divide_by_temperature=self.normalize_detailed_balance,
            )
            sample_vals *= detailed_balance_factor_correction

        # Evaluate resolution model
        resolution_vals = self.resolution_components.evaluate(
            self._energy_grid.energy_dense_centered
        )

        # Convolution
        convolved = fftconvolve(sample_vals, resolution_vals, mode='same')
        convolved *= self._energy_grid.energy_dense_step  # normalize

        if self.upsample_factor is not None:
            # interpolate back to original energy grid
            convolved = np.interp(
                self.energy.values,
                self._energy_grid.energy_dense,
                convolved,
                left=0.0,
                right=0.0,
            )

        return convolved

convolution()

Calculate the convolution of the sample and resolution models at the values given in energy. Includes detailed balance correction if temperature is provided.

Returns:

Type	Description
ndarray	np.ndarray The convolved values evaluated at energy.

Source code in src/easydynamics/convolution/numerical_convolution.py

def convolution(
    self,
) -> np.ndarray:
    """Calculate the convolution of the sample and resolution models
    at the values given in energy. Includes detailed balance
    correction if temperature is provided.

    Returns:
        np.ndarray
            The convolved values evaluated at energy.
    """

    # Give warnings if peaks are very wide or very narrow
    self._check_width_thresholds(
        model=self.sample_components,
        model_name='sample model',
    )
    self._check_width_thresholds(
        model=self.resolution_components,
        model_name='resolution model',
    )

    # Evaluate sample model. If called via the Convolution class,
    # delta functions are already filtered out.
    sample_vals = self.sample_components.evaluate(
        self._energy_grid.energy_dense - self._energy_grid.energy_even_length_offset
    )

    # Detailed balance correction
    if self.temperature is not None:
        detailed_balance_factor_correction = detailed_balance_factor(
            energy=self._energy_grid.energy_dense,
            temperature=self.temperature,
            energy_unit=self.energy.unit,
            divide_by_temperature=self.normalize_detailed_balance,
        )
        sample_vals *= detailed_balance_factor_correction

    # Evaluate resolution model
    resolution_vals = self.resolution_components.evaluate(
        self._energy_grid.energy_dense_centered
    )

    # Convolution
    convolved = fftconvolve(sample_vals, resolution_vals, mode='same')
    convolved *= self._energy_grid.energy_dense_step  # normalize

    if self.upsample_factor is not None:
        # interpolate back to original energy grid
        convolved = np.interp(
            self.energy.values,
            self._energy_grid.energy_dense,
            convolved,
            left=0.0,
            right=0.0,
        )

    return convolved

numerical_convolution_base

NumericalConvolutionBase

Bases: ConvolutionBase

Base class for numerical convolutions of sample and resolution models. Provides methods to handle upsampling, extension, and detailed balance correction. This base class has no convolution functionality.

energy : np.ndarray or scipp.Variable 1D array of energy values where the convolution is evaluated. sample_components : ComponentCollection or ModelComponent The components to be convolved. resolution_components : ComponentCollection or ModelComponent The resolution components to convolve with. upsample_factor : int, optional The factor by which to upsample the input data before convolution. Default is 5. extension_factor : float, optional The factor by which to extend the input data range before convolution. Default is 0.2. temperature : Parameter, float, or None, optional The temperature to use for detailed balance correction. Default is None. temperature_unit : str or sc.Unit, optional The unit of the temperature parameter. Default is 'K'. energy_unit : str or sc.Unit, optional The unit of the energy. Default is 'meV'. normalize_detailed_balance : bool, optional Whether to normalize the detailed balance correction. Default is True.

Source code in src/easydynamics/convolution/numerical_convolution_base.py

class NumericalConvolutionBase(ConvolutionBase):
    """Base class for numerical convolutions of sample and resolution
    models. Provides methods to handle upsampling, extension, and
    detailed balance correction. This base class has no convolution
    functionality.

    Args:
    energy : np.ndarray or scipp.Variable
        1D array of energy values where the convolution is evaluated.
    sample_components : ComponentCollection or ModelComponent
        The components to be convolved.
    resolution_components : ComponentCollection or ModelComponent
        The resolution components to convolve with.
    upsample_factor : int, optional
        The factor by which to upsample the input data
        before convolution. Default is 5.
    extension_factor : float, optional
        The factor by which to extend the input data range
        before convolution. Default is 0.2.
    temperature : Parameter, float, or None, optional
        The temperature to use for detailed balance correction.
        Default is None.
    temperature_unit : str or sc.Unit, optional
        The unit of the temperature parameter. Default is 'K'.
    energy_unit : str or sc.Unit, optional
        The unit of the energy. Default is 'meV'.
    normalize_detailed_balance : bool, optional
        Whether to normalize the detailed balance correction.
        Default is True.
    """

    def __init__(
        self,
        energy: np.ndarray | sc.Variable,
        sample_components: ComponentCollection | ModelComponent,
        resolution_components: ComponentCollection | ModelComponent,
        upsample_factor: Numerical = 5,
        extension_factor: float = 0.2,
        temperature: Parameter | float | None = None,
        temperature_unit: str | sc.Unit = 'K',
        energy_unit: str | sc.Unit = 'meV',
        normalize_detailed_balance: bool = True,
    ):
        super().__init__(
            energy=energy,
            sample_components=sample_components,
            resolution_components=resolution_components,
            energy_unit=energy_unit,
        )

        if temperature is not None and not isinstance(temperature, (Numerical, Parameter)):
            raise TypeError('Temperature must be None, a number or a Parameter.')

        if not isinstance(temperature_unit, (str, sc.Unit)):
            raise TypeError('Temperature_unit must be a string or sc.Unit.')
        self._temperature_unit = temperature_unit
        self._temperature = None
        self.temperature = temperature

        self._normalize_detailed_balance = normalize_detailed_balance

        self._upsample_factor = upsample_factor
        self._extension_factor = extension_factor

        # Create a dense grid to improve accuracy.
        # When upsample_factor>1, we evaluate on this grid and
        # interpolate back to the original values at the end
        self._energy_grid = self._create_energy_grid()

    @ConvolutionBase.energy.setter
    def energy(self, energy: np.ndarray) -> None:
        ConvolutionBase.energy.fset(self, energy)
        # Recreate dense grid when energy is updated
        self._energy_grid = self._create_energy_grid()

    @property
    def upsample_factor(self) -> Numerical:
        """Get the upsample factor."""

        return self._upsample_factor

    @upsample_factor.setter
    def upsample_factor(self, factor: Numerical) -> None:
        """Set the upsample factor and recreate the dense grid."""
        if factor is None:
            self._upsample_factor = factor
            self._energy_grid = self._create_energy_grid()
            return

        if not isinstance(factor, Numerical):
            raise TypeError('Upsample factor must be a numerical value or None.')
        factor = float(factor)
        if factor <= 1.0:
            raise ValueError('Upsample factor must be greater than 1.')

        self._upsample_factor = factor

        # Recreate dense grid when upsample factor is updated
        self._energy_grid = self._create_energy_grid()

    @property
    def extension_factor(self) -> float:
        """Get the extension factor.

        The extension factor determines how much the energy range is
        extended on both sides before convolution.
        0.2 means extending by 20% of the original energy span
        on each side
        """

        return self._extension_factor

    @extension_factor.setter
    def extension_factor(self, factor: Numerical) -> None:
        """
        Set the extension factor and recreate the dense grid.
        The extension factor determines how much the energy range is
        extended on both sides before convolution.
        0.2 means extending by 20% of the original energy span
        on each side.

        Args:
            factor : float
                The new extension factor.

        Raises:
            TypeError: If factor is not a number.
        """
        if not isinstance(factor, Numerical):
            raise TypeError('Extension factor must be a number.')
        if factor < 0.0:
            raise ValueError('Extension factor must be non-negative.')

        self._extension_factor = factor
        # Recreate dense grid when extension factor is updated
        self._energy_grid = self._create_energy_grid()

    @property
    def temperature(self) -> Optional[Parameter]:
        """Get the temperature."""

        return self._temperature

    @temperature.setter
    def temperature(self, temp: Parameter | float | None) -> None:
        """Set the temperature.

        If None, disables detailed balance
        correction and removes the temperature parameter.
        Args:
            temp : Parameter, float, or None
                The temperature to set. The unit will be the same as
                the existing temperature parameter if it exists,
                otherwise 'K'.
        Raises:
            TypeError: If temp is not a float, Parameter, or None.
        """

        if temp is None:
            self._temperature = None
        elif isinstance(temp, Numerical):
            if self._temperature is not None:
                self._temperature.value = float(temp)
            else:
                self._temperature = Parameter(
                    name='temperature',
                    value=float(temp),
                    unit=self._temperature_unit,
                    fixed=True,
                )
        elif isinstance(temp, Parameter):
            self._temperature = temp
        else:
            raise TypeError('Temperature must be None, a float or a Parameter.')

    @property
    def normalize_detailed_balance(self) -> bool:
        """Get whether to normalize the detailed balance factor."""

        return self._normalize_detailed_balance

    @normalize_detailed_balance.setter
    def normalize_detailed_balance(self, normalize: bool) -> None:
        """Set whether to normalize the detailed balance factor.

        If True, the detailed balance factor is divided by temperature.
        Args:
            normalize : bool
                Whether to normalize the detailed balance factor.
        Raises:
            TypeError: If normalize is not a bool.
        """

        if not isinstance(normalize, bool):
            raise TypeError('normalize_detailed_balance must be True or False.')

        self._normalize_detailed_balance = normalize

    def _create_energy_grid(
        self,
    ) -> EnergyGrid:
        """Create a dense grid by upsampling and extending the energy
        array.

        If upsample_factor is None, no upsampling or extension is
        performed.
        This dense grid is used for convolution to improve accuracy.
        Returns:
            EnergyGrid
                The dense grid created by upsampling and extending
                energy.
        The EnergyGrid has the following attributes:
           energy_dense : np.ndarray
               The upsampled and extended energy array.
              energy_dense_centered : np.ndarray
                The centered version of energy_dense
                (used for resolution evaluation).
            energy_dense_step : float
                The spacing of energy_dense
                (used for width checks and normalization).
            energy_span_dense : float
                The total span of energy_dense. (used for width checks).
            energy_even_length_offset : float
                The offset to apply if energy_dense has even length
                (used for convolution alignment).
        """
        if self.upsample_factor is None:
            # Check if the array is uniformly spaced.
            energy_diff = np.diff(self.energy.values)
            is_uniform = np.allclose(energy_diff, energy_diff[0])
            if not is_uniform:
                raise ValueError(
                    'Input array `energy` must be uniformly spaced if upsample_factor is not given.'  # noqa: E501
                )
            energy_dense = self.energy.values

            energy_span_dense = self.energy.values.max() - self.energy.values.min()
        else:
            # Create an extended and upsampled energy grid
            energy_min, energy_max = self.energy.values.min(), self.energy.values.max()
            energy_span_original = energy_max - energy_min
            extra = self.extension_factor / 2 * energy_span_original
            extended_min = energy_min - extra
            extended_max = energy_max + extra
            num_points = round(len(self.energy.values) * self.upsample_factor)
            energy_dense = np.linspace(extended_min, extended_max, num_points)
            energy_span_dense = extended_max - extended_min

        if len(energy_dense) < 2:
            raise ValueError('Energy array must have at least two points.')
        energy_dense_step = energy_dense[1] - energy_dense[0]

        # Handle offset for even length of energy_dense in convolution.
        # The convolution of two arrays of length N is of length 2N-1.
        #  When using 'same' mode, only the central N points are kept,
        # so the output has the same length as the input.
        # However, if N is even, the center falls between two points,
        # leading to a half-bin offset.
        # For example, if N=4, the convolution has length 7, and when we
        # select the 4 central points we either get
        # indices [2,3,4,5] or [1,2,3,4], both of which are offset by
        # 0.5*dx from the true center at index 3.5.
        if len(energy_dense) % 2 == 0:
            energy_even_length_offset = -0.5 * energy_dense_step
        else:
            energy_even_length_offset = 0.0

        # Handle the case when energy_dense is not symmetric around 0.
        # The resolution is still centered around zero (or close to it),
        # so it needs to be evaluated there.
        if not np.isclose(energy_dense.mean(), 0.0):
            energy_dense_centered = np.linspace(
                -0.5 * energy_span_dense, 0.5 * energy_span_dense, len(energy_dense)
            )
        else:
            energy_dense_centered = energy_dense

        energy_grid = EnergyGrid(
            energy_dense=energy_dense,
            energy_dense_centered=energy_dense_centered,
            energy_dense_step=energy_dense_step,
            energy_span_dense=energy_span_dense,
            energy_even_length_offset=energy_even_length_offset,
        )

        return energy_grid

    def _check_width_thresholds(
        self,
        model: ComponentCollection | ModelComponent,
        model_name: str,
    ) -> None:
        """Helper function to check and warn if components are wide
        compared to the span of the data, or narrow compared to the
        spacing.

        In both cases, the convolution accuracy may be compromised.
        Args:
            model : ComponentCollection or ModelComponent
                The model to check.
            model_name : str
                A string indicating whether the model is a
                'sample model' or 'resolution model' for
                warning messages.
        returns:
            None
        warns:
            UserWarning
                If the component widths are not appropriate for the data
                span or bin spacing.
        """

        # Handle ComponentCollection or ModelComponent
        if isinstance(model, ComponentCollection):
            components = model.components
        else:
            components = [model]  # Treat single ModelComponent as a list

        for comp in components:
            if hasattr(comp, 'width'):
                if comp.width.value > LARGE_WIDTH_THRESHOLD * self._energy_grid.energy_span_dense:
                    warnings.warn(
                        f"The width of the {model_name} component '{comp.unique_name}' \
                            ({comp.width.value}) is large compared to the span of the input "
                        f'array ({self._energy_grid.energy_span_dense}). \
                            This may lead to inaccuracies in the convolution. \
                                Increase extension_factor to improve accuracy.',
                        UserWarning,
                    )
                if comp.width.value < SMALL_WIDTH_THRESHOLD * self._energy_grid.energy_dense_step:
                    warnings.warn(
                        f"The width of the {model_name} component '{comp.unique_name}' \
                            ({comp.width.value}) is small compared to the spacing of the input "
                        f'array ({self._energy_grid.energy_dense_step}). \
                            This may lead to inaccuracies in the convolution. \
                                Increase upsample_factor to improve accuracy.',
                        UserWarning,
                    )

    def __repr__(self) -> str:
        return (
            f'{self.__class__.__name__}('
            f'energy=array of shape {self.energy.values.shape},\n '
            f'sample_components={repr(self.sample_components)}, \n'
            f'resolution_components={repr(self.resolution_components)},\n '
            f'energy_unit={self._energy_unit}, '
            f'upsample_factor={self.upsample_factor}, '
            f'extension_factor={self.extension_factor}, '
            f'temperature={self.temperature}, '
            f'normalize_detailed_balance={self.normalize_detailed_balance})'
        )

extension_factor property writable

Get the extension factor.

The extension factor determines how much the energy range is extended on both sides before convolution. 0.2 means extending by 20% of the original energy span on each side

normalize_detailed_balance property writable

Get whether to normalize the detailed balance factor.

temperature property writable

Get the temperature.

upsample_factor property writable

Get the upsample factor.