Implement computation of Reshetikhin--Turaev knot polynomials along with miscellaneous features for Tangle

setup.py

@@ -205,7 +205,8 @@ def run(self):
                     'networkx',
                     'packaging',
                     'snappy_manifolds>=1.4',
-                    'knot_floer_homology>=1.2.2']
+                    'knot_floer_homology>=1.2.2',
+                    'opt_einsum>=3.4.0']
 
 setup( name = 'spherogram',
        version = version,
@@ -214,9 +215,11 @@ def run(self):
        dependency_links = [],
        packages = ['spherogram', 'spherogram.links', 'spherogram.links.bands',
                    'spherogram.links.test', 'spherogram.codecs',
-                   'spherogram.dev', 'spherogram.dev.dev_jennet'],
+                   'spherogram.dev', 'spherogram.dev.dev_jennet',
+                   'spherogram.links.reshetikhin_turaev'],
        package_dir = {'spherogram' : 'spherogram_src', 'spherogram.dev':'dev'},
-       package_data = {'spherogram.links'  :  ['doc.pdf']},
+       package_data = {'spherogram.links'  :  ['doc.pdf'],
+                       'spherogram.links.reshetikhin_turaev': ['R_matrices/*/*']},
        ext_modules = ext_modules,
        cmdclass =  {'clean': SpherogramClean,
                     'test': SpherogramTest,

spherogram_src/__init__.py

@@ -1,6 +1,7 @@
 from .presentations import *
 from .links import *
 from .codecs import *
+from .links.reshetikhin_turaev import DictLaurentPolynomial, RMatrix, colored_links_gould_R_matrices, colored_jones_R_matrices, prefactor_colored_jones
 
 # Make the module version number easily accessible.
 from . import version as _version
@@ -23,4 +24,5 @@ def version():
            # from spherogram.links.tangles:
            'Tangle', 'CapTangle', 'CupTangle', 'RationalTangle',
            'ZeroTangle', 'InfinityTangle', 'MinusOneTangle', 'OneTangle', 'IntegerTangle',
-           'IdentityBraid', 'BraidTangle', 'ComponentTangle', 'join_strands']
+           'IdentityBraid', 'BraidTangle', 'ComponentTangle', 'join_strands',
+           'DictLaurentPolynomial', 'RMatrix', 'colored_links_gould_R_matrices', 'colored_jones_R_matrices', 'prefactor_colored_jones']

spherogram_src/links/invariants.py

@@ -88,7 +88,6 @@ def sage_braid_as_int_word(braid):
     see the documentation for the "sage_link" method for details.
 """
 
-
 class Link(links_base.Link):
     __doc__ = links_base.Link.__doc__ + extra_docstring
 
@@ -331,6 +330,157 @@ def alexander_polynomial(self, multivar=True, v='no', method='default',
         if multivar and factored:  # it's easier to view this way
             return p.factor()
         return p
+
+    def colored_links_gould_polynomial(self, 
+                                       n, 
+                                       sage_output=_within_sage,
+                                       sage_polynomials=False, 
+                                       timed=False):
+        """
+        Computes the n-colored Links--Gould polynomial of a link. 
+        The output is an instance of Sage's LaurentPolynomial if in sage,
+        otherwise a DictLaurentPolynomial.
+
+        Colored Links--Gould polynomials are bivariate, for which we default to
+        using DictLaurentPolynomial during the procedure to reduce RAM consumption. 
+
+        >>> Link('3_1').colored_links_gould_polynomial(1)
+        t^2*q^2 - t*q^3 - t*q + 2*q^2 - t^-1*q^3 + 1 - t^-1*q + t^-2*q^2
+        >>> Link('4_1').colored_links_gould_polynomial(1)
+        t^2 - 3*t*q + 2*q^2 - 3*t*q^-1 + 7 - 3*t^-1*q + 2*q^-2 - 3*t^-1*q^-1 + t^-2
+
+        Mirror image is equal to substituting q with q^-1:
+
+        >>> mtref = Link('3_1').mirror()
+        >>> mtref_LG = mtref.colored_links_gould_polynomial(1, sage_output=False)
+        >>> mtref_LG.change_vars({'q': 'q^-1'})
+        t^2*q^2 - t*q^3 - t*q + 2*q^2 - t^-1*q^3 + 1 - t^-1*q + t^-2*q^2
+
+        The colored Links--Gould polynomial specializes to the square of 
+        the Alexander polynomial:
+
+        >>> tref_LG = Link('3_1').colored_links_gould_polynomial(1, sage_output=False)
+        >>> tref_LG.change_vars({'q': '1'})
+        t^2 - 2*t + 3 - 2*t^-1 + t^-2
+        >>> fig8_LG = Link('4_1').colored_links_gould_polynomial(1, sage_output=False)
+        >>> fig8_LG.change_vars({'q': '1'})
+        t^2 - 6*t + 11 - 6*t^-1 + t^-2
+
+        1-colored Links--Gould polynomial is invariant under mutation:
+
+        >>> conway_LG = Link('11n34').colored_links_gould_polynomial(1)
+        >>> KT_LG = Link('11n42').colored_links_gould_polynomial(1)
+        >>> conway_LG == KT_LG
+        True
+
+        A mutation pair with the same 2-colored Links--Gould polynomial:
+
+        >>> K1 = Link('12n364')
+        >>> K2 = Link('12n365').mirror()
+        >>> K1.colored_links_gould_polynomial(2) == K2.colored_links_gould_polynomial(2)
+        True
+
+        Some higher colored values for the trefoil:
+
+        >>> K = Link('3_1')
+        >>> K.colored_links_gould_polynomial(2) # doctest: +NORMALIZE_WHITESPACE
+        -t^2*q^5 + t*q^6 + t^2*q^4 - t*q^5 + t*q^4 - 2*q^5 + t^-1*q^6 + t^2*q^2 
+        - 2*t*q^3 + 2*q^4 - t^-1*q^5 + t^-1*q^4 - t^-2*q^5 - t*q + 2*q^2 - 
+        2*t^-1*q^3 + t^-2*q^4 + 1 - t^-1*q + t^-2*q^2
+
+        >>> K.colored_links_gould_polynomial(3) # doctest: +NORMALIZE_WHITESPACE
+        -t*q^27 + t^2*q^24 + t*q^25 - t^-1*q^27 - t^2*q^22 + t*q^23 + 2*q^24 + 
+        t^-1*q^25 - t^2*q^20 - 2*t*q^21 - 2*q^22 + t^-1*q^23 + t^-2*q^24 + 
+        t^2*q^18 + 2*t*q^19 - 2*q^20 - 2*t^-1*q^21 - t^-2*q^22 - t^2*q^16 + 
+        t*q^17 + 2*q^18 + 2*t^-1*q^19 - t^-2*q^20 - 2*t*q^15 - 2*q^16 + t^-1*q^17
+        + t^-2*q^18 + t^2*q^12 + t*q^13 - 2*t^-1*q^15 - t^-2*q^16 + 2*q^12 + 
+        t^-1*q^13 - 2*t*q^9 + t^-2*q^12 + t^2*q^6 - 2*t^-1*q^9 + 2*q^6 - t*q^3 + 
+        t^-2*q^6 - t^-1*q^3 + 1
+
+        >>> K.colored_links_gould_polynomial(4) # doctest: +NORMALIZE_WHITESPACE
+        t*q^24 - t^2*q^22 - t*q^23 + t^2*q^21 - t*q^22 + t^-1*q^24 + t^2*q^20 -
+        2*q^22 - t^-1*q^23 + 2*t*q^20 + 2*q^21 - t^-1*q^22 - t^2*q^18 - t*q^19 + 
+        2*q^20 - t^-2*q^22 - 2*t*q^18 + 2*t^-1*q^20 + t^-2*q^21 + t^2*q^16 + 
+        t*q^17 - 2*q^18 - t^-1*q^19 + t^-2*q^20 - t^2*q^15 + 2*t*q^16 - 2*t^-1*q^18 
+        - t^2*q^14 + 2*q^16 + t^-1*q^17 - t^-2*q^18 - 2*t*q^14 - 2*q^15 + 2*t^-1*q^16 
+        + t^2*q^12 + 2*t*q^13 - 2*q^14 + t^-2*q^16 - t^2*q^11 + t*q^12 - 2*t^-1*q^14 
+        - t^-2*q^15 + 2*q^12 + 2*t^-1*q^13 - t^-2*q^14 - 2*t*q^10 - 2*q^11 + t^-1*q^12 
+        + t^2*q^8 + t*q^9 + t^-2*q^12 - 2*t^-1*q^10 - t^-2*q^11 + 2*q^8 + t^-1*q^9 - 
+        2*t*q^6 + t^2*q^4 + t^-2*q^8 - 2*t^-1*q^6 + 2*q^4 - t*q^2 + t^-2*q^4 - t^-1*q^2 
+        + 1
+        """
+        from .reshetikhin_turaev import colored_links_gould_R_matrices, DictLaurentPolynomial
+
+        ans = self.min_long_diagram().reshetikhin_turaev_network(colored_links_gould_R_matrices(n, sage_polynomials=sage_polynomials)).evaluate(timed=timed)
+
+        if sage_output:
+            if not sage_polynomials:
+                ans = (ans[0].to_sage(), ans[1])
+        else:
+            if sage_polynomials:
+                ans = (DictLaurentPolynomial.from_sage(ans[0]), ans[1])
+            else:
+                ans = (ans[0].to_checked(), ans[1])
+
+        if timed:
+            return ans
+        else:
+            return ans[0]
+
+    def colored_jones_polynomial(self, 
+                                 n, 
+                                 sage_output=_within_sage,
+                                 sage_polynomials=_within_sage,
+                                 timed=False):
+        """
+        Computes the n-colored Jones polynomial of a link. 
+        The output is an instance of Sage's PuiseuxSeries if in sage,
+        otherwise a DictLaurentPolynomial.
+
+        The 1-colored Jones polynomial is equal to the usual Jones polynomial.
+        Here we follow the ordinary convention of variables for Jones polynomials,
+        instead of the squared q in jones_polynomial()
+
+        Colored Jones polynomials are univariate, for whom sage's PuiseuxSeries
+        has highly optimized multiplications, hence we default to use sage 
+        during the procedure whenever possible.
+
+        >>> Link('3_1').colored_jones_polynomial(1)
+        -q^-4 + q^-3 + q^-1
+        >>> Link('4_1').colored_jones_polynomial(1)
+        q^-2 - q^-1 + 1 - q + q^2
+        >>> Link('L2a1').colored_jones_polynomial(1)
+        q^(-5/2) + q^(-1/2)
+
+        Some values of higher colored Jones polynomials for the trefoil:
+
+        >>> Link('3_1').colored_jones_polynomial(2)
+        q^-11 - q^-10 - q^-9 + q^-8 - q^-7 + q^-5 + q^-2
+        >>> Link('3_1').colored_jones_polynomial(3) # doctest: +NORMALIZE_WHITESPACE
+        -q^-21 + q^-20 + q^-19 - q^-17 + q^-15 - q^-14 - q^-13 + q^-11 - q^-10 + q^-7 
+        + q^-3
+        >>> Link('3_1').colored_jones_polynomial(4) # doctest: +NORMALIZE_WHITESPACE
+        q^-34 - q^-33 - q^-32 + 2*q^-29 - q^-28 + 2*q^-24 - q^-23 - q^-22 + q^-19 - 
+        q^-18 - q^-17 + q^-14 - q^-13 + q^-9 + q^-4
+        """
+        from .reshetikhin_turaev import colored_jones_R_matrices, prefactor_colored_jones, DictLaurentPolynomial
+
+        ans = self.min_long_diagram().reshetikhin_turaev_network(colored_jones_R_matrices(n, sage_polynomials=sage_polynomials)).evaluate(timed=timed)
+        ans = (ans[0] * prefactor_colored_jones(n, self.writhe(), sage_polynomial=sage_polynomials), ans[1])
+
+        if sage_output:
+            if not sage_polynomials:
+                ans = (ans[0].to_sage(), ans[1])
+        else:
+            if sage_polynomials:
+                ans = (DictLaurentPolynomial.from_sage(ans[0]), ans[1])
+            else:
+                ans = (ans[0].to_checked(), ans[1])
+
+        if timed:
+            return ans
+        else:
+            return ans[0]
 
     def knot_floer_homology(self, prime=2, complex=False):
         """

spherogram_src/links/links_base.py

@@ -893,7 +893,8 @@ def reverse_orientation(self, component_index):
         """
         Reverse the orientation of components specified by component_index.
 
-        component_index: either a single index of component or a list of indices of components
+        component_index: either a single index of component 
+                         or a list of indices of components
 
         >>> L = Link([(4, 0, 5, 3), (0, 6, 1, 5), (6, 2, 7, 1), (2, 4, 3, 7)])
         >>> L
@@ -1370,6 +1371,57 @@ def keep(C):
                         A[a] = B[b]
 
         return type(self)(final_crossings, check_planarity=False)
+
+    def long_diagram(self, cut_at=None):
+        """
+        Returns the long diagram of self obtained 
+        by cutting open the strand specified by cut_at.
+
+        cut_at should be a pair of integers (i, j) 
+        representing the j-th strand of the i-th crossing.
+        If not specified, the first strand of the first 
+        crossing will be chosen by default.
+
+        >>> T = Link('4_1').long_diagram()
+        >>> T.PD_code()
+        ((1, 1), [(0, 5, 1, 6), (4, 1, 5, 2), (2, 8, 3, 7), (6, 4, 7, 3)], [0, 8])
+        """
+        from .tangles import Tangle
+        L = self.copy()
+
+        if cut_at is None:
+            strand = L.crossings[0].crossing_strands()[0]
+        else:
+            i, j = cut_at
+            strand = L.crossings[i].crossing_strands()[j]
+
+        open_strands = [strand, strand.opposite()]
+
+        for c in L.crossings:
+            c._clear()
+
+        return Tangle((1,1), L.crossings, open_strands)
+
+    def min_long_diagram(self):
+        """
+        Return the long diagram of self with the minimal contraction width
+        """
+
+        if not self.crossings:
+            return self.long_diagram()
+
+        min_width = None
+
+        for i in range(len(self.crossings)):
+            entry_indices = [3, 0] if self.crossings[i].sign == 1 else [0, 1]
+            for j in entry_indices:
+                diagram = self.long_diagram(cut_at=(i,j))
+                width = diagram.contraction_width()[0]
+                if min_width is None or width < min_width:
+                    ans = diagram
+                    min_width = width
+
+        return ans  
 
     def __len__(self):
         return len(self.crossings)