MAINT: provenance for resolution smearing validation

validation/scripts/generate_smeared_data.ipynb

@@ -0,0 +1,201 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "47acb03d-06fe-41eb-a3e6-e48672462c9c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from functools import lru_cache\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import scipy\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import test_discovery\n",
+    "from refnx.reflect import abeles"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "39f3ca00-7aa7-4a40-8098-25373ecf135c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# test systems to calculate smeared reflectivity for\n",
+    "unpolarised_systems = [\"test4\", \"test5\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c3a8d5a4-c925-4caa-9e1d-1ecfcf11e6ec",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# the function for calculating unsmeared reflectivity\n",
+    "# It shouldn't matter which one we use, as long as all the packages agree\n",
+    "fkernel = abeles\n",
+    "\n",
+    "# Integrate smearing over (-_INTLIMIT, _INTLIMIT)\n",
+    "_INTLIMIT = 3.5\n",
+    "\n",
+    "# Gaussian Quadrature order\n",
+    "QUAD_ORDER = 10001"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0e9a9f73-c73f-4b5a-b7ba-50bd4bb10951",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@lru_cache(maxsize=128)\n",
+    "def gauss_legendre(n):\n",
+    "    \"\"\"\n",
+    "    Calculate gaussian quadrature abscissae and weights\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    n : int\n",
+    "        Gaussian quadrature order.\n",
+    "    Returns\n",
+    "    -------\n",
+    "    (x, w) : tuple\n",
+    "        The abscissae and weights for Gauss Legendre integration.\n",
+    "    \"\"\"\n",
+    "    return scipy.special.p_roots(n)\n",
+    "\n",
+    "\n",
+    "def smeared_kernel_pointwise(qvals, w, dqvals, quad_order=17, threads=-1):\n",
+    "    \"\"\"\n",
+    "    Resolution smearing that uses fixed order Gaussian quadrature integration\n",
+    "    for the convolution.\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    qvals : array-like\n",
+    "        The Q values for evaluation\n",
+    "    w : array-like\n",
+    "        The uniform slab model parameters in 'layer' form.\n",
+    "    dqvals : array-like\n",
+    "        dQ values corresponding to each value in `qvals`. Each dqval is the\n",
+    "        FWHM of a Gaussian approximation to the resolution kernel.\n",
+    "    quad-order : int, optional\n",
+    "        Specify the order of the Gaussian quadrature integration for the\n",
+    "        convolution.\n",
+    "    threads: int, optional\n",
+    "        Specifies the number of threads for parallel calculation. This\n",
+    "        option is only applicable if you are using the ``_creflect``\n",
+    "        module. The option is ignored if using the pure python calculator,\n",
+    "        ``_reflect``. If `threads == -1` then all available processors are\n",
+    "        used.\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    reflectivity : np.ndarray\n",
+    "        The smeared reflectivity\n",
+    "    \"\"\"\n",
+    "    # get the gauss-legendre weights and abscissae\n",
+    "    abscissa, weights = gauss_legendre(quad_order)\n",
+    "\n",
+    "    # get the normal distribution at that point\n",
+    "    prefactor = 1.0 / np.sqrt(2 * np.pi)\n",
+    "\n",
+    "    def gauss(x):\n",
+    "        return np.exp(-0.5 * x * x)\n",
+    "\n",
+    "    gaussvals = prefactor * gauss(abscissa * _INTLIMIT)\n",
+    "\n",
+    "    va = qvals - _INTLIMIT * dqvals\n",
+    "    vb = qvals + _INTLIMIT * dqvals\n",
+    "\n",
+    "    assert np.all(va > 0)\n",
+    "\n",
+    "    va = va[..., np.newaxis]\n",
+    "    vb = vb[..., np.newaxis]\n",
+    "\n",
+    "    qvals_for_res = (abscissa[np.newaxis, :] * (vb - va) + vb + va) / 2.0\n",
+    "    smeared_rvals = fkernel(qvals_for_res, w, threads=threads)\n",
+    "\n",
+    "    smeared_rvals *= np.atleast_2d(gaussvals * weights)\n",
+    "    return np.sum(smeared_rvals, -1) * _INTLIMIT"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "aba06376-c1e4-4957-ae3d-4b823be0836d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# grab the layers and data for the systems that need resolution smearing (re-)applied\n",
+    "need_to_smear = {\n",
+    "    td[0].split(\".txt\")[0]: td\n",
+    "    for td in test_discovery.get_test_data()\n",
+    "    if td[0].split(\".txt\")[0] in unpolarised_systems\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d8cf1305-f980-4933-903d-9b5a0a09be8a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for name, td in need_to_smear.items():\n",
+    "    layers, data = td[1:3]\n",
+    "    qvals = data[:, 0]\n",
+    "    dqvals = data[:, 3]\n",
+    "    r_existing = data[:, 1]\n",
+    "\n",
+    "    # calculate smeared data\n",
+    "    r_smeared = smeared_kernel_pointwise(qvals, layers, dqvals, quad_order=QUAD_ORDER)\n",
+    "\n",
+    "    try:\n",
+    "        np.testing.assert_allclose(r_existing, r_smeared, rtol=5e-3)\n",
+    "    except AssertionError as e:\n",
+    "        print(name)\n",
+    "        raise e\n",
+    "\n",
+    "    data[:, 1] = r_smeared\n",
+    "\n",
+    "    pth = Path(\"..\", \"test\", \"unpolarised\", \"data\", f\"{name}.dat\")\n",
+    "    np.savetxt(pth, data)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c22b6e2c-ecbb-4eb8-9fc5-2402680a11ec",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}