diff --git a/docs/docs/tutorials/index.md b/docs/docs/tutorials/index.md index 23a36adc5..10e87adb5 100644 --- a/docs/docs/tutorials/index.md +++ b/docs/docs/tutorials/index.md @@ -53,6 +53,9 @@ tutorials. - [DeltaLorentz](delta_lorentz.ipynb) – Learn how to create and use a model with a Delta function and a Lorentzian that have a shared Debye-Waller-like Q-dependence. +- [Mittag-Leffler relaxation](mittag_leffler.ipynb) – Learn how to use + the multiscale relaxation model of Hassani et al. (2022), both as a + single lineshape and as a diffusion model across Q. - [Sample model](sample_model.ipynb) – Learn how to create a model of the scattering from your sample including model components and diffusion models. diff --git a/docs/docs/tutorials/mittag_leffler.ipynb b/docs/docs/tutorials/mittag_leffler.ipynb new file mode 100644 index 000000000..69b84864d --- /dev/null +++ b/docs/docs/tutorials/mittag_leffler.ipynb @@ -0,0 +1,561 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "mittag-leffler-0", + "metadata": {}, + "source": [ + "# Mittag-Leffler relaxation\n", + "\n", + "The multiscale relaxation model of\n", + "\n", + "> A. N. Hassani, L. Haris, M. Appel, T. Seydel, A. M. Stadler and G. R. Kneller,\n", + "> *Multiscale relaxation dynamics and diffusion of myelin basic protein in solution studied by\n", + "> quasielastic neutron scattering*, J. Chem. Phys. **156**, 025102 (2022),\n", + "> [doi:10.1063/5.0077100](https://doi.org/10.1063/5.0077100).\n", + "\n", + "Internal protein dynamics does not relax exponentially. It relaxes with a power-law tail, which\n", + "reflects a broad distribution of relaxation rates rather than a single one. The paper models that\n", + "with a Mittag-Leffler (ML) relaxation function (their Eq. 28)\n", + "\n", + "$$ \\phi_{ML}(t) = E_\\alpha\\left(-(|t|/\\tau_R)^\\alpha\\right), \\qquad 0 < \\alpha \\le 1 $$\n", + "\n", + "which decays as $(t/\\tau_R)^{-\\alpha}/\\Gamma(1-\\alpha)$ at long times and becomes a simple\n", + "exponential at $\\alpha = 1$. For a protein in solution this internal relaxation is multiplied by\n", + "the global translational diffusion of the whole molecule (their Eqs. 38-39),\n", + "\n", + "$$ F^{(+)}(t) = e^{-\\epsilon|t|}\\left(EISF + (1 - EISF)\\,\\phi_{ML}(|t|)\\right),\n", + "\\qquad \\epsilon = \\hbar D q^2 $$\n", + "\n", + "and the measured spectrum is its Fourier transform, their Eq. (41):\n", + "\n", + "$$ S^{(+)}(x) = EISF\\,\\frac{1}{\\pi}\\frac{\\epsilon}{x^2+\\epsilon^2}\n", + "+ (1 - EISF)\\,\\tilde{\\phi}^{(\\epsilon)}_{ML}(|x|). $$\n", + "\n", + "Note that the elastic term is a **Lorentzian of half width $\\epsilon$**, not a delta function:\n", + "global diffusion broadens it.\n", + "\n", + "EasyDynamics provides this at two levels.\n", + "\n", + "- `DiffusionDampedMittagLeffler` is the quasi-elastic lineshape\n", + " $\\tilde{\\phi}^{(\\epsilon)}_{ML}$ of their Eq. (42), a `ModelComponent` with a free `damping`\n", + " that you can fit however you like.\n", + "- `MittagLefflerDiffusion` is the full Eq. (41) across $Q$, a `DiffusionModel` in which\n", + " `damping` is constrained to $\\hbar D q^2$ with a single global $D$.\n", + "\n", + "## Two notes on Eq. (42)\n", + "\n", + "The component evaluates Eq. (42) with two deliberate differences from the printed equation.\n", + "\n", + "1. **The printed equation is written for $\\tau_R = 1$.** You can see this without any derivation:\n", + " it forms $(\\omega^2+\\epsilon^2)^\\alpha + 1$, and you cannot add a pure number to something\n", + " carrying dimensions of frequency$^{2\\alpha}$ unless the frequency is already dimensionless. The\n", + " component evaluates the printed expression at the reduced variables $\\omega/\\Gamma$ and\n", + " $\\epsilon/\\Gamma$ and divides by $\\Gamma = \\hbar/\\tau_R$, which restores a general relaxation\n", + " time while leaving the printed expression itself untouched. $\\Gamma$ is the component's\n", + " `width`.\n", + "2. **The $1/\\pi$ prefactor is missing from the printed equation.** Applying the paper's own\n", + " Eq. (40), $\\tilde{\\phi}^{(\\epsilon)} = \\frac{1}{\\pi}\\Re\\{\\hat{\\phi}(\\epsilon+i\\omega)\\}$, to\n", + " the ML Laplace transform of their Eq. (32) reproduces Eq. (42) exactly, prefactor included.\n", + " Without it the lineshape integrates to $\\pi$ rather than to 1, which would be inconsistent with\n", + " the normalised elastic term in Eq. (41). It is included here, so `scale` is the integrated area\n", + " of the profile." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-1", + "metadata": {}, + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.optimize import curve_fit\n", + "\n", + "import easydynamics as edyn\n", + "\n", + "%matplotlib widget" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-2", + "metadata": {}, + "source": [ + "## The component\n", + "\n", + "`DiffusionDampedMittagLeffler` has four parameters: `scale` (the integrated area), `alpha` (the\n", + "form parameter $\\alpha$), `width` (the relaxation rate $\\Gamma = \\hbar/\\tau_R$, as an energy) and\n", + "`damping` (the diffusion damping $\\epsilon$, also as an energy).\n", + "\n", + "The form parameter is what makes this more than a Lorentzian. At $\\alpha = 1$ the relaxation is\n", + "exponential and the profile is exactly a Lorentzian of half width at half maximum\n", + "`width + damping`. Lowering $\\alpha$ broadens the distribution of relaxation rates, which\n", + "sharpens the peak and fattens the wings at the same time." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-3", + "metadata": {}, + "source": [ + "x = np.linspace(-2, 2, 801)\n", + "\n", + "plt.figure()\n", + "for alpha in [1.0, 0.8, 0.6, 0.4]:\n", + " mittag_leffler = edyn.DiffusionDampedMittagLeffler(\n", + " display_name=f'alpha={alpha}',\n", + " scale=1.0,\n", + " alpha=alpha,\n", + " width=0.1,\n", + " damping=0.02,\n", + " )\n", + " plt.plot(x, mittag_leffler.evaluate(x), label=mittag_leffler.display_name)\n", + "\n", + "# alpha=1 must be a plain Lorentzian of HWHM width+damping, shown here as a check.\n", + "lorentzian = edyn.Lorentzian(display_name='Lorentzian, HWHM=0.12', area=1.0, width=0.12)\n", + "plt.plot(x, lorentzian.evaluate(x), 'k--', label=lorentzian.display_name)\n", + "\n", + "plt.yscale('log')\n", + "plt.xlabel('Energy (meV)')\n", + "plt.ylabel('Intensity (arb. units)')\n", + "plt.legend()\n", + "plt.show()" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-4", + "metadata": {}, + "source": [ + "`damping` must be strictly positive. It is what keeps the profile finite at zero energy transfer:\n", + "without it the ML spectrum diverges as $|\\omega|^{\\alpha-1}$, as the paper notes below its\n", + "Eq. (34). Increasing it broadens the profile without changing its area." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-5", + "metadata": {}, + "source": [ + "plt.figure()\n", + "for damping in [0.005, 0.02, 0.1, 0.4]:\n", + " mittag_leffler = edyn.DiffusionDampedMittagLeffler(\n", + " scale=1.0, alpha=0.6, width=0.1, damping=damping\n", + " )\n", + " plt.plot(x, mittag_leffler.evaluate(x), label=f'damping={damping} meV')\n", + "\n", + "plt.yscale('log')\n", + "plt.xlabel('Energy (meV)')\n", + "plt.ylabel('Intensity (arb. units)')\n", + "plt.legend()\n", + "plt.show()\n", + "\n", + "# The scale parameter is the integrated area. The wings fall off as slowly as |x|**(-1-alpha),\n", + "# so this check needs a very wide, logarithmically spaced grid to converge.\n", + "half_grid = np.logspace(-5, 9, 30000)\n", + "xx = np.concatenate([-half_grid[::-1], half_grid])\n", + "for alpha in [1.0, 0.6]:\n", + " mittag_leffler = edyn.DiffusionDampedMittagLeffler(\n", + " scale=1.0, alpha=alpha, width=0.1, damping=0.02\n", + " )\n", + " area = np.trapezoid(mittag_leffler.evaluate(xx), xx)\n", + " print(f'alpha={alpha}: area under the curve = {area:.4f}')" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-6", + "metadata": {}, + "source": [ + "## The diffusion model\n", + "\n", + "`MittagLefflerDiffusion` assembles Eq. (41) at every $Q$: an elastic `Lorentzian` of area\n", + "$K\\cdot EISF$ and half width $\\epsilon = \\hbar D q^2$, plus a `DiffusionDampedMittagLeffler` of\n", + "scale $K\\cdot(1-EISF)$ and the same damping $\\epsilon$. The diffusion coefficient $D$ is the\n", + "single global parameter that ties the $Q$ values together, exactly as in\n", + "`BrownianTranslationalDiffusion`.\n", + "\n", + "The numbers below are those of MBP in $D_2O$ buffer at 283 K: $D = 3.31$ Ų/ns from dynamic light\n", + "scattering (Table I of the paper), and a $Q$ range matching IN16B's\n", + "0.8 Å⁻¹ < $q$ < 1.8 Å⁻¹." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-7", + "metadata": {}, + "source": [ + "Q = np.linspace(0.8, 1.8, 6)\n", + "\n", + "# hbar in meV*ps, used to turn a relaxation time in ps into a rate in meV.\n", + "HBAR_MEV_PS = 0.6582119569\n", + "\n", + "model = edyn.MittagLefflerDiffusion(\n", + " display_name='MBP in D2O, 283 K',\n", + " scale=1.0,\n", + " diffusion_coefficient=3.31e-11, # 3.31 A^2/ns, from DLS\n", + " A_0=0.02,\n", + " relaxation_rate=HBAR_MEV_PS / 50.0, # tau_R = 50 ps\n", + " alpha=0.85,\n", + " Q=Q,\n", + ")\n", + "print(model)\n", + "print(model.get_component_collections()[0].list_component_names())" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "code", + "id": "mittag-leffler-8", + "metadata": {}, + "source": [ + "energy = np.linspace(-0.15, 0.15, 1201)\n", + "collections = model.get_component_collections()\n", + "cmap = plt.cm.viridis\n", + "\n", + "plt.figure()\n", + "for i, collection in enumerate(collections):\n", + " plt.plot(\n", + " energy,\n", + " collection.evaluate(energy),\n", + " color=cmap(i / (len(collections) - 1)),\n", + " label=f'Q={Q[i]:.1f} 1/A',\n", + " )\n", + "plt.yscale('log')\n", + "plt.xlabel('Energy (meV)')\n", + "plt.ylabel('S(Q, E) (arb. units)')\n", + "plt.legend()\n", + "plt.show()\n", + "\n", + "# The two components at a single Q: the elastic Lorentzian and the Mittag-Leffler term.\n", + "plt.figure()\n", + "collection = collections[-1]\n", + "for name in collection.list_component_names():\n", + " plt.plot(energy, collection.evaluate_component(energy, name), label=name)\n", + "plt.plot(energy, collection.evaluate(energy), 'k--', label='Sum')\n", + "plt.yscale('log')\n", + "plt.xlabel('Energy (meV)')\n", + "plt.ylabel('S(Q, E) (arb. units)')\n", + "plt.title(f'Q = {Q[-1]:.1f} 1/A')\n", + "plt.legend()\n", + "plt.show()" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-9", + "metadata": {}, + "source": [ + "The damping $\\epsilon = \\hbar D q^2$ is shared by both components by construction, and it is what\n", + "`calculate_width` returns. At these $Q$ values it is a few µeV, comparable with IN16B's 3.5 µeV\n", + "resolution — which is exactly the difficulty the paper describes: the slow relaxation modes\n", + "overlap with the elastic line and with global diffusion." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-10", + "metadata": {}, + "source": [ + "Q_fine = np.linspace(0.1, 2.0, 101)\n", + "\n", + "plt.figure()\n", + "plt.plot(Q_fine, model.calculate_width(Q_fine) * 1e3)\n", + "plt.xlabel('Q (1/A)')\n", + "plt.ylabel('epsilon = hbar*D*Q^2 (ueV)')\n", + "plt.axhline(3.5, color='k', linestyle='dashed', label='IN16B resolution, 3.5 ueV')\n", + "plt.legend()\n", + "plt.show()\n", + "\n", + "for collection, epsilon in zip(collections, model.calculate_width(), strict=True):\n", + " print(\n", + " f'epsilon = {epsilon * 1e3:6.3f} ueV | '\n", + " f'Lorentzian width = {collection[0].width.value * 1e3:6.3f} ueV | '\n", + " f'ML damping = {collection[1].damping.value * 1e3:6.3f} ueV'\n", + " )" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-11", + "metadata": {}, + "source": [ + "### Parameters that vary with Q\n", + "\n", + "In the paper only $D$ is global. The relaxation time $\\tau$, the form parameter $\\alpha$ and the\n", + "EISF are all fitted **independently at each $q$** (their Fig. 5). Pass `allow_Q_variation` to get\n", + "that: each $Q$ then owns its own parameters, while $D$ and the overall scale stay shared." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-12", + "metadata": {}, + "source": [ + "model_per_Q = edyn.MittagLefflerDiffusion(\n", + " display_name='MBP in D2O, 283 K, per-Q parameters',\n", + " scale=1.0,\n", + " diffusion_coefficient=3.31e-11,\n", + " A_0=0.02,\n", + " relaxation_rate=HBAR_MEV_PS / 50.0,\n", + " alpha=0.85,\n", + " allow_Q_variation={'A_0': True, 'relaxation_rate': True, 'alpha': True},\n", + " Q=Q,\n", + ")\n", + "\n", + "print('Global parameters :', [p.name for p in model_per_Q.get_global_variables()])\n", + "print('Free parameters :', len(model_per_Q.get_free_parameters()))" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-13", + "metadata": {}, + "source": [ + "The paper does not tabulate the fitted $\\tau(q)$, $\\alpha(q)$ and $EISF(q)$ — they appear only as\n", + "plots in its Fig. 5. The values used below are therefore **indicative**, read off that figure for\n", + "MBP in $D_2O$ at 283 K, and are here only to show the machinery and the trends the paper\n", + "describes: $\\tau$ falls with $q$ because higher $q$ probes more localised and faster motion, while\n", + "$\\alpha$ rises towards 1 because localised motion relaxes more nearly exponentially." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-14", + "metadata": {}, + "source": [ + "# Indicative values read off Fig. 5 of the paper; not published numbers.\n", + "tau_ps = np.array([110.0, 75.0, 55.0, 42.0, 33.0, 27.0])\n", + "alpha_per_Q = np.array([0.80, 0.83, 0.86, 0.88, 0.90, 0.92])\n", + "eisf_per_Q = np.array([0.04, 0.03, 0.02, 0.02, 0.01, 0.01])\n", + "\n", + "for i in range(len(Q)):\n", + " model_per_Q._relaxation_rate_list[i].value = HBAR_MEV_PS / tau_ps[i] # ruff: ignore[private-member-access]\n", + " model_per_Q._alpha_list[i].value = alpha_per_Q[i] # ruff: ignore[private-member-access]\n", + " model_per_Q._A_0_list[i].value = eisf_per_Q[i] # ruff: ignore[private-member-access]\n", + "\n", + "fig, axes = plt.subplots(3, 1, sharex=True, figsize=(5, 7))\n", + "axes[0].plot(Q, HBAR_MEV_PS / model_per_Q.calculate_relaxation_rate(), 'o-')\n", + "axes[0].set_ylabel('tau (ps)')\n", + "axes[0].set_ylim(0, 130)\n", + "axes[1].plot(Q, model_per_Q.calculate_alpha(), 'o-')\n", + "axes[1].set_ylabel('alpha')\n", + "axes[1].set_ylim(0, 1.2)\n", + "axes[2].plot(Q, model_per_Q.calculate_EISF(), 'o-')\n", + "axes[2].set_ylabel('EISF')\n", + "axes[2].set_ylim(0, 0.2)\n", + "axes[2].set_xlabel('Q (1/A)')\n", + "fig.suptitle('Layout of Fig. 5, with indicative values')\n", + "plt.show()" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-15", + "metadata": {}, + "source": [ + "## Comparing with the paper\n", + "\n", + "Most of the figures in the paper cannot be reproduced from the paper alone.\n", + "\n", + "- **Figs. 1 and 2** show measured spectra and intermediate scattering functions. They need the raw\n", + " IN16B data, which lives at [doi.ill.fr/10.5291/ILL-DATA.8-04-874](http://doi.ill.fr/10.5291/ILL-DATA.8-04-874).\n", + "- **Fig. 4** shows circular dichroism results, which are not tabulated.\n", + "- **Figs. 5 and 6** are built from the per-$q$ fit parameters, which are only ever shown as plots.\n", + " The layout of Fig. 5 is reproduced above and Fig. 6 below, but with indicative parameters rather\n", + " than the published ones.\n", + "\n", + "Two figures *are* reproducible from what the paper publishes, and both are done below from the\n", + "library's own code: the **Arrhenius plot of Fig. 3** from the diffusion coefficients in Table I,\n", + "and the **energy barrier distribution of Fig. 7** from the closed-form Eq. (48)." + ] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-16", + "metadata": {}, + "source": [ + "### Fig. 3, left panel: Arrhenius plot of the diffusion coefficient\n", + "\n", + "Table I of the paper lists the diffusion coefficients measured by dynamic light scattering for\n", + "both solvents. Fitting $D(T) = D_0 e^{-\\Delta G / k_B T}$ (their Eq. 44) to them recovers the\n", + "activation energies quoted in the text. MBP coagulates at 323 K in pure $D_2O$, so the paper\n", + "replaces the measured 0.31 Ų/ns there by the value extrapolated from lower temperatures, and so\n", + "do we." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-17", + "metadata": {}, + "source": [ + "# Table I of Hassani et al. (2022). D_DLS in A^2/ns.\n", + "temperature = np.array([283.0, 293.0, 303.0, 313.0, 323.0])\n", + "D_D2O = np.array([3.31, 5.72, 7.49, 9.85, 11.94]) # 323 K: extrapolated, MBP coagulates\n", + "D_TFE = np.array([2.05, 3.17, 4.42, 6.32, 8.15])\n", + "\n", + "BOLTZMANN_KCAL = 1.987204259e-3 # kcal/(mol*K)\n", + "\n", + "\n", + "def arrhenius(T: np.ndarray, D_0: float, delta_G: float) -> np.ndarray:\n", + " \"\"\"Eq. (44): D(T) = D_0 exp(-dG / kB T), with dG in kcal/mol.\"\"\"\n", + " return D_0 * np.exp(-delta_G / (BOLTZMANN_KCAL * T))\n", + "\n", + "\n", + "plt.figure()\n", + "T_fine = np.linspace(280.0, 326.0, 200)\n", + "for label, D, colour, quoted in [\n", + " ('MBP D2O', D_D2O, 'tab:blue', 5.10),\n", + " ('MBP 30% TFE', D_TFE, 'tab:orange', 6.05),\n", + "]:\n", + " (D_0, delta_G), _ = curve_fit(arrhenius, temperature, D, p0=[1e5, 5.0])\n", + " plt.plot(1000.0 / temperature, D, 'o', color=colour, label=label)\n", + " plt.plot(1000.0 / T_fine, arrhenius(T_fine, D_0, delta_G), '-', color=colour)\n", + " print(f'{label:12s}: fitted dG = {delta_G:.2f} kcal/mol (paper quotes {quoted:.2f})')\n", + "\n", + "plt.yscale('log')\n", + "plt.xlabel('1000/T (1/K)')\n", + "plt.ylabel('D (A^2/ns)')\n", + "plt.legend()\n", + "plt.title('Fig. 3, left panel')\n", + "plt.show()" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-18", + "metadata": {}, + "source": [ + "Both activation energies come out at the values quoted in the paper, 5.10 and 6.05 kcal/mol. Note\n", + "that Eq. (44) has to be fitted to $D$ itself rather than to $\\ln D$: least squares on the\n", + "logarithm weights the low-temperature points much more heavily and gives 5.69 and 6.28 kcal/mol\n", + "instead." + ] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-19", + "metadata": {}, + "source": [ + "### Fig. 7, right panel: the distribution of energy barriers\n", + "\n", + "Reading the ML relaxation as diffusion in Zwanzig's rough harmonic potential, the paper derives a\n", + "distribution of dimensionless barrier heights $h = \\Delta E / k_B T$ (their Eq. 48),\n", + "\n", + "$$ P_{ML}(h) = \\frac{2h\\sin(\\pi\\alpha)}\n", + "{\\pi\\left(e^{-\\alpha h^2} + e^{\\alpha h^2} + 2\\cos(\\pi\\alpha)\\right)}. $$\n", + "\n", + "This is a closed form in $h$ and $\\alpha$ only, so it reproduces exactly. As $\\alpha \\to 1$ it\n", + "collapses onto $\\delta(h)$ — a smooth potential with no barriers — and as $\\alpha \\to 0$ it\n", + "broadens to include arbitrarily high ones. That is the physical reading of the form parameter:\n", + "**a smaller $\\alpha$ means a rougher energy landscape.**" + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-20", + "metadata": {}, + "source": [ + "h = np.linspace(0.0, 10.0, 201)\n", + "alpha_grid = np.linspace(0.02, 1.0, 100)\n", + "\n", + "barrier_model = edyn.MittagLefflerDiffusion(Q=np.array([1.0]))\n", + "distribution = np.empty((alpha_grid.size, h.size))\n", + "for i, alpha in enumerate(alpha_grid):\n", + " barrier_model.alpha = alpha\n", + " distribution[i] = barrier_model.calculate_energy_barrier_distribution(h)[0]\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(projection='3d')\n", + "H, A = np.meshgrid(h, alpha_grid)\n", + "ax.plot_surface(H, A, distribution, cmap='Blues_r', linewidth=0, antialiased=True)\n", + "ax.set_xlabel('h')\n", + "ax.set_ylabel('alpha')\n", + "ax.set_zlabel('P(h)')\n", + "ax.set_zlim(0, 1)\n", + "ax.set_title('Fig. 7, right panel')\n", + "plt.show()" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "mittag-leffler-21", + "metadata": {}, + "source": [ + "### Fig. 6: relaxation rate spectra\n", + "\n", + "The same model has a closed-form distribution of relaxation *rates* (their Eq. 37),\n", + "\n", + "$$ p_{ML}(\\lambda) = \\frac{\\sin(\\pi\\alpha)}\n", + "{\\pi\\lambda\\left((\\lambda\\tau_R)^{-\\alpha} + (\\lambda\\tau_R)^{\\alpha} + 2\\cos(\\pi\\alpha)\\right)}, $$\n", + "\n", + "which is what the paper plots against $q$ in its Fig. 6. It depends on the fitted $\\tau(q)$ and\n", + "$\\alpha(q)$, so the surface below uses the indicative values set above rather than the published\n", + "ones. The trend the paper reports is nevertheless the one that shows: as $q$ grows the peak moves\n", + "to higher rates and the distribution narrows, because $\\alpha$ approaches 1." + ] + }, + { + "cell_type": "code", + "id": "mittag-leffler-22", + "metadata": {}, + "source": [ + "# lambda in 1/ps, converted to the model's meV via hbar.\n", + "rate_per_ps = np.linspace(0.001, 0.15, 200)\n", + "spectrum = model_per_Q.calculate_relaxation_rate_spectrum(rate_per_ps * HBAR_MEV_PS)\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(projection='3d')\n", + "RATE, QQ = np.meshgrid(rate_per_ps, Q)\n", + "ax.plot_surface(RATE, QQ, spectrum, cmap='autumn', linewidth=0, antialiased=True)\n", + "ax.set_xlabel('lambda (1/ps)')\n", + "ax.set_ylabel('Q (1/A)')\n", + "ax.set_zlabel('p(lambda)')\n", + "ax.set_title('Layout of Fig. 6, with indicative values')\n", + "plt.show()" + ], + "execution_count": null, + "outputs": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "default", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/docs/user-guide/concept.md b/docs/docs/user-guide/concept.md index 17014b924..d4cc2c030 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`, +`DiffusionDampedMittagLeffler` 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/docs/mkdocs.yml b/docs/mkdocs.yml index 6575f0c95..b8659848e 100644 --- a/docs/mkdocs.yml +++ b/docs/mkdocs.yml @@ -212,6 +212,7 @@ nav: - Detailed balance: tutorials/detailed_balance.ipynb - DeltaLorentz: tutorials/delta_lorentz.ipynb - Diffusion model: tutorials/diffusion_model.ipynb + - Mittag-Leffler relaxation: tutorials/mittag_leffler.ipynb - Sample model: tutorials/sample_model.ipynb - Instrument model: tutorials/instrument_model.ipynb - Experiment: tutorials/experiment.ipynb diff --git a/src/easydynamics/__init__.py b/src/easydynamics/__init__.py index de7391225..b3ce8f1a4 100644 --- a/src/easydynamics/__init__.py +++ b/src/easydynamics/__init__.py @@ -29,12 +29,14 @@ from easydynamics.sample_model import DampedHarmonicOscillator from easydynamics.sample_model import DeltaFunction from easydynamics.sample_model import DeltaLorentz +from easydynamics.sample_model import DiffusionDampedMittagLeffler from easydynamics.sample_model import Exponential from easydynamics.sample_model import ExpressionComponent from easydynamics.sample_model import Gaussian from easydynamics.sample_model import InstrumentModel from easydynamics.sample_model import JumpTranslationalDiffusion from easydynamics.sample_model import Lorentzian +from easydynamics.sample_model import MittagLefflerDiffusion from easydynamics.sample_model import Polynomial from easydynamics.sample_model import ResolutionModel from easydynamics.sample_model import SampleModel @@ -62,6 +64,7 @@ 'DeltaFunction', 'DeltaLorentz', 'DetailedBalanceSettings', + 'DiffusionDampedMittagLeffler', 'EasyDynamicsBase', 'EasyDynamicsModelBase', 'Experiment', @@ -72,6 +75,7 @@ 'InstrumentModel', 'JumpTranslationalDiffusion', 'Lorentzian', + 'MittagLefflerDiffusion', 'MultiQPosteriorSampler', 'ParameterAnalysis', 'ParameterLabels', diff --git a/src/easydynamics/sample_model/__init__.py b/src/easydynamics/sample_model/__init__.py index 07f434a18..775cb3d23 100644 --- a/src/easydynamics/sample_model/__init__.py +++ b/src/easydynamics/sample_model/__init__.py @@ -7,6 +7,9 @@ DampedHarmonicOscillator, ) from easydynamics.sample_model.components.delta_function import DeltaFunction +from easydynamics.sample_model.components.diffusion_damped_mittag_leffler import ( + DiffusionDampedMittagLeffler, +) from easydynamics.sample_model.components.exponential import Exponential from easydynamics.sample_model.components.expression_component import ExpressionComponent from easydynamics.sample_model.components.gaussian import Gaussian @@ -20,6 +23,9 @@ from easydynamics.sample_model.diffusion_model.jump_translational_diffusion import ( JumpTranslationalDiffusion, ) +from easydynamics.sample_model.diffusion_model.mittag_leffler_diffusion import ( + MittagLefflerDiffusion, +) from easydynamics.sample_model.instrument_model import InstrumentModel from easydynamics.sample_model.resolution_model import ResolutionModel from easydynamics.sample_model.sample_model import SampleModel @@ -31,12 +37,14 @@ 'DampedHarmonicOscillator', 'DeltaFunction', 'DeltaLorentz', + 'DiffusionDampedMittagLeffler', 'Exponential', 'ExpressionComponent', 'Gaussian', 'InstrumentModel', 'JumpTranslationalDiffusion', 'Lorentzian', + 'MittagLefflerDiffusion', 'Polynomial', 'ResolutionModel', 'SampleModel', diff --git a/src/easydynamics/sample_model/components/__init__.py b/src/easydynamics/sample_model/components/__init__.py index 940f41a7e..951b6cd9b 100644 --- a/src/easydynamics/sample_model/components/__init__.py +++ b/src/easydynamics/sample_model/components/__init__.py @@ -5,6 +5,9 @@ DampedHarmonicOscillator, ) from easydynamics.sample_model.components.delta_function import DeltaFunction +from easydynamics.sample_model.components.diffusion_damped_mittag_leffler import ( + DiffusionDampedMittagLeffler, +) from easydynamics.sample_model.components.exponential import Exponential from easydynamics.sample_model.components.expression_component import ExpressionComponent from easydynamics.sample_model.components.gaussian import Gaussian @@ -15,6 +18,7 @@ __all__ = [ 'DampedHarmonicOscillator', 'DeltaFunction', + 'DiffusionDampedMittagLeffler', 'Exponential', 'ExpressionComponent', 'Gaussian', diff --git a/src/easydynamics/sample_model/components/diffusion_damped_mittag_leffler.py b/src/easydynamics/sample_model/components/diffusion_damped_mittag_leffler.py new file mode 100644 index 000000000..9e04a60c7 --- /dev/null +++ b/src/easydynamics/sample_model/components/diffusion_damped_mittag_leffler.py @@ -0,0 +1,448 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +from __future__ import annotations + +import numpy as np +import scipp as sc +from easyscience.variable import Parameter + +from easydynamics.sample_model.components.mixins import CreateParametersMixin +from easydynamics.sample_model.components.model_component import ModelComponent +from easydynamics.utils.utils import Numeric + +MINIMUM_ALPHA = 1e-10 # To avoid a vanishing form parameter +MAXIMUM_ALPHA = 1.0 # Above 1 the Mittag-Leffler function is no longer a relaxation function + + +class DiffusionDampedMittagLeffler(CreateParametersMixin, ModelComponent): + r""" + Model of the diffusion-damped Mittag-Leffler relaxation spectrum, Eq. (42) of Hassani et al. + + This is the spectral (energy-domain) form of the multiscale relaxation model used for + intrinsically disordered proteins in + + A. N. Hassani, L. Haris, M. Appel, T. Seydel, A. M. Stadler and G. R. Kneller, *Multiscale + relaxation dynamics and diffusion of myelin basic protein in solution studied by quasielastic + neutron scattering*, J. Chem. Phys. **156**, 025102 (2022), + [doi:10.1063/5.0077100](https://doi.org/10.1063/5.0077100). + + In the time domain the internal dynamics is described by the Mittag-Leffler (ML) relaxation + function $\phi_{ML}(t) = E_\alpha(-(|t|/\tau_R)^\alpha)$ (their Eq. 28), damped by global + translational diffusion, $e^{-\epsilon |t|}$ with $\epsilon = D q^2$ (their Eqs. 38-39). The + Fourier transform of that product is the "generalised Lorentzian" of their Eq. (42), + + $$ \tilde{\phi}^{(\epsilon)}_{ML}(\omega) = \frac{1}{\pi} + \frac{\epsilon(\omega^2+\epsilon^2)^{\alpha/2} + + \omega\sin(\alpha\arg(\epsilon+i|\omega|)) + + \epsilon\cos(\alpha\arg(\epsilon+i|\omega|))} + {(\omega^2+\epsilon^2)\left(\left((\omega^2+\epsilon^2)^{\alpha}+1\right) + (\omega^2+\epsilon^2)^{-\alpha/2} + 2\cos(\alpha\arg(\epsilon+i|\omega|))\right)} $$ + + which this component evaluates as + + $$ I(x) = \frac{A}{\Gamma} \tilde{\phi}^{(\epsilon/\Gamma)}_{ML}\left(\frac{|x|}{\Gamma}\right) + $$ + + where $A$ is the scale factor (``scale``), $\alpha$ is the form parameter (``alpha``), $\Gamma + = \hbar/\tau_R$ is the ML relaxation rate expressed as an energy (``width``) and $\epsilon = + \hbar D q^2$ is the diffusion damping, also expressed as an energy (``damping``). scale has + unit = x_unit * y_unit; width and damping have unit = x_unit; alpha is dimensionless. + + Two remarks on the relation to the printed Eq. (42): + + - Eq. (42) is written for $\tau_R = 1$, i.e. for $\omega$ and $\epsilon$ measured in units of + the relaxation rate. Rescaling $\omega \to \omega/\Gamma$ and $\epsilon \to \epsilon/\Gamma$ + and dividing by $\Gamma$ restores a general relaxation time while leaving the printed + expression intact; that is the form implemented here. + - The $1/\pi$ prefactor is absent from the printed Eq. (42), but is required for the lineshape + to be normalised the way Eq. (41) of the same paper assumes, i.e. + $\int d\omega\, \tilde{\phi}^{(\epsilon)}_{ML}(\omega) = 1$. It is included here, so + ``scale`` is the integrated area of the profile, and the $\alpha \to 1$ limit is a Lorentzian + of area ``scale`` and half width at half maximum ``width + damping``. + + The profile is symmetric about $x = 0$, so there is no center parameter. Because a strictly + positive ``damping`` is enforced, the profile stays regular at $x = 0$, where the undamped ML + spectrum would diverge as $|\omega|^{\alpha - 1}$. + + Examples + -------- + **Creating a diffusion-damped Mittag-Leffler component** + + ```python + import numpy as np + import easydynamics as edyn + + ml = edyn.DiffusionDampedMittagLeffler(scale=1.0, alpha=0.8, width=0.05, damping=0.01) + x = np.linspace(-0.5, 0.5, 200) + values = ml.evaluate(x) + ``` + + **Modifying parameters after construction** + + ```python + import easydynamics as edyn + + ml = edyn.DiffusionDampedMittagLeffler(name='MBP internal dynamics') + ml.scale = 0.9 + ml.alpha = 0.75 + ml.width = 0.02 + ml.damping = 0.005 + ``` + """ + + def __init__( + self, + scale: Numeric = 1.0, + alpha: Numeric = 1.0, + width: Numeric = 1.0, + damping: Numeric = 1.0, + x_unit: str | sc.Unit = 'meV', + y_unit: str | sc.Unit = 'dimensionless', + name: str = 'DiffusionDampedMittagLeffler', + display_name: str | None = None, + unique_name: str | None = None, + ) -> None: + """ + Initialize the diffusion-damped Mittag-Leffler component. + + Parameters + ---------- + scale : Numeric, default=1.0 + Scale factor in front of the normalised profile, i.e. the integrated area of the + profile. Unit is ``x_unit * y_unit``. Must be non-negative. + alpha : Numeric, default=1.0 + Form parameter of the Mittag-Leffler relaxation function. Must lie in (0, 1]. alpha=1 + gives exponential relaxation, and hence a Lorentzian profile, while smaller values give + a broader relaxation rate spectrum. + width : Numeric, default=1.0 + Mittag-Leffler relaxation rate ``hbar / tau_R`` in x_unit. Must be strictly positive. + damping : Numeric, default=1.0 + Diffusion damping ``hbar * D * q**2`` in x_unit. Must be strictly positive; it is what + keeps the profile finite at x=0. + x_unit : str | sc.Unit, default='meV' + Unit of the x-axis. width and damping are stored in this unit. scale_unit = x_unit * + y_unit. + y_unit : str | sc.Unit, default='dimensionless' + Unit of the y-axis (output). + name : str, default='DiffusionDampedMittagLeffler' + 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._scale = self._create_scale_parameter( + scale=scale, name=name, x_unit=self.x_unit, y_unit=self.y_unit + ) + self._alpha = self._create_alpha_parameter(alpha=alpha, name=name) + # These methods live in CreateParametersMixin + self._width = self._create_width_parameter(width=width, name=name, x_unit=self.x_unit) + self._damping = self._create_width_parameter( + width=damping, name=name, param_name='damping', x_unit=self.x_unit + ) + + @staticmethod + def _create_scale_parameter( + scale: Numeric, + name: str, + x_unit: str | sc.Unit, + y_unit: str | sc.Unit, + ) -> Parameter: + """ + Create the scale Parameter with unit = x_unit * y_unit. + + Parameters + ---------- + scale : Numeric + Initial scale value. Must be non-negative. + name : str + Base name used to label the Parameter (``name + ' scale'``). + x_unit : str | sc.Unit + X-axis unit. The resulting scale unit is ``x_unit * y_unit``. + y_unit : str | sc.Unit + Y-axis unit. The resulting scale unit is ``x_unit * y_unit``. + + Returns + ------- + Parameter + Configured scale Parameter with ``unit = x_unit * y_unit`` and ``min = 0``. + + Raises + ------ + TypeError + If *scale* is not a numeric type. + ValueError + If *scale* is not finite, or is negative. + """ + if not isinstance(scale, Numeric): + raise TypeError('scale must be a number.') + if not np.isfinite(scale): + raise ValueError('scale must be a finite number.') + if float(scale) < 0: + raise ValueError('scale must be non-negative.') + return Parameter( + name=name + ' scale', + value=float(scale), + unit=str(sc.Unit(x_unit) * sc.Unit(y_unit)), + min=0.0, + ) + + @staticmethod + def _create_alpha_parameter(alpha: Numeric, name: str) -> Parameter: + """ + Create the dimensionless form parameter alpha, bounded to (0, 1]. + + Parameters + ---------- + alpha : Numeric + Initial value of the form parameter. + name : str + Base name used to label the Parameter (``name + ' alpha'``). + + Returns + ------- + Parameter + Configured alpha Parameter with ``unit = 'dimensionless'``. + + Raises + ------ + TypeError + If *alpha* is not a numeric type. + ValueError + If *alpha* is not finite, or does not lie in (0, 1]. + """ + if not isinstance(alpha, Numeric): + raise TypeError('alpha must be a number.') + if not np.isfinite(alpha): + raise ValueError('alpha must be a finite number.') + if not MINIMUM_ALPHA <= float(alpha) <= MAXIMUM_ALPHA: + raise ValueError('alpha must be greater than zero and at most one.') + return Parameter( + name=name + ' alpha', + value=float(alpha), + unit='dimensionless', + min=MINIMUM_ALPHA, + max=MAXIMUM_ALPHA, + ) + + @property + def scale(self) -> Parameter: + """ + Get the scale parameter. + + Returns + ------- + Parameter + The scale Parameter with unit ``x_unit * y_unit``. It is the integrated area of the + profile. + """ + return self._scale + + @scale.setter + def scale(self, value: Numeric) -> None: + """ + Parameters + ---------- + value : Numeric + New scale value (in current scale 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 scale parameter's bounds, + e.g. a negative value against its ``min=0``. + """ + self._set_bounded_parameter_value(self._scale, value, 'scale') + + @property + def alpha(self) -> Parameter: + """ + Get the form parameter of the Mittag-Leffler relaxation function. + + Returns + ------- + Parameter + The dimensionless alpha Parameter, bounded to (0, 1]. + """ + return self._alpha + + @alpha.setter + def alpha(self, value: Numeric) -> None: + """ + Parameters + ---------- + value : Numeric + New form parameter. Must lie in (0, 1]. + + Notes + ----- + A ``TypeError`` propagates from the shared value setter if *value* is not a numeric type, + and a ``ValueError`` propagates from it if *value* falls outside (0, 1]. + """ + self._set_bounded_parameter_value(self._alpha, value, 'alpha') + + @property + def width(self) -> Parameter: + """ + Get the width parameter (Mittag-Leffler relaxation rate). + + Returns + ------- + Parameter + The relaxation rate ``hbar / tau_R`` Parameter with unit ``x_unit``. + """ + return self._width + + @width.setter + def width(self, value: Numeric) -> None: + """ + Parameters + ---------- + value : Numeric + New relaxation rate in x_unit. Must be strictly positive. + + Raises + ------ + TypeError + If *value* is not a numeric type. + ValueError + If *value* is not positive, or violates the width parameter's bounds. + """ + if not isinstance(value, Numeric): + raise TypeError('width must be a number') + if float(value) <= 0: + raise ValueError('width must be positive') + self._set_bounded_parameter_value(self._width, value, 'width') + + @property + def damping(self) -> Parameter: + """ + Get the damping parameter (the diffusion damping ``hbar * D * q**2``). + + Returns + ------- + Parameter + The diffusion damping Parameter with unit ``x_unit``. + """ + return self._damping + + @damping.setter + def damping(self, value: Numeric) -> None: + """ + Parameters + ---------- + value : Numeric + New diffusion damping in x_unit. Must be strictly positive. + + Raises + ------ + TypeError + If *value* is not a numeric type. + ValueError + If *value* is not positive, or violates the damping parameter's bounds. + """ + if not isinstance(value, Numeric): + raise TypeError('damping must be a number') + if float(value) <= 0: + raise ValueError('damping must be positive') + self._set_bounded_parameter_value(self._damping, value, 'damping') + + def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndarray: + r""" + Evaluate the diffusion-damped Mittag-Leffler spectrum at x_vals. + + Eq. (42) of Hassani et al. is evaluated in the reduced variables $\omega = |x|/\Gamma$ and + $\epsilon' = \epsilon/\Gamma$, and the result divided by $\Gamma$, which turns the printed + $\tau_R = 1$ expression into the general one. 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 Mittag-Leffler spectrum values at x_vals. + """ + width = self._resolve_param_value(self._width, eval_unit) + damping = self._resolve_param_value(self._damping, eval_unit) + scale = self._resolve_param_value(self._scale, self._eval_area_unit(eval_unit)) + alpha = self._alpha.value + + # Reduced (tau_R = 1) variables, so that Eq. (42) applies verbatim. + omega = np.abs(x_vals) / width + epsilon = damping / width + + modulus_squared = omega**2 + epsilon**2 + # arg(epsilon + i|omega|); damping > 0 is enforced, so modulus_squared is never zero. + phase = alpha * np.arctan2(omega, epsilon) + numerator = ( + epsilon * modulus_squared ** (alpha / 2) + + omega * np.sin(phase) + + epsilon * np.cos(phase) + ) + denominator = modulus_squared * ( + (modulus_squared**alpha + 1) * modulus_squared ** (-alpha / 2) + 2 * np.cos(phase) + ) + return scale * numerator / (np.pi * width * denominator) + + def convert_x_unit(self, new_x_unit: str | sc.Unit) -> None: + """ + Convert x-axis parameters (width, damping) and scale to new_x_unit. + + The dimensionless alpha is unaffected. + + 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._width, self._damping], + area_param=self._scale, + ) + + def convert_y_unit(self, new_y_unit: str | sc.Unit) -> None: + """ + Convert the y-axis (output) unit by rescaling the scale parameter. + + The scale 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._scale) + + def __repr__(self) -> str: + """ + Return a string representation of the diffusion-damped Mittag-Leffler component. + + Returns + ------- + str + A string representation of the diffusion-damped Mittag-Leffler component. + """ + 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' scale = {self.scale},\n' + f' alpha = {self.alpha},\n' + f' width = {self.width},\n' + f' damping = {self.damping})' + ) diff --git a/src/easydynamics/sample_model/diffusion_model/__init__.py b/src/easydynamics/sample_model/diffusion_model/__init__.py index 778abc129..b8c079391 100644 --- a/src/easydynamics/sample_model/diffusion_model/__init__.py +++ b/src/easydynamics/sample_model/diffusion_model/__init__.py @@ -8,9 +8,13 @@ from easydynamics.sample_model.diffusion_model.jump_translational_diffusion import ( JumpTranslationalDiffusion, ) +from easydynamics.sample_model.diffusion_model.mittag_leffler_diffusion import ( + MittagLefflerDiffusion, +) __all__ = [ 'BrownianTranslationalDiffusion', 'DeltaLorentz', 'JumpTranslationalDiffusion', + 'MittagLefflerDiffusion', ] diff --git a/src/easydynamics/sample_model/diffusion_model/mittag_leffler_diffusion.py b/src/easydynamics/sample_model/diffusion_model/mittag_leffler_diffusion.py new file mode 100644 index 000000000..b324ef82b --- /dev/null +++ b/src/easydynamics/sample_model/diffusion_model/mittag_leffler_diffusion.py @@ -0,0 +1,1330 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + + +import numpy as np +import scipp as sc +from easyscience.variable import DescriptorNumber +from easyscience.variable import Parameter + +from easydynamics.sample_model.component_collection import ComponentCollection +from easydynamics.sample_model.components import DiffusionDampedMittagLeffler +from easydynamics.sample_model.components import Lorentzian +from easydynamics.sample_model.components.diffusion_damped_mittag_leffler import MAXIMUM_ALPHA +from easydynamics.sample_model.components.diffusion_damped_mittag_leffler import MINIMUM_ALPHA +from easydynamics.sample_model.diffusion_model.diffusion_model_base import DiffusionModelBase +from easydynamics.utils.fit_target import FitTarget +from easydynamics.utils.utils import CANONICAL_Q_UNIT +from easydynamics.utils.utils import Numeric +from easydynamics.utils.utils import Q_type +from easydynamics.utils.utils import angstrom +from easydynamics.utils.utils import convert_parameter_unit +from easydynamics.utils.utils import hbar +from easydynamics.utils.utils import verify_Q_index + +MINIMUM_WIDTH = 1e-10 # To avoid division by zero + + +class MittagLefflerDiffusion(DiffusionModelBase): + r""" + Multiscale relaxation model of Hassani et al.: Mittag-Leffler internal dynamics damped by + global translational diffusion. + + A. N. Hassani, L. Haris, M. Appel, T. Seydel, A. M. Stadler and G. R. Kneller, *Multiscale + relaxation dynamics and diffusion of myelin basic protein in solution studied by quasielastic + neutron scattering*, J. Chem. Phys. **156**, 025102 (2022), + [doi:10.1063/5.0077100](https://doi.org/10.1063/5.0077100). + + The intermediate scattering function is written as an internal relaxation, described by a + Mittag-Leffler function, times a global diffusion factor (their Eqs. 38-39), + + $$ F^{(+)}(t) = K e^{-\epsilon |t|} \left( EISF + (1 - EISF) E_\alpha(-(|t|/\tau_R)^\alpha) + \right), \qquad \epsilon = \hbar D q^2 $$ + + whose Fourier transform is their Eq. (41), + + $$ S^{(+)}(x) = K \left[ EISF \frac{1}{\pi} \frac{\epsilon}{x^2 + \epsilon^2} + + (1 - EISF) \tilde{\phi}^{(\epsilon)}_{ML}(|x|) \right]. $$ + + This model builds that sum at every Q as a + [`ComponentCollection`][easydynamics.sample_model.ComponentCollection] of two components: + + - a [`Lorentzian`][easydynamics.sample_model.Lorentzian] carrying the elastic term, with area + $K \cdot EISF$ and half width at half maximum $\epsilon = \hbar D q^2$. Note that the elastic + line is a Lorentzian and not a delta function, because global diffusion broadens it. + - a + [`DiffusionDampedMittagLeffler`][easydynamics.sample_model.DiffusionDampedMittagLeffler] + carrying the quasi-elastic term, with scale $K \cdot (1 - EISF)$, damping $\epsilon = \hbar D + q^2$, width $\hbar/\tau_R$ and form parameter $\alpha$. + + The diffusion coefficient $D$ is global: it is the single parameter that ties the Q values + together, exactly as in + [`BrownianTranslationalDiffusion`][easydynamics.sample_model.BrownianTranslationalDiffusion]. + In the paper $\tau_R$, $\alpha$ and the EISF are instead fitted independently at every Q (their + Fig. 5); pass ``allow_Q_variation={'A_0': True, 'relaxation_rate': True, 'alpha': True}`` to + reproduce that. Q is assumed to be in 1/angstrom and $D$ in m^2/s. + + Examples + -------- + **Creating a MittagLefflerDiffusion model with the paper's per-Q parameters** + + ```python + import numpy as np + import easydynamics as edyn + + Q = np.linspace(0.8, 1.8, 6) + model = edyn.MittagLefflerDiffusion( + scale=1.0, + diffusion_coefficient=3.3e-11, + A_0=0.05, + relaxation_rate=0.02, + alpha=0.85, + allow_Q_variation={'A_0': True, 'relaxation_rate': True, 'alpha': True}, + Q=Q, + ) + component_collections = model.get_component_collections() + ``` + + See also the tutorials. + """ + + def __init__( + self, + scale: Numeric = 1.0, + diffusion_coefficient: Numeric = 1.0, + A_0: Numeric = 0.0, + relaxation_rate: Numeric = 1.0, + alpha: Numeric = 1.0, + allow_Q_variation: dict | None = None, + Q: Q_type | None = None, + x_unit: str | sc.Unit = 'meV', + y_unit: str | sc.Unit = 'dimensionless', + name: str = 'MittagLefflerDiffusion', + display_name: str | None = None, + lorentzian_name: str = 'Elastic Lorentzian', + lorentzian_display_name: str | None = None, + mittag_leffler_name: str = 'Mittag-Leffler', + mittag_leffler_display_name: str | None = None, + unique_name: str | None = None, + ) -> None: + """ + Initialize a new MittagLefflerDiffusion model. + + Parameters + ---------- + scale : Numeric, default=1.0 + Scale factor K for the model. Must be non-negative. Its unit is ``x_unit * y_unit``. + diffusion_coefficient : Numeric, default=1.0 + Global translational diffusion coefficient D in m^2/s. Sets the diffusion damping + ``epsilon = hbar * D * Q**2`` shared by both components. Must be non-negative. + A_0 : Numeric, default=0.0 + Elastic incoherent structure factor (EISF), the elastic fraction of the intensity. Must + lie in [0, 1]. The paper finds it close to zero over most of its Q range. + relaxation_rate : Numeric, default=1.0 + Mittag-Leffler relaxation rate ``hbar / tau_R`` in x_unit. Must be strictly positive. + alpha : Numeric, default=1.0 + Form parameter of the Mittag-Leffler relaxation function. Must lie in (0, 1]; alpha=1 + reduces the quasi-elastic term to a Lorentzian. + allow_Q_variation : dict | None, default=None + Dict describing which of ``'A_0'``, ``'relaxation_rate'`` and ``'alpha'`` are free at + every Q instead of shared, with boolean values. If None, none of them vary with Q. The + paper's analysis corresponds to all three being True. + Q : Q_type | None, default=None + Q values for the model in 1/angstrom. If None, Q is not set. + x_unit : str | sc.Unit, default='meV' + Unit of the x-axis (energy). Must be convertible to meV. + y_unit : str | sc.Unit, default='dimensionless' + Unit of the model output (intensity). Determines scale.unit = x_unit * y_unit. + name : str, default='MittagLefflerDiffusion' + Name of the diffusion model. + display_name : str | None, default=None + Display name of the diffusion model. + lorentzian_name : str, default='Elastic Lorentzian' + Name of the elastic Lorentzian component. + lorentzian_display_name : str | None, default=None + Display name of the elastic Lorentzian component. If None, it falls back to + *lorentzian_name*. + mittag_leffler_name : str, default='Mittag-Leffler' + Name of the Mittag-Leffler component. + mittag_leffler_display_name : str | None, default=None + Display name of the Mittag-Leffler component. If None, it falls back to + *mittag_leffler_name*. + unique_name : str | None, default=None + Unique name of the diffusion model. If None, a unique name will be generated. + + Raises + ------ + TypeError + If mittag_leffler_name is not a string, or if mittag_leffler_display_name is not a + string or None. + """ + super().__init__( + scale=scale, + x_unit=x_unit, + y_unit=y_unit, + Q=Q, + lorentzian_name=lorentzian_name, + lorentzian_display_name=lorentzian_display_name, + name=name, + display_name=display_name, + unique_name=unique_name, + ) + + # -------------------------------------------------------------- + # Parameters + # -------------------------------------------------------------- + self._hbar = hbar + self._angstrom = angstrom + + self._diffusion_coefficient = self._create_diffusion_coefficient_parameter( + diffusion_coefficient + ) + self._A_0, self._A_1 = self._create_A0_A1_parameters(A_0) + self._relaxation_rate = self._create_relaxation_rate_parameter(relaxation_rate) + self._alpha = self._create_alpha_parameter(alpha) + + # -------------------------------------------------------------- + # names + # -------------------------------------------------------------- + if not isinstance(mittag_leffler_name, str): + raise TypeError('mittag_leffler_name must be a string.') + + if mittag_leffler_display_name is None: + mittag_leffler_display_name = mittag_leffler_name + + if not isinstance(mittag_leffler_display_name, str): + raise TypeError('mittag_leffler_display_name must be a string or None.') + + self._mittag_leffler_name = mittag_leffler_name + self._mittag_leffler_display_name = mittag_leffler_display_name + + # -------------------------------------------------------------- + # Q variation + # -------------------------------------------------------------- + self._allow_Q_variation = self._create_Q_variation_dict(allow_Q_variation) + + # create_component_collections creates the per-Q parameter lists itself, so the + # components it builds are backed by the very parameters stored in those lists. + self.create_component_collections() + + # ------------------------------------------------------------------ + # Properties + # ------------------------------------------------------------------ + + @property + def diffusion_coefficient(self) -> Parameter: + """ + Get the global diffusion coefficient parameter D. + + Returns + ------- + Parameter + Diffusion coefficient D in m^2/s. + """ + return self._diffusion_coefficient + + @diffusion_coefficient.setter + def diffusion_coefficient(self, diffusion_coefficient: Numeric) -> None: + """ + Set the global diffusion coefficient parameter D. + + Parameters + ---------- + diffusion_coefficient : Numeric + The new value for D in m^2/s. + + Raises + ------ + TypeError + If diffusion_coefficient is not a number. + ValueError + If diffusion_coefficient is negative. + """ + if not isinstance(diffusion_coefficient, Numeric): + raise TypeError('diffusion_coefficient must be a number.') + if float(diffusion_coefficient) < 0: + raise ValueError('diffusion_coefficient must be non-negative.') + self._diffusion_coefficient.value = float(diffusion_coefficient) + + @property + def A_0(self) -> Parameter: + """ + Get the elastic fraction parameter A_0 (the EISF). + + Returns + ------- + Parameter + The dimensionless A_0 parameter, bounded to [0, 1]. + """ + return self._A_0 + + @A_0.setter + def A_0(self, A_0: Numeric) -> None: + """ + Set the elastic fraction parameter A_0. + + Parameters + ---------- + A_0 : Numeric + The new value for A_0. Must lie in [0, 1]. + + Raises + ------ + TypeError + If A_0 is not a number. + ValueError + If A_0 is not between 0 and 1. + """ + if not isinstance(A_0, Numeric): + raise TypeError('A_0 must be a number.') + if float(A_0) < 0 or float(A_0) > 1: + raise ValueError('A_0 must be between 0 and 1.') + self._A_0.value = float(A_0) + + @property + def A_1(self) -> Parameter: + """ + Get the quasi-elastic fraction parameter A_1 = 1 - A_0. + + Returns + ------- + Parameter + The dependent A_1 parameter. + """ + return self._A_1 + + @A_1.setter + def A_1(self, _A_1: Numeric) -> None: + """ + Reject assignment to the dependent A_1 parameter. + + Parameters + ---------- + _A_1 : Numeric + Ignored. + + Raises + ------ + AttributeError + Always; A_1 is derived from A_0 and must be changed through A_0. + """ + raise AttributeError('A_1 is derived from A_0 and cannot be set directly. Set A_0.') + + @property + def relaxation_rate(self) -> Parameter: + """ + Get the Mittag-Leffler relaxation rate parameter hbar / tau_R. + + Returns + ------- + Parameter + The relaxation rate parameter with unit ``x_unit``. + """ + return self._relaxation_rate + + @relaxation_rate.setter + def relaxation_rate(self, relaxation_rate: Numeric) -> None: + """ + Set the Mittag-Leffler relaxation rate parameter. + + Parameters + ---------- + relaxation_rate : Numeric + The new relaxation rate in x_unit. Must be strictly positive. + + Raises + ------ + TypeError + If relaxation_rate is not a number. + ValueError + If relaxation_rate is smaller than the minimum width. + """ + if not isinstance(relaxation_rate, Numeric): + raise TypeError('relaxation_rate must be a number.') + if float(relaxation_rate) < MINIMUM_WIDTH: + raise ValueError(f'relaxation_rate must be at least {MINIMUM_WIDTH}.') + self._relaxation_rate.value = float(relaxation_rate) + + @property + def alpha(self) -> Parameter: + """ + Get the Mittag-Leffler form parameter alpha. + + Returns + ------- + Parameter + The dimensionless alpha parameter, bounded to (0, 1]. + """ + return self._alpha + + @alpha.setter + def alpha(self, alpha: Numeric) -> None: + """ + Set the Mittag-Leffler form parameter alpha. + + Parameters + ---------- + alpha : Numeric + The new form parameter. Must lie in (0, 1]. + + Raises + ------ + TypeError + If alpha is not a number. + ValueError + If alpha does not lie in (0, 1]. + """ + if not isinstance(alpha, Numeric): + raise TypeError('alpha must be a number.') + if not MINIMUM_ALPHA <= float(alpha) <= MAXIMUM_ALPHA: + raise ValueError('alpha must be greater than zero and at most one.') + self._alpha.value = float(alpha) + + @property + def mittag_leffler_name(self) -> str: + """ + Get the name of the Mittag-Leffler component. + + Returns + ------- + str + Name of the Mittag-Leffler component. + """ + return self._mittag_leffler_name + + @mittag_leffler_name.setter + def mittag_leffler_name(self, mittag_leffler_name: str) -> None: + """ + Set the name of the Mittag-Leffler component. + + Parameters + ---------- + mittag_leffler_name : str + The new name for the Mittag-Leffler component. + + Raises + ------ + TypeError + If mittag_leffler_name is not a string. + """ + if not isinstance(mittag_leffler_name, str): + raise TypeError('mittag_leffler_name must be a string.') + self._mittag_leffler_name = mittag_leffler_name + + @property + def mittag_leffler_display_name(self) -> str | None: + """ + Get the display name of the Mittag-Leffler component. + + Returns + ------- + str | None + Display name of the Mittag-Leffler component, or None if not set. + """ + return self._mittag_leffler_display_name + + @mittag_leffler_display_name.setter + def mittag_leffler_display_name(self, mittag_leffler_display_name: str | None) -> None: + """ + Set the display name of the Mittag-Leffler component. + + Parameters + ---------- + mittag_leffler_display_name : str | None + The new display name for the Mittag-Leffler component. + + Raises + ------ + TypeError + If mittag_leffler_display_name is not a string or None. + """ + if not isinstance(mittag_leffler_display_name, (str, type(None))): + raise TypeError('mittag_leffler_display_name must be a string or None.') + self._mittag_leffler_display_name = mittag_leffler_display_name + + # ------------------------------------------------------------------ + # Other methods + # ------------------------------------------------------------------ + + def calculate_width(self, Q: Q_type = None) -> np.ndarray: + """ + Calculate the diffusion damping epsilon = hbar * D * Q**2. + + This is both the half width at half maximum of the elastic Lorentzian and the ``damping`` + of the Mittag-Leffler component; the two share it by construction. + + Parameters + ---------- + Q : Q_type, default=None + Scattering vector in 1/angstrom. If None, the Q stored in the model is used. + + Returns + ------- + np.ndarray + Damping values in the unit of the model (e.g. meV). + """ + Q = self._ensure_Q(Q) + + unit_conversion_factor = self._hbar * self.diffusion_coefficient / (self._angstrom**2) + unit_conversion_factor.convert_unit(self.x_unit) + return Q**2 * unit_conversion_factor.value + + def calculate_relaxation_rate(self, Q: Q_type = None) -> np.ndarray: + """ + Calculate the Mittag-Leffler relaxation rate hbar / tau_R at each Q. + + If the relaxation rate is allowed to vary with Q, the requested Q values are matched + against the Q stored in the model and the corresponding per-Q rates are returned. Otherwise + the shared rate is returned for every Q. + + Parameters + ---------- + Q : Q_type, default=None + Scattering vector in 1/angstrom. If None, the Q stored in the model is used. + + Returns + ------- + np.ndarray + Relaxation rates in the unit of the model (e.g. meV). + + Raises + ------ + ValueError + If Q-variation is enabled but Q has not been set on the model yet. + """ + Q = self._ensure_Q(Q) + + if self._allow_Q_variation['relaxation_rate'] is True: + if not self._relaxation_rate_list: + raise ValueError( + 'Relaxation rate Q-variation list is empty. ' + 'Set Q before calling calculate_relaxation_rate.' + ) + indices = self._match_Q_indices(Q) + return np.array([self._relaxation_rate_list[i].value for i in indices]) + + return self.relaxation_rate.value * np.ones_like(Q) + + def calculate_alpha(self, Q: Q_type = None) -> np.ndarray: + """ + Calculate the Mittag-Leffler form parameter alpha at each Q. + + If alpha is allowed to vary with Q, the requested Q values are matched against the Q stored + in the model and the corresponding per-Q values are returned. Otherwise the shared alpha is + returned for every Q. + + Parameters + ---------- + Q : Q_type, default=None + Scattering vector in 1/angstrom. If None, the Q stored in the model is used. + + Returns + ------- + np.ndarray + Form parameters (dimensionless). + + Raises + ------ + ValueError + If Q-variation is enabled but Q has not been set on the model yet. + """ + Q = self._ensure_Q(Q) + + if self._allow_Q_variation['alpha'] is True: + if not self._alpha_list: + raise ValueError( + 'Alpha Q-variation list is empty. Set Q before calling calculate_alpha.' + ) + indices = self._match_Q_indices(Q) + return np.array([self._alpha_list[i].value for i in indices]) + + return self.alpha.value * np.ones_like(Q) + + def calculate_EISF(self, Q: Q_type = None) -> np.ndarray: + """ + Calculate the Elastic Incoherent Structure Factor (EISF), i.e. A_0. + + Parameters + ---------- + Q : Q_type, default=None + Scattering vector in 1/angstrom. + + Returns + ------- + np.ndarray + EISF values (dimensionless). + """ + Q = self._ensure_Q(Q) + + if self._allow_Q_variation['A_0'] is True: + indices = self._match_Q_indices(Q) + return np.array([self._A_0_list[i].value for i in indices]) + + return self.A_0.value * np.ones_like(Q) + + def calculate_QISF(self, Q: Q_type = None) -> np.ndarray: + """ + Calculate the Quasi-Elastic Incoherent Structure Factor (QISF), i.e. A_1 = 1 - A_0. + + Parameters + ---------- + Q : Q_type, default=None + Scattering vector in 1/angstrom. + + Returns + ------- + np.ndarray + QISF values (dimensionless). + """ + Q = self._ensure_Q(Q) + + if self._allow_Q_variation['A_0'] is True: + indices = self._match_Q_indices(Q) + return np.array([self._A_1_list[i].value for i in indices]) + + return self.A_1.value * np.ones_like(Q) + + def calculate_relaxation_rate_spectrum(self, rate: np.ndarray, Q: Q_type = None) -> np.ndarray: + r""" + Calculate the relaxation rate spectrum p(lambda) of Eq. (37) of Hassani et al. + + $$ p_{ML}(\lambda) = \frac{\sin(\pi\alpha)} {\pi\lambda((\lambda\tau_R)^{-\alpha} + + (\lambda\tau_R)^{\alpha} + 2\cos(\pi\alpha))} $$ + + This is the distribution of exponential relaxation rates whose superposition gives the + Mittag-Leffler relaxation function, and is what the paper plots in its Fig. 6. It is a + property of the internal dynamics alone, so the diffusion damping does not enter. + + Parameters + ---------- + rate : np.ndarray + Relaxation rates lambda at which to evaluate the spectrum, expressed as energies in + x_unit. Must be strictly positive. + Q : Q_type, default=None + Scattering vector in 1/angstrom. If None, the Q stored in the model is used. + + Returns + ------- + np.ndarray + Array of shape ``(len(Q), len(rate))`` holding p(lambda) for each Q, with unit + ``1/x_unit``. + + Raises + ------ + ValueError + If any requested rate is not strictly positive. + """ + rate = np.atleast_1d(np.asarray(rate, dtype=float)) + if np.any(rate <= 0): + raise ValueError('rate must be strictly positive.') + + Q = self._ensure_Q(Q) + alpha = self.calculate_alpha(Q)[:, np.newaxis] + # lambda * tau_R, with tau_R expressed through the relaxation rate hbar/tau_R + reduced_rate = rate[np.newaxis, :] / self.calculate_relaxation_rate(Q)[:, np.newaxis] + + numerator = np.sin(np.pi * alpha) + denominator = ( + np.pi + * rate[np.newaxis, :] + * (reduced_rate**-alpha + reduced_rate**alpha + 2 * np.cos(np.pi * alpha)) + ) + return numerator / denominator + + def calculate_energy_barrier_distribution( + self, barrier_height: np.ndarray, Q: Q_type = None + ) -> np.ndarray: + r""" + Calculate the energy barrier distribution P(h) of Eq. (48) of Hassani et al. + + $$ P_{ML}(h) = \frac{2h\sin(\pi\alpha)} {\pi(e^{-\alpha h^2} + e^{\alpha h^2} + + 2\cos(\pi\alpha))} $$ + + where $h = \Delta E / (k_B T)$ is the dimensionless barrier height of Zwanzig's rough + harmonic potential. This is what the paper plots in the right panel of its Fig. 7. As + $\alpha \to 1$ it collapses onto $\delta(h)$, a smooth potential; as $\alpha \to 0$ it + broadens to include arbitrarily high barriers. + + Parameters + ---------- + barrier_height : np.ndarray + Dimensionless barrier heights h at which to evaluate the distribution. + Q : Q_type, default=None + Scattering vector in 1/angstrom. If None, the Q stored in the model is used. + + Returns + ------- + np.ndarray + Array of shape ``(len(Q), len(barrier_height))`` holding P(h) for each Q. + """ + barrier_height = np.atleast_1d(np.asarray(barrier_height, dtype=float)) + + Q = self._ensure_Q(Q) + alpha = self.calculate_alpha(Q)[:, np.newaxis] + h = barrier_height[np.newaxis, :] + + # Factor exp(alpha * h**2) out of the denominator, so tall barriers underflow smoothly to + # zero instead of overflowing exp() and leaving inf/inf behind. + damping = np.exp(-alpha * h**2) + numerator = 2 * h * np.sin(np.pi * alpha) * damping + denominator = np.pi * (damping**2 + 1 + 2 * np.cos(np.pi * alpha) * damping) + return numerator / denominator + + def create_component_collections(self) -> list[ComponentCollection]: + r""" + Create ComponentCollections for the MittagLefflerDiffusion model at the given Q values. + + Each collection holds the elastic Lorentzian (area $K \cdot EISF$, HWHM $\hbar D q^2$) and + the Mittag-Leffler component (scale $K \cdot (1 - EISF)$, damping $\hbar D q^2$). The per-Q + parameter lists are recreated here so the built components are backed by the very + parameters stored in the lists, keeping ``calculate_*`` in sync with the components. The + created collections are installed on the model, so the returned list is the live one. + + Returns + ------- + list[ComponentCollection] + List of ComponentCollections with a Lorentzian and a Mittag-Leffler component for each + Q value. + """ + if self.Q is None: + self._A_0_list = [] + self._A_1_list = [] + self._relaxation_rate_list = [] + self._alpha_list = [] + self._component_collections = [] + return self._component_collections + + Q = self.Q.values + + if self._allow_Q_variation['A_0'] is True: + self._A_0_list, self._A_1_list = self._create_A0_A1_parameter_lists() + else: + self._A_0_list = [] + self._A_1_list = [] + + if self._allow_Q_variation['relaxation_rate'] is True: + self._relaxation_rate_list = self._create_relaxation_rate_parameter_list() + else: + self._relaxation_rate_list = [] + + if self._allow_Q_variation['alpha'] is True: + self._alpha_list = self._create_alpha_parameter_list() + else: + self._alpha_list = [] + + component_collection_list = [None] * len(Q) + for i, Q_value in enumerate(Q): + component_collection_list[i] = ComponentCollection( + name=f'{self.name}_Q{Q_value:.2f}', + display_name=f'{self.display_name}_Q{Q_value:.2f}', + x_unit=self.x_unit, + y_unit=self.y_unit, + ) + + # easyscience propagates inf bounds through arithmetic, producing inf/inf=nan + # as a transient intermediate. Python's min/max ignore nan so the final bounds + # are correct; suppress the spurious numpy RuntimeWarning. + with np.errstate(invalid='ignore'): + component_collection_list[i].append_component( + self._create_lorentzian_component(Q_value, i) + ) + component_collection_list[i].append_component( + self._create_mittag_leffler_component(Q_value, i) + ) + + self._component_collections = component_collection_list + return self._component_collections + + def get_fit_targets(self) -> list[FitTarget]: + """ + Get the fittable predictions of the MittagLefflerDiffusion model as FitTargets. + + The model predicts three Q-dependent quantities: ``'width'`` (the shared damping ``hbar * D + * Q**2``), ``'area'`` (the Mittag-Leffler weight ``scale * QISF(Q)``) and + ``'elastic_area'`` (the elastic Lorentzian weight ``scale * EISF(Q)``). The base class + implementation is replaced rather than extended, because here it is the *elastic* line that + is a Lorentzian, so the base's ``'area'`` key would point at the wrong component. + + Returns + ------- + list[FitTarget] + The fittable predictions of this model. + """ + return [ + FitTarget( + name='area', + dataset_key=f'{self.mittag_leffler_name} scale', + function=lambda Q, model=self, **_: model.calculate_QISF(Q) * model.scale.value, + label=f'{self.display_name} area', + x_unit=CANONICAL_Q_UNIT, + y_unit=str(self.scale.unit), + ), + FitTarget( + name='width', + dataset_key=f'{self.lorentzian_name} width', + function=lambda Q, model=self, **_: model.calculate_width(Q), + label=f'{self.display_name} width', + x_unit=CANONICAL_Q_UNIT, + y_unit=str(self.x_unit), + ), + FitTarget( + name='elastic_area', + dataset_key=f'{self.lorentzian_name} area', + function=lambda Q, model=self, **_: model.calculate_EISF(Q) * model.scale.value, + label=f'{self.display_name} elastic_area', + x_unit=CANONICAL_Q_UNIT, + y_unit=str(self.scale.unit), + ), + ] + + def get_global_variables(self) -> list[Parameter]: + """ + Get all global variables from the diffusion model. + + Returns + ------- + list[Parameter] + A list of all global variables from the diffusion model. + """ + variables = [self.scale, self.diffusion_coefficient] + + if self._allow_Q_variation['A_0'] is False: + variables.append(self.A_0) + variables.append(self.A_1) + + if self._allow_Q_variation['relaxation_rate'] is False: + variables.append(self.relaxation_rate) + + if self._allow_Q_variation['alpha'] is False: + variables.append(self.alpha) + + return variables + + def get_independent_variables(self, Q_index: int | None = None) -> list[Parameter]: + """ + Get the independent variables from the diffusion model. + + The per-Q relaxation rate and alpha parameters are the components' own parameters, so they + are reached through the component collections; only the per-Q A_0/A_1 pairs, which the + component areas merely depend on, are listed here. + + Parameters + ---------- + Q_index : int | None, default=None + The index of the Q value for which to get the independent variables. If None, + independent variables for all Q values will be included. + + Returns + ------- + list[Parameter] + List of independent variables in the model. + """ + verify_Q_index(Q_index=Q_index, Q=self.Q, allow_none=True) + + variables = [] + if self._allow_Q_variation['A_0'] is True: + if Q_index is None: + variables.extend(self._A_0_list) + variables.extend(self._A_1_list) + else: + variables.append(self._A_0_list[Q_index]) + variables.append(self._A_1_list[Q_index]) + + return variables + + def get_all_variables(self, Q_index: int | None = None) -> list[DescriptorNumber]: + """ + Get a list of all variables (Parameters and Descriptors) in the model. + + Parameters + ---------- + Q_index : int | None, default=None + The index of the Q value for which to get the variables. If None, variables for all Q + values will be included. + + Returns + ------- + list[DescriptorNumber] + List of all variables in the model. + """ + verify_Q_index(Q_index=Q_index, Q=self.Q, allow_none=True) + + variables = self.get_global_variables() + variables.extend(self.get_independent_variables(Q_index=Q_index)) + + if Q_index is None: + for component_collection in self._component_collections: + variables.extend(component_collection.get_all_variables()) + else: + variables.extend(self._component_collections[Q_index].get_all_variables()) + + return variables + + # ------------------------------------------------------------------ + # Private methods for init + # ------------------------------------------------------------------ + + def _create_Q_variation_dict(self, allow_Q_variation: dict | None) -> dict: + """ + Create the allow_Q_variation dict, ensuring it has the correct keys and default values. + + Parameters + ---------- + allow_Q_variation : dict | None + Dict describing whether to allow Q variation of A_0, relaxation_rate and alpha. + + Raises + ------ + TypeError + If allow_Q_variation is not a dict or None. + ValueError + If allow_Q_variation contains unknown keys. + + Returns + ------- + dict + A dict with keys 'A_0', 'relaxation_rate' and 'alpha'. + """ + allow_Q_variation_default = { + 'A_0': False, + 'relaxation_rate': False, + 'alpha': False, + } + allowed_keys = set(allow_Q_variation_default) + + if allow_Q_variation is None: + allow_Q_variation = {} + if not isinstance(allow_Q_variation, dict): + raise TypeError('allow_Q_variation must be a dict or None.') + + unknown_keys = set(allow_Q_variation) - allowed_keys + if unknown_keys: + raise ValueError(f'Unknown keys in allow_Q_variation: {unknown_keys}') + + return {**allow_Q_variation_default, **allow_Q_variation} + + @staticmethod + def _create_diffusion_coefficient_parameter(diffusion_coefficient: Numeric) -> Parameter: + """ + Create the global diffusion coefficient parameter. + + Parameters + ---------- + diffusion_coefficient : Numeric + The value for D in m^2/s. + + Raises + ------ + TypeError + If diffusion_coefficient is not a number. + ValueError + If diffusion_coefficient is negative. + + Returns + ------- + Parameter + The created diffusion coefficient parameter. + """ + if not isinstance(diffusion_coefficient, Numeric): + raise TypeError('diffusion_coefficient must be a number.') + if float(diffusion_coefficient) < 0: + raise ValueError('diffusion_coefficient must be non-negative.') + return Parameter( + name='diffusion_coefficient', + value=float(diffusion_coefficient), + fixed=False, + unit='m**2/s', + min=0.0, + ) + + @staticmethod + def _create_A0_A1_parameters(A_0: Numeric) -> tuple[Parameter, Parameter]: + """ + Create the shared A_0 and A_1 parameters. + + Parameters + ---------- + A_0 : Numeric + The value for the A_0 parameter. + + Raises + ------ + TypeError + If A_0 is not a number. + ValueError + If A_0 is not between 0 and 1. + + Returns + ------- + tuple[Parameter, Parameter] + A tuple containing the A_0 and A_1 parameters. + """ + if not isinstance(A_0, Numeric): + raise TypeError('A_0 must be a number.') + if float(A_0) < 0 or float(A_0) > 1: + raise ValueError('A_0 must be between 0 and 1.') + + A_0_parameter = Parameter(name='A_0', value=float(A_0), fixed=False, min=0.0, max=1.0) + A_1_parameter = Parameter.from_dependency( + name='A_1', + dependency_expression='1 - A_0', + dependency_map={'A_0': A_0_parameter}, + ) + return A_0_parameter, A_1_parameter + + def _create_relaxation_rate_parameter(self, relaxation_rate: Numeric) -> Parameter: + """ + Create the shared relaxation rate parameter. + + Parameters + ---------- + relaxation_rate : Numeric + The value for the relaxation rate in x_unit. + + Raises + ------ + TypeError + If relaxation_rate is not a number. + ValueError + If relaxation_rate is less than the minimum width. + + Returns + ------- + Parameter + The created relaxation rate parameter. + """ + if not isinstance(relaxation_rate, Numeric): + raise TypeError('relaxation_rate must be a number.') + if float(relaxation_rate) < MINIMUM_WIDTH: + raise ValueError(f'relaxation_rate must be at least {MINIMUM_WIDTH}.') + + return Parameter( + name='relaxation_rate', + value=float(relaxation_rate), + fixed=False, + min=MINIMUM_WIDTH, + unit=self.x_unit, + ) + + @staticmethod + def _create_alpha_parameter(alpha: Numeric) -> Parameter: + """ + Create the shared form parameter alpha. + + Parameters + ---------- + alpha : Numeric + The value for the form parameter. + + Raises + ------ + TypeError + If alpha is not a number. + ValueError + If alpha does not lie in (0, 1]. + + Returns + ------- + Parameter + The created alpha parameter. + """ + if not isinstance(alpha, Numeric): + raise TypeError('alpha must be a number.') + if not MINIMUM_ALPHA <= float(alpha) <= MAXIMUM_ALPHA: + raise ValueError('alpha must be greater than zero and at most one.') + + return Parameter( + name='alpha', + value=float(alpha), + fixed=False, + min=MINIMUM_ALPHA, + max=MAXIMUM_ALPHA, + unit='dimensionless', + ) + + def _create_A0_A1_parameter_lists(self) -> tuple[list[Parameter], list[Parameter]]: + """ + Create per-Q A_0 and A_1 parameters, seeded from the shared A_0. + + Returns + ------- + tuple[list[Parameter], list[Parameter]] + The per-Q A_0 parameters and the per-Q A_1 parameters derived from them. + """ + A_0_list = [] + A_1_list = [] + for _ in range(len(self.Q)): + # The per-Q amplitudes carry the model name so they do not collide with other + # models' parameters. The name is the same at every Q on purpose: parameters are + # tracked across Q by name (unique within a Q, shared across Q). + a_0 = Parameter( + name=f'{self.name} A_0', + display_name='A_0', + value=float(self.A_0.value), + fixed=False, + min=0.0, + max=1.0, + ) + a_1 = Parameter.from_dependency( + name=f'{self.name} A_1', + dependency_expression='1 - A_0', + dependency_map={'A_0': a_0}, + ) + A_0_list.append(a_0) + A_1_list.append(a_1) + + return A_0_list, A_1_list + + def _create_relaxation_rate_parameter_list(self) -> list[Parameter]: + """ + Create per-Q relaxation rate parameters, seeded from the shared relaxation rate. + + Returns + ------- + list[Parameter] + The per-Q relaxation rate parameters, named after the Mittag-Leffler component's own + width parameter so they slot straight into it. + """ + return [ + Parameter( + name=f'{self.mittag_leffler_name} width', + value=float(self.relaxation_rate.value), + fixed=False, + min=MINIMUM_WIDTH, + unit=self.x_unit, + ) + for _ in range(len(self.Q)) + ] + + def _create_alpha_parameter_list(self) -> list[Parameter]: + """ + Create per-Q alpha parameters, seeded from the shared alpha. + + Returns + ------- + list[Parameter] + The per-Q form parameters, named after the Mittag-Leffler component's own alpha + parameter so they slot straight into it. + """ + return [ + Parameter( + name=f'{self.mittag_leffler_name} alpha', + value=float(self.alpha.value), + fixed=False, + min=MINIMUM_ALPHA, + max=MAXIMUM_ALPHA, + unit='dimensionless', + ) + for _ in range(len(self.Q)) + ] + + # ------------------------------------------------------------------ + # Private methods + # ------------------------------------------------------------------ + + def _create_lorentzian_component(self, Q_value: float, Q_index: int) -> Lorentzian: + """ + Build the elastic Lorentzian for one Q value. + + Its width is the diffusion damping ``hbar * D * Q**2`` and its area is ``scale * A_0``; + both are always dependent parameters, since the elastic line is entirely determined by the + global diffusion coefficient and the elastic fraction. + + Parameters + ---------- + Q_value : float + Scattering vector in 1/angstrom. + Q_index : int + Index of this Q value in the stored Q. + + Returns + ------- + Lorentzian + The configured elastic Lorentzian component. + """ + component = Lorentzian( + name=self.lorentzian_name, + display_name=self.lorentzian_display_name, + x_unit=self.x_unit, + y_unit=self.y_unit, + ) + component.width.make_dependent_on( + dependency_expression=self._write_damping_dependency_expression(Q_value), + dependency_map=self._write_damping_dependency_map_expression(), + desired_unit=self.x_unit, + ) + component.area.make_dependent_on( + dependency_expression='scale * A_0', + dependency_map=self._write_amplitude_dependency_map_expression(Q_index, elastic=True), + ) + return component + + def _create_mittag_leffler_component( + self, Q_value: float, Q_index: int + ) -> DiffusionDampedMittagLeffler: + """ + Build the Mittag-Leffler component for one Q value. + + Its damping is the diffusion damping ``hbar * D * Q**2`` and its scale is ``scale * A_1``, + both always dependent. Its width and alpha are either the per-Q parameters from the + corresponding lists, or made dependent on the shared parameters. + + Parameters + ---------- + Q_value : float + Scattering vector in 1/angstrom. + Q_index : int + Index of this Q value in the stored Q. + + Returns + ------- + DiffusionDampedMittagLeffler + The configured Mittag-Leffler component. + """ + component = DiffusionDampedMittagLeffler( + name=self.mittag_leffler_name, + display_name=self.mittag_leffler_display_name, + x_unit=self.x_unit, + y_unit=self.y_unit, + ) + + if self._allow_Q_variation['relaxation_rate'] is True: + component._width = self._relaxation_rate_list[Q_index] # ruff: ignore[private-member-access] + else: + component.width.make_dependent_on( + dependency_expression='relaxation_rate', + dependency_map={'relaxation_rate': self.relaxation_rate}, + desired_unit=self.x_unit, + ) + + if self._allow_Q_variation['alpha'] is True: + component._alpha = self._alpha_list[Q_index] # ruff: ignore[private-member-access] + else: + component.alpha.make_dependent_on( + dependency_expression='alpha', + dependency_map={'alpha': self.alpha}, + desired_unit='dimensionless', + ) + + component.damping.make_dependent_on( + dependency_expression=self._write_damping_dependency_expression(Q_value), + dependency_map=self._write_damping_dependency_map_expression(), + desired_unit=self.x_unit, + ) + component.scale.make_dependent_on( + dependency_expression='scale * A_1', + dependency_map=self._write_amplitude_dependency_map_expression(Q_index, elastic=False), + ) + return component + + def _on_Q_change(self) -> None: + """ + Handle changes to the Q values. + + Rebuilds the component collections; the per-Q parameter lists are recreated inside + ``create_component_collections``. + """ + self.create_component_collections() + + def _convert_extra_x_unit_parameters(self, unit_str: str) -> None: + """ + Convert the shared relaxation rate template to the new x-axis unit. + + The per-Q relaxation rate list (when Q-variation is enabled) holds the very Parameter + objects used by the components, so those are converted in place with the collections. + + Parameters + ---------- + unit_str : str + The new x-axis unit as a string. + """ + convert_parameter_unit(self._relaxation_rate, unit_str) + + def _write_damping_dependency_expression(self, Q: float) -> str: + """ + Write the dependency expression for the diffusion damping ``hbar * D * Q**2``. + + Parameters + ---------- + Q : float + Scattering vector in 1/angstrom. + + Raises + ------ + TypeError + If Q is not a float. + + Returns + ------- + str + Dependency expression for the damping. + """ + if not isinstance(Q, (float)): + raise TypeError('Q must be a float.') + + # Q is given as a float, so we need to add the units + return f'hbar * D * {Q}**2 * 1/(angstrom**2)' + + def _write_damping_dependency_map_expression(self) -> dict[str, DescriptorNumber]: + """ + Write the dependency map for the diffusion damping. + + Returns + ------- + dict[str, DescriptorNumber] + Dependency map for the damping. + """ + return { + 'D': self.diffusion_coefficient, + 'hbar': self._hbar, + 'angstrom': self._angstrom, + } + + def _write_amplitude_dependency_map_expression( + self, Q_index: int, elastic: bool + ) -> dict[str, DescriptorNumber]: + """ + Write the dependency map for a component's amplitude. + + Parameters + ---------- + Q_index : int + Index of the Q value, used to pick the per-Q amplitude when Q-variation is enabled. + elastic : bool + True for the elastic Lorentzian's area (``scale * A_0``), False for the Mittag-Leffler + component's scale (``scale * A_1``). + + Returns + ------- + dict[str, DescriptorNumber] + Dependency map for the amplitude. + """ + if self._allow_Q_variation['A_0'] is True: + amplitude = self._A_0_list[Q_index] if elastic else self._A_1_list[Q_index] + else: + amplitude = self.A_0 if elastic else self.A_1 + + return {'scale': self.scale, 'A_0' if elastic else 'A_1': amplitude} + + # ------------------------------------------------------------------ + # dunder methods + # ------------------------------------------------------------------ + + def __repr__(self) -> str: + """ + String representation of the MittagLefflerDiffusion model. + + Returns + ------- + str + String representation of the MittagLefflerDiffusion model. + """ + return ( + f'MittagLefflerDiffusion(display_name={self.display_name}, ' + f'x_unit={self.x_unit}, y_unit={self.y_unit}, \n' + f' diffusion_coefficient={self.diffusion_coefficient}, \n' + f' A_0={self.A_0}, A_1={self.A_1}, \n' + f' relaxation_rate={self.relaxation_rate}, \n' + f' alpha={self.alpha}, \n' + f' scale={self.scale})' + ) diff --git a/tests/unit/easydynamics/sample_model/components/test_diffusion_damped_mittag_leffler.py b/tests/unit/easydynamics/sample_model/components/test_diffusion_damped_mittag_leffler.py new file mode 100644 index 000000000..96acf85c7 --- /dev/null +++ b/tests/unit/easydynamics/sample_model/components/test_diffusion_damped_mittag_leffler.py @@ -0,0 +1,400 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +from copy import copy + +import numpy as np +import pytest +import scipp as sc +from easyscience.variable import Parameter +from scipp import UnitError +from scipy.integrate import quad + +from easydynamics.sample_model import DiffusionDampedMittagLeffler + + +def laplace_reference(x, scale, alpha, width, damping): + r""" + Reference lineshape built from the Laplace transform of the Mittag-Leffler function. + + Eqs. (32) and (40) of Hassani et al. give the spectrum as ``(1/pi) Re{phi_hat(damping + + i|x|)}``, with ``phi_hat(s) = s**(alpha-1)/(s**alpha + width**alpha)``. Evaluating that with + complex arithmetic is independent of the real-valued modulus/argument form of Eq. (42) that the + component implements. + """ + s = damping + 1j * np.abs(np.asarray(x, dtype=float)) + return scale * np.real(s ** (alpha - 1) / (s**alpha + width**alpha)) / np.pi + + +class TestDiffusionDampedMittagLeffler: + @pytest.fixture + def mittag_leffler(self): + return DiffusionDampedMittagLeffler( + name='TestMLName', + display_name='TestML', + scale=2.0, + alpha=0.7, + width=0.3, + damping=0.05, + x_unit='meV', + ) + + ############# + # Construction + ############# + + def test_init_no_inputs(self): + # WHEN THEN + ml = DiffusionDampedMittagLeffler() + + # EXPECT + assert ml.display_name == 'DiffusionDampedMittagLeffler' + assert ml.scale.value == pytest.approx(1.0) + assert ml.alpha.value == pytest.approx(1.0) + assert ml.width.value == pytest.approx(1.0) + assert ml.damping.value == pytest.approx(1.0) + assert ml.x_unit == 'meV' + assert ml.y_unit == 'dimensionless' + + def test_initialization(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN THEN EXPECT + assert mittag_leffler.display_name == 'TestML' + assert mittag_leffler.scale.value == pytest.approx(2.0) + assert mittag_leffler.alpha.value == pytest.approx(0.7) + assert mittag_leffler.width.value == pytest.approx(0.3) + assert mittag_leffler.damping.value == pytest.approx(0.05) + assert mittag_leffler.x_unit == 'meV' + + def test_parameter_units(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN THEN EXPECT scale = x_unit * y_unit, width/damping = x_unit, alpha dimensionless + assert mittag_leffler.scale.unit == 'meV' + assert mittag_leffler.alpha.unit == 'dimensionless' + assert mittag_leffler.width.unit == 'meV' + assert mittag_leffler.damping.unit == 'meV' + + @pytest.mark.parametrize( + 'kwargs, expected_message', + [ + ({'scale': 'invalid'}, 'scale must be a number'), + ({'alpha': 'invalid'}, 'alpha must be a number'), + ({'width': 'invalid'}, 'width must be a number'), + ({'damping': 'invalid'}, 'damping 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): + DiffusionDampedMittagLeffler(**kwargs) + + @pytest.mark.parametrize( + 'kwargs, expected_message', + [ + ({'scale': -1.0}, 'scale must be non-negative'), + ({'scale': np.nan}, 'scale must be a finite number'), + ({'alpha': 0.0}, 'alpha must be greater than zero and at most one'), + ({'alpha': 1.5}, 'alpha must be greater than zero and at most one'), + ({'alpha': np.nan}, 'alpha must be a finite number'), + ({'width': -0.6}, 'must be greater than zero'), + ({'damping': -0.6}, 'must be greater than zero'), + ], + ) + def test_input_value_validation_raises(self, kwargs, expected_message): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match=expected_message): + DiffusionDampedMittagLeffler(**kwargs) + + ############# + # Property setters + ############# + + @pytest.mark.parametrize( + 'prop, valid_value', + [ + ('scale', 3.0), + ('alpha', 0.5), + ('width', 0.7), + ('damping', 0.2), + ], + ) + def test_property_setters( + self, + mittag_leffler: DiffusionDampedMittagLeffler, + prop, + valid_value, + ): + # WHEN THEN: set a valid value + setattr(mittag_leffler, prop, valid_value) + + # EXPECT + assert getattr(mittag_leffler, prop).value == valid_value + + # WHEN: set an invalid value — THEN EXPECT + with pytest.raises(TypeError, match='must be a number'): + setattr(mittag_leffler, prop, 'invalid') + + @pytest.mark.parametrize( + 'prop, invalid_value, expected_message', + [ + ('scale', -1.0, 'violates the parameter bounds'), + ('alpha', 1.5, 'violates the parameter bounds'), + ('alpha', -0.1, 'violates the parameter bounds'), + ('width', -0.5, 'width must be positive'), + ('width', 1e-12, 'violates the parameter bounds'), + ('damping', -0.5, 'damping must be positive'), + ('damping', 1e-12, 'violates the parameter bounds'), + ], + ) + def test_setters_out_of_bounds_raise( + self, + mittag_leffler: DiffusionDampedMittagLeffler, + prop, + invalid_value, + expected_message, + ): + # WHEN + original = getattr(mittag_leffler, prop).value + + # THEN EXPECT the assignment raises instead of being silently clamped + with pytest.raises(ValueError, match=expected_message): + setattr(mittag_leffler, prop, invalid_value) + assert getattr(mittag_leffler, prop).value == pytest.approx(original) + + def test_get_all_parameters(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN THEN + params = mittag_leffler.get_all_parameters() + + # EXPECT + assert len(params) == 4 + assert all(isinstance(param, Parameter) for param in params) + expected_names = { + 'TestMLName scale', + 'TestMLName alpha', + 'TestMLName width', + 'TestMLName damping', + } + assert {param.name for param in params} == expected_names + + ############# + # Evaluation + ############# + + @pytest.mark.parametrize('alpha', [0.3, 0.5, 0.8, 1.0]) + @pytest.mark.parametrize('damping', [0.01, 0.5]) + def test_evaluate_matches_laplace_reference(self, alpha, damping): + # WHEN Eq. (42) is the real-valued form of (1/pi) Re{phi_hat(damping + i|omega|)} + ml = DiffusionDampedMittagLeffler(scale=2.0, alpha=alpha, width=0.4, damping=damping) + x = np.linspace(-5.0, 5.0, 401) + + # THEN + result = ml.evaluate(x) + + # EXPECT + expected = laplace_reference(x, scale=2.0, alpha=alpha, width=0.4, damping=damping) + np.testing.assert_allclose(result, expected, rtol=1e-10, atol=1e-14) + + def test_evaluate_reproduces_printed_equation_42(self): + # WHEN width=1 reduces the model to Eq. (42) as printed (tau_R = 1), up to the 1/pi that + # the printed equation omits but Eq. (41) requires + alpha, epsilon = 0.6, 0.2 + ml = DiffusionDampedMittagLeffler(scale=1.0, alpha=alpha, width=1.0, damping=epsilon) + omega = np.linspace(-4.0, 4.0, 201) + + # THEN + result = ml.evaluate(omega) + + # EXPECT: Eq. (42) transcribed term by term + abs_omega = np.abs(omega) + square = abs_omega**2 + epsilon**2 + arg = alpha * np.angle(epsilon + 1j * abs_omega) + numerator = ( + epsilon * square ** (alpha / 2) + abs_omega * np.sin(arg) + epsilon * np.cos(arg) + ) + denominator = square * ((square**alpha + 1) * square ** (-alpha / 2) + 2 * np.cos(arg)) + np.testing.assert_allclose(result, numerator / denominator / np.pi, rtol=1e-10) + + def test_evaluate_alpha_one_is_lorentzian(self): + # WHEN alpha=1 the Mittag-Leffler function reduces to a simple exponential, so the + # spectrum is a Lorentzian of HWHM width+damping (Eq. 31 of the paper) + ml = DiffusionDampedMittagLeffler(scale=2.0, alpha=1.0, width=0.3, damping=0.1) + x = np.linspace(-4.0, 4.0, 301) + + # THEN + result = ml.evaluate(x) + + # EXPECT + hwhm = 0.3 + 0.1 + expected = 2.0 * hwhm / np.pi / (x**2 + hwhm**2) + np.testing.assert_allclose(result, expected, rtol=1e-10) + + def test_evaluate_is_symmetric(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN + x = np.linspace(0.0, 3.0, 101) + + # THEN + positive = mittag_leffler.evaluate(x) + negative = mittag_leffler.evaluate(-x) + + # EXPECT + np.testing.assert_allclose(positive, negative, rtol=1e-12) + + def test_evaluate_finite_and_maximal_at_zero( + self, mittag_leffler: DiffusionDampedMittagLeffler + ): + # WHEN a strictly positive damping keeps the otherwise singular spectrum regular at x=0 + x = np.linspace(-1.0, 1.0, 201) + + # THEN + values = mittag_leffler.evaluate(x) + + # EXPECT + assert np.all(np.isfinite(values)) + assert np.argmax(values) == 100 + assert mittag_leffler.evaluate(0.0)[0] == pytest.approx(np.max(values)) + + @pytest.mark.parametrize('alpha', [0.4, 0.7, 1.0]) + def test_scale_is_the_integrated_area(self, alpha): + # WHEN THEN + ml = DiffusionDampedMittagLeffler(scale=2.5, alpha=alpha, width=0.3, damping=0.05) + integral, _ = quad(lambda x: ml.evaluate(x)[0], -np.inf, np.inf, limit=400) + + # EXPECT + assert integral == pytest.approx(2.5, rel=1e-4) + + def test_smaller_alpha_broadens_the_wings(self): + # WHEN two profiles differ only in alpha, at fixed scale + narrow = DiffusionDampedMittagLeffler(alpha=1.0, width=0.3, damping=0.05) + broad = DiffusionDampedMittagLeffler(alpha=0.5, width=0.3, damping=0.05) + + # THEN + far_out = np.array([20.0]) + + # EXPECT the sub-exponential relaxation leaves more intensity in the wings + assert broad.evaluate(far_out)[0] > narrow.evaluate(far_out)[0] + + def test_evaluate_scipp_input_converts_units(self): + # WHEN + ml = DiffusionDampedMittagLeffler(scale=1.0, alpha=0.7, width=0.3, damping=0.05) + x_mev = np.linspace(-2.0, 2.0, 51) + x_microev = sc.array(dims=['x'], values=x_mev * 1e3, unit='microeV') + + # THEN + from_mev = ml.evaluate(x_mev) + from_microev = ml.evaluate(x_microev) + + # EXPECT the same dimensionless intensities: scale, width and damping are all resolved to + # the unit of x, so the shape only depends on the ratios between them + np.testing.assert_allclose(from_microev, from_mev, rtol=1e-10) + + def test_evaluate_scipp_output(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN + x = np.linspace(-2.0, 2.0, 50) + + # THEN + result = mittag_leffler.evaluate(x, output='scipp') + + # EXPECT + assert isinstance(result, sc.Variable) + assert result.unit == sc.Unit('dimensionless') + np.testing.assert_allclose(result.values, mittag_leffler.evaluate(x, output='numpy')) + + ############# + # Unit conversion + ############# + + def test_convert_x_unit(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN THEN + mittag_leffler.convert_x_unit('microeV') + + # EXPECT + assert mittag_leffler.x_unit == 'microeV' + assert mittag_leffler.scale.value == pytest.approx(2.0 * 1e3) + assert mittag_leffler.width.value == pytest.approx(0.3 * 1e3) + assert mittag_leffler.damping.value == pytest.approx(0.05 * 1e3) + # EXPECT the dimensionless form parameter is untouched + assert mittag_leffler.alpha.value == pytest.approx(0.7) + assert mittag_leffler.alpha.unit == 'dimensionless' + + def test_convert_x_unit_invalid_type_raises( + self, mittag_leffler: DiffusionDampedMittagLeffler + ): + # WHEN THEN EXPECT + with pytest.raises(TypeError, match=r'x_unit must be a string or sc\.Unit'): + mittag_leffler.convert_x_unit(123) + + def test_convert_x_unit_rollback_on_failure( + self, mittag_leffler: DiffusionDampedMittagLeffler + ): + # WHEN THEN + with pytest.raises(UnitError): + mittag_leffler.convert_x_unit('m') + + # EXPECT: state rolled back + assert mittag_leffler.x_unit == 'meV' + assert mittag_leffler.scale.value == pytest.approx(2.0) + assert mittag_leffler.width.value == pytest.approx(0.3) + assert mittag_leffler.damping.value == pytest.approx(0.05) + + def test_convert_y_unit(self): + # WHEN: x_unit='meV', y_unit='1/meV' → scale_unit='dimensionless' + ml = DiffusionDampedMittagLeffler(scale=1.0, width=0.3, damping=0.05, y_unit='1/meV') + + # THEN + ml.convert_y_unit('1/eV') + + # EXPECT + assert ml.y_unit == '1/eV' + assert ml.scale.value == pytest.approx(1e3) + + def test_convert_y_unit_invalid_type_raises( + self, mittag_leffler: DiffusionDampedMittagLeffler + ): + # WHEN THEN EXPECT + with pytest.raises(TypeError): + mittag_leffler.convert_y_unit(123) + + def test_convert_y_unit_rollback_on_failure( + self, mittag_leffler: DiffusionDampedMittagLeffler + ): + # WHEN THEN + with pytest.raises(UnitError): + mittag_leffler.convert_y_unit('K') + + # EXPECT: state rolled back + assert mittag_leffler.y_unit == 'dimensionless' + assert mittag_leffler.scale.value == pytest.approx(2.0) + + def test_y_unit_setter_raises(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN THEN EXPECT + with pytest.raises(AttributeError): + mittag_leffler.y_unit = '1/meV' + + ############# + # Copy and repr + ############# + + def test_copy(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN THEN + ml_copy = copy(mittag_leffler) + + # EXPECT + assert ml_copy is not mittag_leffler + assert ml_copy.display_name == mittag_leffler.display_name + for prop in ('scale', 'alpha', 'width', 'damping'): + assert getattr(ml_copy, prop).value == getattr(mittag_leffler, prop).value + assert getattr(ml_copy, prop).fixed == getattr(mittag_leffler, prop).fixed + assert ml_copy.x_unit == mittag_leffler.x_unit + + def test_repr(self, mittag_leffler: DiffusionDampedMittagLeffler): + # WHEN THEN + repr_str = repr(mittag_leffler) + + # EXPECT + assert 'DiffusionDampedMittagLeffler' in repr_str + assert 'name = TestMLName' in repr_str + assert 'x_unit = meV' in repr_str + assert 'scale =' in repr_str + assert 'alpha =' in repr_str + assert 'width =' in repr_str + assert 'damping =' in repr_str diff --git a/tests/unit/easydynamics/sample_model/diffusion_model/test_mittag_leffler_diffusion.py b/tests/unit/easydynamics/sample_model/diffusion_model/test_mittag_leffler_diffusion.py new file mode 100644 index 000000000..95b056ab6 --- /dev/null +++ b/tests/unit/easydynamics/sample_model/diffusion_model/test_mittag_leffler_diffusion.py @@ -0,0 +1,479 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +import numpy as np +import pytest +import scipp as sc +from easyscience.variable import Parameter + +from easydynamics.sample_model import DiffusionDampedMittagLeffler +from easydynamics.sample_model import Lorentzian +from easydynamics.sample_model import MittagLefflerDiffusion +from easydynamics.utils.utils import angstrom +from easydynamics.utils.utils import hbar + +ALL_Q_VARIATION = {'A_0': True, 'relaxation_rate': True, 'alpha': True} + + +class TestMittagLefflerDiffusion: + @pytest.fixture + def model(self): + return MittagLefflerDiffusion( + scale=2.0, + diffusion_coefficient=3.3e-11, + A_0=0.05, + relaxation_rate=0.02, + alpha=0.85, + ) + + @pytest.fixture + def model_with_Q(self): + return MittagLefflerDiffusion( + scale=2.0, + diffusion_coefficient=3.3e-11, + A_0=0.05, + relaxation_rate=0.02, + alpha=0.85, + allow_Q_variation=ALL_Q_VARIATION, + Q=np.linspace(0.8, 1.8, 5), + ) + + @pytest.fixture + def model_with_Q_no_variation(self): + return MittagLefflerDiffusion( + scale=2.0, + diffusion_coefficient=3.3e-11, + A_0=0.05, + relaxation_rate=0.02, + alpha=0.85, + Q=np.linspace(0.8, 1.8, 5), + ) + + ############# + # Construction + ############# + + def test_init_default(self): + # WHEN THEN + model = MittagLefflerDiffusion() + + # EXPECT + assert model.scale.value == pytest.approx(1.0) + assert model.diffusion_coefficient.value == pytest.approx(1.0) + assert model.A_0.value == pytest.approx(0.0) + assert model.A_1.value == pytest.approx(1.0) + assert model.relaxation_rate.value == pytest.approx(1.0) + assert model.alpha.value == pytest.approx(1.0) + assert model.x_unit == 'meV' + assert model.get_component_collections() == [] + + def test_init_with_Q_builds_two_components_per_Q(self, model_with_Q): + # WHEN THEN + collections = model_with_Q.get_component_collections() + + # EXPECT one elastic Lorentzian and one Mittag-Leffler component at every Q + assert len(collections) == 5 + for collection in collections: + assert collection.list_component_names() == ['Elastic Lorentzian', 'Mittag-Leffler'] + assert isinstance(collection[0], Lorentzian) + assert isinstance(collection[1], DiffusionDampedMittagLeffler) + + def test_parameter_units(self, model): + # WHEN THEN EXPECT + assert model.scale.unit == 'meV' + assert model.diffusion_coefficient.unit == 'm^2/s' + assert model.relaxation_rate.unit == 'meV' + assert model.alpha.unit == 'dimensionless' + + @pytest.mark.parametrize( + 'kwargs, expected_exception, expected_message', + [ + ({'diffusion_coefficient': 'invalid'}, TypeError, 'must be a number'), + ({'diffusion_coefficient': -1.0}, ValueError, 'must be non-negative'), + ({'A_0': 'invalid'}, TypeError, 'A_0 must be a number'), + ({'A_0': 1.5}, ValueError, 'A_0 must be between 0 and 1'), + ({'relaxation_rate': 'invalid'}, TypeError, 'must be a number'), + ({'relaxation_rate': -1.0}, ValueError, 'relaxation_rate must be at least'), + ({'alpha': 'invalid'}, TypeError, 'alpha must be a number'), + ({'alpha': 1.5}, ValueError, 'greater than zero and at most one'), + ({'alpha': 0.0}, ValueError, 'greater than zero and at most one'), + ({'mittag_leffler_name': 123}, TypeError, 'mittag_leffler_name must be a string'), + ({'allow_Q_variation': 'invalid'}, TypeError, 'must be a dict or None'), + ({'allow_Q_variation': {'nope': True}}, ValueError, 'Unknown keys'), + ], + ) + def test_input_validation_raises(self, kwargs, expected_exception, expected_message): + # WHEN THEN EXPECT + with pytest.raises(expected_exception, match=expected_message): + MittagLefflerDiffusion(**kwargs) + + @pytest.mark.parametrize( + 'prop, valid_value', + [ + ('diffusion_coefficient', 5e-11), + ('A_0', 0.3), + ('relaxation_rate', 0.05), + ('alpha', 0.6), + ], + ) + def test_setters(self, model, prop, valid_value): + # WHEN THEN + setattr(model, prop, valid_value) + + # EXPECT + assert getattr(model, prop).value == pytest.approx(valid_value) + + # WHEN: an invalid value — THEN EXPECT + with pytest.raises(TypeError, match='must be a number'): + setattr(model, prop, 'invalid') + + def test_A_1_setter_raises(self, model): + # WHEN THEN EXPECT A_1 is derived from A_0 + with pytest.raises(AttributeError, match='derived from A_0'): + model.A_1 = 0.5 + + def test_A_1_follows_A_0(self, model): + # WHEN THEN + model.A_0 = 0.25 + + # EXPECT + assert model.A_1.value == pytest.approx(0.75) + + ############# + # Q-dependent quantities + ############# + + def test_calculate_width_is_hbar_D_Q_squared(self, model_with_Q): + # WHEN THEN + Q = model_with_Q.Q.values + widths = model_with_Q.calculate_width() + + # EXPECT epsilon = hbar * D * Q**2, matching Eq. (39) of Hassani et al. + factor = hbar * model_with_Q.diffusion_coefficient / angstrom**2 + factor.convert_unit('meV') + np.testing.assert_allclose(widths, Q**2 * factor.value, rtol=1e-10) + + def test_calculate_EISF_and_QISF_sum_to_one(self, model_with_Q): + # WHEN THEN + eisf = model_with_Q.calculate_EISF() + qisf = model_with_Q.calculate_QISF() + + # EXPECT + np.testing.assert_allclose(eisf + qisf, np.ones(5), rtol=1e-12) + + def test_calculate_relaxation_rate_and_alpha_shared(self, model_with_Q_no_variation): + # WHEN THEN EXPECT the shared values are returned at every Q + np.testing.assert_allclose( + model_with_Q_no_variation.calculate_relaxation_rate(), np.full(5, 0.02) + ) + np.testing.assert_allclose(model_with_Q_no_variation.calculate_alpha(), np.full(5, 0.85)) + + def test_per_Q_parameters_are_independent(self, model_with_Q): + # WHEN a single Q gets its own alpha, relaxation rate and elastic fraction + model_with_Q._alpha_list[2].value = 0.6 + model_with_Q._relaxation_rate_list[2].value = 0.05 + model_with_Q._A_0_list[2].value = 0.2 + + # THEN + alphas = model_with_Q.calculate_alpha() + rates = model_with_Q.calculate_relaxation_rate() + eisf = model_with_Q.calculate_EISF() + + # EXPECT only that Q changed + np.testing.assert_allclose(alphas, [0.85, 0.85, 0.6, 0.85, 0.85]) + np.testing.assert_allclose(rates, [0.02, 0.02, 0.05, 0.02, 0.02]) + np.testing.assert_allclose(eisf, [0.05, 0.05, 0.2, 0.05, 0.05]) + + def test_per_Q_parameters_back_the_components(self, model_with_Q): + # WHEN + model_with_Q._alpha_list[2].value = 0.6 + model_with_Q._relaxation_rate_list[2].value = 0.05 + model_with_Q._A_0_list[2].value = 0.2 + + # THEN + collection = model_with_Q.get_component_collections()[2] + + # EXPECT the component parameters are the very ones in the per-Q lists + assert collection[1].alpha.value == pytest.approx(0.6) + assert collection[1].width.value == pytest.approx(0.05) + assert collection[0].area.value == pytest.approx(2.0 * 0.2) + assert collection[1].scale.value == pytest.approx(2.0 * 0.8) + + def test_calculate_raises_when_Q_variation_enabled_but_Q_unset(self): + # WHEN + model = MittagLefflerDiffusion(allow_Q_variation=ALL_Q_VARIATION) + + # THEN EXPECT + with pytest.raises(ValueError, match='Q must be provided'): + model.calculate_relaxation_rate() + with pytest.raises(ValueError, match='Q must be provided'): + model.calculate_alpha() + + ############# + # Parameter wiring + ############# + + def test_damping_is_shared_by_both_components(self, model_with_Q): + # WHEN THEN + collections = model_with_Q.get_component_collections() + widths = model_with_Q.calculate_width() + + # EXPECT the elastic Lorentzian's HWHM and the ML component's damping are both epsilon + for i, collection in enumerate(collections): + assert collection[0].width.value == pytest.approx(widths[i]) + assert collection[1].damping.value == pytest.approx(widths[i]) + + def test_global_diffusion_coefficient_drives_every_Q(self, model_with_Q): + # WHEN + collections = model_with_Q.get_component_collections() + before = [collection[1].damping.value for collection in collections] + + # THEN + model_with_Q.diffusion_coefficient = 2 * model_with_Q.diffusion_coefficient.value + + # EXPECT + after = [collection[1].damping.value for collection in collections] + np.testing.assert_allclose(after, np.array(before) * 2, rtol=1e-10) + + def test_shared_alpha_and_rate_drive_every_Q(self, model_with_Q_no_variation): + # WHEN + collections = model_with_Q_no_variation.get_component_collections() + + # THEN + model_with_Q_no_variation.alpha = 0.5 + model_with_Q_no_variation.relaxation_rate = 0.03 + + # EXPECT + for collection in collections: + assert collection[1].alpha.value == pytest.approx(0.5) + assert collection[1].width.value == pytest.approx(0.03) + + def test_scale_and_A_0_drive_the_component_amplitudes(self, model_with_Q_no_variation): + # WHEN + collections = model_with_Q_no_variation.get_component_collections() + + # THEN + model_with_Q_no_variation.scale = 4.0 + model_with_Q_no_variation.A_0 = 0.25 + + # EXPECT area = scale * EISF for the elastic line, scale * (1 - EISF) for the ML term + for collection in collections: + assert collection[0].area.value == pytest.approx(4.0 * 0.25) + assert collection[1].scale.value == pytest.approx(4.0 * 0.75) + + def test_matches_equation_41_by_hand(self, model_with_Q_no_variation): + # WHEN the collection is Eq. (41): EISF * Lorentzian(eps) + (1-EISF) * ML(eps) + collection = model_with_Q_no_variation.get_component_collections()[0] + epsilon = model_with_Q_no_variation.calculate_width()[0] + x = np.linspace(-0.3, 0.3, 201) + + # THEN + result = collection.evaluate(x) + + # EXPECT + elastic = Lorentzian(area=2.0 * 0.05, width=epsilon) + quasi_elastic = DiffusionDampedMittagLeffler( + scale=2.0 * 0.95, alpha=0.85, width=0.02, damping=epsilon + ) + expected = elastic.evaluate(x) + quasi_elastic.evaluate(x) + np.testing.assert_allclose(result, expected, rtol=1e-10) + + ############# + # Variables and fit targets + ############# + + def test_get_global_variables_with_Q_variation(self, model_with_Q): + # WHEN THEN + names = [variable.name for variable in model_with_Q.get_global_variables()] + + # EXPECT only scale and D are global once everything else varies with Q + assert names == ['scale', 'diffusion_coefficient'] + + def test_get_global_variables_without_Q_variation(self, model_with_Q_no_variation): + # WHEN THEN + names = [variable.name for variable in model_with_Q_no_variation.get_global_variables()] + + # EXPECT + assert names == [ + 'scale', + 'diffusion_coefficient', + 'A_0', + 'A_1', + 'relaxation_rate', + 'alpha', + ] + + def test_get_independent_variables_holds_the_per_Q_amplitudes(self, model_with_Q): + # WHEN THEN + for_all_Q = model_with_Q.get_independent_variables() + for_one_Q = model_with_Q.get_independent_variables(Q_index=1) + + # EXPECT one A_0/A_1 pair per Q + assert len(for_all_Q) == 10 + assert len(for_one_Q) == 2 + assert for_one_Q[0] is model_with_Q._A_0_list[1] + + def test_free_parameters_match_the_papers_fit(self, model_with_Q): + # WHEN the paper fits tau, alpha and the EISF per Q against a single global D + free = model_with_Q.get_free_parameters() + + # THEN + names = [parameter.name for parameter in free] + + # EXPECT scale + D + 5 * (A_0, alpha, width) + assert names.count('scale') == 1 + assert names.count('diffusion_coefficient') == 1 + assert names.count('MittagLefflerDiffusion A_0') == 5 + assert names.count('Mittag-Leffler alpha') == 5 + assert names.count('Mittag-Leffler width') == 5 + assert all(isinstance(parameter, Parameter) for parameter in free) + + def test_get_fit_targets(self, model_with_Q): + # WHEN THEN + targets = {target.name: target for target in model_with_Q.get_fit_targets()} + + # EXPECT the quasi-elastic weight points at the Mittag-Leffler component, and the + # Lorentzian keys at the elastic line + assert set(targets) == {'area', 'width', 'elastic_area'} + assert targets['area'].dataset_key == 'Mittag-Leffler scale' + assert targets['width'].dataset_key == 'Elastic Lorentzian width' + assert targets['elastic_area'].dataset_key == 'Elastic Lorentzian area' + + Q = model_with_Q.Q.values + np.testing.assert_allclose( + targets['area'].function(Q), model_with_Q.calculate_QISF(Q) * 2.0 + ) + np.testing.assert_allclose( + targets['elastic_area'].function(Q), model_with_Q.calculate_EISF(Q) * 2.0 + ) + + ############# + # Formulas from the paper + ############# + + def test_relaxation_rate_spectrum_is_normalised(self, model_with_Q_no_variation): + # WHEN Eq. (37) is a distribution over relaxation rates, so it integrates to 1 + rate = np.logspace(-8, 4, 400_001) + + # THEN + spectrum = model_with_Q_no_variation.calculate_relaxation_rate_spectrum(rate) + + # EXPECT + assert spectrum.shape == (5, rate.size) + assert np.trapezoid(spectrum[0], rate) == pytest.approx(1.0, rel=1e-3) + + def test_relaxation_rate_spectrum_alpha_one_is_peaked_at_the_rate(self): + # WHEN alpha=1 the relaxation is exponential, so p(lambda) collapses onto a delta at 1/tau + model = MittagLefflerDiffusion(alpha=1.0, relaxation_rate=0.02, Q=np.array([1.0])) + rate = np.linspace(0.001, 0.1, 2001) + + # THEN + spectrum = model.calculate_relaxation_rate_spectrum(rate) + + # EXPECT + assert rate[np.argmax(spectrum[0])] == pytest.approx(0.02, abs=1e-3) + + def test_relaxation_rate_spectrum_rejects_non_positive_rate(self, model_with_Q): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match='rate must be strictly positive'): + model_with_Q.calculate_relaxation_rate_spectrum(np.array([0.0, 1.0])) + + @pytest.mark.parametrize('alpha', [0.99, 0.7, 0.3]) + def test_energy_barrier_distribution_carries_half_the_rate_spectrum(self, alpha): + # WHEN h >= 0 maps onto lambda * tau_R <= 1 only, which by the lambda -> 1/lambda symmetry + # of Eq. (37) is half of p(lambda) + model = MittagLefflerDiffusion(alpha=alpha, Q=np.array([1.0])) + h = np.linspace(0.0, 30.0, 60_001) + + # THEN + distribution = model.calculate_energy_barrier_distribution(h)[0] + + # EXPECT + assert np.all(np.isfinite(distribution)) + assert np.trapezoid(distribution, h) == pytest.approx(0.5, rel=1e-3) + + def test_energy_barrier_distribution_narrows_as_alpha_goes_to_one(self): + # WHEN Eq. (48) tends to delta(h) as alpha -> 1 and broadens as alpha -> 0 + h = np.linspace(0.0, 10.0, 2001) + smooth = MittagLefflerDiffusion(alpha=0.99, Q=np.array([1.0])) + rough = MittagLefflerDiffusion(alpha=0.3, Q=np.array([1.0])) + + # THEN + smooth_distribution = smooth.calculate_energy_barrier_distribution(h)[0] + rough_distribution = rough.calculate_energy_barrier_distribution(h)[0] + + # EXPECT the rough landscape peaks at, and reaches, much higher barriers + assert h[np.argmax(smooth_distribution)] < h[np.argmax(rough_distribution)] + assert smooth_distribution[-1] < rough_distribution[-1] + + def test_energy_barrier_distribution_vanishes_at_zero_barrier(self, model_with_Q): + # WHEN THEN + distribution = model_with_Q.calculate_energy_barrier_distribution(np.array([0.0])) + + # EXPECT + np.testing.assert_allclose(distribution, np.zeros((5, 1)), atol=1e-15) + + ############# + # Units and Q handling + ############# + + def test_convert_x_unit(self, model_with_Q): + # WHEN + widths_before = model_with_Q.calculate_width() + + # THEN + model_with_Q.convert_x_unit('microeV') + + # EXPECT the shared template, the per-Q rates and the components all follow + assert model_with_Q.x_unit == 'microeV' + assert sc.Unit(model_with_Q.relaxation_rate.unit) == sc.Unit('microeV') + np.testing.assert_allclose(model_with_Q.calculate_width(), widths_before * 1e3, rtol=1e-10) + for collection in model_with_Q.get_component_collections(): + assert collection[1].width.value == pytest.approx(20.0) + assert collection[1].alpha.value == pytest.approx(0.85) + + def test_convert_y_unit_rescales_the_scale(self): + # WHEN scale_unit = x_unit * y_unit = meV/s + model = MittagLefflerDiffusion(scale=2.0, y_unit='1/s', Q=np.linspace(0.8, 1.8, 3)) + + # THEN + model.convert_y_unit('1/ms') + + # EXPECT the scale follows the y-unit, 1/s -> 1/ms being a factor 1e-3 + assert model.y_unit == '1/ms' + assert model.scale.value == pytest.approx(2.0e-3) + assert sc.Unit(model.scale.unit) == sc.Unit('meV/ms') + + def test_on_Q_change_rebuilds_the_per_Q_lists(self, model): + # WHEN + model._allow_Q_variation = dict(ALL_Q_VARIATION) + model.Q = np.linspace(0.5, 1.5, 3) + + # THEN + collections = model.get_component_collections() + + # EXPECT + assert len(collections) == 3 + assert len(model._alpha_list) == 3 + assert len(model._relaxation_rate_list) == 3 + assert collections[1][1].alpha is model._alpha_list[1] + + def test_clear_Q_empties_the_model(self, model_with_Q): + # WHEN THEN + model_with_Q.clear_Q(confirm=True) + + # EXPECT + assert model_with_Q.get_component_collections() == [] + assert model_with_Q._alpha_list == [] + assert model_with_Q._relaxation_rate_list == [] + assert model_with_Q._A_0_list == [] + + def test_repr(self, model): + # WHEN THEN + repr_str = repr(model) + + # EXPECT + assert 'MittagLefflerDiffusion' in repr_str + assert 'diffusion_coefficient=' in repr_str + assert 'relaxation_rate=' in repr_str + assert 'alpha=' in repr_str