diff --git a/validation/scripts/generate_smeared_data.ipynb b/validation/scripts/generate_smeared_data.ipynb new file mode 100644 index 0000000..ce9ef96 --- /dev/null +++ b/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 +} diff --git a/validation/test/unpolarised/data/Untitled.ipynb b/validation/test/unpolarised/data/Untitled.ipynb deleted file mode 100644 index ffb9ceb..0000000 --- a/validation/test/unpolarised/data/Untitled.ipynb +++ /dev/null @@ -1,107 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "layers = np.array([[0, 2.07, 0, 0],\n", - " [100, 3.45, 0.1, 3],\n", - " [200, 5.0, 0.01, 1],\n", - "[0, 6., 0, 5]])" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[0.00e+00, 2.07e+00, 0.00e+00, 0.00e+00],\n", - " [1.00e+02, 3.45e+00, 1.00e-01, 3.00e+00],\n", - " [2.00e+02, 5.00e+00, 1.00e-02, 1.00e+00],\n", - " [0.00e+00, 6.00e+00, 0.00e+00, 5.00e+00]])" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "layers" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "np.savetxt('test0.layers', layers)" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "x = np.linspace(0.005, 0.5, 1001)\n", - "\n", - "data = np.load('refl1d.npy')" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "np.savetxt('test0.txt', np.array([x, data]).T)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "n" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "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.7.3" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -}