Merge pull request #136 from scipp/source-optimize

Optimize source using chopper cascade acceptance diagram

Author: nvaytet

docs/index.md

@@ -63,8 +63,9 @@ sources
 components
 multiple-pulses
 wfm
-dashboard
+optimizations
 ess-instruments
+dashboard
 api-reference/index
 developer/index
 about/index

docs/optimizations.ipynb

@@ -0,0 +1,230 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# Optimizations\n",
+    "\n",
+    "## Optimize for a given chopper cascade\n",
+    "\n",
+    "Most times when we run a `tof` model,\n",
+    "the vast majority of neutrons in a pulse get blocked by the first few choppers in the beam path.\n",
+    "\n",
+    "For example, using the chopper settings for the Odin (ESS) instrument:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import scipp as sc\n",
+    "import matplotlib.pyplot as plt\n",
+    "import tof"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "s1 = tof.Source(facility=\"ess\", neutrons=1_000_000)\n",
+    "beamline = tof.facilities.ess.odin(pulse_skipping=True)\n",
+    "m1 = tof.Model(source=s1, **beamline)\n",
+    "r1 = m1.run()\n",
+    "r1"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3",
+   "metadata": {},
+   "source": [
+    "We can see that out of 1M neutrons that left the source, over 900K were blocked by the first chopper.\n",
+    "In the end, only ~57K make it to the detector.\n",
+    "\n",
+    "This is incredibly wasteful in both memory and compute.\n",
+    "\n",
+    "We can however use the information of the opening and closing times of the choppers to predict which parts of the pulse (in the `birth_time`/`wavelength` space) will make it through and which regions will be blocked.\n",
+    "This is otherwise known as a 'chopper acceptance diagram'.\n",
+    "\n",
+    "This can be visualized by looking at the birth times and wavelengths of the neutrons that made it to the detector,\n",
+    "and compare that to the original distribution of neutrons in the source."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig1 = s1.data.hist(wavelength=300, birth_time=300).plot(norm='log', title=\"Sampled from source\")\n",
+    "fig2 = r1['detector'].data.hist(wavelength=300, birth_time=300).plot(norm='log', title=\"Neutrons that make it to the detector\")\n",
+    "\n",
+    "fig1 + fig2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5",
+   "metadata": {},
+   "source": [
+    "The source has an in-built method to apply the chopper acceptance criteria during the sampling of neutrons;\n",
+    "it is activated via the `optimize_for` argument.\n",
+    "It expects a list of choppers, and only neutrons that would make it through all choppers end up in the source."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Filter out choppers from list of Odin components\n",
+    "choppers = {comp.name: comp for comp in beamline['components'] if isinstance(comp, tof.Chopper)}\n",
+    "\n",
+    "# Create optimized source\n",
+    "s2 = tof.Source(facility=\"ess\", neutrons=1_000_000, optimize_for=choppers)\n",
+    "\n",
+    "m2 = tof.Model(source=s2, **beamline)\n",
+    "r2 = m2.run()\n",
+    "r2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7",
+   "metadata": {},
+   "source": [
+    "We can now see that **all 1M neutrons make to the detector**,\n",
+    "and plotting the birth time/wavelength distribution illustrates the optimization:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig1 = s2.data.hist(wavelength=300, birth_time=300).plot(norm='log', title=\"Sampled from source\")\n",
+    "fig2 = r2['detector'].data.hist(wavelength=300, birth_time=300).plot(norm='log', title=\"Neutrons that make it to the detector\")\n",
+    "\n",
+    "fig1 + fig2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9",
+   "metadata": {},
+   "source": [
+    "It is also clear that the signal recorded at detector is much less noisy due to the improved statistics.\n",
+    "It is also important to note that the overall shape of the data (relative intensities) was not changed by the optimization."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "\n",
+    "r1['detector'].toa.plot(ax=ax[0])\n",
+    "r2['detector'].toa.plot(color='C1', ax=ax[0].twinx())\n",
+    "\n",
+    "r1['detector'].wavelength.plot(ax=ax[1])\n",
+    "r2['detector'].wavelength.plot(color='C1', ax=ax[1].twinx())\n",
+    "\n",
+    "fig.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11",
+   "metadata": {},
+   "source": [
+    "## Pulse skipping\n",
+    "\n",
+    "Some instruments (such as Odin) use a pulse skipping chopper to 'skip' every other pulse, thus allowing to record a wider wavelength range at the detector without having issues where neutrons from successive pulses mix (also known as pulse-overlap).\n",
+    "\n",
+    "In such a setup, when running multiple pulses, all neutrons from every other pulse are rendered useless.\n",
+    "Yet, `tof` naively treats them as normal neutrons and tries to follow them all the way to the detector."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "s1 = tof.Source(facility=\"ess\", neutrons=1_000_000, pulses=4)\n",
+    "m1 = tof.Model(source=s1, **beamline)\n",
+    "r1 = m1.run()\n",
+    "r1.plot(blocked_rays=5000)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "13",
+   "metadata": {},
+   "source": [
+    "To avoid wasting all the neutrons in the second pulse, a simple trick is to override the frequency of the source.\n",
+    "Here we set it to 7 Hz (half of the original 14 Hz), meaning that the second pulse above will not exist at all."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "s2 = tof.Source(facility=\"ess\", neutrons=1_000_000, pulses=2, frequency=sc.scalar(7, unit=\"Hz\"))\n",
+    "m2 = tof.Model(source=s2, **beamline)\n",
+    "r2 = m2.run()\n",
+    "r2.plot(blocked_rays=5000)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15",
+   "metadata": {},
+   "source": [
+    "We have now used half as many neutrons to achieve the same results.\n",
+    "In combination with the `optimize_for` option introduced above, these optimizations can lead to significant speedups."
+   ]
+  }
+ ],
+ "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.12.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}