From 5436730a2e75efbd97674c3153b5cc1fd2767fb3 Mon Sep 17 00:00:00 2001 From: lkdvos <37111893+lkdvos@users.noreply.github.com> Date: Sat, 1 Oct 2022 20:26:17 +0200 Subject: [PATCH] fermions, convenience and docs Various bugfixes, additional documentation, convenience functions and fermion functionality. --- docs/requirements.txt | 2 +- docs/src/conf.py | 4 +- docs/src/man/tensor.rst | 206 ++++++++++++++++++++- src/tensors/FusionTree.m | 55 +++++- src/tensors/Tensor.m | 263 ++++++++++++++++++++++----- src/tensors/charges/AbstractCharge.m | 49 +++++ src/tensors/kernels/AbelianBlock.m | 158 +++------------- src/tensors/kernels/AbstractBlock.m | 29 ++- src/tensors/kernels/MatrixBlock.m | 82 +++++++-- src/tensors/spaces/AbstractSpace.m | 20 +- src/tensors/spaces/CartesianSpace.m | 14 +- src/tensors/spaces/ComplexSpace.m | 2 +- src/tensors/spaces/GradedSpace.m | 103 +++++++++-- src/utility/Options.m | 11 ++ src/utility/linalg/contract.m | 27 +++ src/utility/memsize.m | 2 +- test/TestTensor.m | 28 ++- 17 files changed, 813 insertions(+), 242 deletions(-) diff --git a/docs/requirements.txt b/docs/requirements.txt index 20914e8..dbf71cb 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -6,4 +6,4 @@ nbsphinx==0.8.9 sphinx-gallery==0.10.1 myst-parser==0.17.2 linkify-it-py==2.0.0 - +jinja2==3.0.3 diff --git a/docs/src/conf.py b/docs/src/conf.py index 72377ee..c760c53 100644 --- a/docs/src/conf.py +++ b/docs/src/conf.py @@ -11,7 +11,7 @@ docs_src = os.path.abspath(os.path.dirname(__file__)) # docs source folder docs_root = os.path.abspath(os.path.join(docs_src, '..')) # docs root folder repo_root = os.path.abspath(os.path.join(docs_src, '..', '..')) # repo root folder -GITHUBBASE = 'https://github.com/lkdvos/TensorTrack' # repo link +GITHUBBASE = 'https://github.com/QuantumGhent/TensorTrack' # repo link # -- Project information @@ -67,7 +67,7 @@ # # This is also used if you do content translation via gettext catalogs. # Usually you set "language" from the command line for these cases. -language = None +language = 'en' # The name of the Pygments (syntax highlighting) style to use. pygments_style = 'default' diff --git a/docs/src/man/tensor.rst b/docs/src/man/tensor.rst index de533d0..e5bcc1e 100644 --- a/docs/src/man/tensor.rst +++ b/docs/src/man/tensor.rst @@ -7,7 +7,9 @@ Definitions We start with the definition of a tensor. According to the `Link Wikipedia page ` on tensors, there are 3 equivalent ways of describing a tensor object: * As a multidimensional array. + * As a multilinear map. + * Using tensor products. From a MATLAB perspective, the former approach is the most natural, as such tensors are easily created with several builtin methods. We thus say a tensor of size ``sz`` can be thought of as the array that is created by calling any of the standard functions ``t = zeros(sz)``, ``t = ones(sz)``, ``t = rand(sz)``, ... @@ -38,18 +40,218 @@ In summary, we consider a tensor to be an object which can be represented in eit These definitions are equivalent, and it is always possible to go from one to the other. The main reason for introducing both is that for some algorithms, one is more convenient than the other, and thus this conversion will prove useful. -Creating and accessing tensors ------------------------------- +Creating tensors +---------------- There are several options provided for creating tensors, which are designed to either mimic the existing MATLAB constructors, or to be conveniently generalisable when dealing with symmetries. Almost each of these methods is a wrapper of some sorts around the general static constructor ``Tensor.new``: .. automethod:: src.tensors.Tensor.new :noindex: +In summary, for tensors without any internal symmetry, the easiest way of constructing tensors is by using the form :code:`Tensor.new(fun, [dims...])`, where optionally it is possible to specify arrows for each of the dimensions with :code:`Tensor.new(fun, [dims...], [arrows...])`. +Additionally, it is possible to make the partition of domain and codomain by supplying a rank for the tensor :code:`Tensor.new(fun, [dims...], [arrows...], 'Rank', [r1 r2])`. +While the latter two options do not alter any of the results in the non symmetric case, adding arrows is often useful from a debugging point of view, as it will then become less likely to accidentaly connect tensors that are not compatible. Additionally, the rank determines how the tensor acts as a linear map, which becomes important for decompositions, eigenvalues, etc. +Finally, the most general form is found by explicitly constructing the vector spaces upon which the tensor will act. +For tensors without symmetries, :class:`ComplexSpace` and :class:`CartesianSpace` represent spaces with or without arrows respectively, while :class:`GradedSpace` is used for tensors which have internal symmetries. +These allow for the final form :code:`Tensor.new(fun, codomain, domain)`. +On top of this, several convenience methods exist for commonly used initializers, such as: + +* :code:`Tensor.zeros` +* :code:`Tensor.ones` +* :code:`Tensor.eye` +* :code:`Tensor.rand` +* :code:`Tensor.randn` +* :code:`Tensor.randc` +* :code:`Tensor.randnc` + +Often, it can prove useful to create new tensors based on existing ones, either to ensure that they can be added, subtracted, or contracted. +In these cases, :code:`Tensor.similar` will be the method of choice, which acts as a wrapper to avoid manually having to select spaces from the tensors. + +.. automethod:: src.tensors.Tensor.similar + :noindex: + + +Accessing tensor data +--------------------- + +In general, :class:`Tensor` objects do not allow indexing or slicing operations, as this is not compatible with their internal structure. +Nevertheless, in accordance with the two ways in which a tensor can be represented, it is possible to access the tensor data in two ways. + +First off, when we consider a tensor as a linear operator which maps the domain to the codomain, we can ask for the matrix representation of this map. +For tensors with internal symmetries this will in general be a list of matrices, that represent the block-diagonal part of the tensor, along with the corresponding charge through :code:`matrixblocks(t)`. + +Secondly, when we want to think of a tensor more in terms of an N-dimensional array, we can access these as :code:`tensorblocks(t)`. +For tensors with internal symmetries, this will generate a list of all channels that are not explicitly zero by virtue of the symmetry, along with a representation of these channels, which is called a :class:`FusionTree`. + +Additional details for the symmetric tensors can be found in the :ref:`Symmetries` section. + +As an example, it could prove educational to understand the sizes of the lists and the sizes of the blocks generated by the example code below: + +.. code-block:: matlab + + >> A = Tensor.zeros([2 2 2], 'Rank', [2 1]); + >> Ablocks = matrixblocks(A) + + Ablocks = + + 1×1 cell array + + {4×2 double} + + >> Atensors = tensorblocks(A) + + Atensors = + + 1×1 cell array + + {2×2×2 double} + + >> z2space = GradedSpace.new(Z2(0, 1), [1 1], false); + >> B = Tensor.zeros([z2space z2space], z2space); + >> [Bblocks, Bcharges] = matrixblocks(B) + + Bblocks = + + 1×2 cell array + + {2×1 double} {2×1 double} + + + Bcharges = + + 1×2 Z2: + + logical data: + 0 1 + + >> [Btensors, Btrees] = tensorblocks(B) + + Btensors = + + 4×1 cell array + + {[0]} + {[0]} + {[0]} + {[0]} + + + Btrees = + + (2, 1) Z2 FusionTree array: + + isdual: + 0 0 0 + charges: + + + | + | + + - - | + | + + - + | - | - + + - | - | - + +In the very same way, in order to write data into a tensor, the same two formats can be used. + +First off, :code:`t = fill_matrix(t, blocks)` will take a list of blocks and fill these into the tensor. +This requires the list to be full, and sorted according to the charge, or in other words it has to be of the same shape and order as the output of :code:`matrixblocks`. +If it is only necessary to change some of the blocks, :code:`t = fill_matrix(t, blocks, charges)` additionally passes on an array of charges which specifies which block will be filled. + +Similarly, :code:`t = fill_tensor(t, tensors)` will take a list of N-dimensional arrays and fill these into the tensor, in the same order and shape of the output of :code:`tensorblocks`. +If it is required to only change some of the tensors, an array of :class:`FusionTree` s can be passed in as :code:`t = fill_tensor(t, tensors, trees)` to specify which tensors should be changed. + +.. code-block:: matlab + + >> Ablocks{1} = ones(size(Ablocks{1})); + >> A = fill_matrix(A, Ablocks); + >> Atensors{1} = rand(size(Atensors{1})); + >> A = fill_matrix(A, Atensors); + + >> Bblocks = cellfun(@(x) ones(size(x)), Bblocks, 'UniformOutput', false); + >> B = fill_matrix(B, Bblocks); + >> Bblocks{1} = Bblocks{1} + 1; + >> B = fill_matrix(B, Bblocks(1), Bcharges(1)); + +Additionally, it is also possible to use initializers instead of a list of data. +These initializers should have signature :code:`fun(dims, identifier)`. +For non symmetric tensors, ``identifier`` will be omitted, but for symmetric tensors the matrix case uses charges as ``identifier``, while the tensor case uses fusion trees as ``identifier``. +Again, it is possible to select only some of the blocks through the third argument. + +.. code-block:: matlab + + >> f = @(dims, identifier) ones(dims); + >> A = fill_matrix(A, f); + >> A = fill_tensor(A, f); + + >> g = @(dims, charge) qdim(charge) * ones(dims); + >> B = fill_matrix(B, g); + >> h = @(dims, tree) qdim(f.coupled) * ones(dims); + >> B = fill_tensor(B, h); + +Finally, we mention that for most tensors, it is possible to generate an N-dimensional array representation, at the cost of losing all information about the symmetries. +This can sometimes be useful as a tool for debugging, and can be accessed through :code:`a = double(t)`. + + Index manipulations ------------------- +Once a tensor has been created, it is possible to manipulate the order and partition of the indices through the use of :code:`permute(t, p, r)`. +This methods works similarly to :code:`permute` for arrays, as it requires a permutation vector ``p`` for determining the new order of the tensor legs. +Additionally and optionally, one may specify a rank `r` to determine the partition of the resulting tensor. +In order to only change the partition without permuting indices, :code:`repartition(t, r)` also is implemented. + +.. code-block:: matlab + + >> A = Tensor.rand([1 2 3]) + + A = + + Rank (3, 0) Tensor: + + 1. CartesianSpace of dimension 1 + + 2. CartesianSpace of dimension 2 + + 3. CartesianSpace of dimension 3 + + >> A2 = repartition(A, [1 2]) + + A2 = + + Rank (1, 2) Tensor: + + 1. CartesianSpace of dimension 1 + + 2. CartesianSpace of dimension 2 + + 3. CartesianSpace of dimension 3 + + >> A3 = permute(A, [3 2 1]) + + A3 = + + Rank (3, 0) Tensor: + + 1. CartesianSpace of dimension 3 + + 2. CartesianSpace of dimension 2 + + 3. CartesianSpace of dimension 1 + + >> A3 = permute(A, [3 2 1], [2 1]) + + A3 = + + Rank (2, 1) Tensor: + + 1. CartesianSpace of dimension 3 + + 2. CartesianSpace of dimension 2 + + 3. CartesianSpace of dimension 1 + +.. note:: + + While the partition of tensor indices might seem of no importance for tensors without internal structure, it can still have non-trivial consequences. + This is demonstrated by comparing the ``matrixblocks`` and the ``tensorblocks`` before and after repartitioning. Contractions ------------ diff --git a/src/tensors/FusionTree.m b/src/tensors/FusionTree.m index ba9fb75..176be11 100644 --- a/src/tensors/FusionTree.m +++ b/src/tensors/FusionTree.m @@ -301,6 +301,7 @@ for j = 1:length(ia2) blocks{j} = Rsymbol(abc(j, 1), abc(j, 2), abc(j, 3), inv); end + f.charges = f.charges(order, [2 1 3:end]); f.vertices = f.vertices(order, :); f.isdual(1:2) = f.isdual([2 1]); @@ -773,6 +774,45 @@ [f, p] = sort(f); c = c(:, p); end + + function [c, f] = twist(f, i, inv) + % Compute the coefficients that twist legs. + % + % Arguments + % --------- + % f : :class:`FusionTree` + % tree to repartition. + % + % i : int or logical + % indices of legs to twist. + % + % inv : logical + % flag to determine inverse twisting. + % + % Returns + % ------- + % c : sparse double + % matrix of coefficients that transform input to output trees. + % `f(i) --> c(i,j) * f(j)` + % + % f : :class:`FusionTree` + % twisted trees in canonical form. + + arguments + f + i + inv = false + end + + if istwistless(braidingstyle(f)) + c = speye(length(f)); + return + end + + theta = prod(twist(f.uncoupled(:, i)), 2); + if inv, theta = conj(theta); end + c = spdiags(theta, 0, length(f), length(f)); + end end @@ -781,6 +821,7 @@ function style = braidingstyle(f) style = f.charges.braidingstyle; end + function bool = isallowed(f) if ~hasmultiplicity(fusionstyle(f)) if f.rank(1) == 0 @@ -841,13 +882,16 @@ end function style = fusionstyle(f) - try style = fusionstyle(f.charges); - catch - bla - end end + function [lia, locb] = ismember(f1, f2) + if hasmultiplicity(fusionstyle(f1)) + error('TBA'); + end + + [lia, locb] = ismember(f1.charges, f2.charges, 'rows'); + end function [f, p] = sort(f) % sort - Sort the fusion trees into canonical order. @@ -929,7 +973,8 @@ function varargout = size(f, varargin) if nargin == 1 if nargout < 2 - varargout{1} = [1, size(f.charges, 1)]; + c = f.charges; + varargout{1} = [1, size(c, 1)]; elseif nargout == 2 varargout = {1, size(f.charges, 1)}; else diff --git a/src/tensors/Tensor.m b/src/tensors/Tensor.m index 038909f..5045da2 100644 --- a/src/tensors/Tensor.m +++ b/src/tensors/Tensor.m @@ -108,9 +108,17 @@ if isnumeric(data), data = {data}; end if iscell(data) - t.var = fill_matrix_data(t.var, data, charges); + if isempty(charges) + t.var = fill_matrix_data(t.var, data); + else + t.var = fill_matrix_data(t.var, data, charges); + end else - t.var = fill_matrix_fun(t.var, data, charges); + if isempty(charges) + t.var = fill_matrix_fun(t.var, data); + else + t.var = fill_matrix_fun(t.var, data, charges); + end end end @@ -310,7 +318,8 @@ arguments kwargs.Rank - kwargs.Mode = 'matrix' + kwargs.Mode {mustBeMember(kwargs.Mode, {'matrix', 'tensor'})} ... + = 'matrix' end % Parse special constructors @@ -355,38 +364,41 @@ % Fill tensor if ~isempty(fun) - if strcmp(kwargs.Mode, 'matrix') - t = fill_matrix(t, fun); - elseif strcmp(kwargs.Mode, 'tensor') - t = fill_tensor(t, fun); - else - error('tensors:ArgumentError', 'Unknown options for Mode'); + switch kwargs.Mode + case 'matrix' + t = fill_matrix(t, fun); + case 'tensor' + t = fill_tensor(t, fun); end end end function t = zeros(varargin) - t = Tensor.new(@(dims, charge) zeros(dims), varargin{:}); + t = Tensor.new(@zeros, varargin{:}); end function t = ones(varargin) - t = Tensor.new(@(dims, charge) ones(dims), varargin{:}); + t = Tensor.new(@ones, varargin{:}); end function t = eye(varargin) - t = Tensor.new(@(dims, charge) eye(dims), varargin{:}); + t = Tensor.new(@eye, varargin{:}); end function t = rand(varargin) - t = Tensor.new(@(dims, charge) rand(dims), varargin{:}); + t = Tensor.new(@rand, varargin{:}); + end + + function t = randn(varargin) + t = Tensor.new(@randn, varargin{:}); end function t = randc(varargin) - t = Tensor.new(@(dims, charge) randc(dims), varargin{:}); + t = Tensor.new(@randc, varargin{:}); end function t = randnc(varargin) - t = Tensor.new(@(dims, charge) randnc(dims), varargin{:}); + t = Tensor.new(@randnc, varargin{:}); end end @@ -430,13 +442,61 @@ f = fusiontrees(t.codomain, t.domain); end - function b = tensorblocks(t) + function [b, f] = tensorblocks(t) b = tensorblocks(t.var); + if nargout > 1, f = fusiontrees(t); end end function varargout = matrixblocks(t) [varargout{1:nargout}] = matrixblocks(t.var); end + + function style = braidingstyle(t) + style = braidingstyle(t.codomain, t.domain); + end + + function style = fusionstyle(t) + style = fusionstyle(t.codomain, t.domain); + end + + function tdst = insert_onespace(tsrc, i, dual) + arguments + tsrc + i = nspaces(tsrc) + dual = false + end + + spaces = insertone(space(tsrc), i, dual); + data = matrixblocks(tsrc); + + r = rank(tsrc); + if i <= r(1) + r(1) = r(1) + 1; + else + r(2) = r(2) + 1; + end + tdst = Tensor(spaces(1:r(1)), spaces(r(1)+1:end)'); + tdst.var = fill_matrix_data(tdst.var, data); + end + + function tdst = embed(tsrc, tdst) + + bsrc = tensorblocks(tsrc); + fsrc = fusiontrees(tsrc); + bdst = tensorblocks(tdst); + fdst = fusiontrees(tdst); + + [lia, locb] = ismember(fsrc, fdst); + nsp = nspaces(tdst); + + for i = find(lia).' + sz = min(size(bsrc{i}, 1:nsp), size(bdst{locb(i)}, 1:nsp)); + inds = arrayfun(@(x) 1:x, sz, 'UniformOutput', false); + bdst{locb(i)}(inds{:}) = bsrc{i}(inds{:}); + end + + tdst.var = fill_tensor_data(tdst.var, bdst); + end end @@ -559,6 +619,10 @@ 'Subtraction with scalars only defined for square tensors.'); I = t1(i).eye(t1(i).codomain, t1(i).domain); t1(i).var = t1(i).var - I.var .* t2(i); + elseif iszero(t2(i)) + continue; + elseif iszero(t1(i)) + t1(i) = -t2(i); else assert(isequal(t1(i).domain, t2(i).domain) && ... isequal(t1(i).codomain, t2(i).codomain), 'tensors:SpaceMismatch', ... @@ -590,7 +654,6 @@ t = inv(t1) * t2; end - function t = mrdivide(t1, t2) % Right division of tensors. % @@ -632,35 +695,37 @@ % C : :class:`Tensor` % output tensor. - szA = size(A); - szB = size(B); - assert(szA(2) == szB(1)); - - % Scalar multiplications - if isnumeric(A) || isnumeric(B) + if ~isscalar(A) || ~isscalar(B) + szA = size(A); + szB = size(B); + assert(szA(2) == szB(1)); + + % Tensor multiplications for i = szA(1):-1:1 for j = szB(2):-1:1 - C(i,j) = A(i, 1) .* B(1, j); + C(i,j) = A(i, 1) * B(1, j); for k = 2:szA(2) - C(i,j) = C(i,j) + A(i, k) .* B(k, j); + C(i,j) = C(i,j) + A(i, k) * B(k, j); end end end return end - % Tensor multiplications - for i = szA(1):-1:1 - for j = szB(2):-1:1 - C(i,j) = tensorprod(A(i, 1), B(1, j), ... - rank(A(i, 1), 1) + (1:rank(A(i, 1), 2)), ... - rank(B(1, j), 1):-1:1); - for k = 2:szA(2) - C(i,j) = C(i,j) + tensorprod(A(i, k), B(k, j), ... - rank(A(i, k), 1) + (1:rank(A(i, k), 2)), ... - rank(B(k, j), 1):-1:1); - end - end + if isnumeric(A) || isnumeric(B) + C = A .* B; + return + end + + assert(isequal(A.domain, B.codomain), 'tensors:SpaceMismatch', ... + 'Multiplied spaces incompatible.'); + if ~isempty(A.codomain) || ~isempty(B.domain) + C = Tensor.zeros(A.codomain, B.domain); + C.var = mul(C.var, A.var, B.var); + else + Ablocks = matrixblocks(A.var); + Bblocks = matrixblocks(B.var); + C = horzcat(Ablocks{:}) * vertcat(Bblocks{:}); end end @@ -689,6 +754,7 @@ case {1, 2, Inf} n = 0; for i = 1:numel(t) + if iszero(t(i)), continue; end mblocks = matrixblocks(t(i).var); for j = 1:length(mblocks) n = max(norm(mblocks{j}, p), n); @@ -698,6 +764,7 @@ case 'fro' n = 0; for i = 1:numel(t) + if iszero(t(i)), continue; end [mblocks, mcharges] = matrixblocks(t(i).var); qdims = qdim(mcharges); for j = 1:length(mblocks) @@ -727,6 +794,10 @@ 'Addition with scalars only defined for square tensors.'); I = t1(i).eye(t1(i).codomain, t1(i).domain); t1(i).var = t1(i).var + I.var .* t2(i); + elseif iszero(t2(i)) + continue; + elseif iszero(t1(i)) + t1(i) = t2(i); else assert(isequal(t1(i).domain, t2(i).domain) && ... isequal(t1(i).codomain, t2(i).codomain), 'tensors:SpaceMismatch', ... @@ -904,9 +975,16 @@ med = get(cache, key); if isempty(med) med = struct; + A_ = similar(@(x,charge) uninit(x), A, iA, 'Rank', rA); med.varA = A_.var; - med.mapA = permute(fusiontrees(A), iA, rA); + [med.mapA, f] = permute(fusiontrees(A), iA, rA); + for i = rA(1) + (1:rA(2)) + if ~isdual(space(A_, i)) + [c2, f] = twist(f, i); + med.mapA = med.mapA * c2; + end + end B_ = similar(@(x,charge) uninit(x), B, iB, 'Rank', rB); med.varB = B_.var; @@ -915,7 +993,7 @@ assert(isequal(A_.domain, B_.codomain), 'tensors:SpaceMismatch', ... 'Contracted spaces incompatible.'); if ~isempty(A_.codomain) || ~isempty(B_.domain) - med.C = Tensor(A_.codomain, B_.domain); + med.C = Tensor.zeros(A_.codomain, B_.domain); else med.C = []; end @@ -936,6 +1014,10 @@ end A = permute(A, iA, rA); + for i = rA(1) + (1:rA(2)) + if ~isdual(space(A, i)), A = twist(A, i); end + end +% A = twist(A, [false(1, length(A.codomain)) ~isdual(A.domain')]); B = permute(B, iB, rB); assert(isequal(A.domain, B.codomain), 'tensors:SpaceMismatch', ... @@ -1032,6 +1114,53 @@ error('tensors:TBA', 'This method has not been implemented.'); end + function t = twist(t, i, inv) + % Twist the spaces of a tensor. + % + % Arguments + % --------- + % t : :class:`Tensor` + % input tensor. + % + % i : (1, :) int or logical + % indices to twist. + % + % inv : logical + % flag to indicate inverse twisting. + % + % Returns + % ------- + % t : :class:`Tensor` + % permuted tensor with desired rank. + + arguments + t + i + inv = false + end + + if isempty(i) || ~any(i) || istwistless(braidingstyle(t)) + return + end + + persistent cache + if isempty(cache), cache = LRU; end + + if Options.CacheEnabled() + key = GetMD5({GetMD5_helper(t.codomain), GetMD5_helper(t.domain), i, inv}, ... + 'Array', 'hex'); + med = get(cache, key); + if isempty(med) + med = twist(fusiontrees(t), i, inv); + cache = set(cache, key, med); + end + else + med = twist(fusiontrees(t), i, inv); + end + + t.var = axpby(1, t.var, 0, t.var, 1:nspaces(t), med); + end + function t = uplus(t) end @@ -1493,6 +1622,10 @@ U.var = fill_matrix_data(U.var, Us, dims.charges); S.var = fill_matrix_data(S.var, Ss, dims.charges); V.var = fill_matrix_data(V.var, Vs, dims.charges); + + if nargout <= 1 + U = S; + end end end @@ -1752,6 +1885,10 @@ end end + function bool = iszero(t) + bool = isempty(t.var); + end + function r = cond(t, p) % Condition number with respect to inversion. This is defined as % :math:`||t|| * ||t^{-1}||` in the p-norm. For well conditioned @@ -1874,10 +2011,10 @@ b_vec = vectorize(b); b_sz = size(b_vec); - if ~isempty(x0) + if ~isempty(x0) && ~iszero(x0) x0_vec = vectorize(x0); else - x0_vec = []; + x0_vec = zeros(b_sz); end % Convert input operators to handle vectors @@ -2026,6 +2163,7 @@ 'largestreal', 'smallestreal', 'bothendsreal', ... 'largestimag', 'smallestimag', 'bothendsimag'}), ... 'tensors:ArgumentError', 'Invalid choice of eigenvalue selector.'); + nargoutchk(0, 3); x0_vec = vectorize(x0); sz = size(x0_vec); @@ -2038,15 +2176,35 @@ options.KrylovDim = min(sz(1), options.KrylovDim); - [varargout{1:nargout}] = eigs(A_fun, sz(1), howmany, sigma, ... + [V, D, flag] = eigs(A_fun, sz(1), howmany, sigma, ... 'Tolerance', options.Tol, 'MaxIterations', options.MaxIter, ... 'SubspaceDimension', options.KrylovDim, 'IsFunctionSymmetric', ... - options.IsSymmetric, 'StartVector', x0_vec); - if nargout > 1 + options.IsSymmetric, 'StartVector', x0_vec, ... + 'Display', options.Verbosity == 3); + + if nargout <= 1 + varargout = {D}; + elseif nargout == 2 + for i = howmany:-1:1 + varargout{1}(:, i) = devectorize(V(:, i), x0); + end + varargout{2} = D; + else + varargout{2} = D; for i = howmany:-1:1 - V(:, i) = devectorize(varargout{1}(:, i), x0); + varargout{1}(:, i) = devectorize(V(:, i), x0); end - varargout{1} = V; + varargout{3} = flag; + end + + if flag && options.Verbosity > 0 + warning('eigsolve did not converge.'); + end + + if options.Verbosity > 1 + if flag + fprintf('eigsolve converged.\n'); + end end end @@ -2157,6 +2315,19 @@ a = cell2mat(a_cell); end + + function tdst = desymmetrize(tsrc, mode) + arguments + tsrc + mode {mustBeMember(mode, {'Cartesian', 'Complex'})} = 'Complex' + end + + tdst = repartition(Tensor(double(tsrc)), rank(tsrc)); + if strcmp(mode, 'Complex') + if ~isempty(tsrc.domain), tdst.domain = ComplexSpace(tsrc.domain); end + if ~isempty(tsrc.codomain), tdst.codomain = ComplexSpace(tsrc.codomain); end + end + end end diff --git a/src/tensors/charges/AbstractCharge.m b/src/tensors/charges/AbstractCharge.m index 0cfaf9c..246d74b 100644 --- a/src/tensors/charges/AbstractCharge.m +++ b/src/tensors/charges/AbstractCharge.m @@ -686,6 +686,55 @@ reshape(repmat(X{i}, size(Y, 2), 1), size(X{i}, 1), [])]; end end + + function [lia, locb] = ismember_sorted(a, b) + if isempty(a) || isempty(b) + lia = false(size(a)); + locb = zeros(size(a)); + return + end + + if isscalar(a) || isscalar(b) + if any(size(a) == numel(a)) && any(size(b) == numel(b)) + lia = a == b; + else + lia = a(:) == b(:); + end + + if ~any(lia) + lia = false(size(a)); + locb = zeros(size(a)); + return + end + if ~isscalar(b) + locb = find(lia); + locb = locb(1); + lia = any(lia); + else + locb = double(lia); + end + return + end + + [sortab, indsortab] = sort([a(:); b(:)]); + d = sortab(1:end-1) == sortab(2:end); + ndx1 = indsortab(d); + + if nargout <= 1 + lia = ismember(1:length(a), ndx1); + else + szuA = length(a); + d = find(d); + [lia, locb] = ismember(1:szuA, ndx1); + newd = d(locb(lia)); + where = indsortab(newd+1) - szuA; + locb(lia) = where; + end + lia = reshape(lia, size(a)); + if nargout > 1 + locb = reshape(locb, size(a)); + end + end end methods diff --git a/src/tensors/kernels/AbelianBlock.m b/src/tensors/kernels/AbelianBlock.m index 509a6e4..17a162c 100644 --- a/src/tensors/kernels/AbelianBlock.m +++ b/src/tensors/kernels/AbelianBlock.m @@ -19,31 +19,38 @@ %% General case: addition with permutation % tensor indexing to matrix indexing - rrx = rankrange(X(1).rank); - rry = rankrange(Y(1).rank); + rx = X(1).rank; + ry = Y(1).rank; + rrx = rankrange(rx); + rry = rankrange(ry); p_eff(rry) = rrx(p); + p_eff_1 = p_eff(1:ry(1)); + p_eff_2 = p_eff((1:ry(2)) + ry(1)); % extract small tensor blocks doA = a ~= 1; - vars = cell(size(map, 1), 1); + vars_in = cell(size(map, 1), 1); offset = 0; + + Xrowsizes = {X.rowsizes}; + Xcolsizes = {X.colsizes}; + Xtdims = {X.tdims}; + Xvar = {X.var}; for i = 1:length(X) - rowsz = X(i).rowsizes; - colsz = X(i).colsizes; + rowsz = Xrowsizes{i}; + colsz = Xcolsizes{i}; - var_in = X(i).var; + var_in = Xvar{i}; if doA, var_in = a .* var_in; end - oldtdims = X(i).tdims(:, rrx); - newmdims = [prod(X(i).tdims(:, p_eff(1:Y(1).rank(1))), 2) ... - prod(X(i).tdims(:, p_eff((1:Y(1).rank(2)) + Y(1).rank(1))), 2)]; + oldtdims = Xtdims{i}(:, rrx); + newmdims = [prod(oldtdims(:, p_eff_1), 2) prod(oldtdims(:, p_eff_2), 2)]; for k = 1:length(colsz) - 1 for j = 1:length(rowsz) - 1 ind = j + (k-1) * (length(rowsz) - 1); offset = offset + 1; - - vars{offset} = reshape(permute(reshape(... + vars_in{offset} = reshape(permute(reshape(... var_in(rowsz(j)+1:rowsz(j+1), colsz(k)+1:colsz(k+1)), ... oldtdims(ind, :)), ... p_eff), ... @@ -53,8 +60,11 @@ end % apply map - [row, col] = find(map); - vars(col) = vars(row); + vars_out = cell(size(map, 2), 1); + [rows, cols, vals] = find(map); + for i = 1:length(vals) + vars_out{cols(i)} = vals(i) * vars_in{rows(i)}; + end % inject small tensor blocks if b == 0 @@ -65,13 +75,13 @@ if rows < cols m = cell(rows, 1); for n = 1:rows - m{n} = cat(2, vars{offset + n + ((1:cols)-1) * rows}); + m{n} = cat(2, vars_out{offset + n + ((1:cols)-1) * rows}); end Y(i).var = cat(1, m{:}); else m = cell(cols, 1); for n = 1:cols - m{n} = cat(1, vars{offset + (n-1) * rows + (1:rows)}); + m{n} = cat(1, vars_out{offset + (n-1) * rows + (1:rows)}); end Y(i).var = cat(2, m{:}); end @@ -88,13 +98,13 @@ if rows < cols m = cell(rows, 1); for n = 1:rows - m{n} = cat(2, vars{offset + n + ((1:cols)-1) * rows}); + m{n} = cat(2, vars_out{offset + n + ((1:cols)-1) * rows}); end Y(i).var = Y(i).var + cat(1, m{:}); else m = cell(cols, 1); for n = 1:cols - m{n} = cat(1, vars{offset + (n-1) * rows + (1:rows)}); + m{n} = cat(1, vars_out{offset + (n-1) * rows + (1:rows)}); end Y(i).var = Y(i).var + cat(2, m{:}); end @@ -110,128 +120,18 @@ if rows < cols m = cell(rows, 1); for n = 1:rows - m{n} = cat(2, vars{offset + n + ((1:cols)-1) * rows}); + m{n} = cat(2, vars_out{offset + n + ((1:cols)-1) * rows}); end Y(i).var = b .* Y(i).var + cat(1, m{:}); else m = cell(cols, 1); for n = 1:cols - m{n} = cat(1, vars{offset + (n-1) * rows + (1:rows)}); + m{n} = cat(1, vars_out{offset + (n-1) * rows + (1:rows)}); end Y(i).var = b .* Y(i).var + cat(2, m{:}); end offset = offset + rows * cols; end - -% rrx = rankrange(x(1).rank); -% rry = rankrange(y(1).rank); -% p_eff(rry) = rrx(p); - -% %% extract blocks -% doA = a ~= 1; -% vars = cell(1, size(map, 1)); -% offset = 0; -% for block = x -% rowsz = block.rowsizes; -% colsz = block.colsizes; -% -% if doA -% varA = a .* block.var; -% else -% varA = block.var; -% end -% -% olddims = block.tdims(:, rrx); -% newdims = [prod(olddims(:, p_eff(1:y(1).rank(1))), 2) ... -% prod(olddims(:, p_eff(y(1).rank(1) + (1:y(1).rank(2)))), 2)]; -% for k = 1:length(colsz) - 1 -% for j = 1:length(rowsz) - 1 -% ind = j + (k-1) * (length(rowsz) - 1); -% offset = offset + 1; -% % vars{offset} = varA(rowsz(j)+1:rowsz(j+1), colsz(k)+1:colsz(k+1)); -% % vars{offset} = reshape(vars{offset}, olddims(ind, :)); -% % vars{offset} = permute(vars{offset}, p_eff); -% % vars{offset} = reshape(vars{offset}, newdims(ind, :)); -% -% % less readible but faster -% vars{offset} = reshape(permute(reshape(... -% varA(rowsz(j)+1:rowsz(j+1), colsz(k)+1:colsz(k+1)), ... -% olddims(ind, :)), p_eff), ... -% newdims(ind, :)); -% end -% end -% end -% -% %% apply map -% [row, col] = find(map); -% vars(col) = vars(row); -% -% %% inject blocks -% if b == 0 -% offset = 0; -% for i = 1:length(y) -% rows = length(y(i).rowsizes) - 1; -% cols = length(y(i).colsizes) - 1; -% if rows < cols -% m = cell(rows, 1); -% for n = 1:rows -% m{n} = cat(2, vars{offset + n + ((1:cols)-1) * rows}); -% end -% y(i).var = cat(1, m{:}); -% else -% m = cell(cols, 1); -% for n = 1:cols -% m{n} = cat(1, vars{offset + (n-1) * rows + (1:rows)}); -% end -% y(i).var = cat(2, m{:}); -% end -% offset = offset + rows * cols; -% end -% return -% end -% -% if b == 1 -% offset = 0; -% for i = 1:length(y) -% rows = length(y(i).rowsizes) - 1; -% cols = length(y(i).colsizes) - 1; -% if rows < cols -% m = cell(rows, 1); -% for n = 1:rows -% m{n} = cat(2, vars{offset + n + ((1:cols)-1) * rows}); -% end -% y(i).var = y(i).var + cat(1, m{:}); -% else -% m = cell(cols, 1); -% for n = 1:cols -% m{n} = cat(1, vars{offset + (n-1) * rows + (1:rows)}); -% end -% y(i).var = y(i).var + cat(2, m{:}); -% end -% offset = offset + rows * cols; -% end -% return -% end -% -% offset = 0; -% for i = 1:length(y) -% rows = length(y(i).rowsizes) - 1; -% cols = length(y(i).colsizes) - 1; -% if rows < cols -% m = cell(rows, 1); -% for n = 1:rows -% m{n} = cat(2, vars{offset + n + ((1:cols)-1) * rows}); -% end -% y(i).var = b .* y(i).var + cat(1, m{:}); -% else -% m = cell(cols, 1); -% for n = 1:cols -% m{n} = cat(1, vars{offset + (n-1) * rows + (1:rows)}); -% end -% y(i).var = b .* y(i).var + cat(2, m{:}); -% end -% offset = offset + rows * cols; -% end end function typename = underlyingType(b) diff --git a/src/tensors/kernels/AbstractBlock.m b/src/tensors/kernels/AbstractBlock.m index 128b1b5..5511b03 100644 --- a/src/tensors/kernels/AbstractBlock.m +++ b/src/tensors/kernels/AbstractBlock.m @@ -34,13 +34,30 @@ return end - if braidingstyle(codomain, domain) == BraidingStyle.Abelian && ... - fusionstyle(codomain, domain) == FusionStyle.Unique - X = AbelianBlock(codomain, domain); - return - end + persistent cache + if isempty(cache), cache = LRU; end - X = MatrixBlock(codomain, domain); + if Options.CacheEnabled() + key = GetMD5({codomain, domain}, 'Array', 'hex'); + med = get(cache, key); + if isempty(med) + if braidingstyle(codomain, domain) == BraidingStyle.Abelian && ... + fusionstyle(codomain, domain) == FusionStyle.Unique + med = AbelianBlock(codomain, domain); + else + med = MatrixBlock(codomain, domain); + end + cache = set(cache, key, med); + end + else + if braidingstyle(codomain, domain) == BraidingStyle.Abelian && ... + fusionstyle(codomain, domain) == FusionStyle.Unique + med = AbelianBlock(codomain, domain); + else + med = MatrixBlock(codomain, domain); + end + end + X = med; end end diff --git a/src/tensors/kernels/MatrixBlock.m b/src/tensors/kernels/MatrixBlock.m index a23ce38..d74859b 100644 --- a/src/tensors/kernels/MatrixBlock.m +++ b/src/tensors/kernels/MatrixBlock.m @@ -92,7 +92,7 @@ %% Special case 2: addition without permutation if nargin == 4 || (isempty(p) && isempty(map)) || ... - (all(p == 1:length(p)) && all(X(1).rank == Y(1).rank)) + (all(p == 1:length(p)) && isequal(map, speye(size(map)))) % reduce to scalar multiplication if b == 0 % a ~= 0 -> case 1 Y = X .* a; @@ -129,24 +129,32 @@ %% General case: addition with permutation % tensor indexing to matrix indexing - rrx = rankrange(X(1).rank); - rry = rankrange(Y(1).rank); + rx = X(1).rank; + ry = Y(1).rank; + rrx = rankrange(rx); + rry = rankrange(ry); p_eff(rry) = rrx(p); + p_eff_1 = p_eff(1:ry(1)); + p_eff_2 = p_eff((1:ry(2)) + ry(1)); % extract small tensor blocks doA = a ~= 1; vars_in = cell(size(map, 1), 1); offset = 0; + + Xrowsizes = {X.rowsizes}; + Xcolsizes = {X.colsizes}; + Xtdims = {X.tdims}; + Xvar = {X.var}; for i = 1:length(X) - rowsz = X(i).rowsizes; - colsz = X(i).colsizes; + rowsz = Xrowsizes{i}; + colsz = Xcolsizes{i}; - var_in = X(i).var; + var_in = Xvar{i}; if doA, var_in = a .* var_in; end - oldtdims = X(i).tdims(:, rrx); - newmdims = [prod(X(i).tdims(:, p_eff(1:Y(1).rank(1))), 2) ... - prod(X(i).tdims(:, p_eff((1:Y(1).rank(2)) + Y(1).rank(1))), 2)]; + oldtdims = Xtdims{i}(:, rrx); + newmdims = [prod(oldtdims(:, p_eff_1), 2) prod(oldtdims(:, p_eff_2), 2)]; for k = 1:length(colsz) - 1 for j = 1:length(rowsz) - 1 @@ -173,11 +181,13 @@ end % inject small tensor blocks + Yrowsizes = {Y.rowsizes}; + Ycolsizes = {Y.colsizes}; if b == 0 offset = 0; for i = 1:length(Y) - rows = length(Y(i).rowsizes) - 1; - cols = length(Y(i).colsizes) - 1; + rows = length(Yrowsizes{i}) - 1; + cols = length(Ycolsizes{i}) - 1; if rows < cols m = cell(rows, 1); for n = 1:rows @@ -199,8 +209,8 @@ if b == 1 offset = 0; for i = 1:length(Y) - rows = length(Y(i).rowsizes) - 1; - cols = length(Y(i).colsizes) - 1; + rows = length(Yrowsizes{i}) - 1; + cols = length(Ycolsizes{i}) - 1; if rows < cols m = cell(rows, 1); for n = 1:rows @@ -221,8 +231,8 @@ offset = 0; for i = 1:length(Y) - rows = length(Y(i).rowsizes) - 1; - cols = length(Y(i).colsizes) - 1; + rows = length(Yrowsizes{i}) - 1; + cols = length(Ycolsizes{i}) - 1; if rows < cols m = cell(rows, 1); for n = 1:rows @@ -265,8 +275,22 @@ % C : MatrixBlock % Result of computing C = (A .* a) * (B .* b) - [indA, locA] = ismember([C.charge], [A.charge]); - [indB, locB] = ismember([C.charge], [B.charge]); + Acharge = [A.charge]; + Bcharge = [B.charge]; + Ccharge = [C.charge]; + + if isequal(Acharge, Ccharge) + indA = true(size(A)); + locA = 1:length(A); + else + [indA, locA] = ismember_sorted(Ccharge, Acharge); + end + if isequal(Bcharge, Ccharge) + indB = true(size(B)); + locB = 1:length(B); + else + [indB, locB] = ismember_sorted(Ccharge, Bcharge); + end if nargin == 3 || (a == 1 && b == 1) for i = find(indA & indB) @@ -372,7 +396,7 @@ function b = fill_matrix_fun(b, fun, charges) if nargin < 3 || isempty(charges) for i = 1:length(b) - b(i).var = fun(size(b(i).var), b(i).charge); + b(i).var = fun(size(b(i).var)); end else [lia, locb] = ismember(charges, [b.charge]); @@ -382,6 +406,28 @@ end end + function b = fill_tensor_data(b, data) + ctr = 0; + p = rankrange(b(1).rank); + + for i = 1:length(b) + rowsz = b(i).rowsizes; + colsz = b(i).colsizes; + for k = 1:length(colsz) - 1 + for j = 1:length(rowsz) - 1 + ctr = ctr + 1; + b(i).var(rowsz(j)+1:rowsz(j+1), colsz(k)+1:colsz(k+1)) = ... + reshape(permute(data{ctr}, p), ... + rowsz(j+1)-rowsz(j), colsz(k+1)-colsz(k)); + end + end + end + end + + function b = fill_tensor_fun(b, fun, trees) + + end + function Y = minus(X, Y) % Subtraction of X and Y. % diff --git a/src/tensors/spaces/AbstractSpace.m b/src/tensors/spaces/AbstractSpace.m index 881f97d..c6c071d 100644 --- a/src/tensors/spaces/AbstractSpace.m +++ b/src/tensors/spaces/AbstractSpace.m @@ -3,12 +3,12 @@ properties (Access = protected) dimensions % Specification of the internal dimensions - isdual (1,1) logical = false % Flag to indicate if the space is a dual space + dual (1,1) logical = false % Flag to indicate if the space is a dual space end %% Constructors methods - function spaces = AbstractSpace(dimensions, isdual) + function spaces = AbstractSpace(dimensions, dual) % Construct an array of vector spaces. % % Repeating Arguments @@ -26,20 +26,20 @@ arguments (Repeating) dimensions - isdual + dual end if nargin == 0, return; end for i = length(dimensions):-1:1 spaces(i).dimensions = dimensions{i}; - spaces(i).isdual = isdual{i}; + spaces(i).dual = dual{i}; end end end methods (Static) - function spaces = new(dimensions, isdual) + function spaces = new(dimensions, dual) % Construct a vector space. This is a utility method to be able to access the % constructor of a subclass. % @@ -48,7 +48,7 @@ % dimensions : int or struct % a variable which represents the internal dimension of the space. % - % isdual : logical + % dual : logical % a variable which indicates if a space is dual. % % Returns @@ -63,6 +63,10 @@ %% Structure methods + function f = isdual(space) + f = [space.dual]; + end + function c = charges(spaces) % Compute all charge combinations of a space. If no internal structure is % present, this yields an empty result. @@ -179,7 +183,7 @@ % dual spaces. for i = 1:length(spaces) - spaces(i).isdual = ~spaces(i).isdual; + spaces(i).dual = ~spaces(i).dual; end end @@ -301,7 +305,7 @@ % hashable : cell % data which can be accepted by :func:`GetMD5`. - hashable = {spaces.dimensions, spaces.isdual}; + hashable = {spaces.dimensions, spaces.dual}; end end end diff --git a/src/tensors/spaces/CartesianSpace.m b/src/tensors/spaces/CartesianSpace.m index 1c40aa1..af7ddfc 100644 --- a/src/tensors/spaces/CartesianSpace.m +++ b/src/tensors/spaces/CartesianSpace.m @@ -279,7 +279,7 @@ function disp(spaces) if isscalar(spaces) - fprintf('CartesianSpace of dimension %d:\n', spaces.dimensions); + fprintf('CartesianSpace of dimension %d\n', spaces.dimensions); return end @@ -291,5 +291,17 @@ function disp(spaces) fprintf('\n'); end end + + function s = string(spaces) + if numel(spaces) > 1 + s = strings(size(spaces)); + for i = 1:numel(spaces) + s(i) = string(spaces); + end + return + end + + s = sprintf('CartesianSpace of dimension %d', spaces.dimensions); + end end end diff --git a/src/tensors/spaces/ComplexSpace.m b/src/tensors/spaces/ComplexSpace.m index 04b7da2..4b88c6b 100644 --- a/src/tensors/spaces/ComplexSpace.m +++ b/src/tensors/spaces/ComplexSpace.m @@ -190,7 +190,7 @@ % fused space. space.dimensions = space.dimensions * space2.dimensions; - space.isdual = false; + space.dual = false; end function space = prod(spaces) diff --git a/src/tensors/spaces/GradedSpace.m b/src/tensors/spaces/GradedSpace.m index 916ca67..de1d7ff 100644 --- a/src/tensors/spaces/GradedSpace.m +++ b/src/tensors/spaces/GradedSpace.m @@ -3,7 +3,7 @@ %% Constructors methods - function spaces = GradedSpace(dimensions, isdual) + function spaces = GradedSpace(dimensions, dual) % Construct an array of vector spaces. % % Repeating Arguments @@ -12,7 +12,7 @@ % internal structure of the vector space, with fields 'charges' and % 'degeneracies'. % - % isdual : (1, 1) logical + % dual : (1, 1) logical % flag to denote dual spaces. % % Returns @@ -22,7 +22,7 @@ arguments (Repeating) dimensions (1, 1) struct - isdual (1, 1) logical + dual (1, 1) logical end if nargin == 0 @@ -44,11 +44,17 @@ dimensions{i}.degeneracies = dimensions{i}.degeneracies(I); end - args = [dimensions; isdual]; + args = [dimensions; dual]; end spaces = spaces@AbstractSpace(args{:}); end + + function space = one(spaces) + space = GradedSpace(... + struct('charges', one(spaces(1).dimensions.charges), 'degeneracies', 1), ... + false); + end end methods (Static) @@ -58,9 +64,9 @@ % % Usage % ----- - % spaces = GradedSpace.new(charges, degeneracies, isdual, ...) + % spaces = GradedSpace.new(charges, degeneracies, dual, ...) % - % spaces = GradedSpace.new(dimensions, isdual, ...) + % spaces = GradedSpace.new(dimensions, dual, ...) % % Repeating Arguments % ------------------- @@ -73,7 +79,7 @@ % degeneracies : int % degeneracies for the internal structure of the space. % - % isdual : logical + % dual : logical % a variable which indicates if a space is dual. % % Returns @@ -119,14 +125,14 @@ % c : (:, :) :class:`AbstractCharge` % list of charge combinations, where each row is an entry. - if spaces(1).isdual + if isdual(spaces(1)) c = conj(spaces(1).dimensions.charges); else c = spaces(1).dimensions.charges; end for i = 2:length(spaces) - if spaces(i).isdual + if isdual(spaces(i)) c = combvec(c, conj(spaces(i).dimensions.charges)); else c = combvec(c, spaces(i).dimensions.charges); @@ -193,7 +199,7 @@ args = cell(2, sum(rank)); for i = 1:size(args, 2) args{1, i} = charges(spaces(i)); - args{2, i} = spaces(i).isdual; + args{2, i} = spaces(i).dual; end trees = FusionTree.new(rank, args{:}); @@ -257,6 +263,18 @@ space = GradedSpace(newdimensions, false); end + function spaces = insertone(spaces, i, dual) + arguments + spaces + i = length(spaces) + dual = false + end + + trivialspace = one(spaces); + if dual, trivialspace = conj(trivialspace); end + spaces = [spaces(1:i) trivialspace spaces(i+1:end)]; + end + %% Utility function bools = eq(spaces1, spaces2) @@ -277,7 +295,7 @@ return end - bools = [spaces1.isdual] == [spaces2.isdual]; + bools = [spaces1.dual] == [spaces2.dual]; if isscalar(spaces2) for i = 1:length(spaces1) @@ -298,6 +316,10 @@ end end + function bools = ne(spaces1, spaces2) + bools = ~(spaces1 == spaces2); + end + function bool = ge(space1, space2) bool = le(space2, space1); end @@ -378,20 +400,39 @@ end end + + function space = plus(space, space2) + assert(isscalar(space) && isscalar(space2)); + assert(space.dual == space2.dual); + + c = charges(space2); + d = degeneracies(space2); + [lia, locb] = ismember(c, charges(space)); + for i = find(lia)' + space.dimensions.degeneracies(locb(i)) = d(i) + ... + space.dimensions.degeneracies(locb(i)); + end + [space.dimensions.charges, p] = sort([space.dimensions.charges, ... + c(~lia)]); + space.dimensions.degeneracies = [space.dimensions.degeneracies, ... + d(~lia)]; + space.dimensions.degeneracies = space.dimensions.degeneracies(p); + end end methods function disp(spaces) % Custom display of spaces. if isscalar(spaces) - fprintf('%s GradedSpace of dimension %d:\n', ... - class(spaces.dimensions.charges), dims(spaces)); - title_str = strjust(pad(["isdual", "charges", "degeneracies"]), 'right'); - charge_str = strjust(pad([string(spaces.dimensions.charges) - string(spaces.dimensions.degeneracies)]), 'center'); - fprintf('\t%s:\t%s\n', title_str(1), string(spaces.isdual)); - fprintf('\t%s:\t%s\n', title_str(2), join(charge_str(1, :), char(9))); - fprintf('\t%s:\t%s\n', title_str(3), join(charge_str(2, :), char(9))); + fprintf('%s', string(spaces)); +% fprintf('%s GradedSpace of dimension %d:\n', ... +% class(spaces.dimensions.charges), dims(spaces)); +% title_str = strjust(pad(["dual", "charges", "degeneracies"]), 'right'); +% charge_str = strjust(pad([string(spaces.dimensions.charges) +% string(spaces.dimensions.degeneracies)]), 'center'); +% fprintf('\t%s:\t%s\n', title_str(1), string(spaces.dual)); +% fprintf('\t%s:\t%s\n', title_str(2), join(charge_str(1, :), char(9))); +% fprintf('\t%s:\t%s\n', title_str(3), join(charge_str(2, :), char(9))); return end @@ -405,6 +446,30 @@ function disp(spaces) fprintf('\n'); end end + + function s = string(spaces) + title_str = strjust(pad(["dual", "charges", "degeneracies"]), 'right'); + charge_str = strjust(pad([string(spaces.dimensions.charges) + string(spaces.dimensions.degeneracies)]), 'center'); + s = sprintf(... + '%s GradedSpace of dimension %d:\n\t%s:\t%s\n\t%s:\t%s\n\t%s:\t%s\n', ... + class(spaces.dimensions.charges), dims(spaces), ... + title_str(1), string(spaces.dual), ... + title_str(2), join(charge_str(1, :), char(9)), ... + title_str(3), join(charge_str(2, :), char(9))); + end + + function complexspaces = ComplexSpace(gradedspaces) + d = num2cell(dims(gradedspaces)); + isdual = num2cell([gradedspaces.dual]); + args = [d; isdual]; + complexspaces = ComplexSpace(args{:}); + end + + function cartesianspaces = CartesianSpace(gradedspaces) + d = num2cell(dims(gradedspaces)); + cartesianspaces = CartesianSpace(d{:}); + end end diff --git a/src/utility/Options.m b/src/utility/Options.m index 2a65273..753dc5c 100644 --- a/src/utility/Options.m +++ b/src/utility/Options.m @@ -13,6 +13,17 @@ if isempty(isenabled), isenabled = true; end bool = isenabled; end + + function bool = Debug(bool) + persistent dodebug + if nargin > 0 + dodebug = bool; + return + end + + if isempty(dodebug), dodebug = false; end + bool = dodebug; + end end end diff --git a/src/utility/linalg/contract.m b/src/utility/linalg/contract.m index 7deb1e0..71a227d 100644 --- a/src/utility/linalg/contract.m +++ b/src/utility/linalg/contract.m @@ -81,6 +81,11 @@ % contract last pair [dimA, dimB] = contractinds(ia, ib); + +if Options.Debug + contractcheck(A, ia, ca, B, ib, cb); +end + C = tensorprod(A, B, dimA, dimB, ca, cb, 'NumDimensionsA', length(ia)); ia(dimA) = []; ib(dimB) = []; ic = [ia ib]; @@ -135,6 +140,11 @@ [A, ia, ca] = contracttree(tensors, indices, conjlist, tree{1}); [B, ib, cb] = contracttree(tensors, indices, conjlist, tree{2}); [dimA, dimB] = contractinds(ia, ib); + +if Options.Debug + contractcheck(A, ia, ca, B, ib, cb); +end + C = tensorprod(A, B, dimA, dimB, ca, cb, 'NumDimensionsA', length(ia)); ia(dimA) = []; @@ -142,3 +152,20 @@ ic = [ia ib]; cc = false; end + +function contractcheck(A, ia, ca, B, ib, cb) + +Aspaces = space(A); +if ca, Aspaces = conj(Aspaces); end +Bspaces = space(B); +if cb, Bspaces = conj(Bspaces); end + +[dimA, dimB] = contractinds(ia, ib); + +for i = 1:length(dimA) + assert(Aspaces(dimA(i)) == conj(Bspaces(dimB(i))), 'tensors:SpaceMismatch', ... + 'Invalid index %d:\n\t%s\n\tis incompatible with\n\t%s', ... + ia(dimA(i)), string(Aspaces(dimA(i))), string(Bspaces(dimB(i)))); +end + +end diff --git a/src/utility/memsize.m b/src/utility/memsize.m index ebdc278..c266585 100644 --- a/src/utility/memsize.m +++ b/src/utility/memsize.m @@ -13,7 +13,7 @@ props = properties(in); bytes = 0; for ii = 1:length(props) - bytes = bytes + memsize(in.(thisprop), 'B'); + bytes = bytes + memsize({in.(props{ii})}, 'B'); end end diff --git a/test/TestTensor.m b/test/TestTensor.m index 245716c..21180bb 100644 --- a/test/TestTensor.m +++ b/test/TestTensor.m @@ -138,6 +138,7 @@ function permute_via_conversion(tc, spaces) end function multiplication_via_conversion(tc, spaces) + tc.assumeTrue(istwistless(braidingstyle(spaces))); t1 = Tensor.randnc(spaces(1), spaces(2)); t2 = Tensor.randnc(spaces(2), spaces(3)); @@ -166,12 +167,34 @@ function multiplication_via_conversion(tc, spaces) end function tensorprod_via_conversion(tc, spaces) + tc.assumeTrue(istwistless(braidingstyle(spaces))); t1 = Tensor.randnc([], spaces(1:2)); t2 = Tensor.randnc(spaces(1:2), []); tc.assertTrue(isapprox(tensorprod(t1, t2, [1 2], [2 1]), ... tensorprod(double(t1), double(t2), [1 2], [2 1]), ... 'AbsTol', tc.tol, 'RelTol', tc.tol)); + + A = Tensor.randnc(spaces(1), spaces(1:2)); + C = Tensor.randnc(spaces(1), spaces(1)); + AC = contract(A, [-1 -2 1], C, [1 -3]); + + tc.assertTrue(isapprox(double(AC), contract(double(A), [-1 -2 1], double(C), [1 -3]))); + end + + function contract_order(tc, spaces) + A = Tensor.randnc(spaces(1:2), spaces(1)'); + r = Tensor.randnc(spaces(1)', spaces(1)'); + + args = {A, r, conj(A); [-1 2 1], [1 3], [-2 2 3]}; + + r1 = contract(args{:}, 'Rank', rank(r)); + for p = perms(1:3)' + args2 = args(:, p); + r2 = contract(args2{:}, 'Rank', rank(r)); + tc.assertTrue(isapprox(r1, r2), ... + 'Contraction order should leave result invariant.'); + end end function orthogonalize(tc, spaces) @@ -182,7 +205,7 @@ function orthogonalize(tc, spaces) p2 = [1 5]; tc.assumeTrue(spaces(3) * spaces(4) * spaces(2) >= spaces(1)' * spaces(5)') - for alg = ["qr", "qrpos", "polar", "svd" "ql" "qlpos"] + for alg = ["qr", "qrpos", "polar", "svd", "ql", "qlpos"] [Q, R] = leftorth(t, p1, p2, alg); assertTrue(tc, ... @@ -212,7 +235,7 @@ function orthogonalize(tc, spaces) tc.assumeTrue(spaces(3) * spaces(4) <= spaces(1)' * spaces(2)' * spaces(5)'); p1 = [3 4]; p2 = [2 1 5]; - for alg = ["lq", "lqpos", "polar", "svd" "rq" "rqpos"] + for alg = ["lq", "lqpos", "polar", "svd", "rq", "rqpos"] [L, Q] = rightorth(t, p1, p2, alg); assertTrue(tc, ... @@ -291,4 +314,3 @@ function eigenvalues(tc, spaces) end end end -