Linear Algebra

Basics

Scipp supports basic linear algebra operations in 3 dimensions using vectors of length 3 and 3x3 matrices (to define rotations). Variables containing vectors or matrices are created using scipp.vectors (for an array of vectors) scipp.vector (for a single vector) scipp.matrices (for and array of matrices) scipp.matrix (for a single matrix). Let us consider an example, creating a 1-D variable containg vectors:

import scipp as sc
import scipp.spatial
import numpy as np

vecs = sc.vectors(dims=['x'], unit='m', values=np.arange(2 * 3).reshape(2, 3))
vecs

scipp.Variable (48 Bytes)

• • (x: 2)
    vector3
    m
    [0. 1. 2.], [3. 4. 5.]

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

While vecs has only a single dimension x, the values property exposes the underlying vector value as a 2-D array of float64:

vecs.values

array([[0., 1., 2.],
       [3., 4., 5.]])

Access to individual vector components x, y, and z is provided using the special fields property:

vecs.fields.x

scipp.Variable (16 Bytes out of 48 Bytes)

• • (x: 2)
    float64
    m
    0.0, 3.0

    Values:
    array([0., 3.])

These properties return a view to the underlying data and thus support setting element values as well as in-place modification:

vecs.fields.x = 123.0 * sc.units.m
vecs.fields.y += 0.123 * sc.units.m
vecs.values

array([[123.   ,   1.123,   2.   ],
       [123.   ,   4.123,   5.   ]])

This works in an identical manner for matrices with the properties xx, xy, xz, yx, yy, yz, zx, zy, and zz. The values property allows for use of, e.g., numpy functionality. We may, e.g., compute the inverse using numpy.linalg.inv, which is not supported in scipp yet:

np.random.seed(0)  # matrices are not singular for this seed
mats = sc.spatial.linear_transforms(dims=['x'], values=np.random.rand(2, 3, 3))
inv = sc.spatial.linear_transforms(dims=['x'], values=np.linalg.inv(mats.values))
(mats * inv).values.round(decimals=12)

array([[[ 1., -0.,  0.],
        [-0.,  1.,  0.],
        [-0.,  0.,  1.]],

       [[ 1., -0.,  0.],
        [ 0.,  1., -0.],
        [-0., -0.,  1.]]])

Example

We create some random data and use positions to create a plot with projection='3d' to illustrate some concepts. Image this as an image sensor with 12x12 pixels:

nx = 12
ny = 12
x = sc.linspace(dim='x', start=-0.1, stop=0.1, num=nx, unit='m')
y = sc.linspace(dim='y', start=-0.1, stop=0.1, num=ny, unit='m')
sensor = sc.DataArray(
    data=sc.array(dims=['x','y'], values=np.random.rand(nx, ny)),
    coords={'position':sc.geometry.position(x, y, 0.0*sc.units.m)}
)
sc.plot(sensor, projection="3d", positions="position")

We can use the vector-element access properties in combination with slicing to offset some of the pixels:

sensor.coords['position']['x',5:].fields.x += 0.1 * sc.units.m
sensor.coords['position']['y',5:].fields.y += 0.05 * sc.units.m

sc.plot(sensor, projection="3d", positions="position")

We use sc.spatial.linear_transform to create a rotation matrix (in this case to rotate by 30 deg around the y axis) and apply it to the positions:

rotation = sc.spatial.linear_transform(value=[[ 0.8660254, 0.0000000, 0.5000000],
                            [ 0.0000000, 1.0000000, 0.0000000],
                            [-0.5000000, 0.0000000, 0.8660254]])

sensor.coords['position']['x',5:] = rotation * sensor.coords['position']['x',5:]
sc.plot(sensor, projection="3d", positions="position")

More conveniently, we can create our rotations directly from a rotation vector. Scipp provides the rotations_from_rotvecs function as part of the scipp.spatial module. Let’s reverse our original rotation by applying a -30 deg around the y axis.

from scipp.spatial import rotations_from_rotvecs
rotation_back = rotations_from_rotvecs(rotation_vectors=sc.vector(value=[0, -30.0, 0], unit=sc.units.deg))
sensor.coords['position']['x',5:] = rotation_back * sensor.coords['position']['x',5:]
sc.plot(sensor, projection="3d", positions="position")

We can compound rotations by multipying our rotation variable. Lets make rotate the same again in the same axis to give an overall 60 deg around the y axis, applied to every pixel in our sensor.

rotation = rotation * rotation
sensor.coords['position'] = rotation * sensor.coords['position']
sc.plot(sensor, projection="3d", positions="position")

Scipp also provides vector and scalar products of vectors. Lets determine the surface normal vector for our sensor in the new rotation given the knowledge that all points are coplanar.

a = sensor.coords['position']['x', 0]['y', 0]
b = sensor.coords['position']['x', -1]['y', 0]
c = sensor.coords['position']['x', 0]['y', -1]
norm = sc.cross(c - a, c - b)
norm /= sc.norm(norm)
norm

scipp.Variable (24 Bytes)

• • ()
    vector3
    𝟙
    [ 0.86602541 -0. 0.5 ]

    Values:
    array([ 0.86602541, -0.        ,  0.5       ])

Now lets verify that the normal vector we just calculated corresponds to a rotation of our original normal vector (0, 0, 1) by 60 degrees around y

sc.isclose(rotation * sc.vector(value=[0, 0, 1]), norm)

scipp.Variable (1 Bytes)

• • ()
    bool
    True

    Values:
    array(True)