Computation#

General concepts and mechanisms#

Overview#

Binary operations between data arrays or datasets behave as follows:

Property

Action

coord

compare, abort on mismatch

data

apply operation

mask

combine with or

Dimension matching and transposing#

Operations “align” variables based on their dimension labels. That is, an operation between two variables that have a transposed memory layout behave correctly:

[1]:
import numpy as np

import scipp as sc

rng = np.random.default_rng(12345)
a = sc.array(
    dims=['x', 'y'], values=rng.random((2, 4)), variances=rng.random((2, 4)), unit='m'
)
b = sc.array(
    dims=['y', 'x'], values=rng.random((4, 2)), variances=rng.random((4, 2)), unit='s'
)
a / b
[1]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (384 Bytes)
    • (x: 2, y: 4)
      float64
      m/s
      0.326, 0.432, ..., 3.742, 0.401
      σ = 1.201, 1.347, ..., 18.577, 2.132
      Values:
      array([[0.32595151, 0.43159311, 9.77228585, 1.98839842], [1.19798487, 1.51186315, 3.74187126, 0.40141215]])

      Variances (σ²):
      array([[ 1.44120506, 1.81529253, 1352.12494744, 37.4199719 ], [ 17.24465468, 8.08572378, 345.0888559 , 4.54435748]])

Propagation of uncertainties#

If variables have variances, operations correctly propagate uncertainties (the variances), in contrast to a naive implementation using NumPy:

[2]:
result = a/b
result.values
[2]:
array([[0.32595151, 0.43159311, 9.77228585, 1.98839842],
       [1.19798487, 1.51186315, 3.74187126, 0.40141215]])
[3]:
a.values/np.transpose(b.values)
[3]:
array([[0.32595151, 0.43159311, 9.77228585, 1.98839842],
       [1.19798487, 1.51186315, 3.74187126, 0.40141215]])
[4]:
result.variances
[4]:
array([[   1.44120506,    1.81529253, 1352.12494744,   37.4199719 ],
       [  17.24465468,    8.08572378,  345.0888559 ,    4.54435748]])
[5]:
a.variances/np.transpose(b.variances)
[5]:
array([[2.5251612 , 4.8723756 , 2.70819165, 1.11013763],
       [0.81791708, 0.74070147, 0.73816117, 1.47348507]])

The implementation assumes uncorrelated data and is otherwise based on, e.g., Wikipedia: Propagation of uncertainty. See also Propagation of uncertainties for the concrete equations used for error propagation.

Note:

If an operand with variances is also broadcast in an operation then the resulting values would be correlated. Scipp has no way of tracking or handling such correlations. Subsequent operations that combine values of the result, such as computing the mean, would therefore result in underestimated uncertainties.

To avoid such silent underestimation of uncertainties, operations raise VariancesError when an operand with variances is implicitly or explicitly broadcast in an operations. See our publication Systematic underestimation of uncertainties by widespread neutron-scattering data-reduction software for more background.

Broadcasting#

Missing dimensions in the operands are automatically broadcast. Consider:

[6]:
var_xy = sc.array(dims=['x', 'y'], values=np.arange(6).reshape((2,3)))
print(var_xy.values)
[[0 1 2]
 [3 4 5]]
[7]:
var_y = sc.array(dims=['y'], values=np.arange(3))
print(var_y.values)
[0 1 2]
[8]:
var_xy -= var_y
print(var_xy.values)
[[0 0 0]
 [3 3 3]]

Since var_y did not depend on dimension 'x' it is considered as “constant” along that dimension. That is, the same var_y values are subtracted from all slices of dimension 'x' in var_xy.

Given two variables a and b, we see that broadcasting integrates seamlessly with slicing and transposing:

[9]:
a = sc.array(dims=['x', 'y'], values=rng.random((2, 4)), unit='m')
b = sc.array(dims=['y', 'x'], values=rng.random((4, 2)), unit='s')
a.values
[9]:
array([[0.93198836, 0.72478136, 0.86055132, 0.9293378 ],
       [0.54618601, 0.93767296, 0.49498794, 0.27377318]])
[10]:
a -= a['x', 1]
a.values
[10]:
array([[ 0.38580235, -0.2128916 ,  0.36556338,  0.65556462],
       [ 0.        ,  0.        ,  0.        ,  0.        ]])

Both operands may be broadcast, creating an output with the combination of input dimensions:

[11]:
sc.show(a['x', 1])
sc.show(a['y', 1])
sc.show(a['x', 1] + a['y', 1])
dims=('y',), shape=(4,), unit=m, variances=Falsevalues y
dims=('x',), shape=(2,), unit=m, variances=Falsevalues x
dims=('y', 'x'), shape=(4, 2), unit=m, variances=Falsevalues yx

Note that in-place operations such as += will never change the shape of the left-hand-side. That is only the right-hand-side operation can be broadcast, and the operation fails of a broadcast of the left-hand-side would be required.

Units#

Units are required to be compatible:

[12]:
try:
    a + b
except Exception as e:
    print(str(e))
Cannot add m and s.

Coordinate and name matching#

In operations with datasets, data items are paired based on their names when applying operations to datasets. Operations fail if names do not match:

  • In-place operations such as += accept a right-hand-side operand that omits items that the left-hand-side has. If the right-hand-side contains items that are not in the left-hand-side the operation fails.

  • Non-in-place operations such as + return a new dataset with items from the intersection of the inputs.

Coordinates are compared in operations with datasets or data arrays (or items of datasets). Operations fail if there is any mismatch in aligned (see below) coord values.

[13]:
d1 = sc.Dataset(
    data={'a': sc.array(dims=['x', 'y'], values=rng.random((2, 3))),
          'b': sc.array(dims=['y', 'x'], values=rng.random((3, 2))),
          'c': sc.array(dims=['x', 'y'], values=rng.random((2, 3)))},
    coords={
        'x': sc.arange('x', 2.0, unit='m'),
        'y': sc.arange('y', 3.0, unit='m')})
d2 = sc.Dataset(
    data={'a': sc.array(dims=['x', 'y'], values=rng.random((2, 3))),
          'b': sc.array(dims=['y', 'x'], values=rng.random((3, 2)))},
    coords={
        'x': sc.arange('x', 2.0, unit='m'),
        'y': sc.arange('y', 3.0, unit='m')})
[14]:
d1 += d2
[15]:
d2 += d1
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[15], line 1
----> 1 d2 += d1

KeyError: "Expected 'c' in <scipp.Dataset.keys {a, b}>."
[16]:
d3 = d1 + d2
for name in d3:
    print(name)
a
b
[17]:
d3['a'] -= d3['b']  # transposing
d3['a'] -= d3['x', 1]['b']  # broadcasting
[18]:
d3['a'] -= d3['x', 1:2]['b']  # fail due to coordinate mismatch
---------------------------------------------------------------------------
DatasetError                              Traceback (most recent call last)
Cell In[18], line 1
----> 1 d3['a'] -= d3['x', 1:2]['b']  # fail due to coordinate mismatch

DatasetError: Mismatch in coordinate 'x' in operation 'subtract_equals':
(x: 2)    float64              [m]  [0, 1]
vs
(x: 1)    float64              [m]  [1]

Alignment#

Coordinates are “aligned” be default. This means that they are required to match between operands and if they don’t, an exception is raised.

[19]:
da1 = sc.DataArray(sc.arange('x', 8).fold('x', {'x': 4, 'y': 2}),
                   coords={'x': sc.arange('x', 4), 'y': sc.arange('y', 2)})
da2 = sc.DataArray(sc.arange('x', 4),
                   coords={'x': 10*sc.arange('x', 4), 'y': sc.arange('y', 2)})
[20]:
da1 + da2
---------------------------------------------------------------------------
DatasetError                              Traceback (most recent call last)
Cell In[20], line 1
----> 1 da1 + da2

DatasetError: Mismatch in coordinate 'x' in operation 'add':
(x: 4)      int64  [dimensionless]  [0, 1, 2, 3]
vs
(x: 4)      int64  [dimensionless]  [0, 10, 20, 30]

Coordinates can become “unaligned” in some operations. For example, slicing out a single element makes all coordinates in the sliced dimension unaligned:

[21]:
da1_sliced = da1['x', 1]
print(da1_sliced.coords['x'].aligned, da1_sliced.coords['y'].aligned)
False True

Range-slices preserve alignment:

[22]:
da1['x', 0:1].coords['x'].aligned
[22]:
True

Unaligned coordinates are not required to match in operations. They are instead dropped if there is a mismatch:

[23]:
da2_sliced = da2['x', 1]
da1_sliced + da2_sliced
[23]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (1.03 KB)
    • y: 2
    • y
      (y)
      int64
      𝟙
      0, 1
      Values:
      array([0, 1])
    • (y)
      int64
      𝟙
      3, 4
      Values:
      array([3, 4])

Aligned coordinates take precedence over unaligned if both are present:

[24]:
da2.coords.set_aligned('x', False)
da1 + da2
[24]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (1.36 KB)
    • x: 4
    • y: 2
    • x
      (x)
      int64
      𝟙
      0, 1, 2, 3
      Values:
      array([0, 1, 2, 3])
    • y
      (y)
      int64
      𝟙
      0, 1
      Values:
      array([0, 1])
    • (x, y)
      int64
      𝟙
      0, 1, ..., 9, 10
      Values:
      array([[ 0, 1], [ 3, 4], [ 6, 7], [ 9, 10]])

Note:

sc.identical and similar comparison functions ignore alignment when applied to variables directly. They do, however take it into account when comparing coordinates of data arrays and datasets.

[25]:
da1 = sc.DataArray(sc.arange('x', 4), coords={'x': sc.arange('x', 4)})
da2 = da1.copy()
da2.coords.set_aligned('x', False)
[26]:
sc.identical(da1.coords['x'], da2.coords['x'])
[26]:
True
[27]:
sc.identical(da1, da2)
[27]:
False

Arithmetic#

The arithmetic operations +, -, *, and / and their in-place variants +=, -=, *=, and /= are available for variables, data arrays, and datasets. They can also be combined with slicing.

Trigonometrics#

Trigonometric functions like sin are supported for variables. Units for angles provide a safeguard and ensure correct operation when working with either degree or radian:

[28]:
rad = 3.141593*sc.units.rad
deg = 180.0*sc.units.deg
print(sc.sin(rad))
print(sc.sin(deg))
try:
    rad + deg
except Exception as e:
    print(str(e))
<scipp.Variable> ()    float64  [dimensionless]  -3.4641e-07
<scipp.Variable> ()    float64  [dimensionless]  1.22465e-16
Cannot add rad and deg.

Other#

See the list of free functions for an overview of other available operations.