Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make wavelength bands truly optional #53

Merged
merged 7 commits into from
Jan 23, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 5 additions & 12 deletions docs/examples/sans2d.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -60,10 +60,8 @@
"params[DirectBeamFilename] = 'DIRECT_SANS2D_REAR_34327_4m_8mm_16Feb16.hdf5'\n",
"params[OutFilename] = 'reduced.nxs'\n",
"\n",
"band = sc.linspace('wavelength', 2.0, 16.0, num=2, unit='angstrom')\n",
"params[WavelengthBands] = band\n",
"params[WavelengthBins] = sc.linspace(\n",
" 'wavelength', start=band[0], stop=band[-1], num=141\n",
" 'wavelength', start=2.0, stop=16.0, num=141, unit='angstrom'\n",
")\n",
"\n",
"params[sans.sans2d.LowCountThreshold] = sc.scalar(100, unit='counts')\n",
Expand Down Expand Up @@ -229,6 +227,8 @@
"metadata": {},
"outputs": [],
"source": [
"from esssans.isis import plot_flat_detector_xy\n",
"\n",
"monitors = (\n",
" WavelengthMonitor[SampleRun, Incident],\n",
" WavelengthMonitor[SampleRun, Transmission],\n",
Expand All @@ -243,14 +243,7 @@
"\n",
"display(sc.plot({str(key): results[key] for key in monitors}, norm='log'))\n",
"\n",
"display(\n",
" scn.instrument_view(\n",
" results[MaskedData[SampleRun]].hist(),\n",
" pixel_size=0.0075,\n",
" norm='log',\n",
" camera=pp.graphics.Camera(position=(0, 0, 22)),\n",
" )\n",
")\n",
"display(plot_flat_detector_xy(results[MaskedData[SampleRun]].sum('tof'), norm='log'))\n",
"\n",
"parts = {str(key): results[key] for key in parts}\n",
"parts = {key: val if val.bins is None else val.hist() for key, val in parts.items()}\n",
Expand Down Expand Up @@ -312,7 +305,7 @@
"We can also compute $I(Q)$ inside a set of wavelength bands, instead of using the full wavelength range in one go.\n",
"This is useful for debugging purposes.\n",
"\n",
"To achieve this, we need to turn the `WavelengthBands` parameter into a two-dimensional variable,\n",
"To achieve this, we need to supply the `WavelengthBands` parameter (as a two-dimensional variable),\n",
"representing the wavelength range for each band."
]
},
Expand Down
7 changes: 2 additions & 5 deletions docs/examples/zoom.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -100,10 +100,8 @@
"params[NeXusMonitorName[Incident]] = 'monitor3'\n",
"params[NeXusMonitorName[Transmission]] = 'monitor5'\n",
"\n",
"band = sc.linspace('wavelength', 1.75, 16.5, num=2, unit='angstrom')\n",
"params[WavelengthBands] = band\n",
"params[WavelengthBins] = sc.geomspace(\n",
" 'wavelength', start=band[0], stop=band[-1], num=141\n",
" 'wavelength', start=1.75, stop=16.5, num=141, unit='angstrom'\n",
")\n",
"\n",
"params[QBins] = sc.geomspace(dim='Q', start=0.004, stop=0.8, num=141, unit='1/angstrom')\n",
Expand Down Expand Up @@ -271,8 +269,7 @@
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"pygments_lexer": "ipython3"
}
},
"nbformat": 4,
Expand Down
15 changes: 5 additions & 10 deletions src/esssans/beam_center_finder.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,6 @@
QBins,
SampleRun,
UncertaintyBroadcastMode,
WavelengthBands,
WavelengthBins,
)

Expand Down Expand Up @@ -129,7 +128,7 @@ def _iofq_in_quadrants(
norm: sc.DataArray,
graph: dict,
q_bins: Union[int, sc.Variable],
wavelength_range: sc.Variable,
wavelength_bins: sc.Variable,
transform: sc.Variable,
pixel_shape: sc.DataGroup,
) -> Dict[str, sc.DataArray]:
Expand All @@ -148,8 +147,8 @@ def _iofq_in_quadrants(
Coordinate transformation graph.
q_bins:
Bin edges for Q.
wavelength_range:
The wavelength range to use for computing the intensity as a function of Q.
wavelength_bins:
The binning in wavelength to use for computing the intensity as a function of Q.

Returns
-------
Expand Down Expand Up @@ -177,7 +176,7 @@ def _iofq_in_quadrants(
]
params = {}
params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound
params[WavelengthBands] = wavelength_range
params[WavelengthBins] = wavelength_bins
params[QBins] = q_bins
params[DetectorPixelShape[SampleRun]] = pixel_shape
params[LabFrameTransform[SampleRun]] = transform
Expand Down Expand Up @@ -405,10 +404,6 @@ def beam_center_from_iofq(
com_shift = beam_center_from_center_of_mass(data, graph)
logger.info(f'Initial guess for beam center: {com_shift}')

wavelength_range = sc.concat(
[wavelength_bins.min(), wavelength_bins.max()], dim='wavelength'
)

coords = data.transform_coords(
['cylindrical_x', 'cylindrical_y'], graph=graph
).coords
Expand All @@ -421,7 +416,7 @@ def beam_center_from_iofq(
res = minimize(
_cost,
x0=[com_shift.fields.x.value, com_shift.fields.y.value],
args=(data, norm, graph, q_bins, wavelength_range, transform, pixel_shape),
args=(data, norm, graph, q_bins, wavelength_bins, transform, pixel_shape),
bounds=bounds,
method=minimizer,
tol=tolerance,
Expand Down
21 changes: 12 additions & 9 deletions src/esssans/i_of_q.py
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,7 @@ def resample_direct_beam(

def _process_wavelength_bands(
wavelength_bands: Optional[WavelengthBands],
wavelength_bins: WavelengthBins,
) -> Optional[WavelengthBands]:
"""
Perform some checks and potential reshaping on the wavelength bands.
Expand All @@ -143,7 +144,9 @@ def _process_wavelength_bands(
and an end wavelength.
"""
if wavelength_bands is None:
return wavelength_bands
wavelength_bands = sc.concat(
[wavelength_bins.min(), wavelength_bins.max()], dim='wavelength'
)
if wavelength_bands.ndim == 1:
wavelength_bands = sc.concat(
[wavelength_bands[:-1], wavelength_bands[1:]], dim='x'
Expand All @@ -165,6 +168,7 @@ def _process_wavelength_bands(
def merge_spectra(
data: CleanQ[RunType, IofQPart],
q_bins: QBins,
wavelength_bins: WavelengthBins,
wavelength_bands: Optional[WavelengthBands],
dims_to_keep: Optional[DimsToKeep],
) -> CleanSummedQ[RunType, IofQPart]:
Expand All @@ -180,6 +184,8 @@ def merge_spectra(
A DataArray containing the data that is to be converted to Q.
q_bins:
The binning in Q to be used.
wavelength_bins:
The binning in wavelength to be used.
wavelength_bands:
Defines bands in wavelength that can be used to separate different wavelength
ranges that contribute to different regions in Q space. Note that this needs to
Expand All @@ -198,7 +204,9 @@ def merge_spectra(
if dims_to_keep is not None:
dims_to_reduce -= set(dims_to_keep)

wavelength_bands = _process_wavelength_bands(wavelength_bands)
wavelength_bands = _process_wavelength_bands(
wavelength_bands=wavelength_bands, wavelength_bins=wavelength_bins
)

if data.bins is not None:
out = _events_merge_spectra(
Expand Down Expand Up @@ -231,17 +239,14 @@ def _events_merge_spectra(
data_q: sc.DataArray,
q_bins: Union[int, sc.Variable],
dims_to_reduce: List[str],
wavelength_bands: Optional[sc.Variable] = None,
wavelength_bands: sc.Variable,
) -> sc.DataArray:
"""
Merge spectra of event data
"""
q_all_pixels = data_q.bins.concat(dims_to_reduce)
edges = _to_q_bins(q_bins)
q_binned = q_all_pixels.bin(**edges)
if wavelength_bands is None:
return q_binned

dim = 'wavelength'
wav_binned = q_binned.bin({dim: sc.sort(wavelength_bands.flatten(to=dim), dim)})
# At this point we kind of already have what we need, would be cheapest to just
Expand All @@ -263,14 +268,12 @@ def _dense_merge_spectra(
data_q: sc.DataArray,
q_bins: Union[int, sc.Variable],
dims_to_reduce: List[str],
wavelength_bands: Optional[sc.Variable] = None,
wavelength_bands: sc.Variable,
) -> sc.DataArray:
"""
Merge spectra of dense data
"""
edges = _to_q_bins(q_bins)
if wavelength_bands is None:
return data_q.hist(**edges).sum(dims_to_reduce)
bands = []
band_dim = (set(wavelength_bands.dims) - {'wavelength'}).pop()
for wav_range in sc.collapse(wavelength_bands, keep='wavelength').values():
Expand Down
17 changes: 14 additions & 3 deletions tests/sans2d/reduction_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,10 +34,8 @@

def make_params() -> dict:
params = default_parameters.copy()
band = sc.linspace('wavelength', 2.0, 16.0, num=2, unit='angstrom')
params[WavelengthBands] = band
params[WavelengthBins] = sc.linspace(
'wavelength', start=band[0], stop=band[-1], num=141
'wavelength', start=2.0, stop=16.0, num=141, unit='angstrom'
)
params[WavelengthMask] = sc.DataArray(
data=sc.array(dims=['wavelength'], values=[True]),
Expand Down Expand Up @@ -100,6 +98,19 @@ def test_pipeline_can_compute_background_subtracted_IofQ_in_wavelength_slices():
assert result.sizes['band'] == 10


def test_pipeline_wavelength_bands_is_optional():
params = make_params()
pipeline = sciline.Pipeline(sans2d_providers(), params=params)
noband = pipeline.compute(BackgroundSubtractedIofQ)
with pytest.raises(sciline.UnsatisfiedRequirement):
pipeline.compute(WavelengthBands)
Comment on lines +105 to +106
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I opened scipp/sciline#108 for this.

band = sc.linspace('wavelength', 2.0, 16.0, num=2, unit='angstrom')
pipeline[WavelengthBands] = band
assert sc.identical(band, pipeline.compute(WavelengthBands))
withband = pipeline.compute(BackgroundSubtractedIofQ)
assert sc.identical(noband, withband)


def test_workflow_is_deterministic():
params = make_params()
params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop
Expand Down