Merge pull request #92 from scipp/orso-multi-dataset

docs: example writing multiple datasets to orso file

Author: jokasimr

docs/user-guide/amor/amor-reduction.ipynb

@@ -143,8 +143,32 @@
    "source": [
     "## Computing sample reflectivity\n",
     "\n",
-    "We now compute the sample reflectivity from 3 runs that used different sample rotation angles.\n",
-    "The different rotation angles cover different ranges in $Q$."
+    "We now compute the sample reflectivity from 4 runs that used different sample rotation angles.\n",
+    "The measurements at different rotation angles cover different ranges of $Q$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this tutorial we use some Amor data files we have received.\n",
+    "The file paths to the tutorial files are obtained by calling:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "amor.data.amor_sample_run('608')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "When you encounter `amor.data.amor_sample_run` you should imagining replacing that with a path to your own dataset."
    ]
   },
   {
@@ -154,17 +178,31 @@
    "outputs": [],
    "source": [
     "runs = {\n",
-    "    '608': sc.scalar(0.85, unit='deg'),\n",
-    "    '609': sc.scalar(2.25, unit='deg'),\n",
-    "    '610': sc.scalar(3.65, unit='deg'),\n",
-    "    '611': sc.scalar(5.05, unit='deg'),\n",
+    "    '608': {\n",
+    "        # The sample rotation values in the files are slightly off, so we replace\n",
+    "        # them with corrected values.\n",
+    "        SampleRotation[SampleRun]: sc.scalar(0.85, unit='deg'),\n",
+    "        Filename[SampleRun]: amor.data.amor_sample_run('608'),\n",
+    "    },\n",
+    "    '609': {\n",
+    "        SampleRotation[SampleRun]: sc.scalar(2.25, unit='deg'),\n",
+    "        Filename[SampleRun]: amor.data.amor_sample_run('609'),\n",
+    "    },\n",
+    "    '610': {\n",
+    "        SampleRotation[SampleRun]: sc.scalar(3.65, unit='deg'),\n",
+    "        Filename[SampleRun]: amor.data.amor_sample_run('610'),\n",
+    "    },\n",
+    "    '611': {\n",
+    "        SampleRotation[SampleRun]: sc.scalar(5.05, unit='deg'),\n",
+    "        Filename[SampleRun]: amor.data.amor_sample_run('611'),\n",
+    "    },\n",
     "}\n",
     "\n",
     "reflectivity = {}\n",
-    "for file, angle in runs.items():\n",
-    "    workflow[Filename[SampleRun]] = amor.data.amor_sample_run(file)\n",
-    "    workflow[SampleRotation[SampleRun]] = angle\n",
-    "    reflectivity[file] = workflow.compute(ReflectivityOverQ).hist()\n",
+    "for run_number, params in runs.items():\n",
+    "    workflow[Filename[SampleRun]] = params[Filename[SampleRun]]\n",
+    "    workflow[SampleRotation[SampleRun]] = params[SampleRotation[SampleRun]]\n",
+    "    reflectivity[run_number] = workflow.compute(ReflectivityOverQ).hist()\n",
     "\n",
     "sc.plot(reflectivity, norm='log', vmin=1e-4)"
    ]
@@ -186,7 +224,12 @@
    "source": [
     "from ess.reflectometry.tools import scale_reflectivity_curves_to_overlap\n",
     "\n",
-    "scaled_reflectivity_curves, scale_factors = scale_reflectivity_curves_to_overlap(reflectivity.values())\n",
+    "scaled_reflectivity_curves, scale_factors = scale_reflectivity_curves_to_overlap(\n",
+    "    reflectivity.values(),\n",
+    "    # Optionally specify a Q-interval where the reflectivity is known to be 1.0\n",
+    "    critical_edge_interval=(sc.scalar(0.01, unit='1/angstrom'), sc.scalar(0.014, unit='1/angstrom'))\n",
+    ")\n",
+    "\n",
     "sc.plot(dict(zip(reflectivity.keys(), scaled_reflectivity_curves, strict=True)), norm='log', vmin=1e-5)"
    ]
   },
@@ -236,14 +279,14 @@
    "source": [
     "# Start by computing the `ReflectivityData` for each of the files\n",
     "diagnostics = {}\n",
-    "for file, angle in runs.items():\n",
-    "    workflow[Filename[SampleRun]] = amor.data.amor_sample_run(file)\n",
-    "    workflow[SampleRotation[SampleRun]] = angle\n",
-    "    diagnostics[file] = workflow.compute((ReflectivityData, ThetaBins[SampleRun]))\n",
+    "for run_number, params in runs.items():\n",
+    "    workflow[Filename[SampleRun]] = params[Filename[SampleRun]]\n",
+    "    workflow[SampleRotation[SampleRun]] = params[SampleRotation[SampleRun]]\n",
+    "    diagnostics[run_number] = workflow.compute((ReflectivityData, ThetaBins[SampleRun]))\n",
     "\n",
     "# Scale the results using the scale factors computed earlier\n",
-    "for key, scale_factor in zip(reflectivity.keys(), scale_factors, strict=True):\n",
-    "    diagnostics[key][ReflectivityData] *= scale_factor"
+    "for run_number, scale_factor in zip(reflectivity.keys(), scale_factors, strict=True):\n",
+    "    diagnostics[run_number][ReflectivityData] *= scale_factor"
    ]
   },
   {
@@ -263,8 +306,8 @@
     "from ess.amor.figures import wavelength_theta_figure\n",
     "\n",
     "wavelength_theta_figure(\n",
-    "    [res[ReflectivityData] for res in diagnostics.values()],\n",
-    "    theta_bins=[res[ThetaBins[SampleRun]] for res in diagnostics.values()],\n",
+    "    [result[ReflectivityData] for result in diagnostics.values()],\n",
+    "    theta_bins=[result[ThetaBins[SampleRun]] for result in diagnostics.values()],\n",
     "    q_edges_to_display=(sc.scalar(0.018, unit='1/angstrom'), sc.scalar(0.113, unit='1/angstrom'))\n",
     ")"
    ]
@@ -318,8 +361,8 @@
    "source": [
     "from ess.amor.figures import wavelength_z_figure\n",
     "\n",
-    "workflow[Filename[SampleRun]] = amor.data.amor_sample_run('608')\n",
-    "workflow[SampleRotation[SampleRun]] = runs['608']\n",
+    "workflow[Filename[SampleRun]] = runs['608'][Filename[SampleRun]]\n",
+    "workflow[SampleRotation[SampleRun]] = runs['608'][SampleRotation[SampleRun]]\n",
     "wavelength_z_figure(\n",
     "    workflow.compute(FootprintCorrectedData[SampleRun]),\n",
     "    wavelength_bins=workflow.compute(WavelengthBins),\n",
@@ -430,6 +473,30 @@
     ")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now let's repeat this for all the sample measurements!\n",
+    "To do that we can use an utility in `ess.reflectometry.tools`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ess.reflectometry.tools import orso_datasets_from_measurements\n",
+    "\n",
+    "datasets = orso_datasets_from_measurements(\n",
+    "    workflow,\n",
+    "    runs.values(),\n",
+    "    # Optionally scale the curves to overlap using `scale_reflectivity_curves_to_overlap`\n",
+    "    scale_to_overlap=True\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -444,7 +511,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "iofq_dataset.save('amor_reduced_iofq.ort')"
+    "fileio.orso.save_orso(datasets=datasets, fname='amor_reduced_iofq.ort')"
    ]
   },
   {

src/ess/reflectometry/tools.py

@@ -1,11 +1,17 @@
 # SPDX-License-Identifier: BSD-3-Clause
 # Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
-from collections.abc import Sequence
+from collections.abc import Mapping, Sequence
 from itertools import chain
+from typing import Any
 
 import numpy as np
+import sciline
 import scipp as sc
 import scipy.optimize as opt
+from orsopy.fileio.orso import OrsoDataset
+
+from ess.reflectometry import orso
+from ess.reflectometry.types import ReflectivityOverQ
 
 _STD_TO_FWHM = sc.scalar(2.0) * sc.sqrt(sc.scalar(2.0) * sc.log(sc.scalar(2.0)))
 
@@ -286,3 +292,63 @@ def combine_curves(
         ),
         coords={'Q': q_bin_edges},
     )
+
+
+def orso_datasets_from_measurements(
+    workflow: sciline.Pipeline,
+    runs: Sequence[Mapping[type, Any]],
+    *,
+    scale_to_overlap: bool = True,
+) -> list[OrsoDataset]:
+    '''Produces a list of ORSO datasets containing one
+    reflectivity curve for each of the provided runs.
+    Each entry of :code:`runs` is a mapping of parameters and
+    values needed to produce the dataset.
+
+    Optionally, the reflectivity curves can be scaled to overlap in
+    the regions where they have the same Q-value.
+
+    Parameters
+    -----------
+    workflow:
+        The sciline workflow used to compute `ReflectivityOverQ` for each of the runs.
+
+    runs:
+        The sciline parameters to be used for each run
+
+    scale_to_overlap:
+        If True the curves will be scaled to overlap.
+        Note that the curve of the first run is unscaled and
+        the rest are scaled to match it.
+
+    Returns
+    ---------
+    list of the computed ORSO datasets, containing one reflectivity curve each
+    '''
+    reflectivity_curves = []
+    for parameters in runs:
+        wf = workflow.copy()
+        for name, value in parameters.items():
+            wf[name] = value
+        reflectivity_curves.append(wf.compute(ReflectivityOverQ))
+
+    scale_factors = (
+        scale_reflectivity_curves_to_overlap([r.hist() for r in reflectivity_curves])[1]
+        if scale_to_overlap
+        else (1,) * len(runs)
+    )
+
+    datasets = []
+    for parameters, curve, scale_factor in zip(
+        runs, reflectivity_curves, scale_factors, strict=True
+    ):
+        wf = workflow.copy()
+        for name, value in parameters.items():
+            wf[name] = value
+        wf[ReflectivityOverQ] = scale_factor * curve
+        dataset = wf.compute(orso.OrsoIofQDataset)
+        dataset.info.reduction.corrections = orso.find_corrections(
+            wf.get(orso.OrsoIofQDataset)
+        )
+        datasets.append(dataset)
+    return datasets

tests/tools_test.py

@@ -1,12 +1,21 @@
 # SPDX-License-Identifier: BSD-3-Clause
 # Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
+import numpy as np
 import pytest
+import sciline as sl
 import scipp as sc
 from numpy.testing import assert_allclose as np_assert_allclose
+from orsopy.fileio import Orso, OrsoDataset
 from scipp.testing import assert_allclose
 
-from ess.reflectometry import tools
-from ess.reflectometry.tools import combine_curves, scale_reflectivity_curves_to_overlap
+from ess.reflectometry.orso import OrsoIofQDataset
+from ess.reflectometry.tools import (
+    combine_curves,
+    linlogspace,
+    orso_datasets_from_measurements,
+    scale_reflectivity_curves_to_overlap,
+)
+from ess.reflectometry.types import Filename, ReflectivityOverQ, SampleRun
 
 
 def curve(d, qmin, qmax):
@@ -138,31 +147,31 @@ def test_combined_curves():
 
 
 def test_linlogspace_linear():
-    q_lin = tools.linlogspace(
+    q_lin = linlogspace(
         dim='qz', edges=[0.008, 0.08], scale='linear', num=50, unit='1/angstrom'
     )
     expected = sc.linspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom')
     assert sc.allclose(q_lin, expected)
 
 
 def test_linlogspace_linear_list_input():
-    q_lin = tools.linlogspace(
+    q_lin = linlogspace(
         dim='qz', edges=[0.008, 0.08], unit='1/angstrom', scale=['linear'], num=[50]
     )
     expected = sc.linspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom')
     assert sc.allclose(q_lin, expected)
 
 
 def test_linlogspace_log():
-    q_log = tools.linlogspace(
+    q_log = linlogspace(
         dim='qz', edges=[0.008, 0.08], unit='1/angstrom', scale='log', num=50
     )
     expected = sc.geomspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom')
     assert sc.allclose(q_log, expected)
 
 
 def test_linlogspace_linear_log():
-    q_linlog = tools.linlogspace(
+    q_linlog = linlogspace(
         dim='qz',
         edges=[0.008, 0.03, 0.08],
         unit='1/angstrom',
@@ -176,7 +185,7 @@ def test_linlogspace_linear_log():
 
 
 def test_linlogspace_log_linear():
-    q_loglin = tools.linlogspace(
+    q_loglin = linlogspace(
         dim='qz',
         edges=[0.008, 0.03, 0.08],
         unit='1/angstrom',
@@ -190,7 +199,7 @@ def test_linlogspace_log_linear():
 
 
 def test_linlogspace_linear_log_linear():
-    q_linloglin = tools.linlogspace(
+    q_linloglin = linlogspace(
         dim='qz',
         edges=[0.008, 0.03, 0.08, 0.12],
         unit='1/angstrom',
@@ -206,10 +215,35 @@ def test_linlogspace_linear_log_linear():
 
 def test_linlogspace_bad_input():
     with pytest.raises(ValueError, match="Sizes do not match"):
-        _ = tools.linlogspace(
+        _ = linlogspace(
             dim='qz',
             edges=[0.008, 0.03, 0.08, 0.12],
             unit='1/angstrom',
             scale=['linear', 'log'],
             num=[16, 20],
         )
+
+
+@pytest.mark.filterwarnings("ignore:No suitable")
+def test_orso_datasets_tool():
+    def normalized_ioq(filename: Filename[SampleRun]) -> ReflectivityOverQ:
+        return filename
+
+    def orso_dataset(filename: Filename[SampleRun]) -> OrsoIofQDataset:
+        class Reduction:
+            corrections = []  # noqa: RUF012
+
+        return OrsoDataset(
+            Orso({}, Reduction, [], name=f'{filename}.orso'), np.ones((0, 0))
+        )
+
+    workflow = sl.Pipeline(
+        [normalized_ioq, orso_dataset], params={Filename[SampleRun]: 'default'}
+    )
+    datasets = orso_datasets_from_measurements(
+        workflow,
+        [{}, {Filename[SampleRun]: 'special'}],
+        scale_to_overlap=False,
+    )
+    assert len(datasets) == 2
+    assert tuple(d.info.name for d in datasets) == ('default.orso', 'special.orso')