Utility

This section contains the API documentation for the utility module.

Options and verbosity

class Options

Various package settings.

static CacheEnabled(bool)

Enable cache.

Usage

bool = Options.CacheEnabled(bool) sets the cache enabling to bool.

bool = Options.CacheEnabled(bool) returns the current cache enabling.

static Debug(bool)

Enable cache.

Usage

bool = Options.Debug(bool) sets the debug mode to bool.

bool = Options.Debug(bool) returns the current debug mode.

class Verbosity

Bases: uint8

Verbosity level enumeration class.

Levels:

  • off (0): No information

  • warn (1): Information on failure

  • conv (2): Convergence information

  • iter (3): Information about each iteration

  • detail (4): Detailed information about each iteration

  • diagnostics (255): All possible information

Toolbox

between(x1, x, x2)

Restrict x to lie between x1 and x2

colors(i)

Periodically cycle through standard Matlab colors.

dim2str(sz)

Convert size array to string output.

diracdelta(sz)

Construct delta tensor of given size. Note that all dimensions must have the same length.

memsize(in, unit)

Check memory size of object in prefered unit (‘GB’, ‘MB’, ‘KB’ or ‘B’).

Usage

[bytes, unit] = memsize(in, unit)

randc(varargin)

Generate an array containing pseudorandom complex values, both real and imaginary part drawn from the standard uniform distribution on the open interval (0, 1).

Usage

R = randc([m n])

R = randc(m, n, p, ...)

R = randc(..., classname)

R = randc(..., 'like', Y)

Arguments

  • m, n, p, … (int) – integers defining the size of the output array.

  • classname (char) – datatype of the array, default 'double'.

  • Y (numeric) – create an array of the same class as Y.

Returns

R (numeric) – complex pseudorandom values, real and imaginary part distributed from the uniform distribution.

randnc(varargin)

Generate an array containing pseudorandom complex values, both real and imaginary part drawn from the standard normal distribution.

Usage

R = randnc([m n])

R = randnc(m, n, p, ...)

R = randnc(..., classname)

R = randnc(..., 'like', Y)

Arguments

  • m, n, p, … (int) – integers defining the size of the output array.

  • classname (char) – datatype of the array, default 'double'.

  • Y (numeric) – create an array of the same class as Y.

Returns

R (numeric) – complex pseudorandom values, real and imaginary part distributed from the normal distribution.

simulsort(arrays, kwargs)

Simultaneous sorting of several input arrays.

Sorts the elements such that equal elements of array{i} appear sorted by array{i+1}.

This is achieved by sorting from end to front, making use of the fact that SORT is stable and thus will not mess up the order of later elements when earlier elements are equal.

Usage

[I, array1, array2, ...] = simulsort(array1, array2, ..., kwargs)

Arguments

array1, array2, … – arrays of equal size that need to be sorted. These can be of any type that supports SORT.

Keyword Arguments

  • Dimension (int) – determine the dimension along which to sort. This behaves similarly to SORT, by default: - for vectors, sorts the elements - for matrices, sorts each column - for N-D arrays, sorts along the first non-singleton dimension.

  • Direction (char, ‘ascend’ or ‘descend’) – specify the sorting direction, defaults to ‘ascend’.

Returns

  • I (int) – permutation vector that brings the input arrays into sorted order.

  • array1, array2, … – sorted arrays.

simulsortrows(arrays, kwargs)

Simultaneous sorting of several input arrays by row.

Sorts the rows such that equal rows of array{i} appear sorted by the rows of array{i+1}.

This is achieved by sorting rows from end to front, making use of the fact that SORTROWS is stable and thus will not mess up the order of later rows when earlier rows are equal.

Usage

[I, array1, array2, ...] = simulsortrows(array1, array2, ..., kwargs)

Repeating arguments

array1, array2, … – arrays of equal size that need to be sorted. These can be of any type that supports SORTROWS.

Keyword Arguments

  • Col (int) – vector of indices that specifies the columns used for sorting.

  • Direction (char, ‘ascend’ or ‘descend’) – specify the sorting direction. You can also specify a different direction for each column by using a cell array of ‘ascend’ and ‘descend’ of the same size as Col, such that corresponding elements are sorted ascending or descending.

Returns

  • I (int) – permutation vector that brings the input arrays into rowsorted order.

  • array1, array2, … – rowsorted arrays

simulunique(varargin)

Set unique over multiple arrays.

swapvars(A, B)

Exchange two variables.

Usage

[a, b] = swapvars(a, b) stores the value of a in b, and the value of b in a.

time2str(time, unit, precision)

Returns string version of a duration in seconds.

Usage

timeStr = time2str(time, unit, precision)

Indices

contractinds(ia, ib)

Find the contracted dimensions.

[dimA, dimB] = contractinds(ia, ib) locates the repeated indices in two vectors of integers.

mod1(x, y)

Modulus after division, counting from 1:y.

Arguments

  • x (int) – numerator.

  • y (int) – divisor.

Returns

m (int) – remainder after division, where a value of 0 is replaced with y.

next(i, total)

Give the next index in a cyclic loop.

Usage

j = next(i, total) gives the \(i + 1\)’th index, but loops back to 1 when \(j > \text{total}\).

See also

prev()

prev(i, total)

Give the previous index in a cyclic loop.

Usage

j = prev(i, total) gives the \(i - 1\)’th index, but loops back to total when \(j < 1\).

See also

next()

rankrange(rank)

Convert tensor rank into a contiguous range of dimensions.

traceinds(ia)

Find the traced dimensions.

[dimA1, dimA2] = traceinds(ia) locates the repeated indices in a vectors of integers.

treeindsequence(n)

Compute the number of edge labels needed in a splitting or fusing tree.

Usage

t = treeindsequence(n)

Arguments

n (int) – number of external edges

Returns

t (int) – total number of edges

Example

>> treeindsequence(0:4)

ans =

     0     1     2     4     6
unique1(inds)

Find the elements that appear exactly once.

Usage

inds = unique1(inds) deletes all elements that appear more than once.

Linear algebra

contract(tensors, indices, kwargs)

Compute a tensor network contraction.

Usage

C = contract(t1, idx1, t2, idx2, ...)

C = contract(..., 'Conj', conjlist)

C = contract(..., 'Rank', r)

Repeating Arguments

  • tensors (Tensor) – list of tensors that constitute the vertices of the network.

  • indices ((1, :) int) – list of indices that define the links and contraction order, using ncon-like syntax.

Keyword Arguments

  • Conj ((1, :) logical) – optional list to flag that tensors should be conjugated.

  • Rank ((1, 2) int) – optionally specify the rank of the resulting tensor.

Returns

C (Tensor or numeric) – result of the tensor network contraction.

distance(A, B)

Compute the Euclidian distance between two arrays.

eigsolve(A, v, howmany, sigma, options)

Find a few eigenvalues and eigenvectors of an operator.

Usage

[V, D, flag] = eigsolve(A, v, howmany, sigma, kwargs)

D = eigsolve(A, v, ...)

Arguments

  • A (matrix or function_handle) – A square matrix. A function handle which implements one of the following, depending on sigma:

    • A x, if sigma is 0 or ‘smallestabs’

    • (A - sigma * I) x, if sigma is a nonzero scalar

    • A * x, for all other cases

  • v (vector) – initial guess for the eigenvector. If A is a Tensor, this defaults to a random complex Tensor, for function handles this is a required argument.

  • howmany (int) – amount of eigenvalues and eigenvectors that should be computed. By default this is 1, and this should not be larger than the total dimension of A.

  • sigma (char or numeric) – selector for the eigenvalues, should be either one of the following:

    • ‘largestabs’, ‘lm’: default, eigenvalues of largest magnitude

    • ‘largestreal’, ‘lr’: eigenvalues with largest real part

    • ‘largestimag’, ‘li’: eigenvalues with largest imaginary part.

    • ‘smallestabs’, ‘sm’: default, eigenvalues of smallest magnitude

    • ‘smallestreal’, ‘sr’: eigenvalues with smallest real part

    • ‘smallestimag’, ‘si’: eigenvalues with smallest imaginary part.

    • numeric : eigenvalues closest to sigma.

Keyword Arguments

  • Tol (double) – tolerance of the algorithm.

  • Algorithm (char) – choice of eigensolver algorithm. Currently there is a choice between the use of Matlab’s buitin eigs specified by the identifiers ‘eigs’ or ‘KrylovSchur’, or the use of a custom Arnolid algorithm specified by the identifier ‘Arnoldi’.

  • MaxIter (int) – maximum number of iterations, 100 by default.

  • KrylovDim (int) – number of vectors kept in the Krylov subspace.

  • IsSymmetric (logical) – flag to speed up the algorithm if the operator is symmetric, false by default.

  • Verbosity (Verbosity) – Level of output information, by default nothing is printed if flag is returned, otherwise only warnings are given, defaults to Verbosity.warn.

Returns

  • V ((1, howmany) vector) – vector of eigenvectors.

  • D (numeric) – vector of eigenvalues if only a single output argument is asked, diagonal matrix of eigenvalues otherwise.

  • flag (int) – if flag = 0 then all eigenvalues are converged, otherwise not.

isapprox(A, B, tol)

Verify whether two arrays are approximately equal, based on their Euclidean distance.

Arguments

A, B (numeric) – input arrays of the same size.

Keyword Arguments

  • AbsTol (double) – absolute tolerance.

  • RelTol (double) – relative tolerance.

Returns

bool (logical) – true if A and B are approximately equal.

isisometry(A, side, tol)

Verify whether a matrix is an isometry.

Arguments

  • A (numeric) – input matrix.

  • side (char, ‘left’, ‘right’ or ‘both’) – check if A' * A == I, A * A' == I, or both by default.

Keyword Arguments

  • AbsTol (double) – absolute tolerance

  • RelTol (double) – relative tolerance

Returns

bool (logical) – true if A is an isometry.

isposdef(A)

Verify whether a matrix is positive definite.

leftnull(A, alg, atol)

Compute the left nullspace of a matrix, such that N' * A == 0.

Arguments

  • A (numeric) – input matrix.

  • alg (char, ‘qr’ or ‘svd’) – choice of algorithm, default ‘svd’.

  • atol (double) – absolute tolerance for the null space, defaults to max(size(A)) * eps(max(svd(A))).

Returns

N (numeric) – orthonormal basis for the orthogonal complement of the support of the row space of A.

leftorth(A, alg)

Factorize a matrix into an orthonormal basis Q and remainder R, such that A = Q * R.

Usage

[Q, R] = leftorth(A, alg)

Arguments

  • A (numeric) – input matrix to factorize.

  • alg (char or string) – selection of algorithms for the decomposition:

    • ‘qr’ produces an upper triangular remainder R

    • ‘qrpos’ corrects the diagonal elements of R to be positive.

    • ‘ql’ produces a lower triangular remainder R

    • ‘qlpos’ corrects the diagonal elements of R to be positive.

    • ‘polar’ produces a Hermitian and positive semidefinite R.

    • ‘svd’ uses a singular value decomposition.

Returns

  • Q (numeric) – orthonormal basis matrix

  • R (numeric) – remainder matrix, depends on selected algorithm.

qrpos(A, varargin)

Positive orthogonal-triangular decomposition.

rightnull(A, alg, atol)

Compute the right nullspace of a matrix, such that A * N' == 0.

Arguments

  • A (numeric) – input matrix.

  • alg (char, ‘lq’ or ‘svd’) – choice of algorithm, default ‘svd’.

  • atol (double) – absolute tolerance for the null space, defaults to max(size(A)) * eps(max(svd(A))).

Returns

N (numeric) – orthonormal basis for the orthogonal complement of the support of the column space of A.

rightorth(A, alg)

Factorize a matrix into an orthonormal basis Q and remainder L, such that A = L * Q.

Usage

[R, Q] = rightorth(A, alg)

Arguments

  • A (numeric) – input matrix to factorize.

  • alg (char or string) – selection of algorithms for the decomposition:

    • ‘rq’ produces an upper triangular remainder R

    • ‘rqpos’ corrects the diagonal elements of R to be positive.

    • ‘lq’ produces a lower triangular remainder R

    • ‘lqpos’ corrects the diagonal elements of R to be positive.

    • ‘polar’ produces a Hermitian and positive semidefinite R.

    • ‘svd’ uses a singular value decomposition.

Returns

  • R (numeric) – remainder matrix, depends on selected algorithm.

  • Q (numeric) – orthonormal basis matrix.

tensorprod(A, B, dimA, dimB, cA, cB, options)

Tensor products between two tensors.

Usage

C = tensorprod(A, B, dimA, dimB)

returns the tensor product of tensors A and B. The arguments dimA and dimB are vectors that specify which dimensions to contract in A and B. The size of the output tensor is the size of the uncontracted dimensions of A followed by the size of the uncontracted dimensions of B.

C = tensorprod(A, B)

returns the outer product of tensors A and B. This is equivalent to the previous syntax with dimA = dimB = [].

C = tensorprod(_, NumDimensionsA=ndimsA)

optionally specifies the number of dimensions in tensor A in addition to combat the removal of trailing singleton dimensions.

tensortrace(A, i1, i2)

Compute the (partial) trace of a tensor.

Usage

[C, ic] = tensortrace(A, i1)

traces over the indices that appear twice in i1.

[C, ic] = tensortrace(A, i1, i2)

optionally specifies the output indices’ order.

Permutations

invperm(p)

Compute the inverse of an integer permutation.

Usage

invp = invperm(p)

computes the permutation that satisfies invp(p) = 1:length(p).

isperm(p)

Check if a vector is a permutation.

Usage

bool = isperm(p)

perm2swap(p)

Convert a permutation into a sequence of swaps on neighbouring indices.

Arguments

p (int) – permutation vector.

Returns

s (int) – list of swaps that compose into the permutation vector, where i indicates a swap between indices i and i+1.

Sparse Arrays

class SparseArray(varargin)

Class for multi-dimensional sparse arrays.

Limited to arrays with a total number of elements of at most 2^48-1.

SparseArray(varargin)

Create a sparse array.

Usage

a = SparseArray(subs, vals, sz)

uses the rows of subs and vals to generate a sparse array a of size sz = [m1 m2 ... mn]. subs is a p x n array specifying the subscripts of the nonzero values to be inserted into a. The k-th row of subs specifies the subscripts for the k-th value in vals. The argument vals may be scalar, in which case it is expanded to be the same length as subs, i.e., it is equivalent to vals * (p, 1). In the case of duplicate subscripts in subs, the corresponding values are added.

a = SparseArray

Empty constructor.

a = SparseArray(b)

Copies/converts b if it is a SparseArray, a dense array or a sparse matrix.

a = SparseArray(b, sz)

Copies/converts b if it is a SparseArray, a dense array or a sparse matrix, and sets the size of a to sz

Example

>> subs = [1 1 1; 1 1 3; 2 2 2; 4 4 4; 1 1 1; 1 1 1];
>> vals = [0.5; 1.5; 2.5; 3.5; 4.5; 5.5];
>> sz = [4 4 4];
>> a = SparseArray(subs, vals, sz) %<-- sparse 4x4x4
abs(a)

Absolute value.

Arguments

a (SparseArray) – input array.

Returns

b (SparseArray) – output array.

chop(a, tol)

Set all nonzero values in SparseArray who’s absolute value is below a given threshold to zero.

Arguments

  • a (SparseArray) – input array.

  • tol (double , optional) – threshold tolerance for absolute values of entries, defaults to 1e-15.

Returns

b (SparseArray) – sparse array with entries of absolute value below tol set to zero.

conj(a)

Complex conjugate.

Arguments

a (SparseArray) – input array.

Returns

b (SparseArray) – output array.

ctranspose(a)

Complex conjugate transpose.

Only defined for 2D sparse arrays.

Usage

b = ctranspose(a)

b = a'

Arguments

a (SparseArray) – input array.

Returns

b (SparseArray) – output array.

eig(a, varargin)

Find eigenvalues and eigenvectors of a 2D sparse array.

See also

Documentation for builtin Matlab eig.

eigs(a, varargin)

Find a few eigenvalues and eigenvectors of a 2D sparse array.

See also

Documentation for builtin Matlab eigs.

find(a)

Find subscripts of nonzero elements in a sparse array.

Usage

idx = find(a)

[subs, idx, vals] = find(a)

Arguments

a (SparseArray) – input array

Returns

  • subs ((:, :) int) – subscripts of nonzero array entries.

  • idx ((:, 1) int) – linear indices of nonzero array entries.

  • var ((:, 1) double) – values of nonzero array entries.

full(a)

Convert a sparse array to a dense array.

Arguments

a (SparseArray) – input array.

Returns

d (double) – dense output array.

groupind(a, g)

Group (block) specified sets of contiguous indices.

Arguments

  • a (SparseArray) – input array.

  • g ((1, :) int) – list of number of contiguous indices to be grouped in each index of the output tensor.

Returns

b (SparseArray) – output array with grouped indices.

Example

>> a = SparseArray.random([2, 3, 4, 5, 6], .1);
>> b = groupind(a, [3, 2]); %<-- sparse 24x30
imag(a)

Complex imaginary part of sparse array.

Arguments

a (SparseArray) – input array.

Returns

b (SparseArray) – output array with real entries corresponding to the imaginary part of the entries of a.

isnumeric(~)

Determine whether input is numeric.

Arguments

a (SparseArray) – input array.

Returns

bool (logical) – defaults to true for SparseArray.

isscalar(~)

Determine whether input is scalar.

Arguments

a (SparseArray) – input array.

Returns

bool (logical) – defaults to false for SparseArray.

isrow(a)

Determine whether input is row vector.

Arguments

a (SparseArray) – input array.

Returns

bool (logical)

ldivide(a, b)

Elementwise left division for sparse arrays.

Usage

ldivide(a, b)

a .b

where a or b is a SparseArray. a`and :code:`b must have the same size, unless one is a scalar.

Arguments

a, b (SparseArray or double) – input arrays to be divided.

Returns

c (SparseArray) – elementwise left division of b by a.

minus(a, b)

Elementwise subtraction for sparse arrays.

Usage

minus(a, b)

a - b

where a`or :code:`b is a SparseArray. a`and :code:`b must have the same size, unless one of the is scalar. Scalar can be subtracted from a sparse array of any size, resulting in a dense array.

Arguments

a, b (SparseArray or double) – intput arrays.

Returns

c (SparseArray or double) – output array.

Example

>> a = SparseArray.random([4 3 2], .1);
>> b = SparseArray.random([4 3 2], .1);
>> a - b %<-- sparse
>> a - 5 %<-- dense
>> a - 0 %<-- dense
>> a - full(a) %<-- dense
mrdivide(a, b)

Matrix right division for sparse arrays.

Usage

mrdivide(a, b)

a / b

Arguments

  • a (SparseArray) – intput array.

  • b (double) – scalar to divide by.

Returns

c (SparseArray) – output array.

Example

>> a = SparseArray.random([4 3 2], .1);
>> a / 3 %<-- sparse
mtimes(a, b)

Matrix multiplication for 2D sparse arrays.

Usage

mtimes(a, b)

a * b

where a or b is a SparseArray.

See also

mtimes

ndims(t)

Number of dimensions of a sparse array.

Note that unlike the builtin Matlab behavior <https://mathworks.com/help/matlab/ref/double.ndims.html> trailing singleton dimensions are not ignored.

Example

>> a = SparseArray.random([4 3 1], .1);
>> ndims(a) %<-- returns 3
nnz(a)

Number of nonzero elements in sparse array.

nonzeros(a)

Returns a full column vector of the nonzero elements of a.

norm(a)

Frobenius norm of a sparse array.

numel(a)

Number of elements in a sparse array.

outer(a, b)

Outer product of two sparse arrays.

Warning

Not kron.

permute(a, order)

Permute sparse array dimensions.

plus(a, b)

Elementwise addition for sparse arrays.

Usage

plus(a, b)

a + b

where a and b must have the same size, unless one is a scalar. A scalar can be added to a sparse array of any size.

Arguments

a, b (SparseArray or double) – input arrays.

Returns

c (SparseArray) – output array.

Example

>> a = SparseArray.random([4 3 2], .1);
>> a + b %<-- sparse
>> a + 5 %<-- dense
>> a + 0 %<-- dense
>> a + full(a) %<-- dense
power(a, b)

Elementwise power for sparse array.

qr(a, varargin)

Orthogonal-triangular decomposition for a sparse array.

See also

Documentation for builtin Matlab qr.

qrpos(a, varargin)

Positive orthogonal-triangular decomposition for a sparse array.

rdivide(a, b)

Elementwise right division for sparse arrays.

Usage

rdivide(a, b)

a ./ b

where a or b is a SparseArray. a`and :code:`b must have the same size, unless one is a scalar.

Arguments

a, b (SparseArray or double) – input arrays to be divided.

Returns

c (SparseArray) – elementwise left division of a by b.

Example

>> a = SparseArray.random([4 3 2], .1);
>> a ./ 5 %<-- sparse
>> 5 ./ a %<-- dense
>> a ./ full(a) %<-- sparse
>> full(a) ./ a %<-- dense
real(a)

Complex real part of sparse array.

Arguments

a (SparseArray) – input array.

Returns

b (SparseArray) – output array with real entries corresponding to the real part of the entries of a.

reshape(a, varargin)

Reshape sparse array.

sign(a)

Signum function.

sparse(a)

Convert 2D sparse array to a sparse matrix.

spmatrix(a)

Convert 2D sparse array to a sparse matrix.

See also

sparse

squeeze(a)

Remove singleton dimensions from a sparse array.

Usage

b = squeeze(a)

returns a sparse array b with the same elements as a but with all the singleton dimensions removed.

Example

>> squeeze(SparseArray.random([2, 1, 3], 0.5)) %<-- returns a 2 x 3 SparseArray
>> squeeze(SparseArray([1, 1, 1], 1, [1, 1, 1])) %<-- returns a scalar
subsasgn(a, s, rhs)

Subscripted assignment for sparse array.

subsref(a, s)

Subscripted reference for a sparse array.

Usage

a_sub = a(i)

linear indexing for SparseArray with only one non-singleton dimension, returns a scalar.

a_sub = a(i1, i2, ..., iN)

where each in is an integer index, returns a scalar.

a_sub = a(R1, R2, ..., RN)

where each Rn is either a colon, an array of integers representing a slice in this dimension or a logical array representing the logical indexing of this dimension. Returns a sparse array.

a_sub = a(S)

where S is a 1 x n array of integer subscripts (for dimensions that are indexed) or zeroes (for dimensions that are not indexed). Returns a sparse array.

Example

>> a = SparseArray([4, 4, 4; 2, 2, 1; 2, 3, 2], [3; 5; 1], [4, 4, 4]);
>> a(1, 2, 1) %<-- returns zero
>> a(4, 4, 4) %<-- returns 3
>> a(2, :, :) %<-- returns a 1 x 4 x 4 sparse array
>> a([0, 4, 0]) %<-- returns a 4 x 1 x 4 sparse array
>> a([4, 4, 0]) %<-- returns a 1 x 1 x 4 sparse array
sum(a)

Sum of all elements of a SparseArray.

svd(a, varargin)

Singular value decomposition.

svds(a, varargin)

Find a few singular values and vectors.

See also

Documentation for builtin Matlab svds.

times(a, b)

Array multiplication for sparse tensors.

Usage

c = times(a, b)

a .* b

where a or b is a sparse array. a and b must have the same size, unless one is a scalar. A scalar can be be multiplied by a sparse array of any size.

Example

>> a = SparseArray.random([4 3 2], .1);
>> b = SparseArray.random([4 3 2], .1);
>> a .* b %<-- sparse
>> a .* 5 %<-- sparse
>> a .* 0 %<-- sparse
>> a .* full(a) %<-- sparse
transpose(a)

Transpose.

Only defined for 2D sparse arrays.

Usage

b = transpose(a)

b = a.'

Arguments

a (SparseArray) – input array.

Returns

b (SparseArray) – output array.

uminus(a)

Unary minus.

Usage

b = uminus(a)

b = -a

Arguments

a (SparseArray) – input array.

Returns

b (SparseArray) – output array.

uplus(a)

Unary plus.

Usage

b = uplus(a)

b = +a

Arguments

a (SparseArray) – input array.

Returns

b (SparseArray) – output array.

static delta(numinds, inddim)

Create delta- (ghz-) array with given number of indices and index dimension.

Arguments

  • numinds (int) – number of indices of delta-array.

  • inddim (int) – dimension of each index of delta-array.

Returns

a (SparseArray) – output delta-array.

static random(sz, density)

Create a random complex sparse array.

Arguments

  • sz ((1, :) int) – size of the sparse array

  • density (double) – density of nonzero elements (0 < density < 1)

Validations

mustBeEqualLength(a, b)

Validate that the inputs are of equal length.

mustBeEqualLength(a, b) throws an error if length(a) ~= length(b)

Note

Supported by all classes that define these methods: length

mustBeSorted(A)

Validate that input is sorted.

mustBeSorted(A) throws an error if A is not sorted.

Note

Supported by all classes that define these methods: issorted

mustBeUnique(A)

Validate that there are no repeated values.

mustBeUnique(A) throws an error if A contains repeated values.

Note

Supported all classes that define these methods: unique, length

Wigner

Wigner3j(j1, j2, j3, m1, m2, m3, ifs, ifcb)

Computes the Wigner \(3j\) symbols

\[\begin{split}W = \begin{pmatrix} j_1 & j_2 & j_3 \\ m_1 & m_2 & m_3 \end{pmatrix}\end{split}\]

or the Clebsch-Gordan coefficient

\[W = \left\langle j_1,m_1,j_2,m_2 \middle| j_3,m_3 \right\rangle\]

using the Racah formula.

Usage

W = Wigner3j(j1, j2, j3, m1, m2, m3, ifs, ifcb)

Wigner3j(j1, j2, j3, m1, m2, m3)

Arguments

  • j1, j2, j3 – set of three angular momenta, must be arrays of the same size

  • m1, m2, m3 – set of three angular momentum projections, must be arrays of the same size as the angular momenta

Optional Arguments

  • ifs – the switch of the computational methods; although the central part of all the computations is based on the Racah formula, some details differ:

    • 0: the symbolic (presumably accurate) computation with the double-precision output

    • -1: the same symbolic computation as ifs=0 but with the symbolic output (accurate square root/rational-type equation, simplified)

    • -2: the same as ifs=-1 but without a final simplification of the symbolic expression, which can be time-consuming sometimes

    • 1: (default, recommended), 2, 3 - numeric double-precision algorithms (see Notes below)

    • 4: the symbolic computation with the double precision output, but cached.

    all other input values of ifs are set to the closest from the above list.

  • ifcb – if exists and is true, switches to computing the Clebsch-Gordan coefficients instead of the \(3j\)-symbols (default false)

Returns

W – the resulting values of the \(3j\)-symbols (ifcb=false) or the Clebsch-Gordan coefficients (ifcb=true) in either numeric double-precision or symbolic form (see the input parameter ifs); array of the same size as \(j_1\).

Notes

most of the estimates below is based on extensive numerical tests within a range of the quantum nubers <=1000

  1. the numeric algorithms (ifs>0) are usually much faster than the symbolic algorithm (ifs<=0)

  2. the accuracy of the symbolic algorithm remains the uttermost possible

  3. the accuracy of the numeric algorithms can worsen for big quantum numbers

  4. all the “numeric” algorithms switch automatically to the symbolic computations as soon as the numeric overflow occurs, thereby improving the accuracy but slowing down the computations for some big quantum numbers

  5. in cases with at least one of the \(j\) quantum numbers <=2, the explicit accurate equations are applied in all the algorithms ensuring the highest possible accuracy

  6. for relatively small quantum numbers (up to ~20) all the algorithms provide approximately equivalent results with a reasonably high accuracy

  7. the algorithm ifs=1 ensures a reasonably good accuracy; the worst registered case was Wigner3j(56,59,51,-2,14,-12,1)=-0.00019235064 (exact Wigner3j(56,59,51,-2,14,-12,0)=-0.00019235487) with the relative inaccuracy of 2.2e-5 (occuring before switching to the symbolic regime)

  8. the algorithm ifs=2 (comparing to ifs=1) proceeds in a wider range of quantum numbers numerically before switching to the symbolic regime (i.e., can work faster with big quantum numbers) but with somewhat lower precision; in our experimental probing we registered the cases with the worst relative inaccuracy of 0.13 for Wigner3j(80.5,85,68.5,-18.5,-15,33.5,2)=-2.781e-06 (exact Wigner3j(80.5,85,68.5,-18.5,-15,33.5,0)=-2.4637e-06) and abolute inaccuracy of 0.0001 for Wigner3j(82.5,67,90.5,4.5,1,-5.5,2)=0.0045102 (exact 0.0046133=0.0046133); in the range of the numeric calculations of the algorithm ifs=1, the algorithm ifs=2 provides the same characteristic accuracy as ifs=1;

  9. the algorithm ifs=3 (comparing to ifs=2) proceeds in a wider range of quantum numbers numerically before switching to the symbolic regime (i.e., can work faster with big quantum numbers) but for some combinations of big quantum numbers, completely wrong results are observed: e.g., Wigner3j(465.5,488.5,498,13.5,-90.5,77,3); Wigner3j(95,96,99,-17,1,16,3); thereby, we do not recommend using the algorithm ifs=3 at all, and only keep it here in view of the completeness.

References

  1. https://en.wikipedia.org/wiki/3-j_symbol

  2. https://en.wikipedia.org/wiki/Clebsch-Gordan_coefficients

  3. Zare, R. N., & Harter, W. G. (1988). Angular momentum: understanding spatial aspects in chemistry and physics. New York, 120.

Wigner6j(j1, j2, j3, j4, j5, j6, ifs, ifcb)

Computes the Wigner \(6j\) symbols

\[\begin{split}W = \begin{Bmatrix} j_1 & j_2 & j_3 \\ j_4 & j_5 & j_6 \end{Bmatrix}\end{split}\]

or the recoupling matrix element

\[W = \left\langle (j_1,j_2)j_3,j_4,j_5 \middle| j_1,(j_2,j_4)j_6,j_5 \right\rangle\]

using the Racah formula.

VERSION of January, 2020

Programmer: V. B. Sovkov, St. Petersburg State University, Shanxi University

Usage

W = Wigner6j(j1, j2, j3, j4, j5, j6, ifs, ifcb)

Wigner6j(j1, j2, j3, j4, j5, j6)

Arguments

j1, j2, j3, j4, j5, j6 – set of six angular momenta, must be arrays of the same size

Optional Arguments

  • ifs – the switch of the computational methods; although the central part of all the computations is based on the Racah formula, some details differ:

    • 0 - the symbolic (presumably accurate) computation with the double-precision output

    • -1 - the same symbolic computation as ifs=0 but with the symbolic output (accurate square root/rational-type equation, simplified)

    • -2 - the same as ifs=-1 but without a final simplification of the symbolic expression, which can be time-consuming sometimes

    • 1, 2 (default, recommended), 3 - numeric double-precision algorithms (see Notes below)

      all other input values of ifs are set to the closest from the above list.

  • ifcb – if exists and is true, switches to computing the coupling matrix elements instead of the \(6j\)-symbols (default false)

Returns

W – the resulting values of the \(6j\)-symbols (ifcb=false) or the coupling matrix elements coefficients (ifcb=true) in either numeric double-precision or symbolic form (see the input parameter ifs); array of the same size as \(j_1\).

Notes

most of the estimates below is based on extensive numerical tests within a range of the quantum nubers <=1000

  1. the numeric algorithms (ifs>0) are usually much faster than the symbolic algorithm (ifs<=0)

  2. the accuracy of the symbolic algorithm remains the uttermost possible

  3. the accuracy of the numeric algorithms can worsen for big quantum numbers

  4. all the “numeric” algorithms switch automatically to the symbolic computations

  5. as soon as the numeric overflow occurs, thereby improving the accuracy but slowing down the computations for some big quantum numbers

  6. in cases with at least one of the \(j\) quantum numbers <=2, the explicit accurate equations are applied in all the algorithms ensuring the highest possible accuracy

  7. for relatively small quantum numbers (up to ~20) all the algorithms provide approximately equivalent results with a reasonably high accuracy

  8. the algorithms ifs=1 and ifs=2 ensure a reasonably good accuracy; the worst registered ca for both algorithms (occuring before switching to the symbolic regime) was Wigner6j(41.5,52,43.5,36,38.5,46)=-0.00029 with the biggest relative inaccuracy of 3.6e-10

  9. the algorithm ifs=2 (comparing to ifs=1) proceeds in a wider range of quantum numbers numerically before switching to the symbolic regime (i.e., can work faster with big quantum numbers); in our experimental probing it proved to ensure a good enough accuracy and was chosen as a default one

  10. the algorithm ifs=3 (comparing to ifs=2) proceeds in a wider range of quantum numbers numerically before switching to the symbolic regime (i.e., can work faster with big quantum numbers) but for some combinations of big quantum numbers, completely wrong results are observed: e.g., Wigner6j(227,230.5,210.5,249.5,235,232,3); thereby, we do not recommend using the algorithm code:ifs=3 at all, and only keep it here in view of the completeness.

References

  1. https://en.wikipedia.org/wiki/6-j_symbol

  2. Zare, R. N., & Harter, W. G. (1988). Angular momentum: understanding spatial aspects in chemistry and physics. New York, 120.