• home
• publications
• talks
• code
• teaching
• other
• blog

correlator.py code

copy this code into a file called correlator.py and run it from the command line using

python correlator -01 12,12 1122

more instructions on usage

#! /usr/bin/python

## correlator version 0.1

## Copyright 2007 Tom Brown

## 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 .

## See http://www.nworbmot.org/code/correlator.html for more usage
## instructions


## permutation: [1,2,3,0] is (0123)

## operator: [['X',2],['X',3],['Y',1],['Y',0]] is tr(XYXY)

## power of N: [0,1,2] is N + 2N^2

## array: [ [power of N, operator], [power of N, operator] ]


'''correlator for tree-level and one-loop correlators and dilatations
see http://www.nworbmot.org/code/correlator.html for more usage
instructions
'''

__version__ = "0.1"
__author__ = "Tom Brown"
__copyright__ = "Copyright 2007 Tom Brown, GNU GPL 3"




import getopt, sys, os



def processOp(operator):

    operator = operator +','

    gibbon = []

    length = 0

    counter = 0

    while(operator.find(',') >=0):
        trace = operator[:operator.find(',')]
        operator = operator[operator.find(',')+1:]

        for x in range(0,len(trace)):
            gibbon.append([trace[x],length+x+1])
            counter = counter +1

        gibbon[length+len(trace)-1][1] = length

        length = length + len(trace)

    sortOp(gibbon)

    return gibbon




def sortOp(operator):

    for x in range(0,len(operator)):
        operator[x].insert(1,x)


    operator.sort()

    for x in range(0,len(operator)):
        for y in range(0,len(operator)):
            if operator[y][1] == operator[x][2]:
                operator[x][2] = y
                break

    for x in range(0,len(operator)):
        operator[x].pop(1)


def sortOpArray(array):

    for thingie in array:
        sortOp(thingie[1])




def fieldContent(operator):

    current = ''

    fields = []

    counter = -1

    for x in range(0,len(operator)):
        if operator[x][0]!=current:
            fields.append([operator[x][0],1])
            current=operator[x][0]
            counter = counter + 1
        else:
            fields[counter][1] = fields[counter][1] + 1

    return fields



# This yield function is complicated (see
# http://docs.python.org/ref/yield.html)
def all_perms(str):
    if len(str) <=1:
        yield str
    else:
        for perm in all_perms(str[1:]):
            for i in range(len(perm)+1):
                yield perm[:i] + str[0:1] + perm[i:]


def listPerms(toPerm):
    list = []

    for perm in all_perms(toPerm):
        list.append(perm)

    return list


def generatePerms(fieldContent):
    counter = 0

    ultimateList = [[]]

    for x in range(0,len(fieldContent)):
        toPerm = []
        for y in range(0,fieldContent[x][1]):
            toPerm.append(y+counter)

        counter = counter + fieldContent[x][1]

        #copy list by value
        oldList = ultimateList[:]

        ultimateList = []

        for perm in oldList:
            for p in listPerms(toPerm):
                ultimateList.append(perm + p)

    return ultimateList



# counts cycles in a perm, e.g. countCycles([1,0,3,4,5,6,2,7])
# return 3
def countCycles(perm):
    counter = 0
    for i in range(0,len(perm)):
        if perm[i]!='done':
            counter = counter + 1
            next = perm[i]
            perm[i] = 'done'
            while(next!=i):
                dummy = perm[next]
                perm[next]= 'done'
                next = dummy

    return counter


# returns perm1 \circ perm2
def multPerms(perm1,perm2):

    perm = []

    for i in range(0,len(perm1)):
        perm.append(perm1[perm2[i]])

    return perm

def inversePerm(perm):

    inverse = [0]*len(perm)
    for i in range(0,len(perm)):
        inverse[perm[i]] = i

    return inverse



def printPerm(perm):

    store  = [0]*len(perm)

    output = ""

    for i in range(0,len(perm)):
        if store[i] == 0:
            output = output + "("

            next = i
            output = output +  str(next)
            store[next] = 1

            while perm[next] != i:

                output = output +  str(perm[next])
                next = perm[next]
                store[next] = 1

            output = output +  ")"

    return output



def computeCorr(operator1,operator2):

    if(fieldContent(operator1)!=fieldContent(operator2)):
        print "cripes, the field content is not the same!"
        sys.exit(0)

    fields = fieldContent(operator1)

    gammas = generatePerms(fields)

    # remove fields
    alpha1 = []
    alpha2 = []

    for i in range(0,len(operator1)):
        alpha1.append(operator1[i][1])
        alpha2.append(operator2[i][1])

    powersOfN = [0]*(len(operator1)+1)


    # multiply alpha1 gamma alpha2 gamma^inverse
    for gamma in gammas:
        perm = multPerms(alpha1,multPerms(gamma,multPerms(alpha2,inversePerm(gamma))))
        noCycles = countCycles(perm)
        powersOfN[noCycles] =  powersOfN[noCycles] + 1


    return powersOfN




# return 0 or 1
def sameOps(operator1,operator2):

    if(fieldContent(operator1)!=fieldContent(operator2)):
        return 0


    fields = fieldContent(operator1)

    gammas = generatePerms(fields)

    # remove fields
    alpha1 = []
    alpha2 = []

    for i in range(0,len(operator1)):
        alpha1.append(operator1[i][1])
        alpha2.append(operator2[i][1])


    result = 0

    for gamma in gammas:
        if alpha1 == multPerms(gamma,multPerms(alpha2,inversePerm(gamma))):
            result = 1
            break

    return result




def printPowersOfN(powersOfN):

    final = ''

    for i in range(0,len(powersOfN)):
        if powersOfN[i]!= 0:
            if i ==0:
                final = final + ' + ' + str(powersOfN[i])
            else:
                final = final + ' + ' + str(powersOfN[i]) +  'N^' + str(i)

    return final[3:]



def multPowersOfN(powersOfN1,powersOfN2):
    result = [0]*len(powersOfN1)*len(powersOfN2)
    for x in range(0,len(powersOfN1)):
        for y in range(0,len(powersOfN2)):
            result[x+y] = result[x+y] + powersOfN1[x]*powersOfN2[y]
    return result



def addPowersOfN(powersOfN1,powersOfN2):
    if len(powersOfN1) > len(powersOfN2):
        result = powersOfN1[:]
        for x in range(0,len(powersOfN2)):
            result[x] = result[x] + powersOfN2[x]
    else:
        result = powersOfN2[:]
        for x in range(0,len(powersOfN1)):
            result[x] = result[x] + powersOfN1[x]
    return result



def printOperator(operator):
    finalString = ''
    for x in range(0, len(operator)):
        if len(operator[x]) == 2:
            dummy = 1
            start = x
            next = x
            while dummy:
                stringToAdd = operator[next][0]
                if stringToAdd[0] == 's':
                    stringToAdd = stringToAdd[1:]
                finalString = finalString + stringToAdd
                operator[next].append("done")
                next = operator[next][1]
                if next == start:
                    dummy = 0
            finalString = finalString + ','

    for x in range(0, len(operator)):
        if operator[x][2] == "done":
            operator[x].pop()

    # remove final comma
    finalString = finalString[:-1]

    return finalString

# finds the first derivative 
def actOperator(operator,factor,result):
    derivative = 0
    for x in range(0, len(operator)):
        if operator[x][0][0] == "d":
            derivative = x
            break
    if derivative ==0:
        addArray([[factor,convertOperator(operator)]],result)
    else:
        for x in range(0, len(operator)):
            if operator[x][0][0] == operator[derivative][0][1]:
                applyDerivative(operator,derivative,x,factor,result)

def applyDerivative(operator,derivative,target,factor,result):
    newOperator = []
    for x in range(0,len(operator)):
        if x == derivative:
            # this is not seen when we print the operator
            newOperator.append(["n",0,0])
        elif x == target:
            newOperator.append(["n",0,0])
        else:
            if operator[x][1] ==derivative:
                if operator[target][1] == target:
                    newOperator.append([operator[x][0],operator[derivative][1]])
                else:
                    newOperator.append([operator[x][0],operator[target][1]])
            elif operator[x][1] ==target:
                newOperator.append([operator[x][0],operator[derivative][1]])
            else:
                newOperator.append(operator[x])
    if operator[derivative][1] == target:
        factor = multPowersOfN([0,1],factor)
    if operator[target][1] == derivative:
        factor = multPowersOfN([0,1],factor)


    actOperator(newOperator,factor,result)



def printArray(array):
    finalString = ''
    for thingie in array:
        finalString = finalString + printPowersOfN(thingie[0]) + ' ' + printOperator(thingie[1]) + ' + '

    if len(array)==0:
        finalString ='nothing   '

    return finalString[:-2]



##adds array 1 to array2
def addArray(array1,array2):


    for x in range(0,len(array1)):
        dummy = 0
        for y  in range(0,len(array2)):
            if sameOps(array1[x][1],array2[y][1]):
                dummy = 1
                array2[y][0] = addPowersOfN(array1[x][0],array2[y][0])
                if array2[y][0] == [0]*len(array2[y][0]):
                    array2.pop(y)
                break
        if dummy == 0:
            array2.append(array1[x])


def convertOperator(operator):
    string = printOperator(operator)
    resultOp = processOp(string)

    return resultOp


def dilatationSingle(factor,operator):

    #"d" signifies a derivative; "s" stops the derivative acting on this
    #letter


    dilatationOperator = [
        [["sX",1],["sY",2],["dX",3],["dY",0]],
        [["sX",1],["sY",2],["dY",3],["dX",0]],
        [["sY",1],["sX",2],["dY",3],["dX",0]],
        [["sY",1],["sX",2],["dX",3],["dY",0]]
        ]


    fixedOperator = []

    for i in range(0,len(operator)):
        fixedOperator.append([operator[i][0],operator[i][1] + 4])

    result = []
    for x in range(0,4):
        operator = fixedOperator[:]
        actOperator(dilatationOperator[x] + operator,multPowersOfN(factor,[(-1)**x]), result)

    return result



def dilatationArray(array):
    result = []
    for thingie in array:
        addArray(dilatationSingle(thingie[0],thingie[1]),result)

    return result


def computeCorrArray(array1, array2):

    result = [0]
    for opPlusFactor1 in array1:
        for opPlusFactor2 in array2:
            result = addPowersOfN(result,multPowersOfN(computeCorr(opPlusFactor1[1],opPlusFactor2[1]),multPowersOfN(opPlusFactor1[0],opPlusFactor2[0])))

    return result




def convertToBinary(number,digits):

    result = [0]*digits

    for i in range(0,digits):
        power = digits-i-1

        result[power] = number/(2**power)

        number = number -  (2**power)*(number/(2**power))

    return result



def antisymmetrise(operator):

    print "after applying antisymmetrisation on",

    print printOperator(operator),

    print "we get"

    counter = 0

    locations = []


    # clumsy loop to give locations list of where antisymmetrisers are
    for x in range(0,len(operator)):

        dummy = 0

        for gibbon in locations:
            if gibbon[0] == x or gibbon[1] ==x:
                dummy = 1
                break

        if dummy == 1:
            continue

        if operator[x][0] in ['0','1','2','3','4','5','6','7','8','9']:
            counter = counter + 1
            locations.append([x,0])
            dummy = 0
            for y in range(x+1,len(operator)):
                if operator[y][0] == operator[x][0]:
                    locations[counter-1][1] = y
                    dummy = 1-dummy
            if dummy==0:
                print "no matching antisymmetrisers!"
                sys.exit(0)

    array = []

    ## works for counter = 0 too since 2**counter = 1
    ## this loop put in the antisymmetrised X's and Y's
    for i in range(0,2**counter):
        binary = convertToBinary(i,counter)

        # newOperator = operator[:] doesn't work here because the
        # elements of operator[:] are pointers to lists rather than
        # immutable objects such as integers or strings; see Python
        # Cookbook recipe 4.5

        newOperator = []

        for x in range(0,len(operator)):
            newOperator.append([operator[x][0],operator[x][1]])

        sumBin = 0

        for j in range(0,counter):
            if binary[j] == 0:
                newOperator[locations[j][0]][0] = 'X'
                newOperator[locations[j][1]][0] = 'Y'
            else:
                newOperator[locations[j][0]][0] = 'Y'
                newOperator[locations[j][1]][0] = 'X'

            sumBin = sumBin + binary[j]

        sortOp(newOperator)

        addArray([[[(-1)**sumBin],newOperator]],array)



    print printArray(array)

    return array



def multArrayByPowersOfN(powersOfN,array):

    for thingie in array:
        thingie[0] = multPowersOfN(powersOfN,thingie[0])


def antisymmetriseArray(array):

    result = []

    for thingie in array:
        arrayToAdd = antisymmetrise(thingie[1])
        multArrayByPowersOfN(thingie[0],arrayToAdd)
        addArray(arrayToAdd,result)

    return result


def oneLoopArray(array1,array2):

    dilOp2 = dilatationArray(array2)

    sortOpArray(dilOp2)

    return computeCorrArray(array1,dilOp2)




if __name__ == "__main__":

    treeOpt = 0
    dilOpt = 0
    oneloopOpt = 0
    sunOpt = 0

    try:
        options, arguments = getopt.gnu_getopt(sys.argv[1:], 
        '01dhs')
    except getopt.error:
        print 'error: you tried to use an unknown option or the argument for an option that requires it was missing; try \'correlator.py -h\' for more information'
        sys.exit(0)

    for o,a in options:
        if o == '-h':
            print __doc__
            sys.exit(0)

        elif o == '-s':
            sunOpt = 1
            print "no SU(N) option yet"
            sys.exit(0)

        elif o == '-0':
            treeOpt = 1

        elif o == '-1':
            oneloopOpt = 1

        elif o == '-d':
            dilOpt = 1

    if len(arguments) == 0:
        print 'you didn\'t specify an operator'
        sys.exit(0)
    elif len(arguments) == 1:
        dilOpt = 1
        if treeOpt + oneloopOpt > 0:
            print 'you didn\'t specify enough operators'
            sys.exit(0)
    elif len(arguments) == 2:
        if treeOpt + oneloopOpt == 0:
            treeOpt = 1
        if dilOpt == 1:
            print 'you should only feed the dilatation operator a single operator'
            sys.exit(0)
    else:
        print 'you gave too many arguments'
        sys.exit(0)

    if dilOpt:
        print "computing the dilatation of " + arguments[0]

        operator = processOp(arguments[0])

        resultArray = dilatationArray(antisymmetrise(operator))

        print 'final result:'

        print printArray(resultArray)



    if treeOpt:
        print "computing the tree-level correlator of " + arguments[0] + " and " + arguments[1]

        operator1 = processOp(arguments[0])

        operator2 = processOp(arguments[1])

        powersOfN = computeCorrArray(antisymmetrise(operator1),antisymmetrise(operator2))

        print "the tree-level correlator gives "

        print printPowersOfN(powersOfN)



    if oneloopOpt:
        print "computing the one-loop correlator of " + arguments[0] + " and " + arguments[1]

        operator1 = processOp(arguments[0])

        operator2 = processOp(arguments[1])

        powersOfN = oneLoopArray(antisymmetrise(operator1),antisymmetrise(operator2))

        print "the one-loop correlator gives "

        print printPowersOfN(powersOfN)