Model component

.gitignore

@@ -3,4 +3,4 @@ examples/INS_example/*
 examples/Anesthetics
 src/easydynamics/__pycache__
 .vscode/*
-**/__pycache__/*
+**/__pycache__/*

examples/component_example.ipynb

@@ -0,0 +1,136 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "64deaa41",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "from easydynamics.sample_model import Gaussian\n",
+    "from easydynamics.sample_model import Lorentzian\n",
+    "from easydynamics.sample_model import DampedHarmonicOscillator\n",
+    "from easydynamics.sample_model import Polynomial\n",
+    "from easydynamics.sample_model import DeltaFunction\n",
+    "\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "\n",
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "784d9e82",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Creating a component\n",
+    "gaussian=Gaussian(name='Gaussian',width=0.5,area=1)\n",
+    "dho = DampedHarmonicOscillator(name='DHO',center=1.0,width=0.3,area=2.0)\n",
+    "lorentzian = Lorentzian(name='Lorentzian',center=-1.0,width=0.2,area=1.0)\n",
+    "polynomial = Polynomial(name='Polynomial',coefficients=[0.1, 0, 0.5])  # y=0.1+0.5*x^2\n",
+    "\n",
+    "x=np.linspace(-2, 2, 100)\n",
+    "\n",
+    "plt.figure()\n",
+    "y=gaussian.evaluate(x)\n",
+    "plt.plot(x, y, label='Gaussian')\n",
+    "y=dho.evaluate(x)\n",
+    "plt.plot(x, y, label='DHO')\n",
+    "y=lorentzian.evaluate(x)\n",
+    "plt.plot(x, y, label='Lorentzian')\n",
+    "y=polynomial.evaluate(x)\n",
+    "plt.plot(x, y, label='Polynomial')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2f57228c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The area under the DHO curve is indeed equal to the area parameter.\n",
+    "xx=np.linspace(-15, 15, 10000)\n",
+    "yy=dho.evaluate(xx)\n",
+    "area= np.trapezoid(yy, xx)\n",
+    "print(f\"Area under DHO curve: {area:.4f}\")\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6c0929ed",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "delta = DeltaFunction(name='Delta', center=0.0, area=1.0)\n",
+    "x1=np.linspace(-2, 2, 100)\n",
+    "y=delta.evaluate(x1)\n",
+    "x2=np.linspace(-2,2,51)\n",
+    "y2=delta.evaluate(x2)\n",
+    "plt.figure()\n",
+    "plt.plot(x1, y, label='Delta Function')\n",
+    "plt.plot(x2, y2, label='Delta Function (coarser)')\n",
+    "plt.legend()\n",
+    "plt.show()\n",
+    "# The area under the Delta function is indeed equal to the area parameter.\n",
+    "xx=np.linspace(-2, 2, 10000)\n",
+    "yy=delta.evaluate(xx)\n",
+    "area= np.trapezoid(y, x1)\n",
+    "print(area)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f44b125a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import scipp as sc\n",
+    "x1=sc.linspace(dim='x', start=-2.0, stop=2.0, num=100, unit='meV')\n",
+    "x2=sc.linspace(dim='x', start=-2.0*1e3, stop=2.0*1e3, num=101, unit='microeV')\n",
+    "\n",
+    "polynomial = Polynomial(name='Polynomial',coefficients=[0.1, 0, 0.5])  # y=0.1+0.5*x^2\n",
+    "y1=polynomial.evaluate(x1)\n",
+    "y2=polynomial.evaluate(x2)\n",
+    "\n",
+    "plt.figure()\n",
+    "plt.plot(x1.values, y1, label='Polynomial meV',color='blue')\n",
+    "plt.plot(x2.values/1000, y2, label='Polynomial microeV',linestyle='dashed',color='orange')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "newdynamics",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

src/easydynamics/sample_model/__init__.py

@@ -0,0 +1,18 @@
+from .components import (
+    DampedHarmonicOscillator,
+    DeltaFunction,
+    Gaussian,
+    Lorentzian,
+    Polynomial,
+    Voigt,
+)
+
+__all__ = [
+    "SampleModel",
+    "Gaussian",
+    "Lorentzian",
+    "Voigt",
+    "DeltaFunction",
+    "DampedHarmonicOscillator",
+    "Polynomial",
+]

src/easydynamics/sample_model/components/__init__.py

@@ -0,0 +1,15 @@
+from .damped_harmonic_oscillator import DampedHarmonicOscillator
+from .delta_function import DeltaFunction
+from .gaussian import Gaussian
+from .lorentzian import Lorentzian
+from .polynomial import Polynomial
+from .voigt import Voigt
+
+__all__ = [
+    "Gaussian",
+    "Lorentzian",
+    "Voigt",
+    "DeltaFunction",
+    "DampedHarmonicOscillator",
+    "Polynomial",
+]

src/easydynamics/sample_model/components/damped_harmonic_oscillator.py

@@ -0,0 +1,181 @@
+from __future__ import annotations
+
+import warnings
+from typing import Union
+
+import numpy as np
+import scipp as sc
+from easyscience.variable import Parameter
+
+from .model_component import ModelComponent
+
+Numeric = Union[float, int]
+
+MINIMUM_WIDTH = 1e-10  # To avoid division by zero
+
+
+class DampedHarmonicOscillator(ModelComponent):
+    """
+    Damped Harmonic Oscillator (DHO). 2*area*center^2*width/pi / ( (x^2 - center^2)^2 + (2*width*x)^2 )
+
+    Args:
+        name (str): Name of the component.
+        center (Int or float): Resonance frequency, approximately the peak position.
+        width (Int or float): Damping constant, approximately the half width at half max (HWHM) of the peaks.
+        area (Int or float): Area under the curve.
+        unit (str or sc.Unit): Unit of the parameters. Defaults to "meV".
+    """
+
+    def __init__(
+        self,
+        name: str = "DHO",
+        center: Numeric = 1.0,
+        width: Numeric = 1.0,
+        area: Numeric = 1.0,
+        unit: Union[str, sc.Unit] = "meV",
+    ):
+        # Validate inputs
+        if not isinstance(area, Numeric):
+            raise TypeError("area must be a number.")
+        area = float(area)
+        if area < 0:
+            warnings.warn(
+                "The area of the Damped Harmonic Oscillator with name {} is negative, which may not be physically meaningful.".format(
+                    name
+                )
+            )
+
+        if not isinstance(center, Numeric):
+            raise TypeError("center must be a number.")
+
+        center = float(center)
+
+        if not isinstance(width, Numeric):
+            raise TypeError("width must be a number.")
+
+        width = float(width)
+        if width <= 0:
+            raise ValueError(
+                "The width of a DampedHarmonicOscillator must be greater than zero."
+            )
+
+        super().__init__(name=name, unit=unit)
+
+        # Create Parameters from floats
+        self._area = Parameter(name=name + " area", value=area, unit=unit)
+        if area > 0:
+            self._area.min = 0.0
+
+        self._center = Parameter(name=name + " center", value=center, unit=unit)
+
+        self._width = Parameter(
+            name=name + " width", value=width, unit=unit, min=MINIMUM_WIDTH
+        )
+
+    @property
+    def area(self) -> Parameter:
+        """Return the area parameter."""
+        return self._area
+
+    @area.setter
+    def area(self, value: Numeric):
+        """Set the area parameter."""
+        if not isinstance(value, Numeric):
+            raise TypeError("area must be a number.")
+        value = float(value)
+        if value < 0:
+            warnings.warn(
+                "The area of the Damped Harmonic Oscillator with name {} is negative, which may not be physically meaningful.".format(
+                    self.name
+                )
+            )
+        self._area.value = float(value)
+
+    @property
+    def center(self) -> Parameter:
+        """Return the center parameter."""
+        return self._center
+
+    @center.setter
+    def center(self, value: Numeric):
+        """Set the center parameter."""
+        if not isinstance(value, Numeric):
+            raise TypeError("center must be a number.")
+        self._center.value = float(value)
+
+    @property
+    def width(self) -> Parameter:
+        """Return the width parameter."""
+        return self._width
+
+    @width.setter
+    def width(self, value: Numeric):
+        """Set the width parameter."""
+        if not isinstance(value, Numeric):
+            raise TypeError("width must be a number.")
+        value = float(value)
+        if value <= 0:
+            raise ValueError(
+                "The width of a DampedHarmonicOscillator must be greater than zero."
+            )
+        self._width.value = value
+
+    def evaluate(
+        self, x: Union[Numeric, list, np.ndarray, sc.Variable, sc.DataArray]
+    ) -> np.ndarray:
+        """Evaluate the Damped Harmonic Oscillator at the given x values.
+        If x is a scipp Variable, the unit of the DHO will be converted to
+        match x. The DHO evaluates to 2*area*center^2*width/pi / ( (x^2 - center^2)^2 + (2*width*x)^2 )"""
+
+        x = self._prepare_x_for_evaluate(x)
+
+        normalization = 2 * self._center.value**2 * self._width.value / np.pi
+        denominator = (x**2 - self._center.value**2) ** 2 + (
+            2
+            * self._width.value
+            * x  # No division by zero here, width>0 enforced in setter
+        ) ** 2
+
+        return self._area.value * normalization / (denominator)
+
+    def get_parameters(self):
+        """
+        Get all parameters from the model component.
+        Returns:
+        List[Parameter]: List of parameters in the component.
+        """
+        return [self._area, self._center, self._width]
+
+    def convert_unit(self, unit: Union[str, sc.Unit]):
+        """
+        Convert the unit of the Parameters in the component.
+
+        Args:
+            unit (str or sc.Unit): The new unit to convert to.
+        """
+
+        self._area.convert_unit(unit)
+        self._center.convert_unit(unit)
+        self._width.convert_unit(unit)
+        self._unit = unit
+
+    def __copy__(self) -> DampedHarmonicOscillator:
+        """
+        Return a deep copy of this component with independent parameters.
+        """
+        name = "copy of " + self.name
+
+        model_copy = DampedHarmonicOscillator(
+            name=name,
+            area=self._area.value,
+            center=self._center.value,
+            width=self._width.value,
+            unit=self._unit,
+        )
+        model_copy._area.fixed = self._area.fixed
+        model_copy._center.fixed = self._center.fixed
+        model_copy._width.fixed = self._width.fixed
+        return model_copy
+
+    def __repr__(self):
+        return f"DampedHarmonicOscillator(name = {self.name}, unit = {self._unit},\n area = {self._area},\n center = {self._center},\n width = {self._width})"