Linear Algebra

Basics

Scipp supports basic linear algebra operations in 3 dimensions using vectors of length 3 and several transformations. Variables containing vectors or transformations are created using

• scipp.vectors (for an array of vectors)

• scipp.vector (for a single vector)

• functions in scipp.spatial (for transformations), e.g.,

  • scipp.spatial.rotations (for and array of rotations)

  • scipp.spatial.linear_transform (for a single linear transform, i.e. rotation + scaling)

Let us consider an example, creating a 1-D variable containing vectors:

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

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

scipp.Variable (360 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 (272 Bytes out of 304 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 three-dimensional scatter plot 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.spatial.as_vectors(x, y, 0.0 * sc.units.m)},
)
pp.scatter3d(sensor, pos="position", pixel_size=0.01, cbar=True)

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

pp.scatter3d(sensor, pos="position", pixel_size=0.01, cbar=True)

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:]
pp.scatter3d(sensor, pos="position", pixel_size=0.01, cbar=True)

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:]
pp.scatter3d(sensor, pos="position", pixel_size=0.01, cbar=True)

We can compound rotations by multiplying 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']
pp.scatter3d(sensor, pos="position", pixel_size=0.01, cbar=True)

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 (336 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 (257 Bytes)

• • ()
    bool
    True

    Values:
    array(True)