ainfmf.lib

version="1.0";
category="Miscellaneous";
info="
LIBRARY:   ainfmf.lib  A-infinity minimal models of MFs
AUTHOR:    Daniel Murfet
KEYWORDS:  matrix factorisation

PROCEDURES:
";

// NOTE: We include a version of matrix.lib which suppresses some
// unnecessary output from the procedure "rowred"
LIB "linalg.lib";
LIB "matrix.lib";
LIB "ring.lib";
LIB "moduli.lib";
	// we use seqsfrominterval, monomialdict

////////////////////////////////////////////////////////////////////
// USAGE GUIDE
//

////////////////////////////////////////////////////////////////////
// TERMINOLOGY
//
// There are three types of interactions in the Feynman rules:
//
//  - the trivalent vertex we call a Y-type
//  - the bivalent vertex with one boson and one fermion we call I-type
//  - the bivalent vertex with two fermions we call U-type

////////////////////////////////////////////////////////////////////
// DATA FORMATS
//
// For us a tree is a connected, acyclic undirected graph in which 
// each vertex is either a leaf or trivalent, and there is a chosen
// leaf called the root.
//
// We present such a tree as a list ((1,2),(3,4)) with the bracketing
// representing the tree in the usual way:
//
//      1     2    3     4
//      \     /    \     /
//       \   /      \   /
//        \ /        \ /
//         \          /
//          \        /
//           \      /
//            \    /
//             \  /
//              \/
//               |
//               |
//
// Since it is convenient, we refer to leaves other than the root as
// "input vertices" and trivalent vertices as "internal vertices".

////////////////////////////////////////////////////////////////////
// enumerateTrees
//
// Returns all trees with n leaves.

proc enumerateTrees(int q)
{
	if( q < 1 )
	{
		print("[enumerateTrees] Number of inputs needs to be at least 2.");
        return();
	}
	
	return( enumerateTrees_r(q, 0) );
}

proc enumerateTrees_r(int q, int offset)
{
	// Recursive worker for enumerateTrees, generates all
	// trees with q input vertices, with labels {offset+1, offset+2, ... }
	
	// By recursion. Associated to each tree with n inputs a partition
	// of q = l + r where there are l \ge 1 leaves coming into the left
	// branch and r \ge 1 leaves coming into the right branch (where "coming
	// into" means that these leaves traverse that branch on their way to
	// the root)
	
	if( q == 1 )
	{
		return(offset+1);
	}
	
	if( q < 1 )
	{
		print("[enumerateTrees] Integer q needs to be at least 2.");
        return();
	}
	
	if( offset < 0 )
	{
		print("[enumerateTrees] Offset needs to be at least zero.");
        return();
	}
	
	list allTrees;
	
	int i,j,k;
	list L, R;
	list tempTree;
	for(i=1; i<=q-1; i++)
	{
		L = enumerateTrees_r(i,offset);
		R = enumerateTrees_r(q-i,offset+i);
		
		for(j=1; j<=size(L); j++)
		{
			for(k=1; k<=size(R); k++)
			{
				tempTree = list(L[j],R[k]);
				allTrees = allTrees + list(tempTree);
			}
		}
	}
	
	return(allTrees);
}

////////////////////////////////////////////////////////////////////
// enumerateLocations
//
// Given a tree T returns a list of locations in the tree. A location
// is a vertex other than the root, or an internal edge. Vertices are
// denoted by the subtree that they represent, while edges are denoted
// by list(S) where S is the subtree which is connected by the edge to
// the rest of the tree.
//
// The root of the tree is not a location.
//
// For example, in the tree ((1,2),(3,4)) we have the following locations:
//
//   1, 2, 3, 4 - input vertices
//   (1,2), (3,4) - internal vertices
//   ((1,2),(3,4)) - internal vertex
//   ((1,2)), ((3,4)) - internal edges
//   
//
//      1     2    3     4
//      \     /    \     /
//       \   /      \   /
//        \ /        \ /
//       (1,2)      (3,4)
//         \          /
//          \        /
//        ((1,2)) ((3,4))
//            \    /
//             \  /
//              \/
//         ((1,2),(3,4))
//               |
//               |
//
// Note that input vertices can be distinguished by the fact that they are
// integers, not lists, and internal edges by the fact that they are lists
// of size 1, while interval vertices are lists of size 2.

proc enumerateLocations(list T)
{
	if( isValidTree(T) != 1 )
	{
        return();
	}
	
	// If T is a single vertex then it only has one location
	if( size(T) == 1 )
	{
		return(T);
	}
	
	// Otherwise we return the union of the locations in each branch
	list L = T[1];
	list R = T[2];
	
	list loc = enumerateLocations(L) + enumerateLocations(R);
	
	// Add the bottom-most vertex
	loc = loc + list(T);
	
	// Plus the internal edge connecting each branch to the root
	// but e.g. the edge connecting L to the tree is only internal
	// if L is not a singleton
	if( size(L) > 1 )
	{
		loc = loc + list(list(L));
	}
	
	if( size(R) > 1 )
	{
		loc = loc + list(list(R));
	}
	
	return(loc);
}

////////////////////////////////////////////////////////////////////
// enumerateInternalProp
//
// Given a tree T and potential W, enumerates the possible internal
// propagations of boson and fermions. This comes in the form of a list,
// each entry of which is one possible propagation pattern. Such a pattern
// is an oriented labelled graph, stored as a list of tuples
//
//    ( particle type, source, target )
//
// where 
//
//   - particle type is either "fermion" or "boson"
//   - source and target are locations in the tree
//
// So a pattern is a list of tuples of this kind, and we return the list L
// of all valid patterns. We are presented W as its components Wcomp, i.e.
// such that W is x(1) * Wcomp[1] + ... + x(n) * Wcomp[n].

proc enumerateInternalProp(list T, list Wcomp)
{
	int n = nvars(basering);
	
	if( isValidTree(T) != 1 )
	{
		return();
	}
	
	// Generate the list of valid monomials
	int i;
	list valid_monoms;
	for(i=1;i<=size(Wcomp);i++)
	{
		valid_monoms[i] = monoms_in(Wcomp[i]);
	}
	
	// First we enumerate the internal edges, inputs and vertices
	list loc = enumerateLocations(T);
	
	list loc_edge;
	list loc_input;
	list loc_vertex;

	for(i=1;i<=size(loc);i++)
	{
		def x = loc[i];
		if( locIsEdge(x) == 1 ){ loc_edge = loc_edge + list(x); }
		if( locIsInput(x) == 1 ){ loc_input = loc_input + list(x); }
		if( locIsVertex(x) == 1 ){ loc_vertex = loc_vertex + list(x); }
	}
	
	list loc_inputoredge = loc_edge + loc_input;
	
	// First we generate a list of all possible ways of linking up I-vertex
	// generated fermions with U-vertices. That is, we
	//
	//   - choose, for each internal edge e an integer t(e), which means
	//   that the interaction eats an x(e)-variable and emits a theta_e
	//   fermion.
	//
	//   - then we choose at which U-vertex v this theta is annihilated (this
	//   must be downstream from e, and the path from e to the root must
	//   enter v ON THE RIGHT.
	//
	//   - and which Y-vertex emitted the boson (this must be upstream from e)
	//
	// so our patterns should contain two edges (i.e. two tuples of the kind
	// discussed above)
	
	list L1;
	
	int i,k,l;
	
	for(i=1;i<=size(loc_edge);i++)
	{
		def e = loc_edge[i];
		list templist;
		
		for(k=1;k<=size(loc_vertex);k++)
		{
			def v = loc_vertex[k];
			
			// We have to make sure that v is downstream from e,
			// and that e is linked to the right part of v
			list edge_aboveright = list(v[2]);
				
			if( locIsAbove(e, edge_aboveright) )
			{
				// And where does the boson originate
				for(l=1;l<=size(loc_inputoredge);l++)
				{
					def u = loc_inputoredge[l];
					// u must be upstream from e, and _different_ from e
					if( locToString(u) != locToString(e) && locIsAbove(u, e) )
					{
						//print("      And Y-interaction at " + locToString(u) );
						list c = list( "fermion", e, v );
						list d = list( "boson", u, e );
						templist = templist + list(list(c,d));
					}
				}
			}
		}
				
		L1[i] = templist;
	}
		
	return(L1);
} 

////////////////////////////////////////////////////////////////////
// tikzTree_locpos

proc tikzTree_locpos(def x, int input_offset, list T)
{
	int level_offset = 5;
	
	if( locIsInput(x) )
	{
		return( list((x-1) * input_offset,0) );
	}
	
	if( locIsVertex(x) )
	{
		list l = tikzTree_locpos(x[1], input_offset, T);
		list r = tikzTree_locpos(x[2], input_offset, T);
		
		number x_coord = number( l[1] + r[1] ) / 2;
		
		number lowest = l[2];
		if( r[2] < l[2] ){ lowest = r[2]; }
		
		number y_coord = lowest - level_offset;

		return( list(x_coord,y_coord) );
	}
	
	if( locIsEdge(x) )
	{
		list loc = enumerateLocations(T);
		// Find the next vertex down
		def next_vertex;
		
		int i;
		for(i=1;i<=size(loc);i++)
		{
			def y = loc[i];
			
			if( locIsVertex(y) )
			{
				if( locIsVertex(y[1]) && locToString(y[1]) == locToString(x[1]) )
				{
					next_vertex = y;
				}
				
				if( locIsVertex(y[2]) && locToString(y[2]) == locToString(x[1]) )
				{
					next_vertex = y;
				}
			}	 
		}
		
		// We are halfway between x[1] and next_vertex
		list above = tikzTree_locpos(x[1], input_offset, T);
		list below = tikzTree_locpos(next_vertex, input_offset, T);
		
		number x_coord = number( above[1] + below[1] ) / 2;
		number y_coord = number( above[2] + below[2] ) / 2;
		
		return(list(x_coord, y_coord));
	}
}

////////////////////////////////////////////////////////////////////
// tikzTree_intprop
//
// Given a tree and a pattern of internal propagations
// outputs a series of TikZ commands to draw the tree
// in LaTeX (to be added 

proc tikzTree_intprop(list T, list prop)
{
	list L_intro = tikzTree_intro();
	list L_core = tikzTree_core(T);
	list L_end = tikzTree_end();
	
	// We just add to L_core our additional instructions
	// prop is a list of tuples of the form
	// 
	//
	//    ( particle type, source, target )
	//
	// where 
	//
	//   - particle type is either "fermion" or "boson"
	//   - source and target are locations in the tree
	
	int i;
	for(i=1;i<=size(prop);i++)
	{
		def s = prop[i][2];
		def t = prop[i][3];
		
		string A = "\\draw" + prop[i][1] + " (" + locToString(s) + ") to (" + locToString(t) + ");";
		
		L_core = L_core + list(A);
	}
	
	return(L_intro + L_core + L_end);
}

////////////////////////////////////////////////////////////////////
// tikzTree
//
// Given a tree, outputs a series of TikZ commands to draw the tree
// in LaTeX. Returns the commands as a list of strings (one string
// = 1 line of TeX).

proc tikzTree_latexpreamble()
{
	list L;
	L = L + list("\\documentclass[english,letter paper,12pt,leqno]{article}");
	L = L + list("\\usepackage{amsmath, amscd, amssymb, mathrsfs, accents, amsfonts,amsthm}");
	L = L + list("\\usepackage{tikz}");

	L = L + list("\\setlength{\\evensidemargin}{0.1in}");
	L = L + list("\\setlength{\\oddsidemargin}{0.1in}");
	L = L + list("\\setlength{\\textwidth}{6.3in}");
	L = L + list("\\setlength{\\topmargin}{0.0in}");
	L = L + list("\\setlength{\\textheight}{8.5in}");
	L = L + list("\\setlength{\\headheight}{0in}");

	L = L + list("\\begin{document}");

	return(L);
}

proc tikzTree_intro()
{
	list L;
	
	L = L + list("\\def\\drawl{\\draw[color=black!10, line width=1pt]}");
	L = L + list("\\def\\drawfermion{\\draw[color=black, line width=1pt]}");
	L = L + list("\\def\\drawboson{\\draw[color=teal, dashed, line width=1pt]}");
	L = L + list("\\begin{center}");
	L = L + list("\\begin{tikzpicture}[scale=0.2,auto]");
	
	return(L);
}

proc tikzTree_core(list T)
{
	list L;
	
	// Some constants
	int input_offset = 10;
	
	// First we enumerate the internal edges, inputs and vertices
	list loc = enumerateLocations(T);
	
	list loc_edge;
	list loc_input;
	list loc_vertex;

	int i;
	for(i=1;i<=size(loc);i++)
	{
		def x = loc[i];
		if( locIsEdge(x) == 1 ){ loc_edge = loc_edge + list(x); }
		if( locIsInput(x) == 1 ){ loc_input = loc_input + list(x); }
		if( locIsVertex(x) == 1 ){ loc_vertex = loc_vertex + list(x); }
	}
	
	int i;
	for(i=1;i<=size(loc_input);i++)
	{
		string s = "\\node (" + string(i) + ") at (" + string( (i-1)*input_offset ) + ",0) {$" + string(i) + "$};";
		L = L + list(s);
	}
	
	for(i=1;i<=size(loc);i++)
	{
		def x = loc[i];
		
		if( locIsInput(x) != 1 )
		{
			list c = tikzTree_locpos(x, input_offset, T);
			
			string s = "\\coordinate (" + locToString(x) + ") at (" + string(c[1]) + "," + string(c[2]) + ");";
			
			L = L + list(s);
			
			// Draw the lines connecting this node to its parents
			if( locIsVertex(x) )
			{
				string c1 = "\\drawl (" + locToString(x) + ") to (" + locToString(x[1]) + ");";
				string c2 = "\\drawl (" + locToString(x) + ") to (" + locToString(x[2]) + ");";
				
				L = L + list(c1);
				L = L + list(c2);
			}
		}
	}
	
	return(L);
}

proc tikzTree_end()
{
	list L;
	L = L + list("\\end{tikzpicture}");
	L = L + list("\\end{center}");
	
	return(L);
}

proc tikzTree(list T)
{
	list L_intro = tikzTree_intro();
	list L_core = tikzTree_core(T);
	list L_end = tikzTree_end();
	
	return(L_intro + L_core + L_end);
}

////////////////////////////////////////////////////////////////////
// subsetsofinterval
//
// Given an integer N >=1 and M >= 1 return the set of all subsets
// of {1,...,N} of size M. So obviously M should be \le N.

proc subsetsofinterval(int N, int M)
{
	// print("(N,M) = " + string(N) + "," + string(M));
	if( M > N )
	{
		print("[subsetsofinterval] should not have M > N");
		return();
	}
	
	if( M == 1 )
	{
		list L;
		int i;
		for(i=1;i<=N;i++)
		{
			L = L + list(list(i));
		}
		return(L);
	}
	
	// So M > 1 and N > 1 and we can divide the subsets
	// into those containing 1 and those not containing 1.
	
	list L;
	list L_0 = subsetsofinterval(N-1, M-1);
	
	// Generate the subsets that contain 1:
	int i,j;
	for(i=1;i<=size(L_0);i++)
	{
		list S = list(1);
		
		for(j=1;j<=size(L_0[i]);j++)
		{
			S = S + list(L_0[i][j] + 1);
		}
		
		L = L + list(S);
	}
	
	// And the subsets that do not contain 1:
	if( M < N )
	{
		list L_1 = subsetsofinterval(N-1, M);
		
		int i,j;
		for(i=1;i<=size(L_1);i++)
		{
			list S;
		
			for(j=1;j<=size(L_1[i]);j++)
			{
				S = S + list(L_1[i][j] + 1);
			}
		
			L = L + list(S);
		}
	}
	
	return(L);
}

proc monoms_in(poly f)
{
	int numVars = nvars(basering);
    poly xprod = 1;
    int i;
    for(i=1; i<=numVars; i++)
    {
        xprod = xprod * var(i);
    }
    matrix koffer = coef(f, xprod);
    
    // the first row of koffer contains the monomials and the second row the coeffs
    // store this as a list of monomial intvecs and the coeffs
    list mm_list;
    for(i=1; i<=ncols(koffer); i++)
    {
    	mm_list = mm_list + list(leadexp(koffer[1,i])); // leadexp extracts intvec
    }
    
    return(mm_list);
}

////////////////////////////////////////////////////////////////////
// prodlist
//
// Given a list L = (L[1], ..., L[k]) of lists, returns the product
// list, i.e. the list L[1] x ... x L[k] whose entries are lists, 
// the first entry from L[1], the second from L[2], etc. 

proc prodlist(list L)
{
	if( size(L) == 1 )
	{
		print("[prodlist] Don't do that.");
		return();
	}
	
	int i;
	for(i=1;i<=size(L);i++)
	{
		//print(string(size(L[i])));
		if( size(L[i]) == 0 )
		{
			print("[prodlist] No empty lists.");
			return();
		}
	}
	
	// Indices are intvecs I of length k, with 1 <= I[i] <= size(L[i]) for
	// each i. We iterate through these in lexicographic order
	
	int k = size(L);
	
	intvec index;
	int i;
	for(i=1;i<=k;i++)
	{
		index[i] = 1;
	}
	
	list retlist;
	while(1)
	{
		// Add the tuple corresponding to index
		list templist;
		//print(index);
		
		for(i=1;i<=k;i++)
		{
			templist[i] = L[i][ index[i] ];
		}
		retlist = retlist + list(templist);
		
		// Increment
		index[k] = index[k]+1;
		
		int curr = k;
		int ex = 0;
		while( index[curr] > size(L[curr]) )
		{
			if( curr == 1 )
			{
				ex = 1;
				break;
			}
			
			index[curr] = 1;
			index[curr-1] = index[curr-1] + 1;
			curr = curr - 1;
		}
		
		if( ex == 1 ){ break; }
	}
		
	return(retlist);
}

////////////////////////////////////////////////////////////////////
// locToString
//

proc locToString(def x)
{
	if( locIsInput(x) )
	{
		return(string(x));	
	}
	
	if( locIsVertex(x) )
	{
		string l = locToString(x[1]);
		string r = locToString(x[2]);
		return( "<" + l + r + ">" );
	}
	
	if( locIsEdge(x) )
	{
		return( "e" + locToString(x[1]) );
	}
}

////////////////////////////////////////////////////////////////////
// locIsAbove
//
// Given two locations x,y returns 1 if x is above y in the tree (i.e.
// the unique path from x to the root passes through y) and 0 otherwise.
// Note that "above" includes "equal".

proc locIsAbove(def x, def y)
{
	if( typeof(x) == typeof(y) )
	{
		if( locToString(x) == locToString(y) )
		{
			return(1);
		}
	}
	
	if( typeof(y) != "list" )
	{
		return(0);
	}
	
	if( size(y) == 1 )
	{
		return(locIsAbove(x,y[1]));
	}
	
	if( size(y) == 2 )
	{
		int l = locIsAbove(x,y[1]);
		int r = locIsAbove(x,y[2]);
		
		if( l == 1 || r == 1 )
		{
			return(1);
		}
		
		list edge_aboveright = list( y[2] );
		list edge_aboveleft = list( y[1] );
		
		if( locToString(edge_aboveright) == locToString(x) ||
			locToString(edge_aboveleft) == locToString(x) )
		{
			return(1);
		}
		
		return(0);
	}
}

////////////////////////////////////////////////////////////////////
// locIsInput

proc locIsInput(def x)
{
	if( typeof(x) == "int" )
	{
		return(1);
	}
	else
	{
		return(0);
	}
}

////////////////////////////////////////////////////////////////////
// locIsEdge

proc locIsEdge(def x)
{
	if( typeof(x) == "list" && size(x) == 1 )
	{
		return(1);
	}
	else
	{
		return(0);
	}
}

////////////////////////////////////////////////////////////////////
// locIsVertex

proc locIsVertex(def x)
{
	if( typeof(x) == "list" && size(x) == 2 )
	{
		return(1);
	}
	else
	{
		return(0);
	}
}

////////////////////////////////////////////////////////////////////
// isValidTree
//

proc isValidTree(list T)
{
	if( size(T) > 2 )
	{
		print("[isValidTree] Badly formatted tree.");
        return(0);
	}
	
	return(1);
}

////////////////////////////////////////////////////////////////////
// listofints
//
// Given N >= 1 returns the list 1,...,N

proc listofints(int n)
{
	list L;
	
	int i;
	for(i=1;i<=n;i++)
	{
		L = L + list(i);
	}
	
	return(L);
}

////////////////////////////////////////////////////////////////////
// checklist
//
// checklist(l,p) gives 1 if p is an element of the list l,
// otherwise it gives 0.
////////////////////////////////////////////////////////////////////

proc checklist(list l, p)
{
    int i;
    for(int i=1; i<=size(l); i++)
    {
        if(p==l[i])
        {
            return(1);
        }
    }
    return(0);
}