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:
[1]:
import scipp as sc
import numpy as np
vecs = sc.vectors(dims=['x'], unit='m', values=np.arange(2 * 3).reshape(2, 3))
vecs
[1]:
- (x: 2)vector_3_float64m[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
:
[2]:
vecs.values
[2]:
array([[0., 1., 2.],
[3., 4., 5.]])
Access to individual vector components x
, y
, and z
is provided using the special fields
property:
[3]:
vecs.fields.x
[3]:
- (x: 2)float64m0.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:
[4]:
vecs.fields.x = 123.0 * sc.units.m
vecs.fields.y += 0.123 * sc.units.m
vecs.values
[4]:
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:
[5]:
np.random.seed(0) # matrices are not singular for this seed
mats = sc.matrices(dims=['x'], values=np.random.rand(2, 3, 3))
inv = sc.matrices(dims=['x'], values=np.linalg.inv(mats.values))
(mats * inv).values.round(decimals=12)
[5]:
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:
[6]:
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:
[7]:
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.matrix
to create a rotation matrix (in this case to rotate by 30 deg around the y
axis) and apply it to the positions:
[8]:
rotation = sc.matrix(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 from_rotvec
and as_rotvec
functions as part of the scipp.spatial.transform
module. Let’s reverse our original rotation by applying a -30 deg around the y
axis.
[9]:
from scipp.spatial.transform import from_rotvec
rotation_back = from_rotvec(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.
[10]:
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.
[11]:
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
[11]:
- ()vector_3_float64[ 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
[12]:
from scipp.spatial.transform import as_rotvec
sc.isclose(rotation * sc.vector(value=[0, 0, 1]), norm)
[12]:
- ()boolTrue
Values:
array(True)