{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "# Rearranging and Filtering Binned Data\n", "\n", "## Introduction\n", "\n", "Event filtering refers to the process of removing or extracting a subset of events based on some criterion such as the temperature of the measured sample at the time an event was detected.\n", "Instead of extracting based on a single parameter value or interval, we may also want to rearrange data based on the parameter value, providing quick and convenient access to the parameter-dependence of our data.\n", "Scipp's binned data can be used for both of these purposes.\n", "\n", "The [Quick Reference](#Quick-Reference) below provides a brief overview of the options.\n", "A more detailed walkthrough based on actual data can be found in the [Full example](#Full-example)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Quick Reference\n", "\n", "### Extract events matching parameter value\n", "\n", "Use [label-based indexing on the bins property](../../generated/classes/scipp.Bins.rst#scipp.Bins).\n", "This works similar to regular [label-based indexing](../slicing.ipynb#Label-based-indexing) but operates on the unordered bin contents.\n", "Example:\n", "\n", "```python\n", "param_value = sc.scalar(1.2, unit='m')\n", "filtered = da.bins['param', param_value]\n", "```\n", "\n", "- The output data array has the same dimensions as the input `da`.\n", "- `filtered` contains a *copy* of the filtered events.\n", "\n", "### Extract events falling into a parameter interval\n", "\n", "Use [label-based indexing on the bins property](../../generated/classes/scipp.Bins.rst#scipp.Bins).\n", "This works similar to regular [label-based indexing](../slicing.ipynb#Label-based-indexing) but operates on the unordered bin contents.\n", "Example:\n", "\n", "```python\n", "start = sc.scalar(1.2, unit='m')\n", "stop = sc.scalar(1.3, unit='m')\n", "filtered = da.bins['param', start:stop]\n", "```\n", "\n", "- The output data array has the same dimensions as the input `da`.\n", "`filtered` contains a *copy* of the filtered events.\n", "- Note that as usual the upper bound of the interval (here $1.3~\\text{m}$) is *not* included.\n", "\n", "### Split into bins based on a discrete event parameter\n", "\n", "Use [scipp.group](../../generated/functions/scipp.group.rst).\n", "Example:\n", "\n", "```python\n", "split = da.group('param')\n", "```\n", "\n", "- The output data array has a new dimension `'param'` in addition to the dimensions of the input.\n", "- `split` contains a *copy* of the reordered events.\n", "- Pass an explicit variable to `group` listing desired groups to limit what is included in the output.\n", "\n", "### Split into bins based on a continuous event parameter\n", "\n", "Use [scipp.bin](../../generated/functions/scipp.bin.rst).\n", "Example:\n", "\n", "```python\n", "split = da.bin(param=10)\n", "```\n", "\n", "- The output data array has a new dimension `'param'` in addition to the dimensions of the input.\n", "- `split` contains a *copy* of the reordered events.\n", "- Provide an explicit variable to `bin` to limit the parameter interval that is included in the output, or for fine-grained control over the sub-intervals.\n", "\n", "### Compute derived event parameters for subsequent extracting or splitting\n", "\n", "Use [scipp.transform_coords](../../generated/functions/scipp.transform_coords.rst).\n", "Example:\n", "\n", "```python\n", "da2 = da.transform_coords(derived_param=lambda p1, p2: p1 + p2)\n", "```\n", "\n", "`da2` can now be used with any of the methods for extracting or splitting data described above.\n", "The intermediate variable can also be omitted, and we can directly extract or split the result:\n", "\n", "```python\n", "filtered = da.transform_coords(derived_param=lambda p1, p2: p1 + p2) \\\n", " .bin(derived_param=10)\n", "```\n", "\n", "### Compute derived event parameters from time-series or other metadata\n", "\n", "In practice, events are often tagged with a timestamp, which can be used to lookup parameter values from, e.g., a time-series log given by a data array with a single dimension and a coordinate matching the coordinate name of the timestamps.\n", "Use [scipp.lookup](../../generated/functions/scipp.lookup.rst) with [scipp.transform_coords](../../generated/functions/scipp.transform_coords.rst). Example:\n", "\n", "```python\n", "temperature = da.attrs['sample_temperature'].value # temperature value time-series\n", "interp_temperature = sc.lookup(temperature, mode='previous')\n", "filtered = da.transform_coords(temperature=interp_temperature) \\\n", " .bin(temperature=10)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Full example\n", "\n", "### Input data\n", "\n", "In the following we use neutron diffraction data for a stainless steel tensile bar in a loadframe measured at the [VULCAN Engineering Materials Diffractometer](https://neutrons.ornl.gov/vulcan), kindly provided by the SNS.\n", "Scipp's sample data includes an excerpt from the full dataset:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "import scipp as sc\n", "\n", "dg = sc.data.vulcan_steel_strain_data()\n", "dg" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "da = dg['data']\n", "da" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `dspacing` dimension is the interplanar lattice spacing (the spacing between planes in a crystal), and plotting this data we see two diffraction peaks:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "da.hist().plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Extract time interval\n", "\n", "The [mechanical strain](https://en.wikipedia.org/wiki/Strain_(mechanics)) of the steel sample in the loadframe is recorded in the metadata:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "strain = dg['loadframe.strain']\n", "strain.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the strain drops off for some reason at the end.\n", "We can filter out those events, by extracting the rest as outlined in [Extract events matching parameter value](#Extract-events-matching-parameter-value):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "start = strain.coords['time'][0]\n", "stop = strain.coords['time'][np.argmax(strain.values)]\n", "da = da.bins['time', start:stop]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**Note**\n", " \n", "The above is just a concise way of binning into a single time interval and squeezing the time dimension from the result.\n", "\n", "If *multiple* intervals are to be extracted then the mechanism based on `start` and `stop` values becomes highly inefficient, as every time `da.bins['param', start:stop]` is called *all* of the events have to be processed.\n", "Instead prefer using `da.bin(param=param_bin_edges)` and slice the result using regular positional (or label-based) indexing.\n", "Similarly, prefer using `da.group('param')` to extract based on multiple discrete values.\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "### Filter bad pulses\n", "\n", "In the previous example we directly used an existing event-coordinate (`da.bins.coords['time']`) for selecting the desired subset of data.\n", "In many practical cases such a coordinate may not be available yet and needs to be computed as a preparatory step.\n", "Scipp facilitates this using [scipp.transform_coords](../../generated/functions/scipp.transform_coords.rst) and [scipp.lookup](../../generated/functions/scipp.lookup.rst).\n", "When the desired event-coordinate can be computed directly from existing coordinates then `transform_coords` can do the job on its own.\n", "In other cases, such as the following example, we combine it with `lookup` to, e.g., map timestamps to corresponding sensor readings.\n", "\n", "Our data stores the so called *proton charge*, the total charge of protons per pulse (which produced the neutrons scattered off the sample):" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "proton_charge = dg['proton_charge']\n", "proton_charge.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some pulses have a very low proton charge which may indicate a problem with the source, so we may want to remove events that were produced from these pulses.\n", "We can use `lookup` to define the following \"interpolation function\", marking any pulse as \"good\" if it has more than 90% of the mean proton charge:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "good_pulse = sc.lookup(proton_charge > 0.9 * proton_charge.mean(), mode='previous')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`transform_coords` can utilize this interpolation function to compute a new coordinate (`good_pulse`, with `True` and `False` values) from the `da.bins.coords['time']` coordinate.\n", "We used `mode='previous'` above, so an event's `good_pulse` value will be defined by the *previous* pulse, i.e., the one that produced the neutron event.\n", "See the documentation of [scipp.lookup](../../generated/functions/scipp.lookup.rst) for a full list of available options.\n", "\n", "The return value of `transform_coords` can then be used to index the `bins` property, here to extract only the events that have `good_pulse=True`, i.e., were created by a proton pulse that fulfilled the above criterion:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "da = da.transform_coords(good_pulse=good_pulse).bins['good_pulse', sc.index(True)]\n", "da" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rearrange data based on strain\n", "\n", "As outlined in [Split into bins based on a continuous event parameter](#Split-into-bins-based-on-a-continuous-event-parameter) we can rearrange data based on the current strain value at the time a neutron event was recorded:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interpolate_strain = sc.lookup(strain, mode='previous')\n", "filtered = da.transform_coords(strain=interpolate_strain).bin(strain=100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can histogram and plot the result, but the figure is not very illuminating.\n", "This illustrates a common problem, and we will show below how to address it:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "filtered.hist().transpose().plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The problem we run into with the above figure is that we have a lot more data (events) at zero strain than for the other strain values.\n", "We should therefore *normalize* the result to the incident flux.\n", "In this simplified example we can use the proton charge to accomplish this.\n", "We compute the integrated proton charge for a given strain value:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "proton_charge = dg['proton_charge']\n", "charge_per_time_interval = proton_charge.bin(time=strain.coords['time'])\n", "charge_per_time_interval.coords['strain'] = strain.data\n", "charge_per_strain_value = charge_per_time_interval.bin(\n", " strain=filtered.coords['strain']\n", ").hist()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then normalize our data and obtain a more meaningful plot.\n", "It is now clearly visible how one of the diffraction peaks splits into two under increasing mechanical strain on the sample:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "normalized = filtered / charge_per_strain_value\n", "normalized.hist(strain=30, dspacing=300).plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we have reordered our data by strain, data for a specific strain value can now be *cheaply* extracted using positional or label-based indexing.\n", "This returns a view of the events, i.e., not a copy:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "normalized['strain', 80].hist(dspacing=400).plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In practice it can be useful to integrate strain ranges for comparison in a 1-D plot.\n", "Here we simply histogram with a coarser strain binning:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "coarse = normalized.hist(strain=6, dspacing=200)\n", "strains = [coarse['strain', sc.scalar(x)] for x in [0.0, 0.3, 0.6]]\n", "lines = {f\"strain={strain.coords['strain'].values}\": strain for strain in strains}\n", "sc.plot(lines, norm='log')" ] } ], "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.9.15" } }, "nbformat": 4, "nbformat_minor": 4 }