fermions, convenience and docs

Various bugfixes, additional documentation, convenience functions and fermion functionality.

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

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'

docs/src/man/tensor.rst

@@ -7,7 +7,9 @@ Definitions
 We start with the definition of a tensor. According to the `Link Wikipedia page <https://en.wikipedia.org/wiki/Tensor#Definition>` 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 
 ------------

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