macromax.backend package
This package defines a general interface to do parallel operations efficiently on different architectures,
and it provides specific implementations that are used by the macromax.Solution
class.
The specific implementation that is used can be:
- macromax.backend.config(*args, **kwargs)[source]
Configure a specific back-end, overriding the automatic detection. E.g. use from macromax import backend followed by:
backend.config(type='numpy') backend.config(dict(type='numpy')) backend.config(type='torch', device='cpu') backend.config(dict(type='torch', device='cpu')) backend.config(dict(type='torch', device='cuda')) backend.config(dict(type='torch', device='cuda')) backend.config([dict(type='torch', device='cuda'), dict(type='numpy')])
Default back-ends can be specified in the file backend_config.json in the current folder. This configuration file should contain a list of potential back-ends in [JSON format](https://en.wikipedia.org/wiki/JSON), e.g. ```json [
{“type”: “torch”, “device”: “cuda”}, {“type”: “numpy”}, {“type”: “torch”, “device”: “cpu”}
]
The first back-end that loads correctly will be used.
- macromax.backend.load(nb_pol_dims, grid, dtype, config_list=None)[source]
Load the default or the backend specified using backend.config(). This configuration file should contain a list of potential back-ends in [JSON format](https://en.wikipedia.org/wiki/JSON), e.g. ```json [
{“type”: “torch”, “device”: “cuda”}, {“type”: “numpy”}, {“type”: “torch”, “device”: “cpu”}
]
The first back-end that loads correctly will be used.
- type nb_pol_dims:
- param nb_pol_dims:
The number of polarization dimensions: 1 for scalar, 3 for vectorial calculations.
- type grid:
- param grid:
The uniformly-spaced Cartesian calculation grid as a Grid object.
- type dtype:
- param dtype:
The scalar data type. E.g. np.complex64 or np.complex128.
- type config_list:
- param config_list:
List of alternative backend configurations.
- rtype:
- return:
A BackEnd object to start the calculation with.
- class macromax.backend.BackEnd(nb_dims, grid, hardware_dtype=<class 'numpy.complex128'>)[source]
Bases:
ABC
A class that provides methods to work with arrays of matrices or block-diagonal matrices, represented as ndarrays, where the first two dimensions are those of the matrix, and the final dimensions are the coordinates over which the operations are parallelized and the Fourier transforms are applied.
- __init__(nb_dims, grid, hardware_dtype=<class 'numpy.complex128'>)[source]
Construct object to handle parallel operations on square matrices of nb_rows x nb_rows elements. The matrices refer to points in space on a uniform plaid grid.
- property vector_length: int
The shape of the square matrix that transforms a single vector in the set. This is a pair of identical integer numbers.
- property ft_axes: tuple
The integer indices of the spatial dimensions, i.e. on which to Fourier Transform.
- property vectorial: bool
A boolean indicating if this object represents a vector space (as opposed to a scalar space).
- property hardware_dtype
The scalar data type that is processed by this back-end.
- property numpy_dtype
The equivalent hardware data type in numpy
- abstract allocate_array(shape=None, dtype=None, fill_value=None)[source]
Allocates a new vector array of shape grid.shape and word-aligned for efficient calculations.
- Parameters:
shape (
Union
[Complex
,Sequence
,ndarray
,None
]) – (optional) The shape of the returned array.dtype (
Union
[Type
[float64
],Type
[float32
],None
]) – (optional) The scalar data type of the elements in the array. This may be converted to an equivalent, back-end, specific data type.fill_value (
Optional
[Complex
]) – (optional) A scalar value to pre-populate the array. The default (None) does not pre-populate the array, leaving it in a random state.
- Return type:
- Returns:
A reference to the array.
- assign(arr, out)[source]
Assign the values of one array to those of another. The dtype and shape is broadcasted to match that of the output array.
- assign_exact(arr, out)[source]
Assign the values of one array to those of another. The dtype and shape must match that of the output array.
- sign(arr)[source]
Returns an array with values of -1 where arr is negative, 0 where arr is 0, and 1 where arr is positive.
- static expand_dims(arr, axis)[source]
Inserts a new singleton axis at the indicated position, thus increasing ndim by 1.
- Return type:
- property eye: ndarray
Returns an identity tensor that can be multiplied using singleton expansion. This can be useful for scalar additions or subtractions.
- Returns:
an array with the number of dimensions matching that of the BackEnds’s data set.
- allclose(arr, other=0.0)[source]
Returns True if all elements in arr are close to other.
- Return type:
- abstract ft(arr)[source]
Calculates the discrete Fourier transform over the spatial dimensions of E. The computational complexity is that of a Fast Fourier Transform:
O(N.log(N))
.
- abstract ift(arr)[source]
Calculates the inverse Fourier transform over the spatial dimensions of E. The computational complexity is that of a Fast Fourier Transform:
O(N.log(N))
. The scaling is so thatE == self.ift(self.ft(E))
- convolve(operation_ft, arr)[source]
Apply an operator in Fourier space. This is used to perform a cyclic FFT convolution which overwrites its input argument arr.
- Parameters:
- Return type:
- Returns:
the convolved input array.
- is_matrix(arr)[source]
Checks if an ndarray is a matrix as defined by this parallel_ops_column object.
- to_matrix_field(arr)[source]
Converts the input to an array of the full number of dimensions: len(self.matrix_shape) + len(self.grid.shape), and dtype. The size of each dimensions must match that of that set for the
BackEnd
or be 1 (assumes broadcasting). For electric fields in 3-space, self.matrix_shape == (N, N) == (3, 3) The first (left-most) dimensions of the output are either1x1: The identity matrix for a scalar field, as sound waves or isotropic permittivity.
Nx1: A vector for a vector field, as the electric field.
NxN: A matrix for a matrix field, as anisotropic permittivity
None is interpreted as 0. Singleton dimensions are added on the left so that all dimensions are present. Inputs with 1xN are transposed (not conjugate) to Nx1 vectors.
- Parameters:
arr (
Union
[Complex
,Sequence
,ndarray
]) – The input can be scalar, which assumes that its value is assumed to be repeated for all space. The value can be a one-dimensional vector, in which case the vector is assumed to be repeated for all space.- Return type:
- Returns:
An array with ndim == len(self.matrix_shape) + len(self.grid.shape) and with each non-singleton dimension matching those of the nb_rows and data_shape.
- subtract(left_term, right_term)[source]
Point-wise difference of A and B.
- Parameters:
- Return type:
- Returns:
The point-wise difference of both sets of matrices. Singleton dimensions are expanded.
- mul(left_factor, right_factor, out=None)[source]
Point-wise matrix multiplication of A and B. Overwrites right_factor!
- Parameters:
left_factor (
Union
[Complex
,Sequence
,ndarray
]) – The left matrix array, must start with dimensions n x mright_factor (
Union
[Complex
,Sequence
,ndarray
]) – The right matrix array, must have matching or singleton dimensions to those of A, bar the first two dimensions. In case of missing dimensions, singletons are assumed. The first dimensions must be m x p. Where the m matches that of the left hand matrix unless both m and p are 1 or both n and m are 1, in which case the scaled identity is assumed.out (
Optional
[ndarray
]) – (optional) The destination array for the results.
- Return type:
- Returns:
An array of matrix products with all but the first two dimensions broadcast as needed.
- abstract ldivide(denominator, numerator=1.0)[source]
Parallel matrix left division, A^{-1}B, on the final two dimensions of A and B result_lm = A_kl B_km
A and B must have have all but the final dimension identical or singletons. B defaults to the identity matrix.
- curl(field_array)[source]
Calculates the curl of a vector E with the final dimension the vector dimension. The input argument may be overwritten!
- curl_ft(field_array_ft)[source]
Calculates the Fourier transform of the curl of a Fourier transformed E with the final dimension the vector dimension. The final dimension of the output will always be of length 3; however, the input length may be shorter, in which case the missing values are assumed to be zero. The first dimension of the input array corresponds to the first element in the final dimension, if it exists, the second dimension corresponds to the second element etc.
- cross(A, B)[source]
Calculates the cross product of vector arrays A and B.
- Parameters:
- Return type:
- Returns:
A vector array of dimensions
[vector_length, 1, *data_shape]
containing the cross product A x B in the first dimension and the other dimensions remain on the same axes.
- outer(A, B)[source]
Calculates the Dyadic product of vector arrays A and B.
- Parameters:
- Return type:
- Returns:
A matrix array of dimensions
[vector_length, vector_length, *data_shape]
containing the dyadic product \(A \otimes B\) in the first two dimensions and the other dimensions remain on the same axes.
- div(field_array)[source]
Calculate the divergence of the input vector or tensor field E. The input argument may be overwritten!
- div_ft(field_array_ft)[source]
Calculate the Fourier transform of the divergence of the pre-Fourier transformed input electric field.
- Parameters:
field_array_ft (
Union
[Complex
,Sequence
,ndarray
]) – The input array representing the field pre-Fourier-transformed in all spatial dimensions. The input is of the shape[m, n, x, y, z, ...]
- Return type:
- Returns:
The Fourier transform of the divergence of the field in the shape
[n, 1, x, y, z]
.
- transversal_projection(field_array)[source]
Projects vector arrays onto their transverse component. The input argument may be overwritten!
- longitudinal_projection(field_array)[source]
Projects vector arrays onto their longitudinal component. The input argument may be overwritten!
- transversal_projection_ft(field_array_ft)[source]
Projects the Fourier transform of a vector E array onto its transversal component.
- longitudinal_projection_ft(field_array_ft)[source]
Projects the Fourier transform of a vector array onto its longitudinal component. Overwrites self.array_ft_input!
- property k2: ndarray
Helper def for calculation of the Fourier transform of the Green def
- Returns:
\(|k|^2\) for the specified sample grid and output shape
- mat3_eigh(arr)[source]
Calculates the eigenvalues of the 3x3 Hermitian matrices represented by A and returns a new array of 3-vectors, one for each matrix in A and of the same dimensions, baring the second dimension. When the first two dimensions are 3x1 or 1x3, a diagonal matrix is assumed. When the first two dimensions are singletons (1x1), a constant diagonal matrix is assumed and only one eigenvalue is returned. Returns an array of one dimension less: 3 x data_shape. With the exception of the first dimension, the shape is maintained.
Before substituting this for
`numpy.linalg.eigvalsh`
, note that this implementation is about twice as fast as the current numpy implementation for 3x3 Hermitian matrix fields. The difference is even greater for the PyTorch implementation. The implementation below can be made more efficient by limiting it to Hermitian matrices and thus real eigenvalues. In the end, only the maximal eigenvalue is required, so an iterative power iteration may be even faster, perhaps after applying a Given’s rotation or Householder reflection to make Hermitian tridiagonal.- Parameters:
arr (
Union
[Complex
,Sequence
,ndarray
]) – The set of 3x3 input matrices for which the eigenvalues are requested. This must be an ndarray with the first two dimensions of size 3.- Return type:
- Returns:
The set of eigenvalue-triples contained in an ndarray with its first dimension of size 3, and the remaining dimensions equal to all but the first two input dimensions.
- calc_roots_of_low_order_polynomial(C)[source]
Calculates the (complex) roots of polynomials up to order 3 in parallel. The coefficients of the polynomials are in the first dimension of C and repeated in the following dimensions, one for each polynomial to determine the roots of. The coefficients are input for low to high order in each column. In other words, the polynomial is:
np.sum(C * (x**range(C.size))) == 0
This method is used by mat3_eigh, but only for polynomials with real roots.
Submodules
macromax.backend.numpy module
The module providing the pure-python NumPy back-end implementation.
- class macromax.backend.numpy.BackEndNumpy(nb_dims, grid, hardware_dtype=<class 'numpy.complex128'>)[source]
Bases:
BackEnd
A class that provides methods to work with arrays of matrices or block-diagonal matrices, represented as ndarrays, where the first two dimensions are those of the matrix, and the final dimensions are the coordinates over which the operations are parallelized and the Fourier transforms are applied.
- __init__(nb_dims, grid, hardware_dtype=<class 'numpy.complex128'>)[source]
Construct object to handle parallel operations on square matrices of nb_rows x nb_rows elements. The matrices refer to points in space on a uniform plaid grid.
- allocate_array(shape=None, dtype=None, fill_value=None)[source]
Allocates a new vector array of shape grid.shape and word-aligned for efficient calculations.
- Return type:
- ft(arr)[source]
Calculates the discrete Fourier transform over the spatial dimensions of E. The computational complexity is that of a Fast Fourier Transform:
O(N\log(N))
.
- ift(arr)[source]
Calculates the inverse Fourier transform over the spatial dimensions of E. The computational complexity is that of a Fast Fourier Transform:
O(N\log(N))
. The scaling is so thatE == self.ift(self.ft(E))
- abs(arr)
Returns the absolute value (magnitude) of the elements in the input.
- adjoint(mat)
Transposes the elements of individual matrices with complex conjugation.
- asnumpy(arr)
Convert the internal array (or tensor) presentation to a numpy.ndarray.
- assign(arr, out)
Assign the values of one array to those of another. The dtype and shape is broadcasted to match that of the output array.
- assign_exact(arr, out)
Assign the values of one array to those of another. The dtype and shape must match that of the output array.
- astype(arr, dtype=None)
- calc_roots_of_low_order_polynomial(C)
Calculates the (complex) roots of polynomials up to order 3 in parallel. The coefficients of the polynomials are in the first dimension of C and repeated in the following dimensions, one for each polynomial to determine the roots of. The coefficients are input for low to high order in each column. In other words, the polynomial is:
np.sum(C * (x**range(C.size))) == 0
This method is used by mat3_eigh, but only for polynomials with real roots.
- static clear_cache()
- conj(arr)
Returns the conjugate of the elements in the input.
- convolve(operation_ft, arr)
Apply an operator in Fourier space. This is used to perform a cyclic FFT convolution which overwrites its input argument arr.
- Parameters:
- Return type:
- Returns:
the convolved input array.
- cross(A, B)
Calculates the cross product of vector arrays A and B.
- Parameters:
- Return type:
- Returns:
A vector array of dimensions
[vector_length, 1, *data_shape]
containing the cross product A x B in the first dimension and the other dimensions remain on the same axes.
- curl(field_array)
Calculates the curl of a vector E with the final dimension the vector dimension. The input argument may be overwritten!
- curl_ft(field_array_ft)
Calculates the Fourier transform of the curl of a Fourier transformed E with the final dimension the vector dimension. The final dimension of the output will always be of length 3; however, the input length may be shorter, in which case the missing values are assumed to be zero. The first dimension of the input array corresponds to the first element in the final dimension, if it exists, the second dimension corresponds to the second element etc.
- div(field_array)
Calculate the divergence of the input vector or tensor field E. The input argument may be overwritten!
- div_ft(field_array_ft)
Calculate the Fourier transform of the divergence of the pre-Fourier transformed input electric field.
- Parameters:
field_array_ft (
Union
[Complex
,Sequence
,ndarray
]) – The input array representing the field pre-Fourier-transformed in all spatial dimensions. The input is of the shape[m, n, x, y, z, ...]
- Return type:
- Returns:
The Fourier transform of the divergence of the field in the shape
[n, 1, x, y, z]
.
- static evaluate_polynomial(C, X)
Evaluates the polynomial P at X for testing.
- static expand_dims(arr, axis)
Inserts a new singleton axis at the indicated position, thus increasing ndim by 1.
- Return type:
- property eye: ndarray
Returns an identity tensor that can be multiplied using singleton expansion. This can be useful for scalar additions or subtractions.
- Returns:
an array with the number of dimensions matching that of the BackEnds’s data set.
- property ft_axes: tuple
The integer indices of the spatial dimensions, i.e. on which to Fourier Transform.
- property hardware_dtype
The scalar data type that is processed by this back-end.
- inv(mat)
Inverts the set of input matrices M.
- is_matrix(arr)
Checks if an ndarray is a matrix as defined by this parallel_ops_column object.
- is_scalar(arr)
Tests if A represents a scalar field (as opposed to a vector field).
- is_vector(arr)
Tests if A represents a vector field.
- property k2: ndarray
Helper def for calculation of the Fourier transform of the Green def
- Returns:
\(|k|^2\) for the specified sample grid and output shape
- ldivide(denominator, numerator=1.0)[source]
Parallel matrix left division, A^{-1}B, on the final two dimensions of A and B result_lm = A_kl B_km
A and B must have have all but the final dimension identical or singletons. B defaults to the identity matrix.
- longitudinal_projection(field_array)
Projects vector arrays onto their longitudinal component. The input argument may be overwritten!
- longitudinal_projection_ft(field_array_ft)
Projects the Fourier transform of a vector array onto its longitudinal component. Overwrites self.array_ft_input!
- mat3_eigh(arr)
Calculates the eigenvalues of the 3x3 Hermitian matrices represented by A and returns a new array of 3-vectors, one for each matrix in A and of the same dimensions, baring the second dimension. When the first two dimensions are 3x1 or 1x3, a diagonal matrix is assumed. When the first two dimensions are singletons (1x1), a constant diagonal matrix is assumed and only one eigenvalue is returned. Returns an array of one dimension less: 3 x data_shape. With the exception of the first dimension, the shape is maintained.
Before substituting this for
`numpy.linalg.eigvalsh`
, note that this implementation is about twice as fast as the current numpy implementation for 3x3 Hermitian matrix fields. The difference is even greater for the PyTorch implementation. The implementation below can be made more efficient by limiting it to Hermitian matrices and thus real eigenvalues. In the end, only the maximal eigenvalue is required, so an iterative power iteration may be even faster, perhaps after applying a Given’s rotation or Householder reflection to make Hermitian tridiagonal.- Parameters:
arr (
Union
[Complex
,Sequence
,ndarray
]) – The set of 3x3 input matrices for which the eigenvalues are requested. This must be an ndarray with the first two dimensions of size 3.- Return type:
- Returns:
The set of eigenvalue-triples contained in an ndarray with its first dimension of size 3, and the remaining dimensions equal to all but the first two input dimensions.
- mul(left_factor, right_factor, out=None)
Point-wise matrix multiplication of A and B. Overwrites right_factor!
- Parameters:
left_factor (
Union
[Complex
,Sequence
,ndarray
]) – The left matrix array, must start with dimensions n x mright_factor (
Union
[Complex
,Sequence
,ndarray
]) – The right matrix array, must have matching or singleton dimensions to those of A, bar the first two dimensions. In case of missing dimensions, singletons are assumed. The first dimensions must be m x p. Where the m matches that of the left hand matrix unless both m and p are 1 or both n and m are 1, in which case the scaled identity is assumed.out (
Optional
[ndarray
]) – (optional) The destination array for the results.
- Return type:
- Returns:
An array of matrix products with all but the first two dimensions broadcast as needed.
- property numpy_dtype
The equivalent hardware data type in numpy
- outer(A, B)
Calculates the Dyadic product of vector arrays A and B.
- Parameters:
- Return type:
- Returns:
A matrix array of dimensions
[vector_length, vector_length, *data_shape]
containing the dyadic product \(A \otimes B\) in the first two dimensions and the other dimensions remain on the same axes.
- sign(arr)
Returns an array with values of -1 where arr is negative, 0 where arr is 0, and 1 where arr is positive.
- subtract(left_term, right_term)
Point-wise difference of A and B.
- Parameters:
- Return type:
- Returns:
The point-wise difference of both sets of matrices. Singleton dimensions are expanded.
- to_matrix_field(arr)
Converts the input to an array of the full number of dimensions: len(self.matrix_shape) + len(self.grid.shape), and dtype. The size of each dimensions must match that of that set for the
BackEnd
or be 1 (assumes broadcasting). For electric fields in 3-space, self.matrix_shape == (N, N) == (3, 3) The first (left-most) dimensions of the output are either1x1: The identity matrix for a scalar field, as sound waves or isotropic permittivity.
Nx1: A vector for a vector field, as the electric field.
NxN: A matrix for a matrix field, as anisotropic permittivity
None is interpreted as 0. Singleton dimensions are added on the left so that all dimensions are present. Inputs with 1xN are transposed (not conjugate) to Nx1 vectors.
- Parameters:
arr (
Union
[Complex
,Sequence
,ndarray
]) – The input can be scalar, which assumes that its value is assumed to be repeated for all space. The value can be a one-dimensional vector, in which case the vector is assumed to be repeated for all space.- Return type:
- Returns:
An array with ndim == len(self.matrix_shape) + len(self.grid.shape) and with each non-singleton dimension matching those of the nb_rows and data_shape.
- transversal_projection(field_array)
Projects vector arrays onto their transverse component. The input argument may be overwritten!
- transversal_projection_ft(field_array_ft)
Projects the Fourier transform of a vector E array onto its transversal component.
macromax.backend.tensorflow module
The module providing the TensorFlow back-end implementation.
- class macromax.backend.tensorflow.BackEndTensorFlow(nb_dims, grid, hardware_dtype=tensorflow.complex128, device=None, address=None)[source]
Bases:
BackEnd
A class that provides methods to work with arrays of matrices or block-diagonal matrices, represented as ndarrays, where the first two dimensions are those of the matrix, and the final dimensions are the coordinates over which the operations are parallelized and the Fourier transforms are applied.
- __init__(nb_dims, grid, hardware_dtype=tensorflow.complex128, device=None, address=None)[source]
Construct object to handle parallel operations on square matrices of nb_rows x nb_rows elements. The matrices refer to points in space on a uniform plaid grid.
- Parameters:
nb_dims (
int
) – The number of rows and columns in each matrix. 1 for scalar operations, 3 for polarizationgrid (
Grid
) – The grid that defines the position of the matrices.hardware_dtype – (optional) The data type to use for operations.
device (
Optional
[str
]) – (optional) ‘cpu’, ‘gpu’, or ‘tpu’ to indicate where the calculation will happen.
- property numpy_dtype
The equivalent hardware data type in numpy
- astype(arr, dtype=None)[source]
As necessary, convert the ndarray arr to the type dtype.
- Return type:
Tensor
- assign(arr, out)[source]
Assign the values of one array to those of another. The dtype and shape is broadcasted to match that of the output array.
- Parameters:
arr – The values that need to be assigned to another array.
out – (optional) The target array, which must be of the same shape.
- Return type:
Tensor
- Returns:
The target array or a newly allocated array with the same values as those in arr.
- assign_exact(arr, out)[source]
Assign the values of one array to those of another. The dtype and shape must match that of the output array.
- Parameters:
arr – The values that need to be assigned to another array.
out (
Tensor
) – (optional) The target array, which must be of the same shape.
- Return type:
Tensor
- Returns:
The target array or a newly allocated array with the same values as those in arr.
- allocate_array(shape=None, dtype=None, fill_value=None)[source]
Allocates a new vector array of shape grid.shape and word-aligned for efficient calculations.
- Return type:
Tensor
- sign(arr)[source]
Returns an array with values of -1 where arr is negative, 0 where arr is 0, and 1 where arr is positive.
- swapaxes(arr, ax_from, ax_to)[source]
Transpose (permute) two axes of an ndarray.
- Return type:
Tensor
- static expand_dims(arr, axis)[source]
Inserts a new singleton axis at the indicated position, thus increasing ndim by 1.
- Return type:
Tensor
- abs(arr)[source]
Returns the absolute value (magnitude) of the elements in the input.
- Parameters:
arr – The input array.
- Return type:
Tensor
- Returns:
np.abs(arr) or equivalent
- conj(arr)[source]
Returns the conjugate of the elements in the input.
- Parameters:
arr – The input array.
- Return type:
Tensor
- Returns:
arr.conj or equivalent
- allclose(arr, other=0.0)[source]
Returns True if all elements in arr are close to other.
- Return type:
- ft(arr)[source]
Calculates the discrete Fourier transform over the spatial dimensions of E. The computational complexity is that of a Fast Fourier Transform:
O(N\log(N))
.
- ift(arr)[source]
Calculates the inverse Fourier transform over the spatial dimensions of E. The computational complexity is that of a Fast Fourier Transform:
O(N\log(N))
. The scaling is so thatE == self.ift(self.ft(E))
- convolve(operation_ft, arr)
Apply an operator in Fourier space. This is used to perform a cyclic FFT convolution which overwrites its input argument arr.
- Parameters:
- Return type:
Tensor
- Returns:
the convolved input array.
- subtract(left_term, right_term)[source]
Point-wise difference of A and B.
- Parameters:
left_term (
Union
[Complex
,Sequence
,ndarray
,Tensor
]) – The left matrix array, must start with dimensions n x mright_term (
Union
[Complex
,Sequence
,ndarray
,Tensor
]) – The right matrix array, must have matching or singleton dimensions to those of A. In case of missing dimensions, singletons are assumed.
- Return type:
- Returns:
The point-wise difference of both sets of matrices. Singleton dimensions are expanded.
- mul(left_factor, right_factor, out=None)
Point-wise matrix multiplication of A and B. Overwrites right_factor!
- Parameters:
left_factor (
Union
[Complex
,Sequence
,ndarray
,Tensor
]) – The left matrix array, must start with dimensions n x mright_factor (
Union
[Complex
,Sequence
,ndarray
,Tensor
]) – The right matrix array, must have matching or singleton dimensions to those of A, bar the first two dimensions. In case of missing dimensions, singletons are assumed. The first dimensions must be m x p. Where the m matches that of the left hand matrix unless both m and p are 1 or both n and m are 1, in which case the scaled identity is assumed.out (
Optional
[Tensor
]) – (optional) The destination array for the results.
- Return type:
Tensor
- Returns:
An array of matrix products with all but the first two dimensions broadcast as needed.
- ldivide(denominator, numerator=1.0)[source]
Parallel matrix left division, A^{-1}B, on the final two dimensions of A and B result_lm = A_kl B_km
A and B must have have all but the final dimension identical or singletons. B defaults to the identity matrix.
- longitudinal_projection_ft(field_array_ft)[source]
Projects the Fourier transform of a vector array onto its longitudinal component. Overwrites self.array_ft_input!
- transversal_projection_ft(field_array_ft)[source]
Projects the Fourier transform of a vector E array onto its transversal component.
- curl_ft(field_array_ft)[source]
Calculates the Fourier transform of the curl of a Fourier transformed E with the final dimension the vector dimension. The final dimension of the output will always be of length 3; however, the input length may be shorter, in which case the missing values are assumed to be zero. The first dimension of the input array corresponds to the first element in the final dimension, if it exists, the second dimension corresponds to the second element etc.
- calc_roots_of_low_order_polynomial(C)
Calculates the (complex) roots of polynomials up to order 3 in parallel. The coefficients of the polynomials are in the first dimension of C and repeated in the following dimensions, one for each polynomial to determine the roots of. The coefficients are input for low to high order in each column. In other words, the polynomial is:
np.sum(C * (x**range(C.size))) == 0
This method is used by mat3_eigh, but only for polynomials with real roots.
- static clear_cache()
- cross(A, B)
Calculates the cross product of vector arrays A and B.
- Parameters:
- Return type:
- Returns:
A vector array of dimensions
[vector_length, 1, *data_shape]
containing the cross product A x B in the first dimension and the other dimensions remain on the same axes.
- curl(field_array)
Calculates the curl of a vector E with the final dimension the vector dimension. The input argument may be overwritten!
- div_ft(field_array_ft)
Calculate the Fourier transform of the divergence of the pre-Fourier transformed input electric field.
- Parameters:
field_array_ft (
Union
[Complex
,Sequence
,ndarray
]) – The input array representing the field pre-Fourier-transformed in all spatial dimensions. The input is of the shape[m, n, x, y, z, ...]
- Return type:
- Returns:
The Fourier transform of the divergence of the field in the shape
[n, 1, x, y, z]
.
- static evaluate_polynomial(C, X)
Evaluates the polynomial P at X for testing.
- property eye: ndarray
Returns an identity tensor that can be multiplied using singleton expansion. This can be useful for scalar additions or subtractions.
- Returns:
an array with the number of dimensions matching that of the BackEnds’s data set.
- property ft_axes: tuple
The integer indices of the spatial dimensions, i.e. on which to Fourier Transform.
- property hardware_dtype
The scalar data type that is processed by this back-end.
- inv(mat)
Inverts the set of input matrices M.
- is_matrix(arr)
Checks if an ndarray is a matrix as defined by this parallel_ops_column object.
- is_vector(arr)
Tests if A represents a vector field.
- property k2: ndarray
Helper def for calculation of the Fourier transform of the Green def
- Returns:
\(|k|^2\) for the specified sample grid and output shape
- longitudinal_projection(field_array)
Projects vector arrays onto their longitudinal component. The input argument may be overwritten!
- mat3_eigh(arr)[source]
Calculates the eigenvalues of the 3x3 Hermitian matrices represented by A and returns a new array of 3-vectors, one for each matrix in A and of the same dimensions, baring the second dimension. When the first two dimensions are 3x1 or 1x3, a diagonal matrix is assumed. When the first two dimensions are singletons (1x1), a constant diagonal matrix is assumed and only one eigenvalue is returned. Returns an array of one dimension less: 3 x data_shape. With the exception of the first dimension, the shape is maintained.
- Parameters:
arr (
Union
[Complex
,Sequence
,ndarray
,Tensor
]) – The set of 3x3 input matrices for which the eigenvalues are requested. This must be an ndarray with the first two dimensions of size 3.- Return type:
Tensor
- Returns:
The set of eigenvalue-triples contained in an ndarray with its first dimension of size 3, and the remaining dimensions equal to all but the first two input dimensions.
- outer(A, B)
Calculates the Dyadic product of vector arrays A and B.
- Parameters:
- Return type:
- Returns:
A matrix array of dimensions
[vector_length, vector_length, *data_shape]
containing the dyadic product \(A \otimes B\) in the first two dimensions and the other dimensions remain on the same axes.
- to_matrix_field(arr)
Converts the input to an array of the full number of dimensions: len(self.matrix_shape) + len(self.grid.shape), and dtype. The size of each dimensions must match that of that set for the
BackEnd
or be 1 (assumes broadcasting). For electric fields in 3-space, self.matrix_shape == (N, N) == (3, 3) The first (left-most) dimensions of the output are either1x1: The identity matrix for a scalar field, as sound waves or isotropic permittivity.
Nx1: A vector for a vector field, as the electric field.
NxN: A matrix for a matrix field, as anisotropic permittivity
None is interpreted as 0. Singleton dimensions are added on the left so that all dimensions are present. Inputs with 1xN are transposed (not conjugate) to Nx1 vectors.
- Parameters:
arr (
Union
[Complex
,Sequence
,ndarray
]) – The input can be scalar, which assumes that its value is assumed to be repeated for all space. The value can be a one-dimensional vector, in which case the vector is assumed to be repeated for all space.- Return type:
- Returns:
An array with ndim == len(self.matrix_shape) + len(self.grid.shape) and with each non-singleton dimension matching those of the nb_rows and data_shape.
- transversal_projection(field_array)
Projects vector arrays onto their transverse component. The input argument may be overwritten!
macromax.backend.torch module
The module providing the PyTorch back-end implementation.
- class macromax.backend.torch.BackEndTorch(nb_dims, grid, hardware_dtype=torch.complex128, device=None)[source]
Bases:
BackEnd
A class that provides methods to work with arrays of matrices or block-diagonal matrices, represented as ndarrays, where the first two dimensions are those of the matrix, and the final dimensions are the coordinates over which the operations are parallelized and the Fourier transforms are applied.
- __init__(nb_dims, grid, hardware_dtype=torch.complex128, device=None)[source]
Construct object to handle parallel operations on square matrices of nb_rows x nb_rows elements. The matrices refer to points in space on a uniform plaid grid.
- Parameters:
nb_dims (
int
) – The number of rows and columns in each matrix. 1 for scalar operations, 3 for polarizationgrid (
Grid
) – The grid that defines the position of the matrices.hardware_dtype – (optional) The data type to use for operations.
device (
Optional
[str
]) – (optional) ‘cuda’ or ‘cpu’, to indicate where the calculation will happen.
- property numpy_dtype
The equivalent hardware data type in numpy
- astype(arr, dtype=None)[source]
As necessary, convert the ndarray arr to the type dtype.
- Return type:
Tensor
- assign(arr, out)[source]
Assign the values of one array to those of another. The dtype and shape is broadcasted to match that of the output array.
- Parameters:
arr – The values that need to be assigned to another array.
out – (optional) The target array, which must be of the same shape.
- Return type:
Tensor
- Returns:
The target array or a newly allocated array with the same values as those in arr.
- assign_exact(arr, out)[source]
Assign the values of one array to those of another. The dtype and shape must match that of the output array.
- Parameters:
arr – The values that need to be assigned to another array.
out – (optional) The target array, which must be of the same shape.
- Return type:
Tensor
- Returns:
The target array or a newly allocated array with the same values as those in arr.
- allocate_array(shape=None, dtype=None, fill_value=None)[source]
Allocates a new vector array of shape grid.shape and word-aligned for efficient calculations.
- Return type:
Tensor
- sign(arr)[source]
Returns an array with values of -1 where arr is negative, 0 where arr is 0, and 1 where arr is positive.
- swapaxes(arr, ax_from, ax_to)[source]
Transpose (permute) two axes of an ndarray.
- Return type:
Tensor
- static expand_dims(arr, axis)[source]
Inserts a new singleton axis at the indicated position, thus increasing ndim by 1.
- Return type:
Tensor
- abs(arr)[source]
Returns the absolute value (magnitude) of the elements in the input.
- Parameters:
arr – The input array.
- Return type:
Tensor
- Returns:
np.abs(arr) or equivalent
- conj(arr)[source]
Returns the conjugate of the elements in the input.
- Parameters:
arr – The input array.
- Return type:
Tensor
- Returns:
arr.conj or equivalent
- allclose(arr, other=0.0)[source]
Returns True if all elements in arr are close to other.
- Return type:
- ft(arr)[source]
Calculates the discrete Fourier transform over the spatial dimensions of E. The computational complexity is that of a Fast Fourier Transform:
O(N\log(N))
.
- ift(arr)[source]
Calculates the inverse Fourier transform over the spatial dimensions of E. The computational complexity is that of a Fast Fourier Transform:
O(N\log(N))
. The scaling is so thatE == self.ift(self.ft(E))
- mul(left_factor, right_factor, out=None)[source]
Point-wise matrix multiplication of A and B. Overwrites right_factor!
- Parameters:
left_factor (
Union
[Complex
,Sequence
,ndarray
,Tensor
]) – The left matrix array, must start with dimensions n x mright_factor (
Union
[Complex
,Sequence
,ndarray
,Tensor
]) – The right matrix array, must have matching or singleton dimensions to those of A, bar the first two dimensions. In case of missing dimensions, singletons are assumed. The first dimensions must be m x p. Where the m matches that of the left hand matrix unless both m and p are 1 or both n and m are 1, in which case the scaled identity is assumed.out (
Optional
[Tensor
]) – (optional) The destination array for the results.
- Return type:
Tensor
- Returns:
An array of matrix products with all but the first two dimensions broadcast as needed.
- ldivide(denominator, numerator=1.0)[source]
Parallel matrix left division, A^{-1}B, on the final two dimensions of A and B result_lm = A_kl B_km
A and B must have have all but the final dimension identical or singletons. B defaults to the identity matrix.
- calc_roots_of_low_order_polynomial(C)
Calculates the (complex) roots of polynomials up to order 3 in parallel. The coefficients of the polynomials are in the first dimension of C and repeated in the following dimensions, one for each polynomial to determine the roots of. The coefficients are input for low to high order in each column. In other words, the polynomial is:
np.sum(C * (x**range(C.size))) == 0
This method is used by mat3_eigh, but only for polynomials with real roots.
- convolve(operation_ft, arr)
Apply an operator in Fourier space. This is used to perform a cyclic FFT convolution which overwrites its input argument arr.
- Parameters:
- Return type:
- Returns:
the convolved input array.
- cross(A, B)
Calculates the cross product of vector arrays A and B.
- Parameters:
- Return type:
- Returns:
A vector array of dimensions
[vector_length, 1, *data_shape]
containing the cross product A x B in the first dimension and the other dimensions remain on the same axes.
- curl(field_array)
Calculates the curl of a vector E with the final dimension the vector dimension. The input argument may be overwritten!
- curl_ft(field_array_ft)
Calculates the Fourier transform of the curl of a Fourier transformed E with the final dimension the vector dimension. The final dimension of the output will always be of length 3; however, the input length may be shorter, in which case the missing values are assumed to be zero. The first dimension of the input array corresponds to the first element in the final dimension, if it exists, the second dimension corresponds to the second element etc.
- div(field_array)
Calculate the divergence of the input vector or tensor field E. The input argument may be overwritten!
- div_ft(field_array_ft)
Calculate the Fourier transform of the divergence of the pre-Fourier transformed input electric field.
- Parameters:
field_array_ft (
Union
[Complex
,Sequence
,ndarray
]) – The input array representing the field pre-Fourier-transformed in all spatial dimensions. The input is of the shape[m, n, x, y, z, ...]
- Return type:
- Returns:
The Fourier transform of the divergence of the field in the shape
[n, 1, x, y, z]
.
- static evaluate_polynomial(C, X)
Evaluates the polynomial P at X for testing.
- property eye: ndarray
Returns an identity tensor that can be multiplied using singleton expansion. This can be useful for scalar additions or subtractions.
- Returns:
an array with the number of dimensions matching that of the BackEnds’s data set.
- property ft_axes: tuple
The integer indices of the spatial dimensions, i.e. on which to Fourier Transform.
- property hardware_dtype
The scalar data type that is processed by this back-end.
- inv(mat)
Inverts the set of input matrices M.
- is_matrix(arr)
Checks if an ndarray is a matrix as defined by this parallel_ops_column object.
- is_scalar(arr)
Tests if A represents a scalar field (as opposed to a vector field).
- is_vector(arr)
Tests if A represents a vector field.
- property k2: ndarray
Helper def for calculation of the Fourier transform of the Green def
- Returns:
\(|k|^2\) for the specified sample grid and output shape
- longitudinal_projection(field_array)
Projects vector arrays onto their longitudinal component. The input argument may be overwritten!
- longitudinal_projection_ft(field_array_ft)
Projects the Fourier transform of a vector array onto its longitudinal component. Overwrites self.array_ft_input!
- mat3_eigh(arr)
Calculates the eigenvalues of the 3x3 Hermitian matrices represented by A and returns a new array of 3-vectors, one for each matrix in A and of the same dimensions, baring the second dimension. When the first two dimensions are 3x1 or 1x3, a diagonal matrix is assumed. When the first two dimensions are singletons (1x1), a constant diagonal matrix is assumed and only one eigenvalue is returned. Returns an array of one dimension less: 3 x data_shape. With the exception of the first dimension, the shape is maintained.
Before substituting this for
`numpy.linalg.eigvalsh`
, note that this implementation is about twice as fast as the current numpy implementation for 3x3 Hermitian matrix fields. The difference is even greater for the PyTorch implementation. The implementation below can be made more efficient by limiting it to Hermitian matrices and thus real eigenvalues. In the end, only the maximal eigenvalue is required, so an iterative power iteration may be even faster, perhaps after applying a Given’s rotation or Householder reflection to make Hermitian tridiagonal.- Parameters:
arr (
Union
[Complex
,Sequence
,ndarray
]) – The set of 3x3 input matrices for which the eigenvalues are requested. This must be an ndarray with the first two dimensions of size 3.- Return type:
- Returns:
The set of eigenvalue-triples contained in an ndarray with its first dimension of size 3, and the remaining dimensions equal to all but the first two input dimensions.
- outer(A, B)
Calculates the Dyadic product of vector arrays A and B.
- Parameters:
- Return type:
- Returns:
A matrix array of dimensions
[vector_length, vector_length, *data_shape]
containing the dyadic product \(A \otimes B\) in the first two dimensions and the other dimensions remain on the same axes.
- subtract(left_term, right_term)
Point-wise difference of A and B.
- Parameters:
- Return type:
- Returns:
The point-wise difference of both sets of matrices. Singleton dimensions are expanded.
- to_matrix_field(arr)
Converts the input to an array of the full number of dimensions: len(self.matrix_shape) + len(self.grid.shape), and dtype. The size of each dimensions must match that of that set for the
BackEnd
or be 1 (assumes broadcasting). For electric fields in 3-space, self.matrix_shape == (N, N) == (3, 3) The first (left-most) dimensions of the output are either1x1: The identity matrix for a scalar field, as sound waves or isotropic permittivity.
Nx1: A vector for a vector field, as the electric field.
NxN: A matrix for a matrix field, as anisotropic permittivity
None is interpreted as 0. Singleton dimensions are added on the left so that all dimensions are present. Inputs with 1xN are transposed (not conjugate) to Nx1 vectors.
- Parameters:
arr (
Union
[Complex
,Sequence
,ndarray
]) – The input can be scalar, which assumes that its value is assumed to be repeated for all space. The value can be a one-dimensional vector, in which case the vector is assumed to be repeated for all space.- Return type:
- Returns:
An array with ndim == len(self.matrix_shape) + len(self.grid.shape) and with each non-singleton dimension matching those of the nb_rows and data_shape.
- transversal_projection(field_array)
Projects vector arrays onto their transverse component. The input argument may be overwritten!
- transversal_projection_ft(field_array_ft)
Projects the Fourier transform of a vector E array onto its transversal component.