Trying an alternative approach to dealing with the algebra of qft computationally based on wick contraction. We can easily find an algorithm that finds all possible pairs. Then we just reduce it all accordingly. Some problems: The combinatorics are going to freak out pretty fast.
I think my use of functional stuff like maps and lambdas is intensely unclarifying the code.
import numpy as np
#Needs to return a list of lists of pairs
if len(mylist) == 2:
return [[(mylist, mylist)]]
#popoff the first element
element1 = mylist
pairs = 
for i in range(1,len(mylist)):
element2 = mylist[i]
#get all the pairs expluidng the already picked pair
subpairs = pair(mylist[1:i] + mylist[i+1:])
#Put the picked pair back in x is a list of pairs.
pairs = pairs + map(lambda x: [(element1,element2)] + x ,subpairs)
#For Fermionic pairing, it might be conveneitn to extend pairs to (element1,element2,+-1)
#With the sign depending on wheter i is even or odd
pairing = pair([1,2,3,4])
#Here's an idea, I can use a and adag as objects with indices
#and return functions (two point green's functions)
# or could return matrcies that are discretized green's functions.
if twoguys == 'a' and twoguys == 'adag':
#Confusing, but what
pairings = pair(['a','a','a','adag','adag','adag'])
#could abstract out multiplication as a function passed in. Then if contraction returns functions, I could
#return a new function that multiplies the inner functions
#or could replace multiply with np.dot if returning matrices
return reduce(lambda acc, val: acc * val,map(contract,pairing))