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

combine with or

attr

typically ignored or dropped

In the special case of in-place operations such as += or *= Scipp preserves existing attributes and ignores attributes of the right-hand-side.

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]:

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.

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


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.

Coords are compared in operations with datasets or data arrays (or items of datasets). Operations fail if there is any mismatch in coord or label 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'], values=rng.random(2)),
'd': sc.scalar(value=1.0)},
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]:

try:
d2 += d1
except Exception as e:
print(str(e))

"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
try:
d3['a'] -= d3['x', 1:2]['b'] # fail due to coordinate mismatch
except Exception as e:
print(str(e))

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


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:

[18]:

rad = 3.141593*sc.units.rad
deg = 180.0*sc.units.deg
print(sc.sin(deg))
try:

<scipp.Variable> ()    float64  [dimensionless]  -3.4641e-07