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:

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

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:

result = a / b
result.values

array([[0.32595151, 0.43159311, 9.77228585, 1.98839842],
       [1.19798487, 1.51186315, 3.74187126, 0.40141215]])

a.values / np.transpose(b.values)

array([[0.32595151, 0.43159311, 9.77228585, 1.98839842],
       [1.19798487, 1.51186315, 3.74187126, 0.40141215]])

result.variances

array([[   1.44120506,    1.81529253, 1352.12494744,   37.4199719 ],
       [  17.24465468,    8.08572378,  345.0888559 ,    4.54435748]])

a.variances / np.transpose(b.variances)

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:

var_xy = sc.array(dims=['x', 'y'], values=np.arange(6).reshape((2, 3)))
print(var_xy.values)

[[0 1 2]
 [3 4 5]]

var_y = sc.array(dims=['y'], values=np.arange(3))
print(var_y.values)

[0 1 2]

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:

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

array([[0.93198836, 0.72478136, 0.86055132, 0.9293378 ],
       [0.54618601, 0.93767296, 0.49498794, 0.27377318]])

a -= a['x', 1]
a.values

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:

sc.show(a['x', 1])
sc.show(a['y', 1])
sc.show(a['x', 1] + a['y', 1])

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:

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.

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')},
)

d1 += d2

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}>."

d3 = d1 + d2
for name in d3:
    print(name)

a
b

d3['a'] -= d3['b']  # transposing
d3['a'] -= d3['x', 1]['b']  # broadcasting

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.

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)}
)

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:

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

False True

Range-slices preserve alignment:

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

True

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

da2_sliced = da2['x', 1]
da1_sliced + da2_sliced

scipp.DataArray (1.03 KB)

• Dimensions:
  • y: 2
• Coordinates: (1)
  • y
    (y)
    int64
    𝟙
    0, 1

    Values:
    array([0, 1])

• Data: (1)
  • (y)
    int64
    𝟙
    3, 4

    Values:
    array([3, 4])

Aligned coordinates take precedence over unaligned if both are present:

da2.coords.set_aligned('x', False)
da1 + da2

scipp.DataArray (1.36 KB)

• Dimensions:
  • x: 4
  • y: 2
• Coordinates: (2)
  • x
    (x)
    int64
    𝟙
    0, 1, 2, 3

    Values:
    array([0, 1, 2, 3])

  • y
    (y)
    int64
    𝟙
    0, 1

    Values:
    array([0, 1])

• Data: (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.

da1 = sc.DataArray(sc.arange('x', 4), coords={'x': sc.arange('x', 4)})
da2 = da1.copy()
da2.coords.set_aligned('x', False)

sc.identical(da1.coords['x'], da2.coords['x'])

True

sc.identical(da1, da2)

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:

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.