diff --git a/docs/docs/tutorials/components.ipynb b/docs/docs/tutorials/components.ipynb index 0fa5a8884..1b4af94d6 100644 --- a/docs/docs/tutorials/components.ipynb +++ b/docs/docs/tutorials/components.ipynb @@ -7,7 +7,7 @@ "source": [ "# Components\n", "\n", - "Components are the basic ingredients for all models. Currently, the available components are Gaussian, Lorentzian, Voigt (the convolution of a Gaussian with a Lorentzian), delta function, damped harmonic oscillator and polynomial. This notebooks shows how to use the components. \n", + "Components are the basic ingredients for all models. Currently, the available components are Gaussian, Lorentzian, Voigt (the convolution of a Gaussian with a Lorentzian), delta function, damped harmonic oscillator, stretched exponential and polynomial. This notebooks shows how to use the components. \n", "\n", "Note in particular that a Gaussian, Lorentzian, Voigt or delta function where the center has not been given will be centered at 0." ] @@ -203,11 +203,36 @@ "plt.legend()\n", "plt.show()" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e09b85d0", + "metadata": {}, + "outputs": [], + "source": [ + "stretched_lorz = edyn.StretchedExponential(area=1.0, relaxation_time=10.0, beta=1.0)\n", + "lorentzian = edyn.Lorentzian(center=0.0, width=0.6582 / 10, area=1.0)\n", + "stretched = edyn.StretchedExponential(area=1.0, relaxation_time=10.0, beta=0.5)\n", + "\n", + "x = np.linspace(-1, 1, 1001)\n", + "y = stretched_lorz.evaluate(x)\n", + "y2 = lorentzian.evaluate(x)\n", + "y3 = stretched.evaluate(x)\n", + "plt.figure()\n", + "plt.plot(x, y, label='Stretched exponential (β=1)')\n", + "plt.plot(x, y2, label='Lorentzian', linestyle='dotted')\n", + "plt.plot(x, y3, label='Stretched exponential (β=0.5)')\n", + "plt.xlabel('Energy (meV)')\n", + "plt.ylabel('Intensity (arb. units)')\n", + "plt.legend()\n", + "plt.show()" + ] } ], "metadata": { "kernelspec": { - "display_name": "default", + "display_name": "DEFAULT", "language": "python", "name": "python3" }, diff --git a/docs/docs/user-guide/concept.md b/docs/docs/user-guide/concept.md index 17014b924..2aad691a0 100644 --- a/docs/docs/user-guide/concept.md +++ b/docs/docs/user-guide/concept.md @@ -19,8 +19,8 @@ a list of `ComponentCollection`s (one for each Q) describing the model of the measured data. A `ComponentCollection` is essentially a list of `ModelComponents`. A `ModelComponent` can be any of `Gaussian`, `Lorentzian`, `Voigt` (the convolution of a `Gaussian` and -`Lorentzian`), `DeltaFunction`, `DampedHarmonicOscillator` and -`Polynomium`. +`Lorentzian`), `DeltaFunction`, `DampedHarmonicOscillator`, +`StretchedExponential` and `Polynomium`. Each `ModelComponent` has a number of `Parameter`s. The `Gaussian`, for example, has `area`, `center` and `width`. Each of these `Parameter`s diff --git a/pixi.lock b/pixi.lock index d844f8856..0bf41847a 100644 --- a/pixi.lock +++ b/pixi.lock @@ -229,6 +229,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/a4/81502f486f01db95bc8320646a8a12511f5e556cb63d5e224d91816605c4/trove_classifiers-2026.6.1.19-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/c4/bc41eb19b0fd0db868f4132920879019318d80cc522ad8f2bca4611af808/scipy-1.18.0-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/80/44/e002bad11c7c9dc293141395bb2652f2c45a3dcac737c8385a1088ccbafe/format_docstring-0.4.0-py3-none-any.whl @@ -280,7 +281,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/c7/a0/5ff05d1919ca249508012cad89f08fdc6cfbdaa15b41651c5fe6dffaf1d3/dfo_ls-1.6.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d1/d6/3965ed04c63042e047cb6a3e6ed1a63a35087b6a609aa3a15ed8ac56c221/colorama-0.4.6-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl @@ -516,6 +516,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/a4/81502f486f01db95bc8320646a8a12511f5e556cb63d5e224d91816605c4/trove_classifiers-2026.6.1.19-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/80/44/e002bad11c7c9dc293141395bb2652f2c45a3dcac737c8385a1088ccbafe/format_docstring-0.4.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl @@ -564,7 +565,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/c7/da/32c752228ae345f489e3a42499d817b6c3996da7e8a3bc7a04fc806b243b/pillow-12.3.0-cp314-cp314-macosx_11_0_arm64.whl - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d1/d6/3965ed04c63042e047cb6a3e6ed1a63a35087b6a609aa3a15ed8ac56c221/colorama-0.4.6-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl @@ -782,6 +782,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7d/c2/57f54b03d0f22d4044b8afb9ca0e184f8b1afd57b4f735c2fa70883dc601/contourpy-1.3.3-cp314-cp314-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/80/44/e002bad11c7c9dc293141395bb2652f2c45a3dcac737c8385a1088ccbafe/format_docstring-0.4.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl @@ -831,7 +832,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ce/04/d719a0a36930ecc8dfc801ff340f9dcfc4223f8ca5d39d06b4020032fff8/matplotlib-3.11.1-cp314-cp314-win_amd64.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cf/52/6daa2ee9d95e5c98b8128f8df91eb692eb423ab274b8cf08db52152fad26/yarl-1.24.5-cp314-cp314-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl @@ -1075,6 +1075,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7d/a2/c4d99299e9ce7fad561f8bb56babbbbdd3bb6b4fbd7c0ec674c1dbdd2cc5/chardet-7.6.0-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/80/44/e002bad11c7c9dc293141395bb2652f2c45a3dcac737c8385a1088ccbafe/format_docstring-0.4.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl @@ -1126,7 +1127,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cc/8f/ec6289987824b29529d0dfda0d74a07cec60e54b9c92f3c9da4c0ac732de/contourpy-1.3.3-cp312-cp312-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d1/d6/3965ed04c63042e047cb6a3e6ed1a63a35087b6a609aa3a15ed8ac56c221/colorama-0.4.6-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl @@ -1360,6 +1360,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/a4/81502f486f01db95bc8320646a8a12511f5e556cb63d5e224d91816605c4/trove_classifiers-2026.6.1.19-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/80/44/e002bad11c7c9dc293141395bb2652f2c45a3dcac737c8385a1088ccbafe/format_docstring-0.4.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl @@ -1406,7 +1407,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/c7/a0/5ff05d1919ca249508012cad89f08fdc6cfbdaa15b41651c5fe6dffaf1d3/dfo_ls-1.6.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d1/d6/3965ed04c63042e047cb6a3e6ed1a63a35087b6a609aa3a15ed8ac56c221/colorama-0.4.6-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl @@ -1629,6 +1629,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/a4/81502f486f01db95bc8320646a8a12511f5e556cb63d5e224d91816605c4/trove_classifiers-2026.6.1.19-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/b9/87fea2769fe1c47c1b5b01d8310772c9d1a85d485de7cf386ef7a3332b02/numpy-2.5.2-cp312-cp312-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/80/31/0b2517913687895f5904325c2069d6a3b78f66cc641a86a2baf75a05dcbb/multidict-6.7.1-cp312-cp312-win_amd64.whl @@ -1680,7 +1681,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/c7/a0/5ff05d1919ca249508012cad89f08fdc6cfbdaa15b41651c5fe6dffaf1d3/dfo_ls-1.6.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d6/54/da572c98c0b77626a91b5d3b89f0231d8bff5125c225420908632f8b342d/pymdown_extensions-11.0.1-py3-none-any.whl @@ -1914,6 +1914,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/a4/81502f486f01db95bc8320646a8a12511f5e556cb63d5e224d91816605c4/trove_classifiers-2026.6.1.19-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/c4/bc41eb19b0fd0db868f4132920879019318d80cc522ad8f2bca4611af808/scipy-1.18.0-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/80/44/e002bad11c7c9dc293141395bb2652f2c45a3dcac737c8385a1088ccbafe/format_docstring-0.4.0-py3-none-any.whl @@ -1965,7 +1966,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/c7/a0/5ff05d1919ca249508012cad89f08fdc6cfbdaa15b41651c5fe6dffaf1d3/dfo_ls-1.6.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d1/d6/3965ed04c63042e047cb6a3e6ed1a63a35087b6a609aa3a15ed8ac56c221/colorama-0.4.6-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl @@ -2201,6 +2201,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/a4/81502f486f01db95bc8320646a8a12511f5e556cb63d5e224d91816605c4/trove_classifiers-2026.6.1.19-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/80/44/e002bad11c7c9dc293141395bb2652f2c45a3dcac737c8385a1088ccbafe/format_docstring-0.4.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl @@ -2249,7 +2250,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/c7/da/32c752228ae345f489e3a42499d817b6c3996da7e8a3bc7a04fc806b243b/pillow-12.3.0-cp314-cp314-macosx_11_0_arm64.whl - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d1/d6/3965ed04c63042e047cb6a3e6ed1a63a35087b6a609aa3a15ed8ac56c221/colorama-0.4.6-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl @@ -2467,6 +2467,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/7c/c6/76ee9dacedcd8c67d8fa53dd975613733bdd28242a4c41518ff1c8aeaa64/jupytext-1.19.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7d/c2/57f54b03d0f22d4044b8afb9ca0e184f8b1afd57b4f735c2fa70883dc601/contourpy-1.3.3-cp314-cp314-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/bf/3adcb9b3091b36de729dad91c107179c8c7c51adb2b08c31177bb540bef1/copier-9.17.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/80/44/e002bad11c7c9dc293141395bb2652f2c45a3dcac737c8385a1088ccbafe/format_docstring-0.4.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl @@ -2516,7 +2517,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/ca/31/d4e37e9e550c2b92a9cbc2e4d0b7420a27224968580b5a447f420847c975/pytest_xdist-3.8.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ce/04/d719a0a36930ecc8dfc801ff340f9dcfc4223f8ca5d39d06b4020032fff8/matplotlib-3.11.1-cp314-cp314-win_amd64.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cf/52/6daa2ee9d95e5c98b8128f8df91eb692eb423ab274b8cf08db52152fad26/yarl-1.24.5-cp314-cp314-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/d2/f0/834e479e47e499b6478e807fb57b31cc2db696c4db30557bb6f5aea4a90b/mando-0.7.1-py2.py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/08/c2409cb01d5368dcfedcbaffa7d044cc8957d57a9d0855244a5eb4709d30/funcy-2.0-py2.py3-none-any.whl @@ -2719,6 +2719,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/6b/be/92dd42844fe8a78c2c4a87f8078b9263dcc20aabe86b8420302a6fabaf4a/scipp-26.8.0-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/71/43/1947f06babed6b3f1d7f38b0c767f52df66bfb2bc10b468c4a7de9eceff2/aiohappyeyeballs-2.7.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7f/c4/bc41eb19b0fd0db868f4132920879019318d80cc522ad8f2bca4611af808/scipy-1.18.0-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/88/39/799be3f2f0f38cc727ee3b4f1445fe6d5e4133064ec2e4115069418a5bb6/cloudpickle-3.1.2-py3-none-any.whl @@ -2734,7 +2735,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/ab/b5/36c712098e6191d1b4e349304ef73a8d06aed77e56ceaac8c0a306c7bda1/jupyterlab_widgets-3.0.16-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/c7/99/461bd36dbdfac6c1c53efa370bd55a83227542d0d118f1677dbf1a3dacd5/numpy-2.5.2-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/c7/a0/5ff05d1919ca249508012cad89f08fdc6cfbdaa15b41651c5fe6dffaf1d3/dfo_ls-1.6.5-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d5/b7/1da684a04175473fa4cddbf9a2f572e79514c3fd27a74597f43057d4f3da/aiohttp-3.14.3-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/d9/be/7e6bf4088d003432e9a511656b90e3ec2abf3ff54a6057d2fa6e8ecfcbf1/easydynamics-0.9.2-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/e6/90/90a65e6ae1b6e66183b48874d32509fc306c994beecc7924a6fa3d9f8955/easyscience-2.5.1-py3-none-any.whl @@ -2913,6 +2913,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/6b/01/f804208061b504894546fddc479f6075e0f00dfe88cb1703dafaa8c3c67e/python_engineio-4.13.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/71/43/1947f06babed6b3f1d7f38b0c767f52df66bfb2bc10b468c4a7de9eceff2/aiohappyeyeballs-2.7.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/85/ed/0357a015892fd68058bf2d39d3fd1958e459b997a7db30aaa6aaa434ae96/aiohttp-3.14.3-cp314-cp314-macosx_11_0_arm64.whl - pypi: https://files.pythonhosted.org/packages/88/39/799be3f2f0f38cc727ee3b4f1445fe6d5e4133064ec2e4115069418a5bb6/cloudpickle-3.1.2-py3-none-any.whl @@ -2928,7 +2929,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/b9/ee/d08226fc858044355983a6e5b94f08ff6f3969e0a2b160a4a89f0ddb3445/numpy-2.5.2-cp314-cp314-macosx_14_0_arm64.whl - pypi: https://files.pythonhosted.org/packages/c7/a0/5ff05d1919ca249508012cad89f08fdc6cfbdaa15b41651c5fe6dffaf1d3/dfo_ls-1.6.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/c7/da/32c752228ae345f489e3a42499d817b6c3996da7e8a3bc7a04fc806b243b/pillow-12.3.0-cp314-cp314-macosx_11_0_arm64.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/d9/be/7e6bf4088d003432e9a511656b90e3ec2abf3ff54a6057d2fa6e8ecfcbf1/easydynamics-0.9.2-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/e6/90/90a65e6ae1b6e66183b48874d32509fc306c994beecc7924a6fa3d9f8955/easyscience-2.5.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/e7/05/c19819d5e3d95294a6f5947fb9b9629efb316b96de511b418c53d245aae6/cycler-0.12.1-py3-none-any.whl @@ -3099,6 +3099,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/71/43/1947f06babed6b3f1d7f38b0c767f52df66bfb2bc10b468c4a7de9eceff2/aiohappyeyeballs-2.7.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/7d/c2/57f54b03d0f22d4044b8afb9ca0e184f8b1afd57b4f735c2fa70883dc601/contourpy-1.3.3-cp314-cp314-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/7e/85/a5bfaebfd305ac18b57b0854d74e37e586809061a91fda62f0bd50c8518e/narwhals-2.24.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/82/f5/2f77f0bc663c13371d1c00ab8e550e2c9b11fec3c63ebf12a8336c0f534e/bumps-1.0.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/88/39/799be3f2f0f38cc727ee3b4f1445fe6d5e4133064ec2e4115069418a5bb6/cloudpickle-3.1.2-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/8a/a1/8d812e53a5da1687abb10445275d41a8b13adb781bbf7196ddbcf8d88505/lazy_loader-0.5-py3-none-any.whl @@ -3115,7 +3116,6 @@ environments: - pypi: https://files.pythonhosted.org/packages/c3/d4/98078064ccc76b45cb0f6c002452011e93c4bd26f6850344f0951cc1fe89/fonttools-4.63.0-cp314-cp314-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/c7/a0/5ff05d1919ca249508012cad89f08fdc6cfbdaa15b41651c5fe6dffaf1d3/dfo_ls-1.6.5-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ce/04/d719a0a36930ecc8dfc801ff340f9dcfc4223f8ca5d39d06b4020032fff8/matplotlib-3.11.1-cp314-cp314-win_amd64.whl - - pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/cf/52/6daa2ee9d95e5c98b8128f8df91eb692eb423ab274b8cf08db52152fad26/yarl-1.24.5-cp314-cp314-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/d9/be/7e6bf4088d003432e9a511656b90e3ec2abf3ff54a6057d2fa6e8ecfcbf1/easydynamics-0.9.2-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/e0/bf/52f25716bbe93745595800f36fb17b73711f14da59ed0bb2eba141bc9f0f/multidict-6.7.1-cp314-cp314-win_amd64.whl @@ -7546,7 +7546,7 @@ packages: - matplotlib - numpy - pixi-kernel - - plopp + - plopp>=26.9.0 - pooch - scipp - scipy @@ -9352,6 +9352,36 @@ packages: - sqlparse>=0.5.5 ; extra == 'sql' - sqlframe>=3.22.0,!=3.39.3 ; extra == 'sqlframe' requires_python: '>=3.10' +- pypi: https://files.pythonhosted.org/packages/7f/7e/8139f0faa2ff7aee2759300061e4b7e5e661775bcd5bf70810d0793af48a/plopp-26.9.0-py3-none-any.whl + name: plopp + version: 26.9.0 + sha256: 6fdff3744b81dce65a30092491bb219e69cbc573ad65c8ace0bb0714d86a97d2 + requires_dist: + - lazy-loader>=0.4 + - matplotlib>=3.8 + - scipp>=25.11.0 + - anywidget>=0.9.0 ; extra == 'all' + - ipympl>0.8.4 ; extra == 'all' + - pythreejs>=2.4.1 ; extra == 'all' + - mpltoolbox>=24.6.0 ; extra == 'all' + - ipywidgets>=8.1.5 ; extra == 'all' + - graphviz>=0.20.3 ; extra == 'all' + - graphviz>=0.20.3 ; extra == 'test' + - h5py>=3.12 ; extra == 'test' + - ipympl>=0.8.4 ; extra == 'test' + - ipywidgets>=8.1.5 ; extra == 'test' + - ipykernel>=6.26,<7 ; extra == 'test' + - mpltoolbox>=24.6.0 ; extra == 'test' + - pandas>=2.2.2 ; extra == 'test' + - plotly>=5.15.0 ; extra == 'test' + - pooch>=1.5 ; extra == 'test' + - pyarrow>=13.0.0 ; extra == 'test' + - pytest>=8.0 ; extra == 'test' + - pythreejs>=2.4.1 ; extra == 'test' + - scipy>=1.10.0 ; extra == 'test' + - xarray>=2024.5.0 ; extra == 'test' + - anywidget>=0.9.0 ; extra == 'test' + requires_python: '>=3.11' - pypi: https://files.pythonhosted.org/packages/7f/b9/87fea2769fe1c47c1b5b01d8310772c9d1a85d485de7cf386ef7a3332b02/numpy-2.5.2-cp312-cp312-win_amd64.whl name: numpy version: 2.5.2 @@ -10469,36 +10499,6 @@ packages: - pyparsing>=3 - python-dateutil>=2.7 requires_python: '>=3.11' -- pypi: https://files.pythonhosted.org/packages/ce/27/bdaaf32952f052c15c048fab82d971d30f92b63d54b61486f04a241fa994/plopp-26.7.0-py3-none-any.whl - name: plopp - version: 26.7.0 - sha256: 2082beccfc0df47750600ecdd68f5379bfee004a929bfcb0420e2b180b565fb2 - requires_dist: - - lazy-loader>=0.4 - - matplotlib>=3.8 - - scipp>=25.11.0 - - anywidget>=0.9.0 ; extra == 'all' - - ipympl>0.8.4 ; extra == 'all' - - pythreejs>=2.4.1 ; extra == 'all' - - mpltoolbox>=24.6.0 ; extra == 'all' - - ipywidgets>=8.1.5 ; extra == 'all' - - graphviz>=0.20.3 ; extra == 'all' - - graphviz>=0.20.3 ; extra == 'test' - - h5py>=3.12 ; extra == 'test' - - ipympl>=0.8.4 ; extra == 'test' - - ipywidgets>=8.1.5 ; extra == 'test' - - ipykernel>=6.26,<7 ; extra == 'test' - - mpltoolbox>=24.6.0 ; extra == 'test' - - pandas>=2.2.2 ; extra == 'test' - - plotly>=5.15.0 ; extra == 'test' - - pooch>=1.5 ; extra == 'test' - - pyarrow>=13.0.0 ; extra == 'test' - - pytest>=8.0 ; extra == 'test' - - pythreejs>=2.4.1 ; extra == 'test' - - scipy>=1.10.0 ; extra == 'test' - - xarray>=2024.5.0 ; extra == 'test' - - anywidget>=0.9.0 ; extra == 'test' - requires_python: '>=3.11' - pypi: https://files.pythonhosted.org/packages/cf/52/6daa2ee9d95e5c98b8128f8df91eb692eb423ab274b8cf08db52152fad26/yarl-1.24.5-cp314-cp314-win_amd64.whl name: yarl version: 1.24.5 diff --git a/pixi.toml b/pixi.toml index 007abbce4..406b2019b 100644 --- a/pixi.toml +++ b/pixi.toml @@ -3,14 +3,13 @@ ########### [workspace] - # Supported platforms for the lock file (pixi.lock) platforms = [ 'win-64', # Set minimum supported version for glibc to be 2.35 to ensure packages # like `crysfml` that only have wheels for glibc 2.35+ # (manylinux_2_35_x86_64) are used. - #libc = { family = 'glibc', version = '2.35' } + # libc = { family = 'glibc', version = '2.35' } { platform = 'linux-64', glibc = '2.35' }, # Set minimum supported version for macOS to be 14.0 to ensure packages # like `scipp` that only have wheels for macOS 14.0+ (macosx_14_0_arm64) @@ -45,10 +44,10 @@ python = '3.14.*' # editable installations. [feature.dev.dependencies] -nodejs = '*' # Required for Prettier (non-Python formatting) +nodejs = '*' # Required for Prettier (non-Python formatting) jupyterlab = '*' # Jupyter notebooks -ipython = '*' # Interactive Python shell -pixi-kernel = '*' # Pixi Jupyter kernel +ipython = '*' # Interactive Python shell +pixi-kernel = '*' # Pixi Jupyter kernel [feature.dev.pypi-dependencies] pip = '*' @@ -59,8 +58,8 @@ easydynamics = { path = '.', editable = true, extras = ['dev'] } [feature.user.dependencies] jupyterlab = '*' # Jupyter notebooks -ipython = '*' # Interactive Python shell -pixi-kernel = '*' # Pixi Jupyter kernel +ipython = '*' # Interactive Python shell +pixi-kernel = '*' # Pixi Jupyter kernel [feature.user.pypi-dependencies] pip = '*' @@ -71,7 +70,6 @@ easydynamics = '*' ############## [environments] - # The `default` feature is always included in all environments. # Additional features can be specified per environment. @@ -92,7 +90,6 @@ user = { features = ['py-max', 'user'] } ####### [tasks] - ################## # 🧪 Testing Tasks ################## @@ -190,9 +187,9 @@ notebook-exec = { cmd = 'python -m pytest --nbmake docs/docs/tutorials/ --nbmake ] } notebook-prepare = { depends-on = [ - #'notebook-convert', + # 'notebook-convert', 'notebook-strip', - #'notebook-tweak', + # 'notebook-tweak', ] } ######################## @@ -287,7 +284,7 @@ clean-pycache = "find . -type d -name '__pycache__' -prune -exec rm -rf '{}' +" post-install = { depends-on = [ 'npm-config', 'prettier-install', - #'pre-commit-setup', + # 'pre-commit-setup', ] } ########################## diff --git a/prettierrc.toml b/prettierrc.toml index b98c86eb7..d8a956f5e 100644 --- a/prettierrc.toml +++ b/prettierrc.toml @@ -1,22 +1,24 @@ plugins = [ - "prettier-plugin-toml", # use the TOML plugin + 'prettier-plugin-toml', # use the TOML plugin ] -endOfLine = 'lf' # change line endings to LF -proseWrap = 'always' # change wrapping in Markdown files -semi = false # remove semicolons -singleQuote = true # use single quotes instead of double quotes -tabWidth = 2 # change tab width to 2 spaces -useTabs = false # use spaces instead of tabs +endOfLine = 'lf' # change line endings to LF +proseWrap = 'always' # change wrapping in Markdown files +semi = false # remove semicolons +singleQuote = true # use single quotes instead of double quotes +tabWidth = 2 # change tab width to 2 spaces +useTabs = false # use spaces instead of tabs -printWidth = 79 # wrap lines at 79 characters +printWidth = 79 # wrap lines at 79 characters [[overrides]] -files = ["*.md"] +files = ['*.md'] + [overrides.options] -printWidth = 72 # wrap Markdown files at 72 characters +printWidth = 72 # wrap Markdown files at 72 characters [[overrides]] -files = ["*.yml", "*.yaml"] +files = ['*.yml', '*.yaml'] + [overrides.options] -printWidth = 88 # wrap YAML files at 88 characters +printWidth = 88 # wrap YAML files at 88 characters diff --git a/pyproject.toml b/pyproject.toml index c8ff0397c..7f1f68d8b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ [project] name = 'easydynamics' -dynamic = ['version'] # Use versioningit to manage the version +dynamic = ['version'] # Use versioningit to manage the version description = 'QENS data analysis' authors = [{ name = 'EasyScience contributors' }] readme = 'README.md' @@ -23,58 +23,58 @@ classifiers = [ ] requires-python = '>=3.12' dependencies = [ - 'easyscience>=2.5.1', # The base library of the EasyScience framework. 2.5.1 adds fitting.Sampler - 'numpy', # Numerical arrays (used directly throughout the library) - 'scipy', # Numerical routines (convolution, interpolation, special functions) - 'scipp', # Labelled multi-dimensional arrays; backs Experiment data handling - 'h5py', # HDF5 backend for scipp's HDF5 I/O (Experiment.load_hdf5) - 'matplotlib', # Plotting (posterior trace, corner, and predictive plots) - 'pooch', # Data downloader - 'darkdetect', # Detecting dark mode (system-level) - 'plopp', # Plotting library - 'jupyterlab', # Jupyter notebooks - 'pixi-kernel', # Pixi Jupyter kernel - 'ipykernel', # Jupyter kernel (required for running notebooks) - 'ipywidgets', # Widgets (needed for interactive matplotlib backends) - 'ipympl', # Matplotlib Jupyter widget backend (%matplotlib widget) - 'IPython', # Interactive Python shell - 'sympy', # Symbolic mathematics (used for expression components) + 'easyscience>=2.5.1', # The base library of the EasyScience framework. 2.5.1 adds fitting.Sampler + 'numpy', # Numerical arrays (used directly throughout the library) + 'scipy', # Numerical routines (convolution, interpolation, special functions) + 'scipp', # Labelled multi-dimensional arrays; backs Experiment data handling + 'h5py', # HDF5 backend for scipp's HDF5 I/O (Experiment.load_hdf5) + 'matplotlib', # Plotting (posterior trace, corner, and predictive plots) + 'pooch', # Data downloader + 'darkdetect', # Detecting dark mode (system-level) + 'plopp>=26.9.0', # Plotting library + 'jupyterlab', # Jupyter notebooks + 'pixi-kernel', # Pixi Jupyter kernel + 'ipykernel', # Jupyter kernel (required for running notebooks) + 'ipywidgets', # Widgets (needed for interactive matplotlib backends) + 'ipympl', # Matplotlib Jupyter widget backend (%matplotlib widget) + 'IPython', # Interactive Python shell + 'sympy', # Symbolic mathematics (used for expression components) ] [project.optional-dependencies] dev = [ - 'GitPython', # Interact with Git repositories - 'build', # Building the package - 'pre-commit', # Pre-commit hooks - 'jinja2', # Templating - 'nbmake', # Building notebooks - 'nbstripout', # Strip output from notebooks - 'nbqa', # Linting and formatting notebooks - 'pytest', # Testing - 'pytest-cov', # Test coverage - 'pytest-xdist', # Enable parallel testing - 'ruff', # Linting and formatting code - 'radon', # Code complexity and maintainability - 'validate-pyproject[all]', # Validate pyproject.toml - 'versioningit', # Automatic versioning from git tags - 'jupytext', # Jupyter notebook text format support - 'jupyterquiz', # Quizzes in Jupyter notebooks - 'pydoclint', # Docstring linter - 'docstring-parser-fork!=0.0.15', # Parser used by pydoclint. 0.0.15 splits NumPy-style defaults -> DOC105 failures - 'format-docstring', # Docstring formatter - 'docstripy', # Convert docstrings to other formats - 'interrogate', # Docstring coverage checker - 'copier', # Template management - 'mike', # MkDocs: Versioned documentation support - 'mkdocs', # Static site generator - 'mkdocs-material', # Documentation framework on top of MkDocs - 'mkdocs-autorefs', # MkDocs: Auto-references support - 'mkdocs-jupyter', # MkDocs: Jupyter notebook support - 'mkdocs-plugin-inline-svg', # MkDocs: Inline SVG support - 'mkdocs-markdownextradata-plugin', # MkDocs: Markdown extra data support, such as global variables - 'mkdocstrings-python', # MkDocs: Python docstring support - 'pyyaml', # YAML parser - 'spdx-headers', # SPDX license header validation + 'GitPython', # Interact with Git repositories + 'build', # Building the package + 'pre-commit', # Pre-commit hooks + 'jinja2', # Templating + 'nbmake', # Building notebooks + 'nbstripout', # Strip output from notebooks + 'nbqa', # Linting and formatting notebooks + 'pytest', # Testing + 'pytest-cov', # Test coverage + 'pytest-xdist', # Enable parallel testing + 'ruff', # Linting and formatting code + 'radon', # Code complexity and maintainability + 'validate-pyproject[all]', # Validate pyproject.toml + 'versioningit', # Automatic versioning from git tags + 'jupytext', # Jupyter notebook text format support + 'jupyterquiz', # Quizzes in Jupyter notebooks + 'pydoclint', # Docstring linter + 'docstring-parser-fork!=0.0.15', # Parser used by pydoclint. 0.0.15 splits NumPy-style defaults -> DOC105 failures + 'format-docstring', # Docstring formatter + 'docstripy', # Convert docstrings to other formats + 'interrogate', # Docstring coverage checker + 'copier', # Template management + 'mike', # MkDocs: Versioned documentation support + 'mkdocs', # Static site generator + 'mkdocs-material', # Documentation framework on top of MkDocs + 'mkdocs-autorefs', # MkDocs: Auto-references support + 'mkdocs-jupyter', # MkDocs: Jupyter notebook support + 'mkdocs-plugin-inline-svg', # MkDocs: Inline SVG support + 'mkdocs-markdownextradata-plugin', # MkDocs: Markdown extra data support, such as global variables + 'mkdocstrings-python', # MkDocs: Python docstring support + 'pyyaml', # YAML parser + 'spdx-headers', # SPDX license header validation ] [project.urls] @@ -106,7 +106,7 @@ packages = ['src/easydynamics'] allow-direct-references = true [tool.hatch.version] -source = 'versioningit' # Use versioningit to manage the version +source = 'versioningit' # Use versioningit to manage the version ################################ # Configuration for versioningit @@ -122,9 +122,9 @@ source = 'versioningit' # Use versioningit to manage the version # pixi.lock update without any changes to the source code. [tool.versioningit.format] -distance = '{base_version}+dev{distance}' # example: 1.2.3.post4+dev3 -dirty = '{base_version}+dirty{distance}' # example: 0.5.8+dirty3 -distance-dirty = '{base_version}+devdirty{distance}' # example: 0.5.8+devdirty3 +distance = '{base_version}+dev{distance}' # example: 1.2.3.post4+dev3 +dirty = '{base_version}+dirty{distance}' # example: 0.5.8+dirty3 +distance-dirty = '{base_version}+devdirty{distance}' # example: 0.5.8+devdirty3 # Configure how versioningit detects versions from Git # - 'match' ensures it only considers tags starting with 'v' @@ -142,9 +142,9 @@ default-tag = 'v999.0.0' # https://interrogate.readthedocs.io/en/latest/ [tool.interrogate] -fail-under = 0 # Minimum docstring coverage percentage to pass +fail-under = 0 # Minimum docstring coverage percentage to pass verbose = 1 -#exclude = ['src/**/__init__.py'] +# exclude = ['src/**/__init__.py'] ####################################### # Configuration for coverage/pytest-cov @@ -154,13 +154,13 @@ verbose = 1 # https://coverage.readthedocs.io/en/latest/ [tool.coverage.run] -branch = true # Measure branch coverage as well -source = ['src'] # Limit coverage to the source code directory +branch = true # Measure branch coverage as well +source = ['src'] # Limit coverage to the source code directory [tool.coverage.report] show_missing = true # Show missing lines -skip_covered = false # Skip files with 100% coverage in the report -fail_under = 0 # Minimum coverage percentage to pass +skip_covered = false # Skip files with 100% coverage in the report +fail_under = 0 # Minimum coverage percentage to pass ########################## # Configuration for pytest @@ -188,74 +188,74 @@ testpaths = ['tests'] exclude = ['tmp'] indent-width = 4 line-length = 99 # See also `max-line-length` in [tool.ruff.lint.pycodestyle] -preview = true # Enable new rules that are not yet stable, like DOC +preview = true # Enable new rules that are not yet stable, like DOC # Formatting options for Ruff [tool.ruff.format] -docstring-code-format = true # Whether to format code snippets in docstrings -docstring-code-line-length = 99 # Line length for code snippets in docstrings -indent-style = 'space' # PEP 8 recommends using spaces over tabs -quote-style = 'single' # But double quotes in docstrings (PEP 8, PEP 257) +docstring-code-format = true # Whether to format code snippets in docstrings +docstring-code-line-length = 99 # Line length for code snippets in docstrings +indent-style = 'space' # PEP 8 recommends using spaces over tabs +quote-style = 'single' # But double quotes in docstrings (PEP 8, PEP 257) # Linting rules to use with Ruff [tool.ruff.lint] select = [ # Various rules - #'C90', # https://docs.astral.sh/ruff/rules/#mccabe-c90 - #'D', # https://docs.astral.sh/ruff/rules/#pydocstyle-d - 'F', # https://docs.astral.sh/ruff/rules/#pyflakes-f - #'FLY', # https://docs.astral.sh/ruff/rules/#flynt-fly - #'FURB', # https://docs.astral.sh/ruff/rules/#refurb-furb - 'I', # https://docs.astral.sh/ruff/rules/#isort-i - #'N', # https://docs.astral.sh/ruff/rules/#pep8-naming-n - 'NPY', # https://docs.astral.sh/ruff/rules/#numpy-specific-rules-npy - #'PGH', # https://docs.astral.sh/ruff/rules/#pygrep-hooks-pgh - 'PERF', # https://docs.astral.sh/ruff/rules/#perflint-perf + # 'C90', # https://docs.astral.sh/ruff/rules/#mccabe-c90 + # 'D', # https://docs.astral.sh/ruff/rules/#pydocstyle-d + 'F', # https://docs.astral.sh/ruff/rules/#pyflakes-f + # 'FLY', # https://docs.astral.sh/ruff/rules/#flynt-fly + # 'FURB', # https://docs.astral.sh/ruff/rules/#refurb-furb + 'I', # https://docs.astral.sh/ruff/rules/#isort-i + # 'N', # https://docs.astral.sh/ruff/rules/#pep8-naming-n + 'NPY', # https://docs.astral.sh/ruff/rules/#numpy-specific-rules-npy + # 'PGH', # https://docs.astral.sh/ruff/rules/#pygrep-hooks-pgh + 'PERF', # https://docs.astral.sh/ruff/rules/#perflint-perf 'RUF', # https://docs.astral.sh/ruff/rules/#ruff-specific-rules-ruf - #'TRY', # https://docs.astral.sh/ruff/rules/#tryceratops-try # overly restrictive - 'UP', # https://docs.astral.sh/ruff/rules/#pyupgrade-up + # 'TRY', # https://docs.astral.sh/ruff/rules/#tryceratops-try # overly restrictive + 'UP', # https://docs.astral.sh/ruff/rules/#pyupgrade-up # pycodestyle (E, W) rules - 'E', # https://docs.astral.sh/ruff/rules/#error-e - 'W', # https://docs.astral.sh/ruff/rules/#warning-w + 'E', # https://docs.astral.sh/ruff/rules/#error-e + 'W', # https://docs.astral.sh/ruff/rules/#warning-w # Pylint (PL) rules - #'PLC', # https://docs.astral.sh/ruff/rules/#convention-plc - #'PLE', # https://docs.astral.sh/ruff/rules/#error-ple - #'PLR', # https://docs.astral.sh/ruff/rules/#refactor-plr - #'PLW', # https://docs.astral.sh/ruff/rules/#warning-plw # Good to enable + # 'PLC', # https://docs.astral.sh/ruff/rules/#convention-plc + # 'PLE', # https://docs.astral.sh/ruff/rules/#error-ple + # 'PLR', # https://docs.astral.sh/ruff/rules/#refactor-plr + # 'PLW', # https://docs.astral.sh/ruff/rules/#warning-plw # Good to enable # flake8 rules - 'A', # https://docs.astral.sh/ruff/rules/#flake8-builtins-a - 'ANN', # https://docs.astral.sh/ruff/rules/#flake8-annotations-ann - 'ARG', # https://docs.astral.sh/ruff/rules/#flake8-unused-arguments-arg - 'ASYNC', # https://docs.astral.sh/ruff/rules/#flake8-async-async - 'B', # https://docs.astral.sh/ruff/rules/#flake8-bugbear-b - #'BLE', # https://docs.astral.sh/ruff/rules/#flake8-blind-except-ble # enable when base classes have been merged in + 'A', # https://docs.astral.sh/ruff/rules/#flake8-builtins-a + 'ANN', # https://docs.astral.sh/ruff/rules/#flake8-annotations-ann + 'ARG', # https://docs.astral.sh/ruff/rules/#flake8-unused-arguments-arg + 'ASYNC', # https://docs.astral.sh/ruff/rules/#flake8-async-async + 'B', # https://docs.astral.sh/ruff/rules/#flake8-bugbear-b + # 'BLE', # https://docs.astral.sh/ruff/rules/#flake8-blind-except-ble # enable when base classes have been merged in 'C4', # https://docs.astral.sh/ruff/rules/#flake8-comprehensions-c4 - 'COM', # https://docs.astral.sh/ruff/rules/#flake8-commas-com - 'DTZ', # https://docs.astral.sh/ruff/rules/#flake8-datetimez-dtz - #'EM', # https://docs.astral.sh/ruff/rules/#flake8-errmsg-em # Unsure if I want it - requires all error messages to be rewritten - 'FA', # https://docs.astral.sh/ruff/rules/#flake8-future-annotations-fa - #'FBT', # https://docs.astral.sh/ruff/rules/#flake8-boolean-trap-fbt # should eventually be enabled, but may break serialization - 'FIX', # https://docs.astral.sh/ruff/rules/#flake8-fixme-fix - 'G', # https://docs.astral.sh/ruff/rules/#flake8-logging-format-g - 'ICN', # https://docs.astral.sh/ruff/rules/#flake8-import-conventions-icn - 'INP', # https://docs.astral.sh/ruff/rules/#flake8-no-pep420-inp - 'ISC', # https://docs.astral.sh/ruff/rules/#flake8-implicit-str-concat-isc - 'LOG', # https://docs.astral.sh/ruff/rules/#flake8-logging-log - 'PIE', # https://docs.astral.sh/ruff/rules/#flake8-pie-pie - #'PT', # https://docs.astral.sh/ruff/rules/#flake8-pytest-style-pt # Should eventually be enabled, but quite a few errors to fix + 'COM', # https://docs.astral.sh/ruff/rules/#flake8-commas-com + 'DTZ', # https://docs.astral.sh/ruff/rules/#flake8-datetimez-dtz + # 'EM', # https://docs.astral.sh/ruff/rules/#flake8-errmsg-em # Unsure if I want it - requires all error messages to be rewritten + 'FA', # https://docs.astral.sh/ruff/rules/#flake8-future-annotations-fa + # 'FBT', # https://docs.astral.sh/ruff/rules/#flake8-boolean-trap-fbt # should eventually be enabled, but may break serialization + 'FIX', # https://docs.astral.sh/ruff/rules/#flake8-fixme-fix + 'G', # https://docs.astral.sh/ruff/rules/#flake8-logging-format-g + 'ICN', # https://docs.astral.sh/ruff/rules/#flake8-import-conventions-icn + 'INP', # https://docs.astral.sh/ruff/rules/#flake8-no-pep420-inp + 'ISC', # https://docs.astral.sh/ruff/rules/#flake8-implicit-str-concat-isc + 'LOG', # https://docs.astral.sh/ruff/rules/#flake8-logging-log + 'PIE', # https://docs.astral.sh/ruff/rules/#flake8-pie-pie + # 'PT', # https://docs.astral.sh/ruff/rules/#flake8-pytest-style-pt # Should eventually be enabled, but quite a few errors to fix 'PTH', # https://docs.astral.sh/ruff/rules/#flake8-use-pathlib-pth 'PYI', # https://docs.astral.sh/ruff/rules/#flake8-pyi-pyi 'RET', # https://docs.astral.sh/ruff/rules/#flake8-return-ret 'RSE', # https://docs.astral.sh/ruff/rules/#flake8-raise-rse - 'S', # https://docs.astral.sh/ruff/rules/#flake8-bandit-s + 'S', # https://docs.astral.sh/ruff/rules/#flake8-bandit-s 'SIM', # https://docs.astral.sh/ruff/rules/#flake8-simplify-sim 'SLF', # https://docs.astral.sh/ruff/rules/#flake8-self-slf - 'SLOT', # https://docs.astral.sh/ruff/rules/#flake8-slots-slot + 'SLOT', # https://docs.astral.sh/ruff/rules/#flake8-slots-slot 'T20', # https://docs.astral.sh/ruff/rules/#flake8-print-t20 - 'TC', # https://docs.astral.sh/ruff/rules/#flake8-type-checking-tc - 'TD', # https://docs.astral.sh/ruff/rules/#flake8-todos-td + 'TC', # https://docs.astral.sh/ruff/rules/#flake8-type-checking-tc + 'TD', # https://docs.astral.sh/ruff/rules/#flake8-todos-td 'TID', # https://docs.astral.sh/ruff/rules/#flake8-tidy-imports-tid ] @@ -263,28 +263,28 @@ select = [ # Ignore specific rules globally ignore = [ - 'COM812', # https://docs.astral.sh/ruff/rules/missing-trailing-comma/ + 'COM812', # https://docs.astral.sh/ruff/rules/missing-trailing-comma/ # The following is replaced by 'D'/[tool.ruff.lint.pydocstyle] and [tool.pydoclint] 'DOC', # https://docs.astral.sh/ruff/rules/#pydoclint-doc # Disable, as [tool.format_docstring] split one-line docstrings into the canonical multi-line layout - 'D200', # https://docs.astral.sh/ruff/rules/unnecessary-multiline-docstring/ + 'D200', # https://docs.astral.sh/ruff/rules/unnecessary-multiline-docstring/ ] # Ignore specific rules in certain files or directories [tool.ruff.lint.per-file-ignores] '*/__init__.py' = [ - 'F401', # re-exports are intentional in __init__.py + 'F401', # re-exports are intentional in __init__.py ] 'tests/**' = [ - 'ANN', # https://docs.astral.sh/ruff/rules/#flake8-annotations-ann - 'D', # https://docs.astral.sh/ruff/rules/#pydocstyle-d - 'DOC', # https://docs.astral.sh/ruff/rules/#pydoclint-doc - 'INP001', # https://docs.astral.sh/ruff/rules/implicit-namespace-package/ - 'S101', # https://docs.astral.sh/ruff/rules/assert/ - 'SLF', # https://docs.astral.sh/ruff/rules/#flake8-self-slf # may want to eventually enable it, but accessing private methods is quite useful in tests + 'ANN', # https://docs.astral.sh/ruff/rules/#flake8-annotations-ann + 'D', # https://docs.astral.sh/ruff/rules/#pydocstyle-d + 'DOC', # https://docs.astral.sh/ruff/rules/#pydoclint-doc + 'INP001', # https://docs.astral.sh/ruff/rules/implicit-namespace-package/ + 'S101', # https://docs.astral.sh/ruff/rules/assert/ + 'SLF', # https://docs.astral.sh/ruff/rules/#flake8-self-slf # may want to eventually enable it, but accessing private methods is quite useful in tests ] 'docs/**' = [ - 'INP001', # https://docs.astral.sh/ruff/rules/implicit-namespace-package/ - 'T201', # https://docs.astral.sh/ruff/rules/print/ + 'INP001', # https://docs.astral.sh/ruff/rules/implicit-namespace-package/ + 'T201', # https://docs.astral.sh/ruff/rules/print/ ] # Specific options for certain rules @@ -306,7 +306,7 @@ max-complexity = 10 # https://peps.python.org/pep-0008/#maximum-line-length # Use 99 characters as the project-wide maximum for regular code lines. # Use 99 characters for docstrings. -max-line-length = 99 # See also `line-length` in [tool.ruff] +max-line-length = 99 # See also `line-length` in [tool.ruff] max-doc-length = 99 [tool.ruff.lint.pydocstyle] @@ -333,7 +333,7 @@ max-positional-args = 6 # the parameter declarations in the code (in function's signature). [tool.pydoclint] -#exclude = '\.' # Temporarily disable pydoclint until we are ready +# exclude = '\.' # Temporarily disable pydoclint until we are ready style = 'numpy' check-style-mismatch = true check-arg-defaults = true @@ -347,7 +347,7 @@ allow-init-docstring = true # https://github.com/jsh9/format-docstring [tool.format_docstring] -#exclude = '\.' # Temporarily disable format-docstring until we are ready +# exclude = '\.' # Temporarily disable format-docstring until we are ready docstring_style = 'numpy' line_length = 99 fix_rst_backticks = true diff --git a/src/easydynamics/__init__.py b/src/easydynamics/__init__.py index de7391225..a0f9d570d 100644 --- a/src/easydynamics/__init__.py +++ b/src/easydynamics/__init__.py @@ -38,6 +38,7 @@ from easydynamics.sample_model import Polynomial from easydynamics.sample_model import ResolutionModel from easydynamics.sample_model import SampleModel +from easydynamics.sample_model import StretchedExponential from easydynamics.sample_model import Voigt from easydynamics.settings import ConvolutionSettings from easydynamics.settings import DetailedBalanceSettings @@ -81,6 +82,7 @@ 'PosteriorSummary', 'ResolutionModel', 'SampleModel', + 'StretchedExponential', 'Voigt', 'detailed_balance_factor', 'hbar', diff --git a/src/easydynamics/analysis/analysis.py b/src/easydynamics/analysis/analysis.py index f06a818bd..48910fb7f 100644 --- a/src/easydynamics/analysis/analysis.py +++ b/src/easydynamics/analysis/analysis.py @@ -11,7 +11,7 @@ from easyscience.fitting.minimizers.utils import FitResults from easyscience.fitting.multi_fitter import MultiFitter from easyscience.variable import Parameter -from plopp.backends.matplotlib.figure import InteractiveFigure +from plopp.backends.matplotlib.figure import WidgetFigure from scipp import UnitError from easydynamics.analysis.analysis1d import Analysis1d @@ -391,7 +391,7 @@ def plot_data_and_model( plot_residuals: bool = False, energy: sc.Variable | None = None, **kwargs: dict[str, Any], - ) -> InteractiveFigure: + ) -> WidgetFigure: """ Plot the experimental data and the model prediction. @@ -426,8 +426,8 @@ def plot_data_and_model( Returns ------- - InteractiveFigure - A Plopp InteractiveFigure containing the plot of the data and model. + WidgetFigure + A Plopp WidgetFigure containing the plot of the data and model. """ verify_Q_index(Q_index=Q_index, Q=self.Q, allow_none=True) if Q_index is not None: @@ -647,7 +647,7 @@ def plot_parameters( self, names: str | list[str] | None = None, **kwargs: dict[str, Any], - ) -> InteractiveFigure: + ) -> WidgetFigure: """ Plot fitted parameters as a function of Q. @@ -668,8 +668,8 @@ def plot_parameters( Returns ------- - InteractiveFigure - A Plopp InteractiveFigure containing the plot of the parameters. + WidgetFigure + A Plopp WidgetFigure containing the plot of the parameters. """ ds = self.parameters_to_dataset() diff --git a/src/easydynamics/analysis/analysis1d.py b/src/easydynamics/analysis/analysis1d.py index 6ac0b7eb3..0e0acdab2 100644 --- a/src/easydynamics/analysis/analysis1d.py +++ b/src/easydynamics/analysis/analysis1d.py @@ -11,7 +11,7 @@ from easyscience.fitting.minimizers.utils import FitResults from easyscience.variable import DescriptorNumber from easyscience.variable import Parameter -from plopp.backends.matplotlib.figure import InteractiveFigure +from plopp.backends.matplotlib.figure import WidgetFigure from easydynamics.analysis.analysis_base import AnalysisBase from easydynamics.analysis.posterior_labels import ParameterLabels @@ -452,7 +452,7 @@ def plot_data_and_model( plot_residuals: bool = False, energy: sc.Variable | None = None, **kwargs: dict[str, Any], - ) -> InteractiveFigure: + ) -> WidgetFigure: """ Plot the experimental data and the model prediction for the chosen Q index. Optionally also plot the individual components of the model. @@ -476,7 +476,7 @@ def plot_data_and_model( Returns ------- - InteractiveFigure + WidgetFigure A plot of the data and model. """ data_and_model = self.data_and_model_to_datagroup( diff --git a/src/easydynamics/analysis/parameter_analysis.py b/src/easydynamics/analysis/parameter_analysis.py index 1e99f006a..6413a963f 100644 --- a/src/easydynamics/analysis/parameter_analysis.py +++ b/src/easydynamics/analysis/parameter_analysis.py @@ -11,7 +11,7 @@ from easyscience.fitting.multi_fitter import MultiFitter from easyscience.variable import Parameter from matplotlib import rcParams -from plopp.backends.matplotlib.figure import InteractiveFigure +from plopp.backends.matplotlib.figure import WidgetFigure from easydynamics.analysis.analysis import Analysis from easydynamics.analysis.fit_binding import FitBinding @@ -412,9 +412,7 @@ def _chain_parameters(self) -> list[Parameter]: parameters.setdefault(parameter.unique_name, parameter) return list(parameters.values()) - def plot( - self, names: str | list[str] | None = None, **kwargs: dict[str, Any] - ) -> InteractiveFigure: + def plot(self, names: str | list[str] | None = None, **kwargs: dict[str, Any]) -> WidgetFigure: """ Plot the parameters and fit results. @@ -427,7 +425,7 @@ def plot( Returns ------- - InteractiveFigure + WidgetFigure An interactive figure containing the plots of the parameters and fit results. Raises diff --git a/src/easydynamics/analysis/posterior_sampling.py b/src/easydynamics/analysis/posterior_sampling.py index 7d2615938..e1187a9c0 100644 --- a/src/easydynamics/analysis/posterior_sampling.py +++ b/src/easydynamics/analysis/posterior_sampling.py @@ -42,7 +42,7 @@ from easyscience.variable import Parameter from ipywidgets import VBox from matplotlib.figure import Figure - from plopp.backends.matplotlib.figure import InteractiveFigure + from plopp.backends.matplotlib.figure import WidgetFigure from easydynamics.analysis.posterior import BoundsSuggestions @@ -1584,7 +1584,7 @@ def plot_posterior_predictive( credible_interval: float = 68.0, Q_index: int | None = None, **kwargs: dict[str, Any], - ) -> Figure | InteractiveFigure: + ) -> Figure | WidgetFigure: """ Plot the data against the credible band implied by the posterior. @@ -1610,7 +1610,7 @@ def plot_posterior_predictive( Returns ------- - Figure | InteractiveFigure + Figure | WidgetFigure The matplotlib Figure for one Q, or the plopp figure with a Q slider. Raises @@ -1748,7 +1748,7 @@ def _predictive_with_q_slider( n_draws: int, credible_interval: float, **kwargs: dict[str, Any], - ) -> InteractiveFigure: + ) -> WidgetFigure: """ Build the posterior-predictive figure with a Q slider from the per-Q chains. @@ -1770,7 +1770,7 @@ def _predictive_with_q_slider( Returns ------- - InteractiveFigure + WidgetFigure The plopp figure with its Q slider. Raises diff --git a/src/easydynamics/experiment/experiment.py b/src/easydynamics/experiment/experiment.py index 058df9abe..3beeeb737 100644 --- a/src/easydynamics/experiment/experiment.py +++ b/src/easydynamics/experiment/experiment.py @@ -7,7 +7,7 @@ import numpy as np import plopp as pp import scipp as sc -from plopp.backends.matplotlib.figure import InteractiveFigure +from plopp.backends.matplotlib.figure import WidgetFigure from scipp.io import load_hdf5 as sc_load_hdf5 from scipp.io import save_hdf5 as sc_save_hdf5 @@ -440,7 +440,7 @@ def plot_data( slicer: bool = False, transpose_axes: bool = False, **kwargs: dict, - ) -> InteractiveFigure: + ) -> WidgetFigure: """ Plot the dataset using plopp: https://scipp.github.io/plopp/. @@ -456,7 +456,7 @@ def plot_data( Returns ------- - InteractiveFigure + WidgetFigure A plot of the data and model. Raises diff --git a/src/easydynamics/sample_model/__init__.py b/src/easydynamics/sample_model/__init__.py index 07f434a18..97a788aab 100644 --- a/src/easydynamics/sample_model/__init__.py +++ b/src/easydynamics/sample_model/__init__.py @@ -12,6 +12,7 @@ from easydynamics.sample_model.components.gaussian import Gaussian from easydynamics.sample_model.components.lorentzian import Lorentzian from easydynamics.sample_model.components.polynomial import Polynomial +from easydynamics.sample_model.components.stretched_exponential import StretchedExponential from easydynamics.sample_model.components.voigt import Voigt from easydynamics.sample_model.diffusion_model.brownian_translational_diffusion import ( BrownianTranslationalDiffusion, @@ -40,5 +41,6 @@ 'Polynomial', 'ResolutionModel', 'SampleModel', + 'StretchedExponential', 'Voigt', ] diff --git a/src/easydynamics/sample_model/components/__init__.py b/src/easydynamics/sample_model/components/__init__.py index 940f41a7e..62206663d 100644 --- a/src/easydynamics/sample_model/components/__init__.py +++ b/src/easydynamics/sample_model/components/__init__.py @@ -10,6 +10,7 @@ from easydynamics.sample_model.components.gaussian import Gaussian from easydynamics.sample_model.components.lorentzian import Lorentzian from easydynamics.sample_model.components.polynomial import Polynomial +from easydynamics.sample_model.components.stretched_exponential import StretchedExponential from easydynamics.sample_model.components.voigt import Voigt __all__ = [ @@ -20,5 +21,6 @@ 'Gaussian', 'Lorentzian', 'Polynomial', + 'StretchedExponential', 'Voigt', ] diff --git a/src/easydynamics/sample_model/components/stretched_exponential.py b/src/easydynamics/sample_model/components/stretched_exponential.py new file mode 100644 index 000000000..9e33fd5c0 --- /dev/null +++ b/src/easydynamics/sample_model/components/stretched_exponential.py @@ -0,0 +1,784 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +r""" +The stretched exponential (Kohlrausch-Williams-Watts) relaxation, Fourier transformed to energy. + +The model is defined in time, as $I(t) = A e^{-(|t| / \tau)^\beta}$. To evaluate the Fourier +transform, we need to evaluate the following integral: + +$$ G_\beta(w) = \int_0^\infty e^{-u^\beta} \cos(w u) \, du $$ + +The problem is that $\cos(w u)$ oscillates forever at constant amplitude, and $e^{-u^\beta}$ decays +extremely slowly once $\beta$ is small: at $\beta = 0.3$ it is still significant at $u \sim 10^5$. +Evaluating $G_\beta$ therefore means summing an enormous number of nearly cancelling oscillations. +Furthermore, large $w$ makes the oscillation fast, so the steps must be tiny, while small $\beta$ +stretches the tail over more decades, so the range must be huge. + +An FFT of the sampled relaxation cannot span that tail at small $\beta$ and, because $e^{-(t / +\tau)^\beta}$ has a cusp at $t = 0$, converges only as $(\Delta t)^{1 + \beta}$. +``scipy.stats.levy_stable`` is mathematically the same function, but it is 70-100x slower and not +accurate for $\alpha$ close to 1. + +Instead, the integral is evaluated in the complex plane described in :func:`_kww_shape`. + +The four helpers below are module-level functions rather than methods. :func:`_quadrature_nodes` +and :func:`_reduced_hwhm` must be, because ``@lru_cache`` on a method would key on ``self``: every +instance would rebuild its own grid and the cache would keep every component ever created alive. +:func:`_kww_shape` is a pure function of ``(w, beta)`` that touches no component state: area, +center and units are applied afterwards in ``_evaluate_values``, so it does not belong in the +class, and :func:`_hwhm_dependency_expression` only assembles a string from module constants. + +:func:`_reduced_hwhm` is not on any evaluation path. It solves the half width exactly, and the +closed form that ``width`` actually resolves through is fitted to it; it stays here as the +reference that defines `_HWHM_POLY_COEFFS` and as what the tests check that fit against. +""" + +from __future__ import annotations + +import contextlib +from functools import lru_cache +from typing import TYPE_CHECKING + +import numpy as np +from easyscience.variable import DescriptorNumber +from easyscience.variable import Parameter +from scipp import UnitError +from scipy.optimize import brentq + +from easydynamics.sample_model.components.mixins import CreateParametersMixin +from easydynamics.sample_model.components.model_component import ModelComponent +from easydynamics.utils.utils import Numeric +from easydynamics.utils.utils import convert_value_unit +from easydynamics.utils.utils import hbar + +if TYPE_CHECKING: + import scipp as sc + +MINIMUM_RELAXATION_TIME = 1e-10 # ps. Avoids a division by zero in hbar / tau +MINIMUM_BETA = 0.05 # Below this the quadrature grid can no longer resolve the time tail +MAXIMUM_BETA = 2.0 # Above this the transform is no longer positive, so unphysical + +# Tuning of the exp-sinh quadrature used by _kww_shape. The step and the half-range were chosen +# together so the transform reproduces its analytic special cases (beta = 1 and beta = 2) to +# ~1e-11 relative over w in [0, 1e6]; see test_stretched_exponential.py. +_QUAD_STEP = 0.03 +_QUAD_HALF_RANGE = 4.3 +# exp(-u**beta) has fallen below exp(-_QUAD_DECADES) at the far end of the grid. +_QUAD_DECADES = 50.0 +# The quadrature holds one (energies x nodes) temporary, so long energy axes are evaluated in +# blocks to keep that intermediate at a few tens of MB instead of scaling with the axis. +_MAX_BLOCK_ELEMENTS = 1 << 21 +# The reduced half width spans 1.67 at beta = 2 down to ~7e-27 at beta = 0.05, so it is +# bracketed in log10(w). The lower end stays clear of the subnormal range, where the +# 1 / (w sin theta) rescaling inside _kww_shape would overflow. +_HWHM_LOG_BRACKET = (-300.0, 3.0) + +# The width is exposed as a dependent Parameter, and easyscience resolves a dependency from a +# string expression, so the root solved by _reduced_hwhm has to be written in closed form. It +# very nearly is one: beta * ln(HWHM) - ln(beta) is a small, smooth function of beta alone, so +# HWHM = exp((ln(beta) + p(beta)) / beta) with p the polynomial below, in ascending powers. The +# coefficients are a degree 8 least-squares fit over 400 log-spaced beta in [MINIMUM_BETA, +# MAXIMUM_BETA], and reproduce _reduced_hwhm to 1.4e-5 relative for beta >= 0.1. Below that the +# fit degrades to ~1e-2, but so does the root it is fitted to, and both are tens of orders of +# magnitude under any energy grid step by then. Regenerate, or check these against the root they +# are fitted to, with tools/fit_kww_hwhm_coefficients.py. +_HWHM_POLY_COEFFS = ( + -5.5129123891501144e-05, + -0.24257355020449564, + 0.2875849387004671, + -0.04773226324334423, + -0.04385316873071495, + 0.11064865473723488, + -0.09366816058690182, + 0.03431170767881221, + -0.004659825239152873, +) + + +@lru_cache(maxsize=8) +def _quadrature_nodes(half_steps: int) -> tuple[np.ndarray, np.ndarray]: + r""" + Build the exp-sinh quadrature nodes and weights for ``2 * half_steps + 1`` points. + + The exp-sinh (double-exponential) rule maps the trapezoidal rule on the whole real line onto + the half line via $s = e^{(\pi / 2) \sinh t}$, which makes the node density follow the decades + of $s$ rather than its absolute size. The grids are small and reused across evaluations, so + they are cached rather than rebuilt on every call. + + Parameters + ---------- + half_steps : int + Number of steps of size ``_QUAD_STEP`` on each side of $t = 0$. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + The node positions $s$ on the half line and their quadrature weights. + """ + t = np.arange(-half_steps, half_steps + 1) * _QUAD_STEP + nodes = np.exp(0.5 * np.pi * np.sinh(t)) + weights = _QUAD_STEP * nodes * (0.5 * np.pi) * np.cosh(t) + return nodes, weights + + +def _kww_shape(w: np.ndarray, beta: float) -> np.ndarray: + r""" + Evaluate the reduced Fourier cosine transform of a stretched exponential. + + $$ G_\beta(w) = \int_0^\infty e^{-u^\beta} \cos(w u) \, du $$ + + The physical spectrum is $\frac{A}{\pi \Gamma} G_\beta(x / \Gamma)$ with $\Gamma = \hbar / + \tau$, so $w$ is the reduced energy $x / \Gamma$ and everything else is a scale factor. + + The integrand oscillates without decaying, so it is integrated along the rotated ray $u = s + e^{i \theta}$ instead of the real axis. Writing $\cos(w u) = \mathrm{Re}\, e^{i w u}$ makes + the integrand analytic in a wedge around the positive real axis: $e^{-u^\beta}$ keeps decaying + as long as $\arg(u) < \pi / (2 \beta)$, and the arc at infinity vanishes, so by Cauchy's + theorem the ray may be swung up to that angle without changing the value. + + On the tilted ray $e^{i w u}$ becomes $e^{i w s \cos\theta} e^{-w s \sin\theta}$: the + oscillation now decays with increasing $w$. The faster the oscillation, the faster it is + damped, so the number of oscillations before the integrand dies is bounded *independently of* + $w$. + + The angle $\theta = \pi / (4 \beta)$ is half of the $\pi / (2 \beta)$ limit, a deliberate + safety margin; it is capped at $0.45 \pi$ because for $\beta < 0.5$ the formula would exceed + $\pi / 2$ and the ray would cross the imaginary axis, where $e^{-u^\beta}$ grows. + + Writing $e^{-u^\beta + i w u}$ out in real form on that ray, with the extra $\theta$ in the + phase coming from $du = e^{i \theta} ds$, turns the integral into + + $$ G_\beta(w) = \int_0^\infty e^{-P(s)} \cos(Q(s) + \theta) \, ds $$ + + with $P = \cos(\beta \theta) s^\beta + w \sin(\theta) s$ (the decay) and $Q = w \cos(\theta) s + - \sin(\beta \theta) s^\beta$ (what is left of the oscillation), which an exp-sinh rule on a + grid rescaled to the decay length of $P$ integrates to near machine precision for every $w$. + + The rescaling is the other half of the robustness. $P$ has two terms, and whichever dies first + sets the decay length: $\cos(\beta \theta)^{-1 / \beta}$ for the $s^\beta$ term, $1 / (w + \sin\theta)$ for the $w s$ term. Taking the smaller of the two and stretching the grid onto it + means the integrand always falls off around $s = 1$, so a single node layout serves every $w$ + and every $\beta$. + + Accuracy is therefore governed by $w$ and $\beta$ alone. For $\beta \le 1$ the result tracks + the exact power-law tail out to at least $w = 10^8$. For $\beta > 1$ the profile decays fast + enough to reach the floor of double precision, and values below roughly $10^{-11}$ of the peak + are noise. + + Parameters + ---------- + w : np.ndarray + Reduced energies. Only $|w|$ matters: the transform is even in $w$. + beta : float + Stretching exponent, in ``[MINIMUM_BETA, MAXIMUM_BETA]``. + + Returns + ------- + np.ndarray + The transform, clipped at zero. It is non-negative for every $\beta \le 2$; the clip only + removes rounding noise where the true value has underflowed. + """ + abs_w = np.abs(w) + + theta = min(0.45 * np.pi, np.pi / (4.0 * beta)) + cos_rotated = np.cos(beta * theta) + sin_rotated = np.sin(beta * theta) + sin_theta = np.sin(theta) + cos_theta = np.cos(theta) + + # Reach far enough out in s that exp(-cos(beta theta) s**beta) has died: small beta stretches + # the tail over many more decades, so the grid has to follow it. This inverts the node map + # s = exp(pi / 2 sinh(t)) at the s where cos(beta theta) s**beta equals _QUAD_DECADES, i.e. + # where the decaying factor has fallen to exp(-50). + half_range = np.arcsinh((2.0 / np.pi) * np.log(_QUAD_DECADES / cos_rotated) / beta) + half_steps = int(np.ceil(max(half_range, _QUAD_HALF_RANGE) / _QUAD_STEP)) + nodes, weights = _quadrature_nodes(half_steps) + nodes_beta = nodes**beta + + # Rescale the grid onto whichever of the two exponents in P decays first, so the integrand + # always falls off around s = 1 and the same node layout serves every w. beta_scale is where + # the s**beta term dies, w_scale where the w s term does; w_scale is infinite at w = 0, where + # there is no oscillation to damp and beta_scale is the only length in the problem. + beta_scale = cos_rotated ** (-1.0 / beta) + w_scale = np.divide(1.0, abs_w * sin_theta, out=np.full(abs_w.shape, np.inf), where=abs_w > 0) + scale = np.minimum(beta_scale, w_scale) + + out = np.empty(abs_w.shape) + block = max(1, _MAX_BLOCK_ELEMENTS // nodes.size) + for start in range(0, abs_w.size, block): + stop = start + block + block_w = abs_w[start:stop, np.newaxis] + block_scale = scale[start:stop, np.newaxis] + + s = block_scale * nodes + s_beta = block_scale**beta * nodes_beta + decaying = cos_rotated * s_beta + (block_w * sin_theta) * s + oscillating = (block_w * cos_theta) * s - sin_rotated * s_beta + + integral = (np.exp(-decaying) * np.cos(oscillating + theta) * weights).sum(axis=1) + out[start:stop] = integral * scale[start:stop] + + return np.maximum(out, 0.0) + + +@lru_cache(maxsize=128) +def _reduced_hwhm(beta: float) -> float: + r""" + Solve for the half width at half maximum of $G_\beta$ in reduced energy $w = x / \Gamma$. + + $\Gamma = \hbar / \tau$ is the HWHM only at $\beta = 1$. Below that the profile sharpens + dramatically: the peak $G_\beta(0)$ grows while the area stays fixed, so the half width + collapses far faster than $\Gamma$ suggests -- 0.22 $\Gamma$ at $\beta = 0.5$, 2.7e-4 $\Gamma$ + at $\beta = 0.2$, 8e-11 $\Gamma$ at $\beta = 0.1$. There is no closed form, so the crossing of + $G_\beta(w) = \tfrac{1}{2} G_\beta(0)$ is bracketed numerically. + + The search runs in $\log_{10} w$ because the root ranges over more than 25 decades across the + supported $\beta$; a linear bracket could not resolve the small-$\beta$ end. Each solve costs + a few tens of single-point quadratures, under a millisecond, and the result is cached because + it depends on $\beta$ alone. + + The root is well conditioned down to $\beta \approx 0.1$ and increasingly poorly below it: by + $\beta = 0.05$ the profile takes some 60 e-folds of $w$ to fall by half, so $G_\beta$ is nearly + flat in $\log w$ near the crossing and the result is only good to about a percent. That is + immaterial for the one thing the half width is used for -- comparing against the convolution's + energy grid -- because at those $\beta$ it is already tens of orders of magnitude below any + grid step, and the comparison lands the same way either way. + + Parameters + ---------- + beta : float + Stretching exponent, in ``[MINIMUM_BETA, MAXIMUM_BETA]``. + + Returns + ------- + float + The HWHM in units of $\Gamma$. Multiply by $\Gamma$ to get an energy. + """ + half_peak = 0.5 * _kww_shape(np.array([0.0]), beta)[0] + + def offset(log_w: float) -> float: + return _kww_shape(np.array([10.0**log_w]), beta)[0] - half_peak + + return float(10.0 ** brentq(offset, *_HWHM_LOG_BRACKET, xtol=1e-12)) + + +def _hwhm_dependency_expression() -> str: + r""" + Write the half width of :func:`_reduced_hwhm` as an easyscience dependency expression. + + The expression evaluates ``exp((log(beta) + p(beta)) / beta)`` from `_HWHM_POLY_COEFFS`, in + Horner form, and multiplies it by $\hbar / \tau$ to turn the reduced half width into an energy. + + ``exp`` and ``log`` are not defined on DescriptorNumbers, so the reduced half width is built + from ``b.value`` and is a plain float; the surrounding ``* hbar / tau`` is what makes the whole + expression return a DescriptorNumber, as easyscience requires. Referring to ``b.value`` rather + than ``b`` strips beta's uncertainty, which is why the width carries only the relaxation time's + contribution. When exp and log of DescriptorNumbers are supported, this can be rewritten to + propagate beta's uncertainty too. + + Returns + ------- + str + The dependency expression, over the mapped names ``b``, ``hbar`` and ``tau``. + """ + horner = repr(_HWHM_POLY_COEFFS[-1]) + for coefficient in _HWHM_POLY_COEFFS[-2::-1]: + horner = f'({horner}) * b.value + {coefficient!r}' + return f'exp((log(b.value) + ({horner})) / b.value) * hbar / tau' + + +class StretchedExponential(CreateParametersMixin, ModelComponent): + r""" + Model of a stretched exponential (Kohlrausch-Williams-Watts) relaxation, Fourier transformed + from time to energy. + + The model is defined by its intermediate scattering function + + $$ I(t) = A \exp\left[-\left(\frac{|t|}{\tau}\right)^\beta\right] $$ + + where $\tau$ is the relaxation time and $\beta$ the stretching exponent. Here we calculate the + Fourier transform: + + $$ S(x) = \frac{1}{2\pi\hbar} \int I(t)\, e^{-i (x - x_0) t / \hbar} \, \mathrm{d}t = + \frac{A}{\pi \Gamma} \, G_\beta\!\left(\frac{x - x_0}{\Gamma}\right), \qquad \Gamma = + \frac{\hbar}{\tau} $$ + + with $G_\beta(w) = \int_0^\infty e^{-u^\beta} \cos(w u)\, \mathrm{d}u$. $A$ is the area (the + profile integrates to $A$ over $x$), $x_0$ is the center, and $\Gamma$ is the energy scale set + by the relaxation time. area has unit = x_unit * y_unit; center has unit = x_unit; + relaxation_time has unit ps; beta is dimensionless. + + $\beta = 1$ recovers a Lorentzian of HWHM $\Gamma$ and $\beta = 2$ a Gaussian of standard + deviation $\sqrt{2}\,\Gamma$; in between there is no closed form and the transform is evaluated + numerically. $\beta \le 1$ is the physically usual range. $\beta$ is capped at 2 because + $\exp(-|t|^\beta)$ stops being positive definite beyond it, so the transform would go negative. + + Note that the x-axis has to be an energy, since the relaxation time is turned into an energy + scale through $\hbar$. + + Examples + -------- + **Creating a stretched exponential** + + By default the center is fixed at 0 like a Lorentzian:: + ```python + import numpy as np + import easydynamics as edyn + + kww = edyn.StretchedExponential(area=1.0, relaxation_time=5.0, beta=0.7) + x = np.linspace(-2, 2, 100) + values = kww.evaluate(x) + ``` + + **Modifying parameters after construction** + + ```python + import easydynamics as edyn + + kww = edyn.StretchedExponential(area=2.0, relaxation_time=10.0, beta=0.5, name='Polymer') + kww.relaxation_time = 20.0 + kww.beta = 0.6 + ``` + """ + + def __init__( + self, + area: Numeric = 1.0, + center: Numeric | None = None, + relaxation_time: Numeric = 1.0, + beta: Numeric = 1.0, + x_unit: str | sc.Unit = 'meV', + y_unit: str | sc.Unit = 'dimensionless', + name: str = 'StretchedExponential', + display_name: str | None = None, + unique_name: str | None = None, + ) -> None: + r""" + Initialize the StretchedExponential component. + + Parameters + ---------- + area : Numeric, default=1.0 + Integrated area under the transformed profile. Unit is ``x_unit * y_unit``. + center : Numeric | None, default=None + Peak position in x_unit. If None, defaults to 0 and the center parameter is fixed. + relaxation_time : Numeric, default=1.0 + Relaxation time tau in ps. Must be strictly positive. It enters the spectrum as the + energy scale $\Gamma = \hbar / \tau$. + beta : Numeric, default=1.0 + Stretching exponent. Must lie in ``[0.05, 2.0]``; $\beta = 1$ is a Lorentzian and + $\beta = 2$ a Gaussian. + x_unit : str | sc.Unit, default='meV' + Unit of the x-axis. Must be an energy, since $\hbar / \tau$ is converted into it. + center is stored in this unit. area_unit = x_unit * y_unit. + y_unit : str | sc.Unit, default='dimensionless' + Unit of the y-axis (output). + name : str, default='StretchedExponential' + Name of the component. + display_name : str | None, default=None + Display name shown when plotting. Falls back to *name* if None. + unique_name : str | None, default=None + Globally unique identifier. Auto-generated if None. + """ + super().__init__( + x_unit=x_unit, + y_unit=y_unit, + name=name, + display_name=display_name, + unique_name=unique_name, + ) + + self._area = self._create_area_parameter( + area=area, name=name, x_unit=self.x_unit, y_unit=self.y_unit + ) + self._center = self._create_center_parameter( + center=center, name=name, fix_if_none=True, x_unit=self.x_unit + ) + + self._validate_relaxation_time(relaxation_time) + self._relaxation_time = Parameter( + name=name + ' relaxation_time', + value=float(relaxation_time), + unit='ps', + min=MINIMUM_RELAXATION_TIME, + ) + + self._validate_beta(beta) + self._beta = Parameter( + name=name + ' beta', + value=float(beta), + unit='dimensionless', + min=MINIMUM_BETA, + max=MAXIMUM_BETA, + ) + + # hbar has to be in the dependency map, because carrying J*s through the expression is + # what lets the unit algebra turn the reduced half width into an energy. It is a private + # copy rather than the shared ``utils.hbar`` deliberately: make_dependent_on attaches the + # width as an observer of every mapped variable, so mapping the module-level constant + # would append one observer per component to a list that is never emptied, keeping every + # StretchedExponential ever built alive and lengthening each of its notifications. + self._hbar = DescriptorNumber(name=name + ' hbar', value=hbar.value, unit=str(hbar.unit)) + # Built here rather than on first access so that get_all_parameters reports the same set + # whether or not the width has been read yet. A non-energy x_unit has no width to build + # and must still construct: the UnitError is deferred to the first read, as for evaluate. + self._width: Parameter | None = None + with contextlib.suppress(UnitError): + self._width = self._build_width() + + ################################ + # Validation + ################################ + + @staticmethod + def _validate_relaxation_time(value: Numeric) -> None: + """ + Check that a relaxation time is a finite, strictly positive number. + + Parameters + ---------- + value : Numeric + The candidate relaxation time. + + Raises + ------ + TypeError + If *value* is not a numeric type. + ValueError + If *value* is not finite, or is smaller than ``MINIMUM_RELAXATION_TIME``. + """ + if not isinstance(value, Numeric): + raise TypeError('relaxation_time must be a number.') + if not np.isfinite(value): + raise ValueError('relaxation_time must be a finite number.') + if float(value) < MINIMUM_RELAXATION_TIME: + raise ValueError('relaxation_time must be greater than zero.') + + @staticmethod + def _validate_beta(value: Numeric) -> None: + """ + Check that a stretching exponent is a finite number inside the supported range. + + Parameters + ---------- + value : Numeric + The candidate stretching exponent. + + Raises + ------ + TypeError + If *value* is not a numeric type. + ValueError + If *value* is not finite, or lies outside ``[MINIMUM_BETA, MAXIMUM_BETA]``. + """ + if not isinstance(value, Numeric): + raise TypeError('beta must be a number.') + if not np.isfinite(value): + raise ValueError('beta must be a finite number.') + if not MINIMUM_BETA <= float(value) <= MAXIMUM_BETA: + raise ValueError(f'beta must be between {MINIMUM_BETA} and {MAXIMUM_BETA}.') + + ################################ + # Properties + ################################ + + @property + def area(self) -> Parameter: + """ + Get the area parameter. + + Returns + ------- + Parameter + The area Parameter with unit ``x_unit * y_unit``. + """ + return self._area + + @area.setter + def area(self, value: Numeric) -> None: + """ + Parameters + ---------- + value : Numeric + New area value (in current area unit = x_unit * y_unit). + + Notes + ----- + A ``TypeError`` propagates from the shared value setter if *value* is not a numeric type, + and a ``ValueError`` propagates from it if *value* violates the area parameter's bounds + (e.g. a negative value when the area was created non-negative, giving it ``min=0``). + """ + self._set_bounded_parameter_value(self._area, value, 'area') + + @property + def center(self) -> Parameter: + """ + Get the center parameter. + + Returns + ------- + Parameter + The center (x_0) Parameter with unit ``x_unit``. + """ + return self._center + + @center.setter + def center(self, value: Numeric | None) -> None: + """ + Parameters + ---------- + value : Numeric | None + New center value in x_unit. If None, the center is set to 0 and the parameter is + fixed. + + Raises + ------ + TypeError + If *value* is not None and not a numeric type. + """ + if value is None: + value = 0.0 + self._center.fixed = True + if not isinstance(value, Numeric): + raise TypeError('center must be a number') + self._center.value = value + + @property + def relaxation_time(self) -> Parameter: + """ + Get the relaxation time parameter. + + Returns + ------- + Parameter + The relaxation time (tau) Parameter with unit ``ps``. + """ + return self._relaxation_time + + @relaxation_time.setter + def relaxation_time(self, value: Numeric) -> None: + """ + Parameters + ---------- + value : Numeric + New relaxation time in the parameter's current unit. Must be strictly positive. + + Notes + ----- + A ``TypeError`` propagates from the validator if *value* is not a numeric type, and a + ``ValueError`` propagates from it if *value* is not finite or not positive, or from the + shared value setter if *value* violates the parameter's bounds. + """ + self._validate_relaxation_time(value) + self._set_bounded_parameter_value(self._relaxation_time, value, 'relaxation_time') + + @property + def beta(self) -> Parameter: + """ + Get the stretching exponent parameter. + + Returns + ------- + Parameter + The stretching exponent (beta) Parameter, dimensionless. + """ + return self._beta + + @beta.setter + def beta(self, value: Numeric) -> None: + """ + Parameters + ---------- + value : Numeric + New stretching exponent. Must lie in ``[MINIMUM_BETA, MAXIMUM_BETA]``. + + Notes + ----- + A ``TypeError`` propagates from the validator if *value* is not a numeric type, and a + ``ValueError`` propagates from it if *value* is not finite or lies outside the supported + range, or from the shared value setter if *value* violates the parameter's bounds. + """ + self._validate_beta(value) + self._set_bounded_parameter_value(self._beta, value, 'beta') + + def _build_width(self) -> Parameter: + r""" + Build the dependent width Parameter in the component's current x_unit. + + Returns + ------- + Parameter + A Parameter made dependent on :attr:`beta` and :attr:`relaxation_time`. + + Raises + ------ + UnitError + If x_unit is not an energy, so $\hbar / \tau$ cannot be expressed in it. + """ + width = Parameter(name=self.name + ' width', value=1.0, unit=str(self.x_unit)) + try: + # easyscience propagates inf bounds through arithmetic, producing inf/inf=nan as a + # transient intermediate; suppress the spurious numpy warnings as elsewhere. + with np.errstate(invalid='ignore', divide='ignore'): + width.make_dependent_on( + dependency_expression=_hwhm_dependency_expression(), + dependency_map={ + 'b': self._beta, + 'hbar': self._hbar, + 'tau': self._relaxation_time, + }, + desired_unit=str(self.x_unit), + ) + except UnitError as e: + raise UnitError( + f'{self.__class__.__name__} needs an energy x_unit so that hbar / ' + f'relaxation_time can be expressed in it, but got {self.x_unit}.' + ) from e + return width + + @property + def width(self) -> Parameter: + r""" + Get the half width at half maximum of the peak. + + This is a *dependent* Parameter, derived from :attr:`relaxation_time` and :attr:`beta`. It + is accurate to 1.4e-5 relative for $\beta \ge 0.1$. For better accuracy, use + :func:`_reduced_hwhm` directly and multiply by $\hbar / \tau$ in the desired unit. Note + that the uncertainty of the width may be underestimated, because the dependency expression + uses ``b.value`` rather than ``b`` to compute the width, so the uncertainty of beta is not + propagated into the width. + + Returns + ------- + Parameter + The HWHM expressed in the component's own x_unit. + + Notes + ----- + A ``UnitError`` propagates from :meth:`_build_width` if x_unit is not an energy, so $\hbar + / \tau$ cannot be expressed in it. + """ + if self._width is None: + self._width = self._build_width() + return self._width + + ################################ + # Evaluation + ################################ + + def _energy_scale(self, eval_unit: str | None) -> float: + r""" + Get the energy scale $\Gamma = \hbar / \tau$ expressed in the evaluation unit. + + Parameters + ---------- + eval_unit : str | None + The unit x values are expressed in, or None for the component's own x_unit. + + Returns + ------- + float + Gamma in *eval_unit*. + + Raises + ------ + UnitError + If the target unit is not an energy, so $\hbar / \tau$ cannot be expressed in it. + """ + target_unit = eval_unit if eval_unit is not None else self.x_unit + tau_unit = self._relaxation_time.unit + try: + hbar_value = convert_value_unit(hbar.value, hbar.unit, f'({target_unit})*({tau_unit})') + except UnitError as e: + raise UnitError( + f'{self.__class__.__name__} needs an energy x_unit so that hbar / ' + f'relaxation_time can be expressed in it, but got {target_unit}.' + ) from e + return hbar_value / self._relaxation_time.value + + def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndarray: + r""" + Evaluate the transformed stretched exponential at x_vals. + + $$ S(x) = \frac{A}{\pi \Gamma} \, G_\beta\!\left(\frac{x - x_0}{\Gamma}\right), \qquad + \Gamma = \frac{\hbar}{\tau} $$ + + where *A* is ``area``, *x*₀ is ``center``, *tau* is ``relaxation_time``, *beta* is the + stretching exponent, and $G_\beta$ is the reduced transform computed by :func:`_kww_shape`. + Parameters in the model's own units are temporarily converted to eval_unit for the + computation. + + Parameters + ---------- + x_vals : np.ndarray + Raw x values expressed in eval_unit. + eval_unit : str | None + The unit of x_vals. + + Returns + ------- + np.ndarray + Evaluated values at x_vals. + """ + center = self._resolve_param_value(self._center, eval_unit) + area = self._resolve_param_value(self._area, self._eval_area_unit(eval_unit)) + energy_scale = self._energy_scale(eval_unit) + + shape = _kww_shape((x_vals - center) / energy_scale, self._beta.value) + return area / (np.pi * energy_scale) * shape + + ################################ + # Unit conversion + ################################ + + def convert_x_unit(self, new_x_unit: str | sc.Unit) -> None: + r""" + Convert the center and the area to new_x_unit. + + The relaxation time carries a time unit and beta is dimensionless, so neither is affected; + the energy scale $\hbar / \tau$ is re-derived in the new unit on every evaluation. + + Parameters + ---------- + new_x_unit : str | sc.Unit + Target x-axis unit. Must be dimensionally compatible with the current x_unit. + """ + self._convert_x_unit_area_based( + new_x_unit=new_x_unit, + x_params=[self._center], + area_param=self._area, + ) + # The dependency resolves into the unit it was built with, so rebuild it in the new one. + self._width = self._build_width() + + def convert_y_unit(self, new_y_unit: str | sc.Unit) -> None: + """ + Convert the y-axis (output) unit by rescaling the area parameter. + + The area is rescaled from ``x_unit * old_y_unit`` to ``x_unit * new_y_unit``. + + Parameters + ---------- + new_y_unit : str | sc.Unit + Target y-axis unit. + """ + self._convert_y_unit_area_based(new_y_unit=new_y_unit, area_param=self._area) + + def __repr__(self) -> str: + """ + Return a string representation of the StretchedExponential. + + Returns + ------- + str + A string representation of the StretchedExponential. + """ + return ( + f'{self.__class__.__name__}(name = {self.name}, display_name = {self.display_name}, ' + f'x_unit = {self.x_unit}, y_unit = {self.y_unit},\n ' + f' area = {self.area},\n ' + f' center = {self.center},\n ' + f' relaxation_time = {self.relaxation_time},\n ' + f' beta = {self.beta})' + ) diff --git a/src/easydynamics/utils/plotting.py b/src/easydynamics/utils/plotting.py index cb8ee0a98..52360b518 100644 --- a/src/easydynamics/utils/plotting.py +++ b/src/easydynamics/utils/plotting.py @@ -4,7 +4,7 @@ import plopp as pp import scipp as sc -from plopp.backends.matplotlib.figure import InteractiveFigure +from plopp.backends.matplotlib.figure import WidgetFigure from plopp.plotting._slicer import SlicerPlot from plopp.plotting._slicer import _maybe_reduce_dim from plopp.widgets import slice_dims @@ -17,7 +17,7 @@ def slicerplot_with_residuals( keep: list[str] | str | None = None, operation: str = 'sum', **kwargs: object, -) -> InteractiveFigure: +) -> WidgetFigure: """ Create a SlicerPlot with an additional subplot for residuals. @@ -56,7 +56,7 @@ def slicerplot_with_residuals( Returns ------- - InteractiveFigure + WidgetFigure A figure containing the SlicerPlot and the residuals subplot. Raises diff --git a/src/easydynamics/utils/posterior_plotting.py b/src/easydynamics/utils/posterior_plotting.py index d1e45b146..dc930de5a 100644 --- a/src/easydynamics/utils/posterior_plotting.py +++ b/src/easydynamics/utils/posterior_plotting.py @@ -23,7 +23,7 @@ if TYPE_CHECKING: from ipywidgets import VBox from matplotlib.figure import Figure - from plopp.backends.matplotlib.figure import InteractiveFigure + from plopp.backends.matplotlib.figure import WidgetFigure def plot_trace( @@ -719,7 +719,7 @@ def predictive_with_slider( title: str | None = None, credible_interval: float = 68.0, **kwargs: dict[str, Any], -) -> InteractiveFigure: +) -> WidgetFigure: """ Plot per-Q posterior-predictive bands behind a plopp Q slider. @@ -763,7 +763,7 @@ def predictive_with_slider( Returns ------- - InteractiveFigure + WidgetFigure The plopp figure with its Q slider. Raises diff --git a/tests/integration/fitting/test_fitting_with_stretched_exponential.py b/tests/integration/fitting/test_fitting_with_stretched_exponential.py new file mode 100644 index 000000000..fbbac97dc --- /dev/null +++ b/tests/integration/fitting/test_fitting_with_stretched_exponential.py @@ -0,0 +1,112 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +""" +Round-trip fit of a StretchedExponential through Analysis1d. + +The transform of a stretched exponential has no closed form, so it is evaluated by quadrature. +That makes it worth checking not just that the numbers are right at fixed parameters — the unit +tests cover that — but that the profile stays smooth enough in the relaxation time and the +stretching exponent for a least-squares fitter to walk back to the truth from a poor start. +""" + +import numpy as np +import pytest +import scipp as sc + +from easydynamics.analysis.analysis1d import Analysis1d +from easydynamics.experiment import Experiment +from easydynamics.sample_model import InstrumentModel +from easydynamics.sample_model import SampleModel +from easydynamics.sample_model import StretchedExponential + +TRUE_AREA = 1.0 +TRUE_RELAXATION_TIME = 8.0 +TRUE_BETA = 0.6 +NOISE_FRACTION = 0.01 + +TRUTHS = { + 'StretchedExponential area': TRUE_AREA, + 'StretchedExponential relaxation_time': TRUE_RELAXATION_TIME, + 'StretchedExponential beta': TRUE_BETA, +} + + +def build_analysis() -> Analysis1d: + """Fit a badly-started StretchedExponential against noisy data drawn from the true one.""" + energy_values = np.linspace(-1.5, 1.5, 301) + truth = StretchedExponential( + area=TRUE_AREA, + relaxation_time=TRUE_RELAXATION_TIME, + beta=TRUE_BETA, + ) + profile = truth.evaluate(energy_values) + noise = NOISE_FRACTION * profile.max() + observed = profile + np.random.default_rng(0).normal(0.0, noise, size=profile.shape) + + experiment = Experiment( + data=sc.DataArray( + data=sc.array( + dims=['Q', 'energy'], + values=observed[None, :], + variances=np.full_like(observed, noise**2)[None, :], + ), + coords={ + 'Q': sc.array(dims=['Q'], values=[1.0], unit='1/Angstrom'), + 'energy': sc.array(dims=['energy'], values=energy_values, unit='meV'), + }, + ) + ) + + # Start far from the truth: half the relaxation time, an almost unstretched exponent and half + # again the area. The centre is left at its default, fixed at zero. + analysis = Analysis1d( + display_name='StretchedExponentialIntegration', + experiment=experiment, + sample_model=SampleModel( + components=StretchedExponential(area=1.5, relaxation_time=4.0, beta=0.9) + ), + instrument_model=InstrumentModel(), + Q_index=0, + ) + # The energy offset shifts the spectrum exactly as the component's centre does, so leaving + # both free would make the model unidentifiable. + analysis.instrument_model.fix_energy_offset(Q_index=0) + return analysis + + +@pytest.fixture(scope='module') +def fitted_analysis(): + analysis = build_analysis() + analysis.fit() + return analysis + + +class TestFittingWithStretchedExponential: + def test_fit_describes_the_data(self): + # WHEN + analysis = build_analysis() + + # THEN + results = analysis.fit() + + # EXPECT + assert results.success + assert results.reduced_chi2 < 1.5 + + def test_the_free_parameters_are_the_three_shape_parameters(self, fitted_analysis): + # THEN + names = {parameter.name for parameter in fitted_analysis.get_free_parameters()} + + # EXPECT the centre and the energy offset stay fixed + assert names == set(TRUTHS) + + @pytest.mark.parametrize('name', list(TRUTHS)) + def test_fit_recovers_the_true_parameters(self, fitted_analysis, name): + # WHEN + parameter = next(p for p in fitted_analysis.get_free_parameters() if p.name == name) + + # THEN EXPECT the truth within a few standard errors, on a usable uncertainty + assert np.isfinite(parameter.error) + assert parameter.error > 0.0 + assert abs(parameter.value - TRUTHS[name]) < 4 * parameter.error diff --git a/tests/unit/easydynamics/convolution/test_numerical_convolution_base.py b/tests/unit/easydynamics/convolution/test_numerical_convolution_base.py index 3d80f983b..4ce8db631 100644 --- a/tests/unit/easydynamics/convolution/test_numerical_convolution_base.py +++ b/tests/unit/easydynamics/convolution/test_numerical_convolution_base.py @@ -9,6 +9,7 @@ from easydynamics.convolution.energy_grid import EnergyGrid from easydynamics.convolution.numerical_convolution_base import NumericalConvolutionBase from easydynamics.sample_model import Gaussian +from easydynamics.sample_model import StretchedExponential from easydynamics.sample_model import Voigt from easydynamics.sample_model.component_collection import ComponentCollection from easydynamics.settings.convolution_settings import ConvolutionSettings @@ -605,6 +606,30 @@ def test_check_width_small_threshold(self, default_numerical_convolution_base): model_name='ComponentCollection', ) + def test_check_width_small_threshold_for_stretched_exponential( + self, default_numerical_convolution_base + ): + """ + Regression: a stretched exponential at small beta is far narrower than its energy scale + hbar / tau, and used to report that scale as its width. The spike then fell between the + grid points with no warning at all. + """ + # WHEN tau = 5 ps puts hbar / tau at 0.13 meV, well above the grid step, while the true + # half width at beta = 0.2 is 3.5e-5 meV, well below it + narrow_kww = StretchedExponential( + name='ComponentCollection', area=1.0, relaxation_time=5.0, beta=0.2 + ) + + # THEN EXPECT + with pytest.warns( + UserWarning, + match='Increase upsample_factor to improve', + ): + default_numerical_convolution_base._check_width_thresholds( + model=narrow_kww, + model_name='ComponentCollection', + ) + def test_check_width_no_warnings(self, default_numerical_convolution_base): """ Test that _check_width_thresholds does not warn when model diff --git a/tests/unit/easydynamics/sample_model/components/test_stretched_exponential.py b/tests/unit/easydynamics/sample_model/components/test_stretched_exponential.py new file mode 100644 index 000000000..a5943c865 --- /dev/null +++ b/tests/unit/easydynamics/sample_model/components/test_stretched_exponential.py @@ -0,0 +1,693 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +from copy import copy +from itertools import pairwise + +import numpy as np +import pytest +import scipp as sc +from easyscience.variable import Parameter +from scipp import UnitError +from scipy.integrate import simpson +from scipy.special import gamma + +from easydynamics.sample_model import Gaussian +from easydynamics.sample_model import Lorentzian +from easydynamics.sample_model import StretchedExponential +from easydynamics.sample_model.components.stretched_exponential import _kww_shape +from easydynamics.sample_model.components.stretched_exponential import _reduced_hwhm + +# hbar in meV*ps (CODATA), so the tests derive the energy scale independently of the library. +HBAR_MEV_PS = 0.6582119569509066 + + +def sinh_grid(reach: float, n_points: int = 20001) -> np.ndarray: + """A symmetric grid that is dense near zero and reaches *reach*, for the heavy KWW tails.""" + t = np.linspace(-np.arcsinh(reach), np.arcsinh(reach), n_points) + return np.sinh(t) + + +def naive_fft_transform( + area: float, tau: float, beta: float, n_time: int, t_max: float +) -> tuple[np.ndarray, np.ndarray]: + """ + Transform ``I(t) = area * exp(-(t / tau)**beta)`` to energy the obvious way, with an FFT. + + This is the textbook route the implementation deliberately does not take: sample the + relaxation on a uniform time grid and let ``np.fft`` do the cosine transform. It is an + independent check because it shares no code and no idea with the rotated-contour quadrature + in ``_kww_shape`` -- only the physics. + + The profile is even in *t*, so the two-sided transform is twice the one-sided one and + + S(x) = (1 / (pi hbar)) Re [ integral of I(t) exp(-i x t / hbar) dt over t >= 0 ], + + which on the uniform grid is a trapezoid sum (hence the half-weighted t = 0 sample) evaluated + at the FFT frequencies ``x_n = 2 pi hbar n / (n_time dt)``. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + The energy axis in meV and the spectrum on it. + """ + dt = t_max / n_time + time = np.arange(n_time) * dt + intensity = area * np.exp(-((time / tau) ** beta)) + spectrum = (np.real(np.fft.rfft(intensity)) - 0.5 * intensity[0]) * dt / (np.pi * HBAR_MEV_PS) + energy = 2.0 * np.pi * HBAR_MEV_PS * np.arange(spectrum.size) / (n_time * dt) + return energy, spectrum + + +def fft_comparison(beta: float, n_time: int, tau: float = 5.0, area: float = 1.0) -> float: + """Largest relative gap between the component and the FFT, over the resolvable region.""" + # Reach far enough in time that exp(-(t / tau)**beta) has fallen by exp(-40). + energy, reference = naive_fft_transform(area, tau, beta, n_time, tau * 40.0 ** (1.0 / beta)) + inside = (energy > 0.0) & (energy < 2.0) + energy, reference = energy[inside], reference[inside] + + stretched = StretchedExponential(area=area, relaxation_time=tau, beta=beta) + values = stretched.evaluate(energy) + + # Below a thousandth of the peak the FFT reference is dominated by its own truncation error, + # so comparing there would measure the reference rather than the implementation. + resolvable = values > 1e-3 * stretched.evaluate(np.array([0.0]))[0] + return float(np.max(np.abs(values[resolvable] - reference[resolvable]) / values[resolvable])) + + +##################################### +# The Fourier transform: _kww_shape +##################################### + + +@pytest.mark.parametrize( + 'beta, rel', + # The smallest supported beta spreads exp(-u**beta) over so many decades that the grid only + # just reaches its tail, which costs a few digits; everything above it is at machine precision. + [(0.05, 1e-8), (0.1, 1e-11), (0.3, 1e-12), (0.5, 1e-12), (1.0, 1e-12), (2.0, 1e-12)], +) +def test_kww_shape_at_zero(beta, rel): + # WHEN G(0) is the integral of exp(-u**beta), which is gamma(1 + 1/beta) + + # THEN + value = _kww_shape(np.array([0.0]), beta) + + # EXPECT + assert value[0] == pytest.approx(gamma(1.0 + 1.0 / beta), rel=rel) + + +def test_kww_shape_matches_lorentzian_at_beta_one(): + # WHEN beta = 1 the transform of exp(-u) is known in closed form + w = sinh_grid(1e6) + + # THEN + value = _kww_shape(w, 1.0) + + # EXPECT + np.testing.assert_allclose(value, 1.0 / (1.0 + w**2), rtol=1e-9) + + +def test_kww_shape_matches_gaussian_at_beta_two(): + # WHEN beta = 2 the transform of exp(-u**2) is again a Gaussian + w = np.linspace(-12.0, 12.0, 501) + + # THEN + value = _kww_shape(w, 2.0) + + # EXPECT + expected = 0.5 * np.sqrt(np.pi) * np.exp(-(w**2) / 4) + np.testing.assert_allclose(value, expected, atol=1e-14) + + +@pytest.mark.parametrize('beta', [0.3, 0.6, 0.9]) +def test_kww_shape_matches_the_large_w_asymptote(beta): + # WHEN the leading term of the large-w expansion is + # gamma(beta + 1) sin(pi beta / 2) / w**(beta + 1). The next term is smaller by w**-beta, so + # w has to be taken far out before the leading term alone is worth a few digits. + w = np.array([1e9, 1e11]) + + # THEN + value = _kww_shape(w, beta) + + # EXPECT + expected = gamma(beta + 1) * np.sin(np.pi * beta / 2) / w ** (beta + 1) + np.testing.assert_allclose(value, expected, rtol=3e-3) + + +@pytest.mark.parametrize('beta', [0.3, 0.5, 0.8, 1.0, 1.6, 2.0]) +def test_kww_shape_is_even_and_non_negative(beta): + # WHEN + w = sinh_grid(1e6, n_points=2001) + + # THEN + value = _kww_shape(w, beta) + + # EXPECT + np.testing.assert_allclose(value, _kww_shape(-w, beta), rtol=0, atol=0) + assert np.all(value >= 0.0) + + +@pytest.mark.parametrize('beta', [0.3, 0.5, 0.8, 1.0, 1.6, 2.0]) +def test_kww_shape_integrates_to_pi(beta): + # WHEN the transform of a function that is 1 at t = 0 integrates to pi over w + + # THEN + w = sinh_grid(1e8) + integral = simpson(_kww_shape(w, beta), x=w) + + # EXPECT + assert integral == pytest.approx(np.pi, rel=5e-3) + + +def test_kww_shape_is_blocked_without_changing_the_result(): + # WHEN a long energy axis is evaluated in blocks to bound the memory of the quadrature + w = np.linspace(-40.0, 40.0, 997) + reference = _kww_shape(w, 0.6) + + # THEN + blocked = np.concatenate([_kww_shape(chunk, 0.6) for chunk in np.array_split(w, 7)]) + + # EXPECT the block boundaries do not perturb any value + np.testing.assert_array_equal(blocked, reference) + + +class TestStretchedExponential: + @pytest.fixture + def stretched_exponential(self): + return StretchedExponential( + name='StretchedName', + display_name='TestStretched', + area=2.0, + center=0.5, + relaxation_time=3.0, + beta=0.7, + x_unit='meV', + ) + + ############# + # Creation + ############# + + def test_init_no_inputs(self): + # WHEN THEN + stretched = StretchedExponential() + + # EXPECT + assert stretched.display_name == 'StretchedExponential' + assert stretched.area.value == pytest.approx(1.0) + assert stretched.center.value == pytest.approx(0.0) + assert stretched.relaxation_time.value == pytest.approx(1.0) + assert stretched.beta.value == pytest.approx(1.0) + assert stretched.x_unit == 'meV' + assert stretched.y_unit == 'dimensionless' + assert stretched.relaxation_time.unit == 'ps' + assert stretched.beta.unit == 'dimensionless' + assert stretched.center.fixed is True + + def test_initialization(self, stretched_exponential: StretchedExponential): + # WHEN THEN EXPECT + assert stretched_exponential.display_name == 'TestStretched' + assert stretched_exponential.area.value == pytest.approx(2.0) + assert stretched_exponential.center.value == pytest.approx(0.5) + assert stretched_exponential.relaxation_time.value == pytest.approx(3.0) + assert stretched_exponential.beta.value == pytest.approx(0.7) + assert stretched_exponential.center.fixed is False + + @pytest.mark.parametrize( + 'kwargs, expected_message', + [ + ({'area': 'invalid'}, 'area must be a number'), + ({'center': 'invalid'}, 'center must be None or a number'), + ({'relaxation_time': 'invalid'}, 'relaxation_time must be a number'), + ({'beta': 'invalid'}, 'beta must be a number'), + ({'x_unit': 123}, 'unit must be None, a string'), + ({'y_unit': 123}, 'unit must be None, a string'), + ], + ) + def test_input_type_validation_raises(self, kwargs, expected_message): + # WHEN THEN EXPECT + with pytest.raises(TypeError, match=expected_message): + StretchedExponential(**kwargs) + + @pytest.mark.parametrize( + 'kwargs, expected_message', + [ + ({'relaxation_time': 0.0}, 'relaxation_time must be greater than zero'), + ({'relaxation_time': -1.0}, 'relaxation_time must be greater than zero'), + ({'relaxation_time': np.inf}, 'relaxation_time must be a finite number'), + ({'beta': 0.0}, 'beta must be between'), + ({'beta': 2.5}, 'beta must be between'), + ({'beta': np.nan}, 'beta must be a finite number'), + ], + ) + def test_input_value_validation_raises(self, kwargs, expected_message): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match=expected_message): + StretchedExponential(**kwargs) + + def test_negative_area_warns(self): + # WHEN THEN EXPECT + with pytest.warns(UserWarning, match='may not be physically meaningful'): + StretchedExponential(area=-2.0) + + def test_get_all_parameters(self, stretched_exponential: StretchedExponential): + # WHEN THEN + params = stretched_exponential.get_all_parameters() + + # EXPECT + assert all(isinstance(param, Parameter) for param in params) + assert {param.name for param in params} == { + 'StretchedName area', + 'StretchedName center', + 'StretchedName relaxation_time', + 'StretchedName beta', + 'StretchedName width', + } + + def test_copy(self, stretched_exponential: StretchedExponential): + # WHEN THEN + stretched_copy = copy(stretched_exponential) + + # EXPECT + assert stretched_copy is not stretched_exponential + assert stretched_copy.display_name == stretched_exponential.display_name + assert stretched_copy.area.value == stretched_exponential.area.value + assert stretched_copy.center.value == stretched_exponential.center.value + assert stretched_copy.relaxation_time.value == stretched_exponential.relaxation_time.value + assert stretched_copy.beta.value == stretched_exponential.beta.value + assert stretched_copy.x_unit == stretched_exponential.x_unit + + def test_repr(self, stretched_exponential: StretchedExponential): + # WHEN THEN + repr_str = repr(stretched_exponential) + + # EXPECT + assert 'StretchedExponential' in repr_str + assert 'name = StretchedName' in repr_str + assert 'x_unit = meV' in repr_str + assert 'area =' in repr_str + assert 'center =' in repr_str + assert 'relaxation_time =' in repr_str + assert 'beta =' in repr_str + + ############# + # Parameters + ############# + + @pytest.mark.parametrize( + 'prop, valid_value', + [('area', 3.0), ('center', 0.6), ('relaxation_time', 4.0), ('beta', 0.9)], + ) + def test_property_setters( + self, stretched_exponential: StretchedExponential, prop, valid_value + ): + # WHEN: set a valid value + setattr(stretched_exponential, prop, valid_value) + # THEN EXPECT + assert getattr(stretched_exponential, prop).value == valid_value + + # WHEN: set an invalid value — THEN EXPECT + with pytest.raises(TypeError, match=' must be a number'): + setattr(stretched_exponential, prop, 'invalid') + + def test_relaxation_time_must_be_positive(self, stretched_exponential: StretchedExponential): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match='relaxation_time must be greater than zero'): + stretched_exponential.relaxation_time = -1.0 + assert stretched_exponential.relaxation_time.value == pytest.approx(3.0) + + @pytest.mark.parametrize('value', [0.01, 2.5]) + def test_beta_outside_the_supported_range_raises( + self, stretched_exponential: StretchedExponential, value + ): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match='beta must be between'): + stretched_exponential.beta = value + assert stretched_exponential.beta.value == pytest.approx(0.7) + + def test_area_setter_out_of_bounds_raises(self, stretched_exponential: StretchedExponential): + # WHEN the fixture's area was created non-negative, so it carries min=0 + + # THEN EXPECT a negative assignment raises instead of being silently clamped to 0 + with pytest.raises(ValueError, match='violates the parameter bounds'): + stretched_exponential.area = -1.0 + assert stretched_exponential.area.value == pytest.approx(2.0) + + def test_center_is_fixed_if_set_to_None(self, stretched_exponential: StretchedExponential): + # WHEN + assert stretched_exponential.center.fixed is False + + # THEN + stretched_exponential.center = None + + # EXPECT + assert stretched_exponential.center.value == pytest.approx(0.0) + assert stretched_exponential.center.fixed is True + + def test_width_is_the_half_width_at_half_maximum( + self, stretched_exponential: StretchedExponential + ): + # WHEN width is not stored but solved for as the half maximum crossing of the profile + # THEN + width = stretched_exponential.width + + # EXPECT the profile really has fallen to half its peak there, to the accuracy of the + # closed-form fit the dependency evaluates rather than to that of the root itself + center = stretched_exponential.center.value + peak = stretched_exponential.evaluate(np.array([center]))[0] + half = stretched_exponential.evaluate(np.array([center + width.value]))[0] + assert half == pytest.approx(0.5 * peak, rel=1e-4) + assert str(width.unit) == 'meV' + + def test_width_is_below_the_energy_scale_for_stretched_profiles(self): + # WHEN beta < 1 the peak sharpens well beyond Gamma = hbar / tau, so reporting Gamma would + # hide a sub-grid spike from the convolution's width-versus-grid checks + gamma = HBAR_MEV_PS / 5.0 + + # THEN + widths = { + beta: StretchedExponential(relaxation_time=5.0, beta=beta).width.value + for beta in (1.0, 0.5, 0.2) + } + + # EXPECT Gamma only at beta = 1, and orders of magnitude below it further down + assert widths[1.0] == pytest.approx(gamma, rel=2e-5) + assert widths[0.5] == pytest.approx(0.22355 * gamma, rel=1e-4) + assert widths[0.2] == pytest.approx(2.653e-4 * gamma, rel=1e-3) + + def test_width_is_a_dependent_parameter(self, stretched_exponential: StretchedExponential): + # WHEN width is resolved from beta and relaxation_time by a dependency expression + # THEN + width = stretched_exponential.width + + # EXPECT a Parameter, so the usual width attributes resolve, but a dependent one, so a fit + # can never pick it up as a free parameter + assert isinstance(width, Parameter) + assert width.independent is False + assert width not in stretched_exponential.get_fittable_parameters() + + def test_width_matches_the_solved_half_width(self): + # WHEN the dependency evaluates a closed-form fit rather than the root itself + # THEN EXPECT it tracks _reduced_hwhm over the range the fit targets + for beta in (2.0, 1.0, 0.7, 0.5, 0.3, 0.2, 0.1): + stretched = StretchedExponential(relaxation_time=5.0, beta=beta) + expected = _reduced_hwhm(beta) * HBAR_MEV_PS / 5.0 + assert stretched.width.value == pytest.approx(expected, rel=2e-5) + + def test_width_is_recomputed_after_a_raw_beta_write( + self, stretched_exponential: StretchedExponential + ): + # WHEN a minimizer writes the beta Parameter directly, bypassing the component setter + before = stretched_exponential.width.value + + # THEN + stretched_exponential.beta.value = 0.2 + + # EXPECT the dependency re-evaluates: beta is reached as b.value, but it is still a mapped + # variable, so a write to it retriggers the expression + assert stretched_exponential.width.value < before + + def test_width_can_be_chained_onto(self, stretched_exponential: StretchedExponential): + # WHEN another component's width is made to follow this one + lorentzian = Lorentzian(name='Chained', area=1.0, width=0.1) + with np.errstate(invalid='ignore', divide='ignore'): + lorentzian.width.make_dependent_on( + dependency_expression='w * 2', + dependency_map={'w': stretched_exponential.width}, + desired_unit='meV', + ) + + # THEN + stretched_exponential.beta.value = 0.2 + + # EXPECT the change reaches the chained parameter, not just this component + assert lorentzian.width.value == pytest.approx(2 * stretched_exponential.width.value) + + def test_width_tracks_beta(self, stretched_exponential: StretchedExponential): + # WHEN + before = stretched_exponential.width.value + + # THEN + stretched_exponential.beta = 0.3 + + # EXPECT a smaller beta gives a sharper peak, so a smaller half width + assert stretched_exponential.width.value < before + + def test_width_tracks_the_relaxation_time(self, stretched_exponential: StretchedExponential): + # WHEN + before = stretched_exponential.width.value + + # THEN + stretched_exponential.relaxation_time = 6.0 + + # EXPECT doubling the relaxation time halves the width + assert stretched_exponential.width.value == pytest.approx(before / 2.0) + + def test_width_is_read_only(self, stretched_exponential: StretchedExponential): + # WHEN THEN EXPECT the relaxation time is the fittable parameter, not the width + with pytest.raises(AttributeError): + stretched_exponential.width = 0.5 + + def test_width_follows_the_x_unit(self): + # WHEN the component measures energy in microeV + stretched = StretchedExponential(relaxation_time=5.0, beta=1.0, x_unit='ueV') + + # THEN + width = stretched.width + + # EXPECT the half width is expressed in that unit too + assert width.value == pytest.approx(1e3 * HBAR_MEV_PS / 5.0, rel=2e-5) + assert str(width.unit) == str(sc.Unit('ueV')) + + def test_width_follows_an_x_unit_conversion(self): + # WHEN + stretched = StretchedExponential(relaxation_time=5.0, beta=1.0, x_unit='meV') + before = stretched.width.value + + # THEN + stretched.convert_x_unit('ueV') + + # EXPECT the half width is re-derived in the new unit rather than keeping the old one + assert stretched.width.value == pytest.approx(1e3 * before) + assert str(stretched.width.unit) == str(sc.Unit('ueV')) + + def test_width_with_a_non_energy_x_unit_raises(self): + # WHEN THEN EXPECT hbar / relaxation_time cannot be expressed in metres + with pytest.raises(UnitError, match='needs an energy x_unit'): + _ = StretchedExponential(x_unit='m').width + + def test_width_is_the_lorentzian_hwhm_at_beta_one(self): + # WHEN beta = 1 the transform is a Lorentzian of HWHM Gamma = hbar / tau. Gamma is used + # directly rather than via width, so this pins the profile itself and stays independent of + # the closed-form fit the width dependency evaluates. + stretched = StretchedExponential(area=1.0, relaxation_time=4.0, beta=1.0) + + # THEN + lorentzian = Lorentzian(area=1.0, width=HBAR_MEV_PS / 4.0) + + # EXPECT + x = np.linspace(-2.0, 2.0, 101) + np.testing.assert_allclose(stretched.evaluate(x), lorentzian.evaluate(x), rtol=1e-12) + + ############# + # Evaluation + ############# + + def test_evaluate_reduces_to_a_lorentzian_at_beta_one(self): + # WHEN beta = 1, the transform is a Lorentzian of HWHM hbar / tau + stretched = StretchedExponential(area=2.5, relaxation_time=3.0, beta=1.0) + lorentzian = Lorentzian(area=2.5, width=HBAR_MEV_PS / 3.0) + x = np.linspace(-3.0, 3.0, 401) + + # THEN + result = stretched.evaluate(x) + + # EXPECT + np.testing.assert_allclose(result, lorentzian.evaluate(x), rtol=1e-10) + + def test_evaluate_reduces_to_a_gaussian_at_beta_two(self): + # WHEN beta = 2, the transform is a Gaussian of standard deviation sqrt(2) hbar / tau + stretched = StretchedExponential(area=2.5, relaxation_time=3.0, beta=2.0) + gaussian = Gaussian(area=2.5, width=np.sqrt(2) * HBAR_MEV_PS / 3.0) + x = np.linspace(-3.0, 3.0, 401) + + # THEN + result = stretched.evaluate(x) + + # EXPECT + np.testing.assert_allclose(result, gaussian.evaluate(x), atol=1e-13) + + def test_evaluate_peak_height(self, stretched_exponential: StretchedExponential): + # WHEN the peak sits at the center and equals area * gamma(1 + 1/beta) / (pi * Gamma) + energy_scale = HBAR_MEV_PS / 3.0 + + # THEN + result = stretched_exponential.evaluate(np.array([0.5])) + + # EXPECT + expected = 2.0 * gamma(1.0 + 1.0 / 0.7) / (np.pi * energy_scale) + assert result[0] == pytest.approx(expected, rel=1e-10) + + def test_evaluate_is_symmetric_about_the_center( + self, stretched_exponential: StretchedExponential + ): + # WHEN + offset = np.array([0.05, 0.3, 1.7]) + + # THEN + left = stretched_exponential.evaluate(0.5 - offset) + right = stretched_exponential.evaluate(0.5 + offset) + + # EXPECT + np.testing.assert_allclose(left, right, rtol=1e-12) + + @pytest.mark.parametrize('beta', [0.5, 0.7, 1.0, 1.5]) + def test_area_matches_parameter(self, beta): + # WHEN the transform integrates to the area parameter over the whole axis + stretched = StretchedExponential(area=2.0, relaxation_time=5.0, beta=beta) + # The tails are heavy, so integrate on a grid that is dense near the peak and reaches far + x = sinh_grid(1e8) * (HBAR_MEV_PS / 5.0) + + # THEN + numerical_area = simpson(stretched.evaluate(x), x=x) + + # EXPECT + assert numerical_area == pytest.approx(2.0, rel=5e-3) + + def test_relaxation_time_sets_the_width(self): + # WHEN a longer relaxation time means a narrower line + narrow = StretchedExponential(relaxation_time=20.0, beta=0.8) + wide = StretchedExponential(relaxation_time=2.0, beta=0.8) + + # THEN + narrow_peak = narrow.evaluate(np.array([0.0]))[0] + wide_peak = wide.evaluate(np.array([0.0]))[0] + + # EXPECT the peak scales as tau, since the area is conserved + assert narrow_peak == pytest.approx(10.0 * wide_peak, rel=1e-10) + + def test_evaluate_scipp_output(self, stretched_exponential: StretchedExponential): + # WHEN + x = np.linspace(-5, 5, 50) + + # THEN + result = stretched_exponential.evaluate(x, output='scipp') + + # EXPECT + assert isinstance(result, sc.Variable) + assert result.unit == sc.Unit('dimensionless') + np.testing.assert_allclose(result.values, stretched_exponential.evaluate(x)) + + def test_evaluate_with_scipp_x_in_another_energy_unit( + self, stretched_exponential: StretchedExponential + ): + # WHEN x carries microeV rather than the component's meV + x = np.linspace(-2.0, 2.0, 51) + + # THEN + result = stretched_exponential.evaluate( + sc.array(dims=['energy'], values=x * 1e3, unit='microeV') + ) + + # EXPECT the same profile, since the output carries y_unit either way + np.testing.assert_allclose(result, stretched_exponential.evaluate(x), rtol=1e-12) + + def test_evaluate_with_a_non_energy_x_unit_raises(self): + # WHEN hbar / relaxation_time cannot be expressed in the x unit + stretched = StretchedExponential(x_unit='m') + + # THEN EXPECT + with pytest.raises(UnitError, match='needs an energy x_unit'): + stretched.evaluate(np.array([0.0, 1.0])) + + @pytest.mark.parametrize('beta', [0.6, 0.8, 1.0, 1.5, 2.0]) + def test_evaluate_matches_a_naive_fft_transform(self, beta): + # WHEN the same physics is computed the obvious way, by FFT-ing the sampled relaxation + + # THEN + largest_gap = fft_comparison(beta, n_time=2**20) + + # EXPECT the two agree wherever the FFT itself is trustworthy + assert largest_gap < 1e-3 + + def test_the_naive_fft_converges_onto_evaluate(self): + # WHEN the FFT reference is limited by its own time step, not by the implementation. Its + # error comes from the cusp of exp(-(t / tau)**beta) at t = 0, so the trapezoid sum + # converges as dt**(1 + beta) -- refining dt by four should shrink it by 4**1.6 ~ 9. + gaps = [fft_comparison(0.6, n_time=n) for n in (2**16, 2**18, 2**20)] + + # THEN + ratios = [coarse / fine for coarse, fine in pairwise(gaps)] + + # EXPECT the gap closes at the predicted rate, so the FFT is converging onto our values + assert all(ratio > 4.0 for ratio in ratios), (gaps, ratios) + + ################## + # Unit conversion + ################## + + def test_convert_x_unit(self, stretched_exponential: StretchedExponential): + # WHEN THEN + stretched_exponential.convert_x_unit('microeV') + + # EXPECT the time and the exponent are untouched: they carry no x unit + assert stretched_exponential.x_unit == 'microeV' + assert stretched_exponential.area.value == pytest.approx(2.0 * 1e3) + assert stretched_exponential.center.value == pytest.approx(0.5 * 1e3) + assert stretched_exponential.relaxation_time.value == pytest.approx(3.0) + assert stretched_exponential.relaxation_time.unit == 'ps' + assert stretched_exponential.beta.value == pytest.approx(0.7) + + def test_convert_x_unit_keeps_the_profile(self, stretched_exponential: StretchedExponential): + # WHEN + x = np.linspace(-2.0, 2.0, 51) + before = stretched_exponential.evaluate(x) + + # THEN + stretched_exponential.convert_x_unit('microeV') + + # EXPECT the same curve, read off the rescaled axis + np.testing.assert_allclose(stretched_exponential.evaluate(x * 1e3), before, rtol=1e-12) + + def test_convert_x_unit_invalid_type_raises(self, stretched_exponential: StretchedExponential): + # WHEN THEN EXPECT + with pytest.raises(TypeError, match=r'x_unit must be a string or sc\.Unit'): + stretched_exponential.convert_x_unit(123) + + def test_convert_x_unit_rollback_on_failure(self, stretched_exponential: StretchedExponential): + # WHEN THEN + with pytest.raises(UnitError): + stretched_exponential.convert_x_unit('m') + + # EXPECT: state rolled back + assert stretched_exponential.x_unit == 'meV' + assert stretched_exponential.area.value == pytest.approx(2.0) + assert stretched_exponential.center.value == pytest.approx(0.5) + + def test_convert_y_unit(self): + # WHEN: x_unit='meV', y_unit='1/meV' → area_unit='dimensionless' + stretched = StretchedExponential(area=1.0, x_unit='meV', y_unit='1/meV') + + # THEN: convert y_unit to '1/eV' (same dimension, different scale) + stretched.convert_y_unit('1/eV') + + # EXPECT: y_unit updated and area value rescaled (1e3 factor) + assert stretched.y_unit == '1/eV' + assert stretched.area.value == pytest.approx(1e3) + + def test_convert_y_unit_invalid_type_raises(self, stretched_exponential: StretchedExponential): + # WHEN THEN EXPECT + with pytest.raises(TypeError): + stretched_exponential.convert_y_unit(123) + + def test_convert_y_unit_rollback_on_failure(self): + # WHEN + stretched = StretchedExponential(area=1.0, x_unit='meV') + + # THEN + with pytest.raises(UnitError): + stretched.convert_y_unit('K') + + # EXPECT: state rolled back + assert stretched.y_unit == 'dimensionless' + assert stretched.area.value == pytest.approx(1.0) diff --git a/tools/fit_kww_hwhm_coefficients.py b/tools/fit_kww_hwhm_coefficients.py new file mode 100644 index 000000000..3056d2ee9 --- /dev/null +++ b/tools/fit_kww_hwhm_coefficients.py @@ -0,0 +1,193 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause +""" +Refit ``_HWHM_POLY_COEFFS`` in ``easydynamics.sample_model.components.stretched_exponential``. + +``StretchedExponential.width`` is a dependent easyscience ``Parameter``, and easyscience resolves +a dependency from a string expression, so the half width has to be written in closed form. The +half width is not available in closed form: it is the ``w`` solving + + G_beta(w) = 0.5 * G_beta(0), G_beta(w) = integral of exp(-u**beta) cos(w u) du, + +which ``_reduced_hwhm`` brackets numerically. It is very nearly closed form, though, because + + eps(beta) := beta * ln(HWHM(beta)) - ln(beta) + +is small, smooth and O(0.06), so that + + HWHM(beta) = exp((ln(beta) + eps(beta)) / beta) + +with ``eps`` a low-order polynomial. This script fits that polynomial against ``_reduced_hwhm`` +and prints the coefficient block to paste back into the module. + +The fit degrades below beta = 0.1, but so does the root it is fitted to: by beta = 0.05 the +profile takes some 60 e-folds of ``w`` to fall by half, so ``G_beta`` is nearly flat in ``log w`` +near the crossing. Wuttke's reference implementation libkww (arXiv:0911.4796) likewise supports +only ``0.1 <= beta <= 1.9``. Both are far below any energy grid step by then, so the accuracy +that matters is the one reported for ``beta >= 0.1``. + +Usage +----- +``` +pixi run python tools/fit_kww_hwhm_coefficients.py # refit and print the block +pixi run python tools/fit_kww_hwhm_coefficients.py --check # verify the committed values +pixi run python tools/fit_kww_hwhm_coefficients.py --degree 10 --samples 800 +``` +""" + +import argparse +import sys + +import numpy as np +from numpy.polynomial import polynomial as P + +from easydynamics.sample_model.components.stretched_exponential import _HWHM_POLY_COEFFS +from easydynamics.sample_model.components.stretched_exponential import MAXIMUM_BETA +from easydynamics.sample_model.components.stretched_exponential import MINIMUM_BETA +from easydynamics.sample_model.components.stretched_exponential import _reduced_hwhm + +# The committed coefficients were produced with these settings; changing them changes the fit. +DEFAULT_DEGREE = 8 +DEFAULT_SAMPLES = 400 +# Accuracy is quoted over the range the fit is meant to serve, and where libkww also stops. +REPORTING_FLOOR = 0.1 +# Loose enough to absorb a scipy root-finder nudge, tight enough to catch a real regression. +CHECK_TOLERANCE = 5e-5 + + +def sample_grid(samples: int) -> np.ndarray: + """ + Build the beta grid the fit is evaluated on. + + Log spacing matches how the half width actually varies: it spans some 26 decades across the + supported beta, almost all of that below beta = 0.5. + + Parameters + ---------- + samples : int + Number of grid points. + + Returns + ------- + np.ndarray + Log-spaced beta values over ``[MINIMUM_BETA, MAXIMUM_BETA]``. + """ + return np.exp(np.linspace(np.log(MINIMUM_BETA), np.log(MAXIMUM_BETA), samples)) + + +def fit_coefficients(beta: np.ndarray, degree: int) -> np.ndarray: + """ + Least-squares fit the eps polynomial against the solved half width. + + Parameters + ---------- + beta : np.ndarray + Beta values to fit over. + degree : int + Polynomial degree. + + Returns + ------- + np.ndarray + Coefficients in ascending powers of beta. + """ + hwhm = np.array([_reduced_hwhm(b) for b in beta]) + eps = beta * np.log(hwhm) - np.log(beta) + return P.polyfit(beta, eps, degree) + + +def relative_error(beta: np.ndarray, coefficients: np.ndarray) -> np.ndarray: + """ + Compare the closed form against the solved half width. + + Parameters + ---------- + beta : np.ndarray + Beta values to evaluate at. + coefficients : np.ndarray + Polynomial coefficients in ascending powers of beta. + + Returns + ------- + np.ndarray + Relative error of the closed form at each beta. + """ + exact = np.array([_reduced_hwhm(b) for b in beta]) + approx = np.exp((np.log(beta) + P.polyval(beta, coefficients)) / beta) + return np.abs(approx / exact - 1.0) + + +def report(beta: np.ndarray, coefficients: np.ndarray) -> float: + """ + Print the accuracy of a coefficient set and return the error over the reporting range. + + Parameters + ---------- + beta : np.ndarray + Beta values to evaluate at. + coefficients : np.ndarray + Polynomial coefficients in ascending powers of beta. + + Returns + ------- + float + Maximum relative error for ``beta >= REPORTING_FLOOR``. + """ + error = relative_error(beta, coefficients) + above_floor = error[beta >= REPORTING_FLOOR].max() + print(f' max relative error, all beta : {error.max():.2e}') + print(f' max relative error, beta >= {REPORTING_FLOOR} : {above_floor:.2e}') + for probe in (2.0, 1.0, 0.5, 0.2, 0.1, MINIMUM_BETA): + single = relative_error(np.array([probe]), coefficients)[0] + print(f' beta = {probe:<5} -> {single:.2e}') + return float(above_floor) + + +def main() -> int: + """ + Refit the coefficients, or check the committed ones. + + Returns + ------- + int + Process exit status. + """ + parser = argparse.ArgumentParser(description=__doc__.splitlines()[1]) + parser.add_argument( + '--check', + action='store_true', + help='verify the committed coefficients instead of refitting', + ) + parser.add_argument('--degree', type=int, default=DEFAULT_DEGREE) + parser.add_argument('--samples', type=int, default=DEFAULT_SAMPLES) + args = parser.parse_args() + + beta = sample_grid(args.samples) + + if args.check: + print(f'Checking committed _HWHM_POLY_COEFFS ({len(_HWHM_POLY_COEFFS) - 1} degree):') + above_floor = report(beta, np.array(_HWHM_POLY_COEFFS)) + if above_floor > CHECK_TOLERANCE: + print(f'\nFAIL: {above_floor:.2e} exceeds the {CHECK_TOLERANCE:.0e} tolerance.') + print('Rerun without --check and paste the new block into the module.') + return 1 + print(f'\nOK: within the {CHECK_TOLERANCE:.0e} tolerance.') + return 0 + + coefficients = fit_coefficients(beta, args.degree) + print( + f'Fitted degree {args.degree} over {args.samples} log-spaced beta in ' + f'[{MINIMUM_BETA}, {MAXIMUM_BETA}]:' + ) + report(beta, coefficients) + + print('\nPaste into stretched_exponential.py:\n') + print('_HWHM_POLY_COEFFS = (') + for coefficient in coefficients: + print(f' {float(coefficient)!r},') + print(')') + return 0 + + +if __name__ == '__main__': + sys.exit(main())