khoca.py

#!/usr/bin/python3
## @file

#
#    khoca.py --- Part of khoca, a knot homology calculator
#
# Copyright (C) 2018 Lukas Lewark <lukas@lewark.de>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

# The file to be run by the user.
# The main work is done by pui.pyx, which is called by this file.

import sys
sys.path.insert(0, './bin/')

from pui import pCalcSubTangleTree, pCalculateHomology, pGlueReduced, pCopyHomology, pFirstFreeIdx, pSaveComplexToFile, pDeleteComplex, pLoadComplexFromFile, pSum, pDual, pDeleteNonIsos, pResetSimplificationsCounter, pReducify, pIniStack, pSetRoot, pPrintHomology, pPrintCompileInfo, pCalculateEquiHomology
import random
import re
import math

sys.path.append('.')
import KrasnerGaussToMyPDLib, pseudoBraidToKrasnerGaussLib, BraidToMyPD


debugging = False
NUM_THREADS = 1 if debugging else 12



## Converts a link diagram given in the Planar-diagram notation to the mypd format.
# @param pd A list of lists of four integers each, such as [[0,1,2,3],[1,4,5,2],[4,0,3,5] for the positive trefoil.
# @return That link diagram in mypd-notation, type "0" for positive, type "1" for negative crossings.
def pd_to_mypd(pd):
    result = list()
    for i in pd:
        if (i[3] == ((i[1] + 1) % (2 * len(pd)))):
            result.append([2,i[0] - 1,i[1] - 1,i[2] - 1,i[3] - 1])
        else:
            result.append([3,i[1] - 1,i[2] - 1,i[3] - 1,i[0] - 1])
    return result

def generate_binary_tree(length, s = 0):
    if (length == 1):
        return s
    return [generate_binary_tree(length / 2, s), generate_binary_tree((length / 2) + (length % 2), s + (length / 2))]

def generate_standard_tree(l):
    if (len(l) == 1):
        return l[0]
    return [generate_standard_tree(l[:-1]), l[len(l) - 1]]

def isBetterSignature(s1, s2):
    s1Copy = list(s1)
    s2Copy = list(s2)
    s1Copy.sort(reverse = True)
    s2Copy.sort(reverse = True)
    if (s1Copy < s2Copy): return True
    if (s1Copy > s2Copy): return False
    return (s1 > s2)

# Using a randomised algorithm, produces a linear tree such that intermediate tangles have small girth
def generate_girth_optimised_tree(pd):
    pd = [[i] + pd[i][1:] for i in range(len(pd))]
    iterations = 1000

    bestOrder = []
    bestSignature = []
    girth = 0
    for k in range(iterations):
        copyPd = list(pd)
        openEnds = set()
        order = []
        signature = []
        while len(copyPd) > 0:
            girthChange = [4 - 2 * len(openEnds & set(i[1:])) for i in copyPd]    
            bestGirthChange = min(girthChange)
            bestCrossings = [i for i in range(len(girthChange)) if girthChange[i] == bestGirthChange]
            newCrossing = bestCrossings[0] if debugging else random.choice(bestCrossings)
            order.append(copyPd[newCrossing][0])
            openEnds = openEnds ^ set(copyPd.pop(newCrossing)[1:])
            signature.append(len(openEnds))
        if (bestOrder == []) or (isBetterSignature(signature, bestSignature)):
            bestOrder = order
            bestSignature = signature
            girth = max(bestSignature)
    return (girth, generate_standard_tree(bestOrder))

def eliminatedoubles(l):
    i = 0
    while (i < len(l)):
        j = i + 1
        found = False
        while (j < len(l)):
            if l[i] == l[j]:
                l.pop(j)
                l.pop(i)
                found = True
                break
            j = j + 1
        if not found:
            i = i + 1
        
## Stores the information which elementary tangles to load and how and in which order to glue them.
# In particular, all index shifts have been done, so that calcSubTangleTree may
# call functions in pythonInterface.cpp without further ado.
# Details: This represents a directed tree with vertices numbered 0,1,....
# Each vertex of number i has at most two children, whose number is saved in
# sons[i] and daughers[i] (-1 meaning no child). types saves whether the tangle at the vertex
# is elementary and should be loaded from a file (then types[i] is 0,1,2,3, see function loadType in pui.pyx),
# or the tangle is obtained by gluing son and daughter (then types[i] is -1).
# boundaryPoints[i] is an orderd list of the tangles' boundary points (in global numbering).
# when gluing e.g. a daughter and son with boundaryPoints 2,3,0,1 and 4,6,3,2, respectively, one
# first takes disjoint union (giving a tangle with boundaryPoints 2,3,0,1,4,6,3,2), and then
# glues the matching boundaryPoints; which match which is saved in daughterGluePoints and sonGluePoints,
# which are 0,0 and 7,5 here (because index shifts are already contained).
class CSchedule:
    def initialise(self, mypd):
        nodes = 2 * len(mypd) - 1
        self.shift = 0
        self.daughters = [-1] * nodes
        self.sons = [-1] * nodes
        self.types = [-1] * nodes
        self.selfGlue = [list() for i in range(nodes)]
        self.boundaryPoints = [list() for i in range(nodes)]
        self.sonGluePoints = [list() for i in range(nodes)]
        self.daughterGluePoints = [list() for i in range(nodes)]
        return nodes + 1

    def output(self):
        sys.stderr.write(str(self.daughters) + str(self.sons) + str(self.types) + str(self.boundaryPoints) + str(self.daughterGluePoints) + str(self.sonGluePoints) + str(self.selfGlue) + "\n")


    def getFromMypdAndTree(self, mypd, tree, idx):
        if (type(tree) == type(int())):
            self.types[idx] = mypd[tree][0]
            self.boundaryPoints[idx] = mypd[tree][1:]
            unitedBoundaryPoints = list(self.boundaryPoints[idx])
            i = 0
            while (i < len(unitedBoundaryPoints)):
                j = i + 1
                found = False
                while (j < len(unitedBoundaryPoints)):
                    if unitedBoundaryPoints[i] == unitedBoundaryPoints[j]:
                        self.selfGlue[idx].append([i,j])
                        unitedBoundaryPoints.pop(j)
                        unitedBoundaryPoints.pop(i)
                        found = True
                        break
                    j = j + 1
                if not found:
                    i = i + 1
            return idx + 1
        else:    
            self.daughters[idx] = idx + 1
            self.sons[idx] = self.getFromMypdAndTree(mypd, tree[0], idx + 1)
            result = self.getFromMypdAndTree(mypd, tree[1], self.sons[idx])
            beforeDeletion = (self.boundaryPoints[self.daughters[idx]] + self.boundaryPoints[self.sons[idx]])
            self.boundaryPoints[idx] = [ x for x in beforeDeletion if beforeDeletion.count(x) == 1]
            shift = len(self.boundaryPoints[self.daughters[idx]])
            daughterBP = list(self.boundaryPoints[self.daughters[idx]])
            eliminatedoubles(daughterBP)
            sonBP = list(self.boundaryPoints[self.sons[idx]])
            eliminatedoubles(sonBP)
            unitedBoundaryPoints = list(daughterBP)
            unitedBoundaryPoints.extend(sonBP)
            i = 0
            while (i < len(unitedBoundaryPoints)):
                j = i + 1
                found = False
                while (j < len(unitedBoundaryPoints)):
                    if unitedBoundaryPoints[i] == unitedBoundaryPoints[j]:
                        self.daughterGluePoints[idx].append(i)
                        self.sonGluePoints[idx].append(j)
                        unitedBoundaryPoints.pop(j)
                        unitedBoundaryPoints.pop(i)
                        found = True
                        break
                    j = j + 1
                if not found:
                    i = i + 1
            return result

def reducePD(mypd):
    x = [item for subl in mypd for item in subl[1:]]
    j = x.index(max(x))
    mypd[int(j / 4)][1 + j % 4] +=1 
    return mypd

def calcIt(mypd, shift):
    global verbose, progress
    (girth, tree) = generate_girth_optimised_tree(mypd)
    pIniStack(stackMod, stackFrobenius, N, girth, verbose)
    pSetRoot(stackRoot)

    s = CSchedule()
    result = s.initialise(mypd)
    s.getFromMypdAndTree(mypd, tree, 0)
    s.shift = shift

    pCalcSubTangleTree(s, NUM_THREADS, progress)
    return result

def getInts(s, signed):
    return [int(j) for j in re.findall( ("-?" if signed else "") + "[0-9]+", s[3:])]

## @var mypd
# The internal format for a closed object (as a link), glued together from smaller tangles (as tangles).
#
# mypd is a list of tangles; an object is a list of integers.
# The first of that integers is the type of the object (e.g. a positive or a
# negative crossing), all following integers are the numbers of the boundary
# points of that object, the order of which matters.
# %Boundary point numbers need not be consecutive; each boundary point number
# must appear exactly twice in the whole mypd.

def isAllowedStackMod(a):
    if not a.isdigit():
        return False
    x = int(a)
    if (x <= 3):
        return True
    for i in range(2, int(math.sqrt(x) + 2)):
        if (x % i == 0):
            return False
    return True

def run_commandline(argv, printCommand, progress): 
        NUM_ARGUMENTS = 3
        global stackMod, stackFrobenius, N, stackRoot

        # Parsing the arguments
        if (len(argv) <= NUM_ARGUMENTS):
            printCommand(HELP_TEXT)
            return 1
        
        if isAllowedStackMod(argv[1]):
            stackMod = int(argv[1])
        else:
            printCommand("First argument must be 0, 1 or a prime.")
            return 1
        
        N = 0
        if argv[2][0] == "e":
            N = int(argv[2][1:])
            l = []
            if (N < 2):
                printCommand("e must be followed by an integer that is at least 2.")
                return 1
        else:
            l = [int(i) for i in re.findall("-?[0-9]+", argv[2])]
            if (len(l) == 0):
                printCommand("Second argument must be a list of at least one integer, or e (for equivariant).")
                return 1
        stackFrobenius = l
        equivariant = (len(stackFrobenius) == 0)
        
        if argv[3].isdigit() or ((argv[3][0] in ["+","-"]) and (argv[3][1:].isdigit())):
            stackRoot = int(argv[3])
        else:
            printCommand("Third argument must be a signed integer.")
            return 1
        
        # doing the work
        shift = 0
        idxTranslator = [0]
        for i in argv[(NUM_ARGUMENTS + 1):]:
            if i[:3].capitalize() == "Nat":
                l = getInts(i, False)
                if (len(l) % 5 != 0) or (len(l) == 0):
                    printCommand("The native code (\"" + i[3:] +  "\") following the command \"nat\" must be a non-empty list of integers whose length is divisible by five.")
                    sys.exit()
        
                mypd = [list(j) for j in zip(*[iter(l)]*5)]
                mypd = reducePD(mypd)
                shift += calcIt(mypd, shift)
                idxTranslator += [shift]
            if i[:7].capitalize() == "Special":
                l = eval(i[7:])
                shift += calcIt(l, shift)
                idxTranslator += [shift]
            elif i[:5].capitalize() == "Gauss":
                l = getInts(i, True)
                mypd = reducePD(KrasnerGaussToMyPDLib.KrasnerGaussToMyPDmain(l))
                shift += calcIt(mypd, shift)
                idxTranslator += [shift]
            elif i[:11].capitalize() == "PseudoBraid":
                m = getInts(i, True)
                l = pseudoBraidToKrasnerGaussLib.pseudoBraidToKrasnerGaussMain(m)
                mypd = reducePD(KrasnerGaussToMyPDLib.KrasnerGaussToMyPDmain(l))
                shift += calcIt(mypd, shift)
                idxTranslator += [shift]
            elif i[:5].capitalize() == "Braid":
                mypd = reducePD(BraidToMyPD.BraidToMyPD(i[5:]))
                shift += calcIt(mypd, shift)
                idxTranslator += [shift]
            elif i[:2].capitalize() == "Pd":
                l = getInts(i, False)
                l = [x - 1 for x in l]
                mypd = reducePD(pd_to_mypd([list(j) for j in zip(*[iter(l)]*4)]))
                print(mypd)
                shift += calcIt(mypd, shift)
                idxTranslator += [shift]
            elif i[:10].capitalize() == "Calcnonred":
                l = getInts(i, False)
                if (len(l) != 1):
                    printCommand("\"Calcnonred\" must be followed by exactly one unsigned integer.")
                    sys.exit()
                idx = idxTranslator[l[0]]
                firstFreeIdx = shift
                shift += 1
                idxTranslator[-1] += 1
                pResetSimplificationsCounter(idx)
                pCopyHomology(idx, firstFreeIdx)
                printCommand("Result:")
                printCommand("Unreduced Homology:")
                pGlueReduced(firstFreeIdx, 1, NUM_THREADS, progress)
                pCalculateHomology(firstFreeIdx, False, NUM_THREADS, printCommand, equivariant, progress)
                pDeleteComplex(firstFreeIdx)
            elif i[:4].capitalize() == "Calc":
                nice = (i[4:8].capitalize() == "Nice")
                l = getInts(i, False)
                if (len(l) != 1):
                    printCommand("\"Calc\" must be followed by exactly one unsigned integer.")
                    sys.exit()
                idx = idxTranslator[l[0]]
                firstFreeIdx = shift
                secondFreeIdx = shift + 1
                shift += 2
                idxTranslator[-1] += 2
                pResetSimplificationsCounter(idx)
                pCopyHomology(idx, firstFreeIdx)
                pCopyHomology(idx, secondFreeIdx)
                printCommand("Result:")
                printCommand("Reduced Homology:")
                pReducify(secondFreeIdx)
                pCalculateHomology(secondFreeIdx, nice, NUM_THREADS, printCommand, equivariant, progress, True)
                printCommand("Unreduced Homology:")
                pGlueReduced(firstFreeIdx, 1, NUM_THREADS, progress)
                pCalculateHomology(firstFreeIdx, nice, NUM_THREADS, printCommand, equivariant, progress, True)
                pDeleteComplex(firstFreeIdx)
                pDeleteComplex(secondFreeIdx)
            elif i[:4].capitalize() == "Save":
                l = re.match("[0-9]+", i[4:])
                if (l == None):
                    printCommand("\"Save\" must be followed by an unsigned integer and a filename.")
                    sys.exit()
                pSaveComplexToFile(idxTranslator[int(l.group(0))], i[(4 + len(l.group(0))):]) 
            elif i[:4].capitalize() == "Load":
                pLoadComplexFromFile(shift, i[4:])
                shift += 1
                idxTranslator += [shift]
            elif i[:4].capitalize() == "Dual":
                l = re.match("[0-9]+", i[4:])
                if (l == None):
                    printCommand("\"Dual\" must be followed by an unsigned integer.")
                    sys.exit()
                pDual(shift, int(l.group(0)))
                shift += 1
                idxTranslator += [shift]
            elif i[:3].capitalize() == "Sum":
                l = getInts(i, False)
                if (len(l) != 2):
                    printCommand("\"Sum\" must be followed by exactly two unsigned integer.")
                    sys.exit()
                pSum(shift, idxTranslator[l[0]], idxTranslator[l[1]], NUM_THREADS, progress)
                shift += 1
                idxTranslator += [shift]
            else:
                printCommand(str(i) + " is not a valid command.")
                sys.exit()


def parseOptions(argv):
    global verbose, progress
    o = [x for x in argv[1:] if (x[0] == '-')]
    for i in o:
        if (i == '-v'):
            verbose = 1
        elif (i == '-p'):
            progress = 1
        elif (i == '-h'):
            print(HELP_TEXT)
            sys.exit()
        else:
            print("Unknown option " + i + ".")
            sys.exit()
    return [argv[0]] + [x for x in argv[1:] if (x[0] != '-')]

HELP_TEXT = "Expecting at least three arguments:\n(1) The coefficient ring; 0 for integers, 1 for rationals, a prime p for the finite field with p elements,\n(2) the vector [a_0, ... a_{N-1}] defining the Frobenius algebra F[X]/(X^N+a_{N-1}X^{N-1}+...+a_0) or e followed by the number N for equivariant,\n(3) a root of the polynomial given in (2).\nAny following argument is interpreted as command. For example,\n./khoca.py 0 0.0 0 braidaBaB calc0\ncomputes integral Khovanov sl(2)-homology of the figure-eight knot.\nRefer to the README file or to http://lewark.de/lukas/khoca.html for more detailed help."
verbose = 0
progress = 0
a = parseOptions(sys.argv)

if (verbose):
    pPrintCompileInfo()
    sys.stderr.write(("Multithreading with " + str(NUM_THREADS) + " threads.\n") if (NUM_THREADS > 1) else "No multithreading.\n")

run_commandline(a, print, progress)