Skip to content

Commit 95175e5

Browse files
henrikjacobsenfysAndrewSazonov
henrikjacobsenfys
and
AndrewSazonov
authored
Model component (#50)
* add coverage.yaml * fix typo * Add toml, gitignore and basic folder structure * add dummy test * add some dummy code to check if code coverage works * set up coverage * Add BSD 3-Clause License Added the BSD 3-Clause License to the project. * Update pyproject.toml for hatchling build system * move dummy code to add source folder * change dummy tests * change import * Enhance coverage report command in CI workflow * Create __init__.py * Add GitHub Actions workflow for unit testing * remove dummies, add sample model and components * Fixing type issues * update typing, remove a method that will be introduced later * More type * Add example * More type fixing * more typing... * Update detailed balancing * Update eample and detailed balance calculation * Update gitignore * Update example * Update example * Remove pycache * Remove features to make pull request smaller * try to fix coverage * Fix some minor issues * Formatting * remove dummy code * Added many tests and updated a few style things * More tests * And a few more tests * rename folder to sample_model * split compoments into different files * update type hinting * Update units, remove unused tests * Split tests into components * Update handling of incorrect units * Update warnings to handle parameter input * Remove Parameters as input :( * first step to making attributes * remove get_parameter * changed class descriptions * Update minimum width to be nonzero * Added minimum for area * Update evaluate expressions * Update copy * Remove unneeded checks in evaluate * Update docstring of evaluate * Update evaluate of Gaussian * Check NaN and inf * Make equations slightly easier to read * Add setters and getters for gaussian * Update tests * Update doc strings * Add getters and setters * Update tests of gaussian * Add WHEN THEN EXPECT to tests * Add tests of the getters and setters * Add tests to DHO * Update tests * And a bit more test * ruff * Update units to allow scipp units * add WHEN THEN EXPECT to DHO * add WHEN THEN EXPECT to DeltaFunction * add WHEN THEN EXPECT to Gaussian * Add WHEN THEN EXPECT to Lorentzian * add more WHEN THEN EXPECT to tests * fixed warnings in tests * Update polynomial to accept units * Add some explanation to polynomial component * rename a few tests * make evaluate a bit easier * Handle sc.DataArray as input * Add some tests and handle some comments * Update delta function * Add checks to prevent setting negative widths * Move unit to ModelComponent * Update polynomial, various small fixes * Update tests * parametrize * Parametrize tests of setters * Delete unneeded tests * Update test of get_parameters * Update polynomial * Add EPSILON to delta function * Update delta function and make x sorted * Update polynomial * Update copy * Final fixes --------- Co-authored-by: Andrew Sazonov <[email protected]>
1 parent d51e07e commit 95175e5

18 files changed

+2712
-1
lines changed

‎.gitignore‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,4 +3,4 @@ examples/INS_example/*
33
examples/Anesthetics
44
src/easydynamics/__pycache__
55
.vscode/*
6-
**/__pycache__/*
6+
**/__pycache__/*

‎examples/component_example.ipynb‎

Lines changed: 136 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,136 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"id": "64deaa41",
7+
"metadata": {},
8+
"outputs": [],
9+
"source": [
10+
"import numpy as np\n",
11+
"\n",
12+
"from easydynamics.sample_model import Gaussian\n",
13+
"from easydynamics.sample_model import Lorentzian\n",
14+
"from easydynamics.sample_model import DampedHarmonicOscillator\n",
15+
"from easydynamics.sample_model import Polynomial\n",
16+
"from easydynamics.sample_model import DeltaFunction\n",
17+
"\n",
18+
"\n",
19+
"import matplotlib.pyplot as plt\n",
20+
"\n",
21+
"\n",
22+
"%matplotlib widget"
23+
]
24+
},
25+
{
26+
"cell_type": "code",
27+
"execution_count": null,
28+
"id": "784d9e82",
29+
"metadata": {},
30+
"outputs": [],
31+
"source": [
32+
"# Creating a component\n",
33+
"gaussian=Gaussian(name='Gaussian',width=0.5,area=1)\n",
34+
"dho = DampedHarmonicOscillator(name='DHO',center=1.0,width=0.3,area=2.0)\n",
35+
"lorentzian = Lorentzian(name='Lorentzian',center=-1.0,width=0.2,area=1.0)\n",
36+
"polynomial = Polynomial(name='Polynomial',coefficients=[0.1, 0, 0.5]) # y=0.1+0.5*x^2\n",
37+
"\n",
38+
"x=np.linspace(-2, 2, 100)\n",
39+
"\n",
40+
"plt.figure()\n",
41+
"y=gaussian.evaluate(x)\n",
42+
"plt.plot(x, y, label='Gaussian')\n",
43+
"y=dho.evaluate(x)\n",
44+
"plt.plot(x, y, label='DHO')\n",
45+
"y=lorentzian.evaluate(x)\n",
46+
"plt.plot(x, y, label='Lorentzian')\n",
47+
"y=polynomial.evaluate(x)\n",
48+
"plt.plot(x, y, label='Polynomial')\n",
49+
"plt.legend()\n",
50+
"plt.show()"
51+
]
52+
},
53+
{
54+
"cell_type": "code",
55+
"execution_count": null,
56+
"id": "2f57228c",
57+
"metadata": {},
58+
"outputs": [],
59+
"source": [
60+
"# The area under the DHO curve is indeed equal to the area parameter.\n",
61+
"xx=np.linspace(-15, 15, 10000)\n",
62+
"yy=dho.evaluate(xx)\n",
63+
"area= np.trapezoid(yy, xx)\n",
64+
"print(f\"Area under DHO curve: {area:.4f}\")\n",
65+
"\n"
66+
]
67+
},
68+
{
69+
"cell_type": "code",
70+
"execution_count": null,
71+
"id": "6c0929ed",
72+
"metadata": {},
73+
"outputs": [],
74+
"source": [
75+
"delta = DeltaFunction(name='Delta', center=0.0, area=1.0)\n",
76+
"x1=np.linspace(-2, 2, 100)\n",
77+
"y=delta.evaluate(x1)\n",
78+
"x2=np.linspace(-2,2,51)\n",
79+
"y2=delta.evaluate(x2)\n",
80+
"plt.figure()\n",
81+
"plt.plot(x1, y, label='Delta Function')\n",
82+
"plt.plot(x2, y2, label='Delta Function (coarser)')\n",
83+
"plt.legend()\n",
84+
"plt.show()\n",
85+
"# The area under the Delta function is indeed equal to the area parameter.\n",
86+
"xx=np.linspace(-2, 2, 10000)\n",
87+
"yy=delta.evaluate(xx)\n",
88+
"area= np.trapezoid(y, x1)\n",
89+
"print(area)\n"
90+
]
91+
},
92+
{
93+
"cell_type": "code",
94+
"execution_count": null,
95+
"id": "f44b125a",
96+
"metadata": {},
97+
"outputs": [],
98+
"source": [
99+
"import scipp as sc\n",
100+
"x1=sc.linspace(dim='x', start=-2.0, stop=2.0, num=100, unit='meV')\n",
101+
"x2=sc.linspace(dim='x', start=-2.0*1e3, stop=2.0*1e3, num=101, unit='microeV')\n",
102+
"\n",
103+
"polynomial = Polynomial(name='Polynomial',coefficients=[0.1, 0, 0.5]) # y=0.1+0.5*x^2\n",
104+
"y1=polynomial.evaluate(x1)\n",
105+
"y2=polynomial.evaluate(x2)\n",
106+
"\n",
107+
"plt.figure()\n",
108+
"plt.plot(x1.values, y1, label='Polynomial meV',color='blue')\n",
109+
"plt.plot(x2.values/1000, y2, label='Polynomial microeV',linestyle='dashed',color='orange')\n",
110+
"plt.legend()\n",
111+
"plt.show()"
112+
]
113+
}
114+
],
115+
"metadata": {
116+
"kernelspec": {
117+
"display_name": "newdynamics",
118+
"language": "python",
119+
"name": "python3"
120+
},
121+
"language_info": {
122+
"codemirror_mode": {
123+
"name": "ipython",
124+
"version": 3
125+
},
126+
"file_extension": ".py",
127+
"mimetype": "text/x-python",
128+
"name": "python",
129+
"nbconvert_exporter": "python",
130+
"pygments_lexer": "ipython3",
131+
"version": "3.11.13"
132+
}
133+
},
134+
"nbformat": 4,
135+
"nbformat_minor": 5
136+
}

‎src/easydynamics/sample_model/__init__.py‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
from .components import (
2+
DampedHarmonicOscillator,
3+
DeltaFunction,
4+
Gaussian,
5+
Lorentzian,
6+
Polynomial,
7+
Voigt,
8+
)
9+
10+
__all__ = [
11+
"SampleModel",
12+
"Gaussian",
13+
"Lorentzian",
14+
"Voigt",
15+
"DeltaFunction",
16+
"DampedHarmonicOscillator",
17+
"Polynomial",
18+
]

‎src/easydynamics/sample_model/components/__init__.py‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
from .damped_harmonic_oscillator import DampedHarmonicOscillator
2+
from .delta_function import DeltaFunction
3+
from .gaussian import Gaussian
4+
from .lorentzian import Lorentzian
5+
from .polynomial import Polynomial
6+
from .voigt import Voigt
7+
8+
__all__ = [
9+
"Gaussian",
10+
"Lorentzian",
11+
"Voigt",
12+
"DeltaFunction",
13+
"DampedHarmonicOscillator",
14+
"Polynomial",
15+
]

‎src/easydynamics/sample_model/components/damped_harmonic_oscillator.py‎

Lines changed: 181 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,181 @@
1+
from __future__ import annotations
2+
3+
import warnings
4+
from typing import Union
5+
6+
import numpy as np
7+
import scipp as sc
8+
from easyscience.variable import Parameter
9+
10+
from .model_component import ModelComponent
11+
12+
Numeric = Union[float, int]
13+
14+
MINIMUM_WIDTH = 1e-10 # To avoid division by zero
15+
16+
17+
class DampedHarmonicOscillator(ModelComponent):
18+
"""
19+
Damped Harmonic Oscillator (DHO). 2*area*center^2*width/pi / ( (x^2 - center^2)^2 + (2*width*x)^2 )
20+
21+
Args:
22+
name (str): Name of the component.
23+
center (Int or float): Resonance frequency, approximately the peak position.
24+
width (Int or float): Damping constant, approximately the half width at half max (HWHM) of the peaks.
25+
area (Int or float): Area under the curve.
26+
unit (str or sc.Unit): Unit of the parameters. Defaults to "meV".
27+
"""
28+
29+
def __init__(
30+
self,
31+
name: str = "DHO",
32+
center: Numeric = 1.0,
33+
width: Numeric = 1.0,
34+
area: Numeric = 1.0,
35+
unit: Union[str, sc.Unit] = "meV",
36+
):
37+
# Validate inputs
38+
if not isinstance(area, Numeric):
39+
raise TypeError("area must be a number.")
40+
area = float(area)
41+
if area < 0:
42+
warnings.warn(
43+
"The area of the Damped Harmonic Oscillator with name {} is negative, which may not be physically meaningful.".format(
44+
name
45+
)
46+
)
47+
48+
if not isinstance(center, Numeric):
49+
raise TypeError("center must be a number.")
50+
51+
center = float(center)
52+
53+
if not isinstance(width, Numeric):
54+
raise TypeError("width must be a number.")
55+
56+
width = float(width)
57+
if width <= 0:
58+
raise ValueError(
59+
"The width of a DampedHarmonicOscillator must be greater than zero."
60+
)
61+
62+
super().__init__(name=name, unit=unit)
63+
64+
# Create Parameters from floats
65+
self._area = Parameter(name=name + " area", value=area, unit=unit)
66+
if area > 0:
67+
self._area.min = 0.0
68+
69+
self._center = Parameter(name=name + " center", value=center, unit=unit)
70+
71+
self._width = Parameter(
72+
name=name + " width", value=width, unit=unit, min=MINIMUM_WIDTH
73+
)
74+
75+
@property
76+
def area(self) -> Parameter:
77+
"""Return the area parameter."""
78+
return self._area
79+
80+
@area.setter
81+
def area(self, value: Numeric):
82+
"""Set the area parameter."""
83+
if not isinstance(value, Numeric):
84+
raise TypeError("area must be a number.")
85+
value = float(value)
86+
if value < 0:
87+
warnings.warn(
88+
"The area of the Damped Harmonic Oscillator with name {} is negative, which may not be physically meaningful.".format(
89+
self.name
90+
)
91+
)
92+
self._area.value = float(value)
93+
94+
@property
95+
def center(self) -> Parameter:
96+
"""Return the center parameter."""
97+
return self._center
98+
99+
@center.setter
100+
def center(self, value: Numeric):
101+
"""Set the center parameter."""
102+
if not isinstance(value, Numeric):
103+
raise TypeError("center must be a number.")
104+
self._center.value = float(value)
105+
106+
@property
107+
def width(self) -> Parameter:
108+
"""Return the width parameter."""
109+
return self._width
110+
111+
@width.setter
112+
def width(self, value: Numeric):
113+
"""Set the width parameter."""
114+
if not isinstance(value, Numeric):
115+
raise TypeError("width must be a number.")
116+
value = float(value)
117+
if value <= 0:
118+
raise ValueError(
119+
"The width of a DampedHarmonicOscillator must be greater than zero."
120+
)
121+
self._width.value = value
122+
123+
def evaluate(
124+
self, x: Union[Numeric, list, np.ndarray, sc.Variable, sc.DataArray]
125+
) -> np.ndarray:
126+
"""Evaluate the Damped Harmonic Oscillator at the given x values.
127+
If x is a scipp Variable, the unit of the DHO will be converted to
128+
match x. The DHO evaluates to 2*area*center^2*width/pi / ( (x^2 - center^2)^2 + (2*width*x)^2 )"""
129+
130+
x = self._prepare_x_for_evaluate(x)
131+
132+
normalization = 2 * self._center.value**2 * self._width.value / np.pi
133+
denominator = (x**2 - self._center.value**2) ** 2 + (
134+
2
135+
* self._width.value
136+
* x # No division by zero here, width>0 enforced in setter
137+
) ** 2
138+
139+
return self._area.value * normalization / (denominator)
140+
141+
def get_parameters(self):
142+
"""
143+
Get all parameters from the model component.
144+
Returns:
145+
List[Parameter]: List of parameters in the component.
146+
"""
147+
return [self._area, self._center, self._width]
148+
149+
def convert_unit(self, unit: Union[str, sc.Unit]):
150+
"""
151+
Convert the unit of the Parameters in the component.
152+
153+
Args:
154+
unit (str or sc.Unit): The new unit to convert to.
155+
"""
156+
157+
self._area.convert_unit(unit)
158+
self._center.convert_unit(unit)
159+
self._width.convert_unit(unit)
160+
self._unit = unit
161+
162+
def __copy__(self) -> DampedHarmonicOscillator:
163+
"""
164+
Return a deep copy of this component with independent parameters.
165+
"""
166+
name = "copy of " + self.name
167+
168+
model_copy = DampedHarmonicOscillator(
169+
name=name,
170+
area=self._area.value,
171+
center=self._center.value,
172+
width=self._width.value,
173+
unit=self._unit,
174+
)
175+
model_copy._area.fixed = self._area.fixed
176+
model_copy._center.fixed = self._center.fixed
177+
model_copy._width.fixed = self._width.fixed
178+
return model_copy
179+
180+
def __repr__(self):
181+
return f"DampedHarmonicOscillator(name = {self.name}, unit = {self._unit},\n area = {self._area},\n center = {self._center},\n width = {self._width})"

0 commit comments

Comments
 (0)