Linear Algebra#

These classes are meant to provide user friendly, yet high-performance representations of Vectors and Tensors in any kind of skew or orthogonal coordinate frame.

The Data Model#

The data model is a combination of the two approaches suggested by NumPy for creating NumPy-compliant classes. All Vector and Tensor classes are custom array containers that wrap a direct subclass of NumPy’s ndarray class (see NumPy custom array containers and subclassing ndarray for the details). The double structure allows the container class to manage and arbitrary array object in the baclground while maintaing a unified interface to work with directly. For instance the TopologyArray class maintains either a NumPy or an Awkward array in the background, depending on the input arguments without affecting the way you interact with an instance.

Arrays#

class sigmaepsilon.math.linalg.meta.ArrayWrapper(*args, cls_params=None, contiguous: bool = True, **kwargs)[source]#

Base frontend class for array-like classes. Use it like if it was a numpy.ndarray instance.

chop(tol: float = 1e-12) → ArrayWrapper[source]#

Sets very small values (in an absolute sense) to zero.

New in version 1.0.5.

Parameters:: tol (float, Optional) – The values whose absolute value is less than this limit are set to zero. Default is 1e-12.
Returns:: The object the call was made upon.
Return type:: ArrayWrapper`

property dim: int#: Returns the dimension of the array.

property minmax: Tuple[float]#: Returns the minimum and maximum values of the array.

to_numpy() → ndarray[source]#: Returns the data as a pure NumPy array.

class sigmaepsilon.math.linalg.sparse.JaggedArray(data: Iterable | None = None, *, cuts: Iterable | None = None, force_numpy: bool = True, **kwargs)[source]#

A numba-jittable class that handles 2d matrices with a variable number of columns per row.

The class is actually an interface to awkward.Array, with some additional features, specific to 2d jagged arrays.

At the moment a JaggedArray can be constructed from * a flattened 2d jagged array through unflatteing (we use Awkward here) * a list of 1d lists * a list of 2d NumPy arrays

Parameters:

data (Iterable) – A 2d numpy array or a list of them.
cuts (Iterable, Optional) – An iterable that tells how to unflatten an 1d array into a 2d jagged shape. Only if ‘data’ is an 1d array. Default is None.
force_numpy (bool, Optional) – Forces dense inputs to be NumPy arrays in the background. Default is True.

Examples

The following defines a dense matrix from a flattened shape:

>>> import numpy as np
>>> from sigmaepsilon.math.linalg.sparse import JaggedArray
>>> data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> JaggedArray(data, cuts=[3, 3, 3])
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

Since we use the Awkward library to deflatten arrays, if we provide the argument force_numpy=False, which is the default behaviour, the data is stored as an Awkward array:

>>> JaggedArray(data, cuts=[3, 3, 3], force_numpy=False)
JaggedArray([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

The next example defines a jagged array. In this case the data is inevitably stored as an Awkward array.

>>> data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
>>> JaggedArray(data, cuts=[3, 3, 4])
JaggedArray([[1, 2, 3], [4, 5, 6], [7, 8, 9, 10]])

>>> JaggedArray([[1, 2], [3, 4, 5]])
JaggedArray([[1, 2], [3, 4, 5]])

>>> JaggedArray([np.eye(2), np.eye(3)])
JaggedArray([[[1, 0], [0, 1]], [[1, 0, 0], [0, 1, 0], [0, 0, 1]]])

See also

csr_matrix, awkward.Array

flatten(return_cuts: bool = False)[source]#: Returns the flattened equivalent of the array.

is_jagged() → bool[source]#: Returns True if the topology is jagged, False otherwise.

property shape#: Returns the shape of the data as a tuple. If the topology is jagged, the second item is an iterable.

New in version 0.0.8.

property size: int#: Number of elements in the array.

to_ak() → Array[source]#: Returns underlying data as an Awkward array.

See also

awkward.Array

to_array() → Array | ndarray[source]#: Returns the underlying data, which either an Awkward or a NumPy array.

to_csr() → csr_matrix[source]#: Returns the topology as a csr_matrix.

See also

csr_matrix

to_list()[source]#: Returns underlying data as lists.

New in version 0.0.8.

to_numpy() → ndarray[source]#: Returns underlying data as a NumPy array. This is only possible for regular topologies.

to_scipy() → csr_matrix[source]#: Returns the array as a sparse SciPy CSR matrix.

See also

scipy.sparse.csr_matrix

unique(*args, **kwargs)[source]#: Returns unique elements, by generalizing the functionality provided by numpy.unique(), see its documentation for the details.

widths() → ndarray[source]#: Returns the number of columns for each row.

Numba-jittable Python class for a sparse matrices in CSR format. The meaning of the input variables is the same as in SciPy, and object creation follows the same pattern.

Parameters:

data (SparseLike) – Contains the non-zero values of the matrix, in the order in which they would be encountered if we walked along the rows left to right and top to bottom. If this is a CSC matrix, the walk happens along the columns. From version 0.0.8, Awkward arrays are also accepted. .. versionmodified:: 0.0.8
indices (numpy.ndarray, Optional) – The indices of the columns (rows) during the walk. Default is None.
indptr (numpy.ndarray, Optional) – Stores row (column) boundaries. Default is None.
shape (Tuple, Optional) – Default is None.

Note

1) At the moment, this class does not support NumPy’s array protocoll. If you want this to be the argument to a numpy function, use the to_scipy() method of this class. 2) The attributed ‘data’, ‘indices’, ‘indptr’ and ‘shape’ are all accessible inside Numba-jitted functions.

Examples

Create from a JaggedArray

>>> import numpy as np
>>> from sigmaepsilon.math.linalg.sparse import JaggedArray, csr_matrix
>>> data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
>>> csr = JaggedArray(data, cuts=[3, 3, 4]).to_csr()
>>> csr
3x4 CSR matrix with 10 nonzero values

You can watch it as a NumPy array

>>> csr.to_numpy()
array([[ 1,  2,  3,  0],
       [ 4,  5,  6,  0],
       [ 7,  8,  9, 10]])

Create from a SciPy sparse matrix

>>> from scipy.sparse import csr_matrix as csr_scipy
>>> csr_scipy((3, 4), dtype=np.int8).toarray()
>>> csr_matrix(scipy_matrix)
3x4 CSR matrix with 0 nonzero values

To create the 10 by 10 identity matrix, do this:

>>> csr_matrix.eye(10)
10x10 CSR matrix with 10 nonzero values

You can access rows and row indices of a CSR matrix in Numba jitted code, even in ‘nopython’ mode:

>>> from numba import jit
>>> @jit(nopython=True)
>>> def numba_nopython(csr: csr_matrix, i: int):
>>>     return csr.row(i), csr.irow(i)
>>> row = np.array([0, 0, 1, 2, 2, 2])
>>> col = np.array([0, 2, 2, 0, 1, 2])
>>> data = np.array([1, 2, 3, 4, 5, 6])
>>> matrix = csr_scipy((data, (row, col)), shape=(3, 3))
>>> matrix.toarray()
array([[1, 0, 2],
       [0, 0, 3],
       [4, 5, 6]], dtype=int32)
>>> csr = csr_matrix(matrix)
>>> numba_nopython(csr, 0)
(array([1., 2.]), array([0, 2]))

static eye(N: int) → csr_matrix[source]#: Returns the NxN identity matrix as a CSR matrix.

irow(i: int = 0) → ndarray[source]#: Returns the colum indices of the values of the i-th row.

New in version 0.0.8.

Note

This method is available inside Numba-jitted functions, even in nopython mode.

row(i: int = 0) → ndarray[source]#

Returns the values of the i-th row.

The behavior was changed in version 0.0.8. After that, the call only returns the data related to the i-th row. For the indices see irow().

Note

This method is available inside Numba-jitted functions, even in nopython mode.

to_numpy() → ndarray[source]#: Returns the matrix as a NumPy array. .. versionadded:: 0.0.8

to_scipy() → csr_matrix[source]#: Returns data as a SciPy object.

Frames#

class sigmaepsilon.math.linalg.ReferenceFrame(axes: ndarray | None = None, *args, name: str | None = None, dim: int | None = None)[source]#

A class for arbitrary reference frames, that facilitates transformation of tensor-like quantities. Instances of this class support NumPy’s functions, universal functions and other standard features of NumPy (see the notes below).

An important feature of the class is that it maintains the property of objectivity of the tensorial quantities that use instances of this class in their representation. Upon transformation of a frame, the components of the associated tensorial quantities transform correspondingly.

Parameters:

axes (numpy.ndarray, Optional) – 2d numpy array of floats specifying cartesian reference frames. Default is None.
dim (int, Optional) – Dimension of the mesh. Default is 3.
name (str, Optional) – The name of the frame. Default is None.

Notes

1) The object is able to handle any reference frame, not just cartesian ones. The only restriction is that the base vectors should be linearly independent. 2) All NumPy universal functions return ReferenceFrame instances. If for a NumPy universal function a ReferenceFrame instance is provided by the parameter ‘out’ (eg. numpy.sqrt(A, out=A)), the instance is modified inplace and all associated tensorial quantities change aggordingly. This also true for operators like +=, -=, @=, etc. that naturally cause inplace modification of the called instance.

Examples

Define a frame and rotate it around axis ‘Z’ with an amount of 180 degrees:

>>> from sigmaepsilon.math.linalg import ReferenceFrame
>>> A = ReferenceFrame(dim=3)
>>> B = A.orient_new('Space', [0, 0, np.pi], 'XYZ')

If we define a tensorial quantity with components in a frame, the components of the quantity change if the frame changes:

>>> from sigmaepsilon.math.linalg import Vector
>>> v = Vector([1., 0, 0], frame=A)
>>> A *= 2.0
>>> v
Array([0.5, 0. , 0. ])

A ReferenceFrame instance can be an argument to a NumPy universal function, but the result is always a simple array object. The following call results in a simple array, leaving the frame (and hence the associated tensors, if there are any) unchanged

>>> A = ReferenceFrame(dim=3)
>>> np.sqrt(A)
Array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

but this one changes the frame (and hence the associated tensors if there are any):

>>> A = ReferenceFrame(dim=3)
>>> np.sqrt(A, out=A)
Array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

Of course, if you want your result to be a ReferenceFrame instance, you can always do this:

>>> A = ReferenceFrame(np.sqrt(ReferenceFrame(dim=3)))

The basis vectors that make up a frame are accessible through the property ‘axes’:

>>> A = ReferenceFrame(dim=3)
>>> A.axes
Array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

If the basis vectors change, so do the tensors associated with the frame:

>>> import numpy as np
>>> v = Vector([1.0, 0, 0], frame=A)
>>> A.axes = np.eye(3) * 2.0
>>> v
Array([0.5, 0. , 0. ])

The same happens if we change the frame of a tensorial quantity:

>>> A = ReferenceFrame(dim=3)
>>> v = Vector([1., 0, 0], frame=A)
>>> v.frame = A * 2.0
>>> v
Array([0.5, 0. , 0. ])

To get the matrix that transforms quantities from frame A to frame B, use the dcm method:

>>> source = ReferenceFrame(dim=3)
>>> target = source.orient_new('Space', [0, 0, 90*np.pi/180],  'XYZ')
>>> DCM = source.dcm(target=target)

or equivalenty

>>> DCM = target.dcm(source=source)

See also

RectangularFrame, CartesianFrame

Gram() → ndarray[source]#: Returns the Gram-matrix of the frame.

property axes: ndarray#

Returns a matrix, where each row (or column) is the component array of a basis vector with respect to the ambient frame.

Return type:: numpy.ndarray

copy(deep: bool = False, name: str | None = None) → ReferenceFrame[source]#: Returns a shallow or deep copy of this object, depending of the argument deepcopy (default is False).

dcm(*, target: ReferenceFrame | None = None, source: ReferenceFrame | None = None) → ndarray[source]#

Returns the direction cosine matrix (DCM) of a transformation from a source (S) to a target (T) frame. The current frame can be the source or the target, depending on the arguments.

If called without arguments, it returns the DCM matrix from the ambient frame to the current frame (S=None, T=self).

If source is not None, then T=self.

If target is not None, then S=self.

Parameters:

source ('ReferenceFrame', Optional) – Source frame. Default is None.
target ('ReferenceFrame', Optional) – Target frame. Default is None.

Returns:

DCM matrix from S to T.

Return type:

numpy.ndarray

Example

>>> from sigmaepsilon.math.linalg import ReferenceFrame
>>> import numpy as np
>>> source = ReferenceFrame(dim=3)
>>> target = source.orient_new('Space', [0, 0, 90*np.pi/180],  'XYZ')
>>> DCM = source.dcm(target=target)
>>> arr_source = np.array([3 ** 0.5 / 2, 0.5, 0])
>>> arr_target = DCM @ arr_source

deepcopy(name: str | None = None) → ReferenceFrame[source]#: Returns a deep copy of the frame.

dual() → ReferenceFrame[source]#: Returns the dual (or reciprocal) frame.

classmethod eye(*args, dim=3, **kwargs) → ReferenceFrame[source]#

Returns a standard orthonormal frame.

Return type:: ReferenceFrame

property is_cartesian#: Returns True if the frame is a cartesian (orthonormal) one.

property is_independent: bool#: Returns True if the base vectors that make up the frame are linearly independent.

property is_rectangular#: Returns True if the frame is a rectangular one.

metric_tensor() → TensorLike[source]#: Returns the metric tensor of the frame.

property name: str#: Returns the name of the frame.

orient(*args, **kwargs) → ReferenceFrame[source]#

Orients the current frame inplace. All arguments are forwarded to rotation_matrix(), see there for the details.

Return type:: ReferenceFrame

Example

Define a standard Cartesian frame and rotate it around axis ‘Z’ with 180 degrees:

>>> A = ReferenceFrame(dim=3)
>>> A.orient('Space', [0, 0, np.pi], 'XYZ')
array([[-1.0000000e+00,  1.2246468e-16,  0.0000000e+00],
       [-1.2246468e-16, -1.0000000e+00,  0.0000000e+00],
       [ 0.0000000e+00,  0.0000000e+00,  1.0000000e+00]])

See also

orient_new(), rotation_matrix()

orient_new(*args, name='', **kwargs) → ReferenceFrame[source]#

Returns a new frame, oriented relative to the called object. All extra positional and keyword arguments are forwarded to rotation_matrix(), see there for the details.

Parameters:: name (str) – Name for the new reference frame.
Returns:: A new ReferenceFrame object.
Return type:: ReferenceFrame

See also

orient(), rotation_matrix()

Example

Define a standard Cartesian frame and rotate it around axis ‘Z’ with 180 degrees:

>>> A = ReferenceFrame(dim=3)
>>> B = A.orient_new('Space', [0, 0, np.pi], 'XYZ')

rotate(*args, inplace: bool = True, **kwargs) → ReferenceFrame[source]#

Alias for orient and orient_new, all extra arguments are forwarded.

New in version 0.0.4.

Parameters:: inplace (bool, Optional) – If True, transformation is carried out on the instance the function call is made upon, otherwise a new frame is returned. Default is True.
Returns:: An exisitng (if inplace is True) or a new ReferenceFrame object.
Return type:: ReferenceFrame

Example

Define a standard Cartesian frame and rotate it around axis ‘Z’ with 180 degrees:

>>> A = ReferenceFrame(dim=3)
>>> A.rotate('Space', [0, 0, np.pi], 'XYZ')
array([[-1.0000000e+00,  1.2246468e-16,  0.0000000e+00],
       [-1.2246468e-16, -1.0000000e+00,  0.0000000e+00],
       [ 0.0000000e+00,  0.0000000e+00,  1.0000000e+00]])

See also

orient(), orient_new()

rotate_new(*args, **kwargs) → ReferenceFrame[source]#

Alias for orient_new, all extra arguments are forwarded.

New in version 0.0.4.

Returns:: A new ReferenceFrame object.
Return type:: ReferenceFrame

See also

orient_new()

show(target: ReferenceFrame | None = None) → ndarray[source]#

Returns the components of the current frame in a target frame. If the target is None, the componants are returned in the ambient frame.

Return type:: numpy.ndarray

transpose(inplace: bool = False) → ReferenceFrame[source]#

Either transposes the array of the frame, or returns a copy of it with the components transposed.

Parameters:: inplace (bool, Optional) – If True, the operation is performed on the instance the call is made upon. Default is False.

Note

The rule of transposition differs from the one implemented in NumPy, as only tensorial axes are being transposed.

volume() → float[source]#: Returns the signed volume of the general parallelepiped described by the base vectors that make up the frame.

class sigmaepsilon.math.linalg.RectangularFrame(*args, assume_rectangular: bool = False, **kwargs)[source]#

A class for rectangular reference frames. The behaviour of a RectanguarFrame instance is similar to that of a ReferenceFrame, with minor differences. One is that the rectangular property is utilized wherever possible, resulting in slightly better performance for some operations. The downside is that operations causing inplace modifications are only permitted if it can be guaranteed, that the result is a valid ReferenceFrame object. This is generally not true for function calls like numpy.sqrt(A, out=A) and an exception is raised. However, operators like +=, -= are fine.

Note

A frame being rectangular only implies that the base vectors that make up the frame are mutually perpendicular, but it has no say on the length of the base vectors. That is, a rectangular frame is not necessarily orthonormal.

Examples

For operators *= and /= the type of the instance is unchanged

>>> from sigmaepsilon.math.linalg import ReferenceFrame
>>> A = RectangularFrame(dim=3)
>>> A *= 2
>>> type(A)
<class 'sigmaepsilon.math.linalg.frame.RectangularFrame'>

but for other operations the type of the object is cast to a more general frame

>>> A += 2
>>> type(A)
<class 'sigmaepsilon.math.linalg.frame.ReferenceFrame'>

See also

ReferenceFrame, CartesianFrame

property is_independent: bool#: Returns True if the base vectors that make up the frame are linearly independent.

property is_rectangular#: Returns True if the frame is a rectangular one.

class sigmaepsilon.math.linalg.CartesianFrame(*args, normalize: bool = False, assume_cartesian: bool = False, **kwargs)[source]#

A class for cartesian (orthonormal) reference frames. Just like the RectangularFrame class, this is similar to the more general ReferenceFrame, but with increased performance and even more limitations.

Examples

For operators *= and /= the type of the instance is cast to RectangularFrame

>>> from sigmaepsilon.math.linalg import ReferenceFrame
>>> A = CartesianFrame(dim=3)
>>> A *= 2
>>> type(A)
<class 'sigmaepsilon.math.linalg.frame.RectangularFrame'>

For other operations, the type of the object is cast to a general one

>>> A += 2
>>> type(A)
<class 'sigmaepsilon.math.linalg.frame.ReferenceFrame'>

The only operation that results in a CartesianFrame object is a rotation:

>>> A = CartesianFrame(dim=3)
>>> A.orient('Space', [0, 0, np.pi], 'XYZ')
array([[-1.0000000e+00,  1.2246468e-16,  0.0000000e+00],
       [-1.2246468e-16, -1.0000000e+00,  0.0000000e+00],
       [ 0.0000000e+00,  0.0000000e+00,  1.0000000e+00]])

>>> type(A)
<class 'sigmaepsilon.math.linalg.frame.CartesianFrame'>

See also

ReferenceFrame, RectangularFrame

Gram() → ndarray[source]#: Returns the Gram-matrix of the frame.

dual() → ReferenceFrame[source]#: Returns the dual (or reciprocal) frame.

property is_cartesian#: Returns True if the frame is a cartesian (orthonormal) one.

property is_rectangular#: Returns True if the frame is a rectangular one.

volume() → float[source]#: Returns the signed volume of the general parallelepiped described by the base vectors that make up the frame.

Vectors#

class sigmaepsilon.math.linalg.Vector(*args, frame: FrameLike | None = None, bulk: bool | None = None, rank: int | None = None, **kwargs)[source]#

Extends NumPy’s ndarray class to handle arrays with associated reference frames. The class also provides a mechanism to transform vectors between different frames. Use it like if it was a numpy.ndarray instance.

All parameters are identical to those of numpy.ndarray, except that this class allows to specify an embedding frame.

Parameters:

args (tuple, Optional) – Positional arguments forwarded to numpy.ndarray.
frame (numpy.ndarray, Optional) – The reference frame the vector is represented by its coordinates.
kwargs (dict, Optional) – Keyword arguments forwarded to numpy.ndarray.

Examples

Import the necessary classes:

>>> import numpy as np
>>> from sigmaepsilon.math.linalg import Vector, ReferenceFrame

Create a default frame in 3d space, and create 2 others, each being rotated with 30 degrees around the third axis.

>>> A = ReferenceFrame(dim=3)
>>> B = A.orient_new('Body', [0, 0, 30*np.pi/180], 'XYZ')
>>> C = B.orient_new('Body', [0, 0, 30*np.pi/180], 'XYZ')

To create a vector in a frame:

>>> vA = Vector([1.0, 1.0, 0.0], frame=A)

To create a vector with a relative transformation to another one:

>>> vB = vA.orient_new('Body', [0, 0, -30*np.pi/180], 'XYZ')

Use the array property to get the componets of a Vector:

>>> vB.array
Array([1.3660254, 0.3660254, 0.       ])

If you want to obtain the components of a vector in a specific target frame C, do this:

>>> vB.show(C)
array([ 1., -1.,  0.])

The reason why the result is represented now as ‘array’ insted of ‘Array’ as in the previous case is that the Vector class is an array container. When you type vB.array, what is returned is a wrapped object, an instance of Array, which is also a class of this library. When you say vB.show(C), a NumPy array is returned. Since the Array class is a direct subclass of NumPy’s ndarray class, it doesn’t really matter and the only difference is the capital first letter.

To create a vector in a target frame C, knowing the components in a source frame A:

>>> vC = Vector(vA.show(C), frame=C)

See also

Tensor, ReferenceFrame

copy(deep: bool = False, name: str | None = None) → Vector[source]#: Returns a shallow or deep copy of this object, depending of the argument deepcopy (default is False).

deepcopy(name: str | None = None) → Vector[source]#: Returns a deep copy of the frame.

dual() → Vector[source]#: Returns the vector described in the dual (or reciprocal) frame.

orient(*args, dcm: ndarray | None = None, **kwargs) → Vector[source]#

Orients the vector inplace. If the transformation is not specified by ‘dcm’, all arguments are forwarded to orient_new.

Parameters:: dcm (numpy.ndarray, Optional) – The DCM matrix of the transformation.
Returns:: The same vector the function is called upon.
Return type:: Vector

See also

orient_new()

orient_new(*args, **kwargs) → Vector[source]#

Returns a transformed version of the instance.

Returns:: A new vector.
Return type:: Vector

See also

orient()

property rank: int#: Returns the tensor rank (or order).

show(target: ReferenceFrame | None = None, *, dcm: ndarray | None = None) → ndarray[source]#

Returns the components in a target frame. If the target is None, the components are returned in the ambient frame.

The transformation can also be specified with a proper DCM matrix.

Parameters:

target (numpy.ndarray, Optional) – Target frame.
dcm (numpy.ndarray, Optional) – The DCM matrix of the transformation.

Returns:

The components of the vector in a specified frame, or the ambient frame, depending on the arguments.

Return type:

numpy.ndarray

Tensors#

class sigmaepsilon.math.linalg.Tensor(*args, frame: FrameLike | None = None, bulk: bool | None = None, rank: int | None = None, **kwargs)[source]#

A class to handle tensors.

Parameters:

args (tuple, Optional) – Positional arguments forwarded to numpy.ndarray.
frame (numpy.ndarray, Optional) – The reference frame the vector is represented by its coordinates.
kwargs (dict, Optional) – Keyword arguments forwarded to numpy.ndarray.

Examples

Import the necessary classes:

>>> from sigmaepsilon.math.linalg import Tensor, ReferenceFrame
>>> from numpy.random import rand

Create a Tensor of order 6 in a frame with random components

>>> frame = ReferenceFrame(dim=3)
>>> array = rand(3, 3, 3, 3, 3, 3)
>>> A = Tensor(array, frame=frame)

Get the tensor in the dual frame:

>>> A_dual = A.dual()

Create an other tensor, in this case a 5th-order one, and calculate their generalized dot product, which is a 9th-order tensor:

>>> from sigmaepsilon.math.linalg import dot
>>> array = rand(3, 3, 3, 3, 3)
>>> B = Tensor(array, frame=frame)
>>> C = dot(A, B, axes=[0, 0])
>>> assert C.rank == (A.rank + B.rank - 2)

See also

Vector, ReferenceFrame

copy(deep: bool = False, name: str | None = None) → Tensor[source]#: Returns a shallow or deep copy of this object, depending of the argument deepcopy (default is False).

deepcopy(name: str | None = None) → Tensor[source]#: Returns a deep copy of the frame.

dual() → Tensor2[source]#: Returns the tensor described in the dual (or reciprocal) frame.

orient(*args, **kwargs) → Tensor[source]#

Orients the vector inplace. All arguments are forwarded to orient_new.

Returns:: The same vector the function is called upon.
Return type:: Vector

See also

orient_new()

orient_new(*args, **kwargs) → Tensor[source]#

Returns a transformed version of the instance.

Returns:: A new vector.
Return type:: Vector

See also

orient()

show(target: ReferenceFrame | None = None, *, dcm: ndarray | None = None) → ndarray[source]#

Returns the components in a target frame. If the target is None, the components are returned in the ambient frame.

The transformation can also be specified with a proper DCM matrix.

Parameters:

target (numpy.ndarray, Optional) – Target frame.
dcm (numpy.ndarray, Optional) – The DCM matrix of the transformation.

Returns:

The components of the tensor in a specified frame, or the ambient frame, depending on the arguments.

Return type:

numpy.ndarray

transform_components(Q: ndarray) → ndarray[source]#: Returns the components of the tensor transformed by the matrix Q.

class sigmaepsilon.math.linalg.Tensor2(*args, frame: FrameLike | None = None, bulk: bool | None = None, rank: int | None = None, **kwargs)[source]#

A class to handle second-order tensors. Some operations have dedicated implementations that provide higher performence utilizing implicit parallelization. Examples for tensors of this class include the metric tensor, or the stress and strain tensors of elasticity.

See also

:class:”~sigmaepsilon.math.linalg.tensor.Tensor2x3”

transform_components(Q: ndarray) → ndarray[source]#: Returns the components of the tensor transformed by the matrix Q.

class sigmaepsilon.math.linalg.Tensor4(*args, frame: FrameLike | None = None, bulk: bool | None = None, rank: int | None = None, **kwargs)[source]#

A class to handle fourth-order tensors. Some operations have dedicated implementations that provide higher performence utilizing implicit parallelization. Examples of this class include the piezo-optical tensor, the elasto-optical tensor, the flexoelectric tensor or the elasticity tensor.

See also

:class:”~sigmaepsilon.math.linalg.tensor.Tensor4x3”

transform_components(dcm: ndarray) → ndarray[source]#: Returns the components of the transformed numerical tensor, based on the provided direction cosine matrix.

Routines#

sigmaepsilon.math.linalg.logical.has_full_column_rank(matrix: ndarray) → bool[source]#: Returns True if the input matrix has full column rank, ie if all its columns are linearly independent.

See also

numpy.linalg.matrix_rank()

sigmaepsilon.math.linalg.logical.has_full_rank(matrix: ndarray) → bool[source]#: Returns True if the input matrix has full rank.

See also

numpy.linalg.matrix_rank()

sigmaepsilon.math.linalg.logical.has_full_row_rank(matrix: ndarray) → bool[source]#: Returns True if the input matrix has full row rank, ie if all its rows are linearly independent.

See also

numpy.linalg.matrix_rank()

sigmaepsilon.math.linalg.logical.is_hermitian(arr: ndarray) → bool[source]#: Returns True if the input is a hermitian array.

sigmaepsilon.math.linalg.logical.is_independent_frame(axes: ndarray, tol: float = 0) → bool[source]#

Returns True if a the base vectors of a frame are linearly independent.

Parameters:: axes (numpy.ndarray) – A matrix where the i-th row is the i-th basis vector.

sigmaepsilon.math.linalg.logical.is_normal_frame(axes: ndarray) → bool[source]#

Returns True if a frame is normal, meaning, that it’s base vectors are all of unit length.

Parameters:: axes (numpy.ndarray) – A matrix where the i-th row is the i-th basis vector.

sigmaepsilon.math.linalg.logical.is_orthonormal_frame(axes: ndarray) → bool[source]#

Returns True if a frame is orthonormal.

Parameters:: axes (numpy.ndarray) – A matrix where the i-th row is the i-th basis vector.

sigmaepsilon.math.linalg.logical.is_pos_def(arr) → bool[source]#: Returns True if the input is positive definite.

sigmaepsilon.math.linalg.logical.is_pos_semidef(arr) → bool[source]#: Returns True if the input is positive semi definite.

sigmaepsilon.math.linalg.logical.is_rectangular_frame(axes: ndarray) → bool[source]#

Returns True if a frame is Cartesian.

Parameters:: axes (numpy.ndarray) – A matrix where the i-th row is the i-th basis vector.

sigmaepsilon.math.linalg.utils.Gram(axes: ndarray) → ndarray[source]#

Returns the Gram matrix of a frame.

Parameters:: axes (numpy.ndarray) – A matrix where the i-th row is the i-th basis vector.

sigmaepsilon.math.linalg.utils.cross(a: TensorLike | ArrayWrapper, b: TensorLike | ArrayWrapper, *args, frame: FrameLike | None = None, **kwargs) → TensorLike | ndarray[source]#

Calculates the cross product of two vectors or one vector and a second order tensor. The behaviour coincides with NumPy when all inputs are arrays and generalizes when they are not, but all inputs must be either all arrays or all tensors of some kind.

Parameters:

*args (Tuple, Optional) – Positional arguments forwarded to NumPy, if all input objects are arrays.
a (TensorLike or ArrayLike) – A tensor or an array.
b (TensorLike or ArrayLike) – A tensor or an array.
frame (FrameLike, Optional) – The target frame of the output. Only if all inputs are TensorLike. If not specified, the returned tensor migh be returned in an arbitrary frame, depending on the inputs. Default is None.
**kwargs (dict, Optional) – Keyword arguments forwarded to numpy.cross. As of NumPy version ‘1.22.4’, there are no keyword arguments for numpy.cross, this is to assure compliance with all future versions of numpy.

Returns:

An 1d or 2d array, or an 1d or 2d tensor, depending on the inputs.

Return type:

numpy.ndarray or TensorLike

References

https://mathworld.wolfram.com/CrossProduct.html

Examples

The cross product of two vectors results in a vector:

>>> from sigmaepsilon.math.linalg import ReferenceFrame, Vector, Tensor2
>>> from sigmaepsilon.math.linalg import cross
>>> frame = ReferenceFrame(np.eye(3))
>>> a = Vector(np.array([1., 0, 0]), frame=frame)
>>> b = Vector(np.array([0, 1., 0]), frame=frame)
>>> cross(a, b)
Array([0., 0., 1.])

The cross product of a second order tensor and a vector result a second order tensor:

>>> A = Tensor2(np.eye(3), frame=frame)
>>> cross(A, b)
Array([[ 0.,  0., -1.],
       [ 0.,  0.,  0.],
       [ 1.,  0.,  0.]])

Returns the dot product (without complex conjugation) of two quantities. The behaviour coincides with NumPy when all inputs are arrays and generalizes when they are not, but all inputs must be either all arrays or all tensors of some kind. The operation for tensors of order 1 and 2 have dedicated implementations, for higher order tensors it generalizes to tensor contraction along specified axes.

Parameters:

a (TensorLike or ArrayLike) – A tensor or an array.
b (TensorLike or ArrayLike) – A tensor or an array.
out (ArrayLike, Optional) – Output argument. This must have the exact kind that would be returned if it was not used. See numpy.dot for the details. Only if all inputs are ArrayLike. Default is None.
frame (FrameLike, Optional) – The target frame of the output. Only if all inputs are TensorLike. If not specified, the returned tensor migh be returned in an arbitrary frame, depending on the inputs. Default is None.
axes (tuple or list, Optional) – The indices along which contraction happens if any of the input tensors have a rank higher than 2. Default is None.

Returns:

An array or a tensor, depending on the inputs.

Return type:

TensorLike or numpy.ndarray or scalar

Notes

For general tensors, the current implementation has an upper limit considering the rank of the input tensors. The sum of the ranks of the input tensors plus the sum of contraction indices must be at most 26.

References

https://mathworld.wolfram.com/DotProduct.html

Examples

When working with NumPy arrays, the behaviour coincides with numpy.dot. To take the dot product of a 2nd order tensor and a vector, use it like this:

>>> from sigmaepsilon.math.linalg import ReferenceFrame, Vector, Tensor2
>>> from sigmaepsilon.math.linalg import dot
>>> frame = ReferenceFrame(np.eye(3))
>>> A = Tensor2(np.eye(3), frame=frame)
>>> v = Vector(np.array([1., 0, 0]), frame=frame)
>>> dot(A, v)
Array([1., 0., 0.])

For general tensors, you have to specify the axes along which contraction happens:

>>> from sigmaepsilon.math.linalg import Tensor
>>> A = Tensor(np.ones((3, 3, 3, 3)), frame=frame)  # a tensor of order 4
>>> B = Tensor(np.ones((3, 3, 3)), frame=frame)  # a tensor of order 3
>>> dot(A, B, axes=(0, 0)).rank
5

sigmaepsilon.math.linalg.utils.dual_frame(axes: ndarray) → ndarray[source]#

Returns the dual frame of the input.

Parameters:: axes (numpy.ndarray) – A matrix where the i-th row is the i-th basis vector.

sigmaepsilon.math.linalg.utils.generalized_inverse(matrix: ndarray) → ndarray[source]#: Returns the generalized inverse of the input matrix, in any of the following cases: 1) The matrix is square and has full rank. In this case the returned matrix is the usual inverse. 2) The matrix has more columns than rows and has full row rank. In this case the generalized right inverse is returned. 3) The matrix has more rows than columns and has full column rank. In this case the generalized left inverse is returned.

sigmaepsilon.math.linalg.utils.generalized_left_inverse(matrix: ndarray) → ndarray[source]#: Returns the generalized left inverse

\begin{equation} \left( \mathbf{A}^{T} \mathbf{A} \right)^{-1} \mathbf{A}^{T} \end{equation}

sigmaepsilon.math.linalg.utils.generalized_right_inverse(matrix: ndarray) → ndarray[source]#: Returns the generalized right inverse

\begin{equation} \mathbf{A}^{T} \left( \mathbf{A} \mathbf{A}^{T} \right)^{-1} \end{equation}

sigmaepsilon.math.linalg.utils.normalize_frame(axes: ndarray) → ndarray[source]#

Returns the frame with normalized base vectors.

Parameters:: axes (numpy.ndarray) – A matrix where the i-th row is the i-th basis vector.

sigmaepsilon.math.linalg.utils.permutation_tensor(dim: int = 3) → ndarray[source]#

Returns the Levi-Civita pseudotensor for N dimensions.

Parameters:: N (int, Optional) – The number of dimensions. Default is 3.

sigmaepsilon.math.linalg.utils.random_pos_semidef_matrix(N) → ndarray[source]#

Returns a random positive semidefinite matrix of shape (N, N).

Example

>>> from sigmaepsilon.math.linalg import random_pos_semidef_matrix, is_pos_semidef
>>> arr = random_pos_semidef_matrix(2)
>>> is_pos_semidef(arr)
True

sigmaepsilon.math.linalg.utils.random_posdef_matrix(N, alpha: float = 1e-12) → ndarray[source]#

Returns a random positive definite matrix of shape (N, N).

All eigenvalues of this matrix are >= alpha.

Example

>>> from sigmaepsilon.math.linalg import random_posdef_matrix, is_pos_def
>>> arr = random_posdef_matrix(2)
>>> is_pos_def(arr)
True

sigmaepsilon.math.linalg.utils.rotation_matrix(rot_type: str, amounts: Iterable, rot_order: str | int = '') → ndarray[source]#

Returns a rotation matrix using the mechanism provided by sympy.physics.vector.ReferenceFrame.orientnew.

Parameters:

rot_type (str) –
The method used to generate the direction cosine matrix. Supported methods are:
- 'Axis': simple rotations about a single common axis
- 'DCM': for setting the direction cosine matrix directly
- 'Body': three successive rotations about new intermediate
  axes, also called “Euler and Tait-Bryan angles”
- 'Space': three successive rotations about the parent
  frames’ unit vectors
- 'Quaternion': rotations defined by four parameters which
  result in a singularity free direction cosine matrix
amounts (Iterable) –
Expressions defining the rotation angles or direction cosine matrix. These must match the rot_type. See examples below for details. The input types are:
- 'Axis': 2-tuple (expr/sym/func, Vector)
- 'DCM': Matrix, shape(3, 3)
- 'Body': 3-tuple of expressions, symbols, or functions
- 'Space': 3-tuple of expressions, symbols, or functions
- 'Quaternion': 4-tuple of expressions, symbols, or
  functions
rot_order (str or int, Optional) – If applicable, the order of the successive of rotations. The string '123' and integer 123 are equivalent, for example. Required for 'Body' and 'Space'.

Returns:

A new ReferenceFrame object.

Return type:

ReferenceFrame

See also

sympy.physics.vector.ReferenceFrame.orientnew()

Example

Define a standard Cartesian frame and rotate it around axis ‘Z’ with 180 degrees:

>>> from sigmaepsilon.math.linalg.utils import rotation_matrix
>>> R = rotation_matrix('Space', [0, 0, np.pi], 'XYZ')

sigmaepsilon.math.linalg.utils.show_vector(dcm: ndarray, arr: ndarray) → ndarray[source]#

Returns the coordinates of a single or multiple vectors in a frame specified by one or several DCM matrices. The function can handle the following scenarios:

a single (1d) vector and a single (2d) dcm matrix (trivial case)

a stack of vectors (2d) and a single (2d) dcm matrix

a stack of fectors (2d) and dcm matrices for each vector in the stack (3d)

New in version 1.0.5.

Parameters:

dcm (numpy.ndarray) – The dcm matrix of the transformation as a 2d or 3d float array.
arr (numpy.ndarray) – 1d or 2d float array of coordinates of a single vector. If it is 2d, then it is assumed that the coordinates of the i-th vector are accessible as arr[i].

Returns:

The new coordinates with the same shape as arr.

Return type:

numpy.ndarray