@ -175,7 +174,7 @@ The defining axioms of the Kauffman polynomial are the following, given a link d
We will later be seeing that the Kauffman polynomial can be defined in a more explicit way, using a recursive definition that is the one we will be using to derive our algorithm.
== Computational Knot Theory
= Computational Knot Theory
The first problem in computational knot theory is to find a good representation for knots and links. There are various common representations in the literature, such as:
@ -193,7 +192,7 @@ The first problem in computational knot theory is to find a good representation
There are also other codes like *Braid representations* and *DT (Dowker-Thistlethwaite) Codes* we will not be using in this project.
=== PD codes
== PD codes
The main source for the section is the article on #link("https://knotinfo.math.indiana.edu/descriptions/pd_notation.html")[PD notation from KnotInfo]. The PD code for a link is generated by labelling each arc of the link with a number. Then we choose a starting point for each component and process each component in order.
@ -286,7 +285,7 @@ Every directed crossing appears only once as an over-crossing so this algorithm
`[(6,1,7,2), (8,3,5,4), (2,5,3,6), (4,7,1,8)]`
}
=== SG Codes
== SG Codes
Gauss originally developed a notation called *Gauss codes* based on labelling each crossing of a knot with a number and keeping track of when we walk an over-crossing or an under-crossing using a sign. This produces a list of numbers where each number appears exactly twice with different signs. This has a few problems like the fact that this doesn't distinguish a knot vs its mirror.
@ -316,7 +315,7 @@ Output: List<List<(Int, Int)>>
Converting one code to the other is not too much work as one just need to first do a labelling step to convert crossing labels and then convert the over/under-strand and left/right-handedness relations between the two notations.
=== Comparison of PD and SG codes
== Comparison of PD and SG codes
We will now see how SG codes are better suited for the manipulations (switching and splicing) we need to do on a link. Indeed the #link("https://katlas.org/wiki/Setup")[KnotTheory Mathematica package] when computing the Kauffman polynomial converts the input link from *PD code* into *SG code* as this is better suited for doing manipulations directly on the crossings.
@ -344,7 +343,7 @@ Finally SG codes are also more "space efficient". Let $N$ be a number of bits to
- PD codes use $approx 4n times N$ bits of information, four numbers for each item.
- SG codes use $approx 2n times (N + 2) + k times ceil(log_2(n))$ bits of information. Each crossing appears twice and we store its id and over/under and handedness information, we also need to store the structure of the list with $k times ceil(log_2(n))$ more bits.
- SG codes use $approx 2n times (N + 1.5) + k times ceil(log_2(n))$ bits of information. Each crossing appears twice and we store its id and over/under and handedness (each pair has the same handedness so we can store it only once per crossing) information, we also need to store the structure of the list with $k times ceil(log_2(n))$ more bits.
So PD codes are simpler and compact to store (and generate from a diagram) but SG codes are more space efficient and easy to manipulate.
@ -586,11 +585,11 @@ So PD codes are simpler and compact to store (and generate from a diagram) but S
#pagebreak()
=== Link reconstruction from code
== Link reconstruction from code
We briefly mention that reconstructing a link from a PD or SG code is not trivial and there are various approaches used by various softwares that can be used for this task.
==== Linear Integer Programming
=== Linear Integer Programming
For example the #link("https://doc.sagemath.org/html/en/reference/knots/sage/knots/link.html")[KnotTheory package in Sage] has a #link("https://doc.sagemath.org/html/en/reference/knots/sage/knots/link.html#sage.knots.link.Link.plot")[`Link.plot()`] method that thats a link that can be constructed using a PD code and then plots it as follows
#[
@ -615,13 +614,13 @@ For example the #link("https://doc.sagemath.org/html/en/reference/knots/sage/kno
Sage internally uses a #link("https://github.com/sagemath/sage/blob/e0cf1e41d419feb9ddbc0e3c54823928a01587dc/src/sage/knots/link.py#L3639")[_mixed integer linear programming_ (MILP)] solver to generate a knot diagram from a PD code. Another library called #link("https://github.com/3-manifolds/Spherogram/blob/725086a1d8c5d1381ff6a70315047efd8e0dac3f/spherogram_src/links/orthogonal.py")[Spherogram] instead uses _network flows_. The problem here is to find an orthogonal presentation for the link with the _minimum number of left and right bends_
#footnote[#link("https://www.sciencedirect.com/science/article/pii/S0096300305002778")[A better heuristic for area-compaction of orthogonal representations]].
==== Planar Graph Embeddings
=== Planar Graph Embeddings
Another approach used by #link("https://knotfol.io/")[KnotFolio] is based on #link("https://en.wikipedia.org/wiki/Tutte_embedding")[Tutte embeddings]. A *Tutte embedding or barycentric embedding* of a simple, 3-vertex-connected, planar graph is a crossing-free straight-line embedding with the properties that the outer face is a convex polygon and that each interior vertex is at the average (or barycenter) of its neighbors' positions.
This condition that every point is the average of its neighbors can be easily expressed as a system of linear equations where some points on a chosen outer face have been fixed. When the graph is planar and 3-vertex-connected the linear system is non degenerate and has a unique solution.
== Algorithm for computing the Kauffman Polynomial
= Computing the Kauffman Polynomial
Let's now recap the main formal algorithm for computing the Kauffman polynomial.
@ -650,6 +649,8 @@ Let's now recap the main formal algorithm for computing the Kauffman polynomial.
)
]
=== First Approach
Let now $K$ be an oriented link with $n$ components so $K = K_1 union dotss union K_n$.
#definition[
@ -669,7 +670,7 @@ Let now $K$ be an oriented link with $n$ components so $K = K_1 union dotss unio
$
]
*Algorithm.* Here follows the algorithm that computes $kL_(K)(a,z)$.
*Closed Form Algorithm.* Here follows the algorithm that computes $kL_(K)(a,z)$.
1. If $K = hat(K)$ is a _standard unknot_ then $kL_K (a, z) colon.eq a^w(K)$
@ -708,13 +709,48 @@ Let now $K$ be an oriented link with $n$ components so $K = K_1 union dotss unio
$
]
=== Second Approach
After a first implementation of the algorithm we noticed that it is not very efficient to compute using this closed form algorithm. In this case the naive implementation of the algorithm of just applying the rules recursively is more efficient and we only use the closed form for handling the base cases and the case of multiple components.
The final algorithm we implemented is the following:
*Algorithm*: We apply the following rules recursively
1. If $K$ is a standard unlink then return $a^w(K)$.
2. If $K$ has groups of components that overlap each other then we can apply the following rules, let $K = K_1 union.sq dotss union.sq K_n$ be these groups of components one overling the other and let $d = (a + a^(-1)) slash z - 1$ then
3. If $K$ is a linked component then we find $hat(K)$, the standard unlink for $K$, and pick the first available crossing $c$ to switch, then we have the following cases based on the crossing over/under type and handedness that can be reduced to the following two cases
$
kL(skein.over) = z (kL(skein.h) + kL(skein.v)) - kL(skein.under)
$
$
kL(skein.under) = z (kL(skein.h) + kL(skein.v)) - kL(skein.over)
$
actually due to the symmetry of the Kauffman polynomial we can just use the same rule for both cases, so we can write
$
kL(K) = z (kL(E_c K) + kL(e_c K)) - kL(S_c K)
$
where $E_c K$ and $e_c K$ have one less crossing and $S_c K$ has the same crossings as $K$ but is one step closer to the standard unlink.
#pagebreak()
== Python Implementation
= Python Implementation
The approach has been a mix of bottom-up and top-down. First we defined a couple of classed `SignedGaussCode` and `PDCode` to work with these codes and easily convert between each other.
=== SG Codes
== SG Codes
We are now going to walk thorough the class that lets use work nicely with *SG codes*. The the classes we are going to use are all _frozen data-classes_ to ensure immutability and enforce a more functional programming style.
@ -725,39 +761,60 @@ We are now going to walk thorough the class that lets use work nicely with *SG c
Let's now walk through the most important methods of this class.
=== Writhe
One of the first important things we need is to compute the *writhe* $w(K)$ of a link, this can easily be done with signed gauss codes as its a list of tuples where the second entry is the crossing sign. Let $L$ be an oriented link with components $C_1, ..., C_k$ each with crossings $c_(i, j)$ with $i = 1, ..., k$ and $j = 1, ..., |C_i|$.
@ -788,7 +845,7 @@ Let's notice that here each crossing appears twice, once as over-crossing and on
)
]
==== Standard Unknot
=== Standard Unknot
The next building block for computing the Kauffman polynomial is detecting and computing the *standard unknot or unlink*. Formally this is done by taking the _planar shadow_ and a directed starting point on each component of the link. Then we can walk along the shadow of each component in order and make each crossing an over-crossing when passing on it on the first time.
if crossing.id not in visited_crossings and crossing.is_under():
return crossing.id
visited_crossings.add(crossing.id)
return False
```
=== Crossing Splices
The splicing code is more involved due to the number of cases to analyze, let's first see formally what we need to do.
@ -934,11 +1004,42 @@ where $C_k = (c_k, s_k)$ as a pair of crossing *id* and *sign*. To apply the spl
)
$
Handling the reversed part is actually more involved than just reversing the list of crossings. We also need to correct all the signs of the crossings to account for the new orientation. To do this we need to flip the crossing sign of all crossing ids that occur in this part we are reversing. Notice this can end up even alter signs of crossings in other components so we need to be careful. The code that handles with reversing is the following
Handling the reversed part is actually more involved than just reversing the list of crossings. We also need to correct all the signs of the crossings to account for the new orientation as showed in @splice-signs-problem.
The code for the *vertical splices* is omitted as the cases are the same as for the horizontal splices with a minor change. All the cases for vertical splices are shown below in the following diagram and we can see that the output lists are the same as in the previous splice case just switched based on the crossing sign.
To do this we need to flip the crossing sign of all crossing ids that occur in this part we are reversing. Notice this can end up even alter signs of crossings in other components so we need to be careful. Let's see for example the code for the negative crossing horizontal splice case
```py
over_crossing_ids = set(c.id for c in l2 if c.is_over())
under_crossing_ids = set(c.id for c in l2 if c.is_under())
c.flip_handedness() if c.id in non_self_crossing_ids else c
for c in strand
]
return SGCode([
*(
update_signs(component)
for i, component in enumerate(self.components)
if i != component_index
),
[
*update_signs(l1),
*update_signs(reversed(l2)),
*update_signs(l3),
],
])
```
The code for the *vertical splices* is omitted as the cases are the same as for the horizontal splices just flipped with respect to handedness. All the cases for vertical splices are shown below in the following diagram and we can see that the output lists are the same as in the previous splice case just switched based on the crossing sign.
#figure(
image("assets/operations-splice-v-cases.png"),
@ -949,88 +1050,469 @@ So the final code is just a conversion of all this cases to list slicing and re-
#pagebreak()
= Appendix
== Code for computing the Kauffman polynomial
== Enhancing SGCode rejoining code
The initial code to do horizontal and vertical splices looked something like the following
The final code for computing the Kauffman polynomial is the following, the main idea of the algorithm was already described previously. The code is implemented in a functional style and uses the `SGCode` class defined above to represent the links. At the top level we define globals variables using the `sympy` library for working with polynomials.
for k, component_ids in enumerate(component_groups):
new_link = link.sublink(component_ids)
if k > 0:
result *= d
result *= kauffman_polynomial(new_link)
return result
```
This is actually wrong as this is not correcting the signs of all crossing from the reversed strand. The corrected code for the second case is not very readable
Given this we can also write a function for computing the normalized Kauffman polynomial $F_K(a, z) = a^(-w(K)) L_K(a, z)$ as follows
```py
over_crossing_ids = set(c.id for c in l2 if c.is_over())
under_crossing_ids = set(c.id for c in l2 if c.is_under())
return (a ** (-link.writhe())) * kauffman_polynomial(link)
```
update_signs = lambda strand: [
c.flip_handedness() if c.id in non_self_crossing_ids else c
for c in strand
]
#pagebreak()
=== Debugging and Optimizations
The code (omitted from above, see the github repository for the full code) has actually been instrumented with some helper function to ease debugging and optimizations.
- Use the `cache` python decorator to memoize all calls to `kauffman_polynomial`.
- Another decorator called `log_input_output` helps with debug printing the traces like the following for the Hopf link
# First we convert to minimal rotated form and only then we relabel,
# this ensures a consistent indexing for the cache.
if 'to_minimal' in optimizations:
link = link.to_minimal()
if 'relabel' in optimizations:
link = link.relabel()
result = func(link)
if 'expand' in optimizations:
result = sympy.expand(result)
return result
return wrapper
return decorator
```
This applies `sympy.expand(...)` to keep the resulting polynomial simplified and prepends each call with the two following optimizations
- `to_minimal`: This "rotates" the list of each component in the SG code, brining each list in _minimal lexicographical order_.
- `relabel`: After minimal rotation we apply a relabelling to increase the chances of hitting the cache decorator.
These two optimization decrease the number of total recursive calls by almost an order of magnitude, for example for the case of knot `12n_888` we have the following optimization lattice with number of calls in blue and time of execution in green.
As we can see the most important optimization is the `to_minimal` one that by its own reduces the total number of calls from $29"k"$ to $8"k"$ but relabelling is still able to help a bit more.
All optimizations are enabled by default in the final program and debugging traces are disabled to not impact the speed.
== Experiments
The main achievement of this project was to check all (normalized) Kauffman polynomials present in the #link("https://knotinfo.math.indiana.edu/")[KnotInfo database]. Using the python package `database_knotinfo` we loaded all the data present in the database and wrote a script called `./check_knotinfo.py` with the following options
Calculate the Kauffman polynomial of a knot or link.
All combinations of `skein-generic` typst function
options:
-h, --help show this help message and exit
--knots Include knots from database
--links Include links from database
-c, --count COUNT Number of knots to test per database
```
#[
#set align(center + horizon)
This uses the `multiprocessing` library to parallelize the computation of the Kauffman polynomial for each call of the Kauffman polynomial and checks all the databases of knots and links in a couple of minutes.
We run this on all knots and links in the database and found the following results:
- For *links*: All 4187 include the Kauffman polynomial and are all correctly computed by our algorithm.
- For *knots*: Out of 12965 knots only 2977 knots include the Kauffman polynomial (they are computed only up to 12 crossing even if there are up to 13 crossings knots in the database) and _we found a mismatch in the computation of the polynomial of $10_125$_.
=== The knot $10_125$
#for kind in ("over", "under") {
[*Kind:* #raw(repr(kind))]
#align(
center,
grid(
columns: 2,
gutter: 1.5em,
image("assets/10_125-crop.png", height: 5cm),
scale(
x: -100%,
image("assets/10_125-crop.png", height: 5cm),
),
),
)
This knot is _chiral_ meaning that it is not equivalent to its mirror. Our algorithm computes the following polynomial for this knot
&+ space z (-4 a^3 - 8 a - 6 / a - 1 / a^3 + 1 / a^5)
+ 3 a^2 + 7 + 3 / a^2
$
Then we noticed that they are related by the substitution $(a mapsto 1 slash a, z mapsto z)$. The Kauffman polynomial has this property that inverts the $a$ when computing the mirror image of a knot, in this sense it is able to distinguish chiral variants of knots.
So we checked using the implementation of the Kauffman polynomial using one in the `KnotTheory` Mathematica package. We used the PD code present in the KnotInfo database and found that the result of Mathematica and of our algorithm *match*. This makes us believe that there is a problem in KnotInfo and that the column for the PD code or for the Kauffman polynomial are not correctly synchronized.
=== Performance Analysis
We also checked the performance of our algorithm on all knots and links in the database up to 12 crossings (these are the ones with the Kauffman polynomial in the database). The results are shown in the following histograms, each bar counts the number of knots that took the amount of time in the relative bin.
let steps = range(13).map(x => (strfmt("{:.1}s", float(x)), x + 0.0, x + 1.0))
let counts = (:)
for t in times-knots {
let (step, _, _) = steps.find(interval => {
let (_, a, b) = interval
return a <= t and t <= b
})
})
).flatten()
)
}
]
if not step in counts {
counts.insert(step, 0)
}
counts.insert(step, counts.at(step) + 1)
}
let ys = counts.values()
let xs = range(ys.len())
lq.diagram(
yaxis: (
ticks: range(0, 700, step: 100),
subticks: none,
),
xaxis: (
ticks: steps
.map(((l, _, _)) => l)
.map(scale.with(80%))
.map(rotate.with(-45deg, reflow: true))
.map(align.with(right))
.enumerate(),
subticks: none,
),
legend: (position: left + top),
margin: 5%,
lq.bar(
xs,
ys,
offset: 0.5,
),
)
},
[
#set text(size: 9pt)
Histogram of knots times
]
),
grid(
rows: 2,
gutter: 1em,
{
import "@preview/lilaq:0.3.0" as lq
import "@preview/oxifmt:0.2.1": strfmt
let max-value = 0
for t in times-links {
if t > max-value {
max-value = t
}
}
let steps = range(10).map(x => (strfmt("{:.1}s", float(x)), x + 0.0, x + 1.0))
let counts = (:)
for t in times-links {
let (step, _, _) = steps.find(interval => {
let (_, a, b) = interval
return a <= t and t <= b
})
if not step in counts {
counts.insert(step, 0)
}
counts.insert(step, counts.at(step) + 1)
}
let ys = counts.values()
let xs = range(ys.len())
lq.diagram(
yaxis: (
ticks: range(0, 1700, step: 200),
subticks: none,
),
xaxis: (
ticks: steps
.map(((l, _, _)) => l)
.map(scale.with(80%))
.map(rotate.with(-45deg, reflow: true))
.map(align.with(right))
.enumerate(),
subticks: none,
),
legend: (position: left + top),
margin: 5%,
lq.bar(
xs,
ys,
offset: 0.5,
),
)
},
[
#set text(size: 9pt)
Histogram of links times
]
),
),
)
We can see that the times are just about of a couple of seconds per knot or link. This could be improved by using more sophisticated caching techniques but we did not implement this as the performance was already good enough for our purposes.
= Conclusion
In this project we implemented the Kauffman polynomial from scratch in Python using signed Gauss codes. We compared the results of our algorithm with the ones present in the KnotInfo database and found an error in the computation of the polynomial of the knot $10_125$. We believe that the error is due to a mismatch between the PD code stored in the database with the corresponding Kauffman polynomial.
// #pagebreak()
// = Appendix
// == Enhancing SGCode rejoining code
// The initial code to do horizontal and vertical splices looked something like the following
// ```py
// ...
// if handedness == HANDED_LEFT:
// return SGCode([
// *(component[:]
// for i, component in enumerate(self.components)
// if i != component_index),
// [
// *l1,
// *l3,
// ],
// l2,
// ])
// else:
// return SGCode([
// *(component[:]
// for i, component in enumerate(self.components)
// if i != component_index),
// [
// *l1,
// *l2[::-1],
// *l3,
// ],
// ])
// ```
// This is actually wrong as this is not correcting the signs of all crossing from the reversed strand. The corrected code for the second case is not very readable
// ```py
// over_crossing_ids = set(c.id for c in l2 if c.is_over())
// under_crossing_ids = set(c.id for c in l2 if c.is_under())
