mfweb.lib

version="1.0";
category="Miscellaneous";
info="
LIBRARY:   mfweb.lib  Compilation of webs of matrix factorisations
AUTHOR:    Nils Carqueville, Daniel Murfet
KEYWORDS:  matrix factorisation

PROCEDURES:
";

LIB "linalg.lib";
LIB "matrix.lib";
LIB "ring.lib";
LIB "blow.lib";
LIB "linkhom.lib";

////////////////////////////////////////////////////////////////////
// USAGE GUIDE
//
// A web is an oriented graph where edges are allowed to begin or end, but not
// both, on some unspecified "boundary". All edges are labelled with a set of variables
// and a potential, and every vertex is assigned a matrix factorisation
// of the sum of all the "outgoing" potential minus the sum of all the incoming
// potentials. There must be at least one variable on each edge, disjoint variables
// on distinct edges, and potentials must have isolated singularities. Throughout
// comments in this file, V denotes the total number of vertices in the web (always >0)
// We implicitly order the vertices and consider them as elements in the set {1,...,V}.
//
// WARNING: We do not allow isolated vertices, with no incoming or outgoing vertices.
// The ordering of the variables in L1 and L2 may matter.

////////////////////////////////////////////////////////////////////
// USAGE GUIDE: parameters
//
// We allow rings of the form (0,r),(variables),dp where r is a parameter.
// If this is an algebraic extension with minimum polynomial "minpoly" then 
// there should also be defined a polynomial minpolyblow, in the first ring variable
// var(1), equal to substituting var(1) for r in minpoly.

/////////////////////////////////////////////////
// DATA TYPE: edge
//
// An edge is a list (S,T,L,W) where S,T are integers in the range 0,...,V, L is a
// nonempty list of variables in the ring, and W is a polynomial in the variables
// present in L. Here S, T are the indices, respectively of the "S"ource and "T"arget
// of the edge, where "0" stands for the boundary.

/////////////////////////////////////////////////
// DATA TYPE: augmented edge
//
// An augmented edge is a list (S,T,L,W,B,C) where (S,T,L,W) is an edge, B is a list
// of size n = size(L), and C is a nxn matrix of polynomials such that for each i
// between 1 and n = |L|, we have
//
//      L[i]^B[i] = C[1,i] * diff(W,L[1]) + ... + C[n,i] * diff(W,L[n]).
//
// The existence of such integers B[i] and polynomials C[i,j] follows from the assumption
// that W defines an isolated singularity, hence the partial derivatives of W generate an
// ideal primary for the ideal generated by all the variables on the edge.
//
// NOTE (8/7/2015): The old code had the C matrix rows and columns confused. The variables
// are indexed by i, and each variable has a corresponding _column_ in C.

/////////////////////////////////////////////////
// TERMINOLOGY
//
// We say that a variable is "adjacent" at a vertex if the variable lies on an edge
// adjacent to the vertex. A variable is "external" if it lies on an edge involving
// the boundary, and "internal" otherwise.
//
// The total potential of a web is the sum of the potentials associated with
// each vertex of the web.

/////////////////////////////////////////////////
// DATA TYPE: web
//
// A web is represented by a tuple (V, L1, L2[,L3]) where V is an integer > 0, giving
// number of vertices in the web, and two or three lists, where
//
//  1. L1 is a list of edges.
//  2. L2 is a list of matrices, of size V, with the matrix in posn i
//  being the differential of a matrix factorisation defined in the 
//  set of variables adjacent with the vertex i, factorising the potential
//  explained in the Usage Guide.
//  3. L3 is an optional list of intvecs, also of size v, making each
//  of the MFs into a graded MF.
//
// NOTE: Distinct edges in the web MUST be assigned disjoint sets of variables!

/////////////////////////////////////////////////
// DATA TYPE: compilation strategy
//
// Compilation of a web means taking the tensor product over all MFs in the
// web and pushing forward to the ring with all internal variables in the web
// deleted. A compilation strategy represents one way of effecting this pushforward
// in practice.
//
// Precisely, a compilation strategy consists of a pair (L1,L2), where
// L1 is a permutation of the set {1,...,V} and L2 is a list
// containing each internal variable in the web precisely once. The idea is that L1 tells
// us the order in which to tensor the MFs in the web together (for example, if
// L1 were 3,1,2 we should first tensor the MF at vertex 3 with the one at vertex 1,
// and then tensor the result with the MF at vertex 2) and L2 tells us the ordering
// on the variables to use in the reduction/pushforward procedure.
//
// NOTE: There is no implied relationship between the variable ordering and the
// vertex ordering, e.g., the variables around the first vertex do not need to
// go first.
//
// IMPORTANT: For a compilation strategy to be valid, each fusion has to involve
// a nonzero number of variables. That is, if L1 = {v1,...,vn} then for each 1 <= i <= n-1
// we are trying to fuse v(i+1) with the result of fusing v1,...,vi. We require that
// there is an edge between one of the v1,...,vi and v(i+1).

/////////////////////////////////////////////////
// Settor and Accessor methods for an edge.
//
// The point of these being that if we change the structure of the lists underlying
// an edge, or web, then we only need change the following routines.

proc edgeSource(list edge){ return(edge[1]); }
proc edgeTarget(list edge){ return(edge[2]); }
proc edgeVariables(list edge){ return(edge[3]); }
proc edgePotential(list edge){ return(edge[4]); }
proc augedgePowers(list augedge){ return(augedge[5]); }
proc augedgeTransform(list augedge){ return(augedge[6]); }
proc setEdgeSource(list edge, int src){ edge[1] = src; return(edge); }
proc setEdgeTarget(list edge, int tar){ edge[2] = tar; return(edge); }
proc setEdgeVariables(list edge, list var){ edge[3] = var; return(edge); }
proc setEdgePotential(list edge, poly W){ edge[4] = W; return(edge); }
proc isEdgeBetweenVertices(list edge, int v1, int v2)
{
    if( edgeSource(edge) == v1 && edgeTarget(edge) == v2 ){ return(1); }
    if( edgeSource(edge) == v2 && edgeTarget(edge) == v1 ){ return(1); }
    return(0);
}
proc edgeVerify(list edge)
{
    if( typeof(edge[1]) != "int" ){ return(0); }
    if( edge[1] < 0 ){ return(0); }
    
    if( typeof(edge[2]) != "int" ){ return(0); }
    if( edge[2] < 0 ){ return(0); }
    
    if( typeof(edge[3]) != "list" ){ return(0); }
    if( size(edge[3]) == 0 ){ return(0); }
    
    int i;
    for(i=1; i<=size(edge[3]); i++)
    {
       // Each entry must be a ring variable
        if( rvar(edge[3][i]) == 0 ){ return(0); }
    }
    
    if( typeof(edge[4]) != "poly" ){ return(0); }
    
    return(1);
}

// NOTE: augmentedEdgeForEdge breaks the setter/getter barrier.

/////////////////////////////////////////////////
// Settor and Accessor methods for a web.

proc webSize(list web){ return(web[1]); }
proc webEdges(list web){ return(web[2]); }
proc webMFs(list web){ return(web[3]); }
proc webMFatVertex(list web, int vertex){ return(web[3][vertex]); }
proc webMFgr(list web){ return(web[4]); }
proc webMFgratVertex(list web, int vertex){ return(web[4][vertex]); }
proc webIsGraded(list web){ int t = (size(web) == 4); return(t); }

/////////////////////////////////////////////////
// Basic methods for vertices

proc vertexVariables(list web, int vertex)
{
    list vars;
    list edgeList = webEdges(web);
    
    int i;
    for(i = 1; i <= size(edgeList); i++)
    {
        list edge = edgeList[i];
        if( edgeSource(edge) == vertex || edgeTarget(edge) == vertex )
        {
            vars = vars + edgeVariables(edge);
        }
        kill edge;
    }
    
    kill edgeList;
    return(vars);
}

/////////////////////////////////////////////////
// vertexIncomingEdges
//
// Return a list (possibly empty) of indices into
// the array of edges giving those edges with target
// the given vertex (possibly the boundary 0).

proc vertexIncomingEdges(list web, int vertex)
{
    list eList;
    list edgeList = webEdges(web);
    
    int i;
    for(i = 1; i <= size(edgeList); i++ )
    {
        list edge = edgeList[i];
        if( edgeTarget(edge) == vertex )
        {
            eList = eList + list(i);
        }
        kill edge;
    }
    
    return(eList);
}

proc vertexOutgoingEdges(list web, int vertex)
{
    list eList;
    list edgeList = webEdges(web);
    
    int i;
    for(i = 1; i <= size(edgeList); i++ )
    {
        list edge = edgeList[i];
        if( edgeSource(edge) == vertex )
        {
            eList = eList + list(i);
        }
        kill edge;
    }
    
    return(eList);
}

proc vertexPotential(list web, int vertex)
{
    poly W;
    list edgeList = webEdges(web);
    
    int i;
    for(i=1; i <= size(edgeList); i++ )
    {
        list edge = edgeList[i];
        int nc;
        
        if( edgeSource(edge) == vertex ){ nc = 1; }
        if( edgeTarget(edge) == vertex ){ nc = -1; }
        
        W = W + nc * edgePotential(edge);
        
        kill edge;
    }
    
    return(W);
}

/////////////////////////////////////////////////
// ringWithoutVars
//
// Given a set of variables in the ring, returns
// the ring with those variables deleted (but otherwise
// preserved, for example parameters). We do NOT change
// the current basering. Note that vars should be a list
// of polynomials, not strings.

proc ringWithOnlyVars(list vars)
{
    list exvars;
    int i;
    for(i=1; i<=nvars(basering); i++)
    {
        if( !checklist(vars,var(i)) )
        {
            exvars = exvars + list(var(i));
        }
    }
    
    def nR = ringWithoutVars(exvars);
    return(nR);
}

proc ringWithoutVars(list vars)
{
    if( size(vars) == nvars(basering) )
    {
        print("[ringWithoutVars] Asked to create a ring with no variables, exiting.");
        return(0);
    }
    
    def SSS = basering;
    list rlist = ringlist(SSS);
    list newVar;
    
    // Go through and form newVar = current variables minus vars
    int i;
    for(i = 1; i <= size(rlist[2]); i++)
    {
        int found = 0;
        int j;
        for(j = 1; j <= size(vars); j++ )
        {
            if( string(vars[j]) == rlist[2][i] )
            {
                found = 1;
                break;
            }
        }
        
        if( !found ){ newVar = newVar + list(rlist[2][i]); }
    }
    
    list newringList = rlist;
    newringList[2] = newVar;

    // Fix the variable weighting
    intvec kk = (1..size(newVar));
    for(i=1; i<=size(newVar); i++)
    {
        kk[i] = 1;
    }
    newringList[3][1][2] = kk;
    kill kk;
    
    int useParameter;
    if( defined(minpolyblow) && defined(r) )
    {
        newringList[1][4][1] = 0;     // set minpoly to zero, put it back to the right value below...
        useParameter = 1;
    }
    
    // Create our new ring
    def nR = ring(newringList);
    kill rlist, newringList;
    
    // Now complete the correct definition of nR by specifying the right minpoly: 
    if( useParameter )
    {
        setring nR;
        
        // Note that fetch uses the position of the ring variables, not the name, so 
        // that fetch(var(1)) is always var(1).
        poly minpolyblow = fetch(SSS,minpolyblow);
        poly z = subst(minpolyblow,var(1),r);
        number nu = leadcoef(z);
        minpoly = nu;
        poly minpolyblow = fetch(SSS,minpolyblow);
        export(minpolyblow);
        // NOTE: The double definition of minpolyblow is because Singular nukes
        // all variables in our ring when we set minpoly
        
        setring SSS;
    }
    
    return(nR);
}


/////////////////////////////////////////////////
// augmentedEdgeForEdge
//
// Given an edge, creates an associated augmented edge

proc augmentedEdgeForEdge(list edge)
{
    int powerCutoff = blowFlags("power_cutoff_jacobi"); 
    int useSanityChecks = blowFlags("sanity_checks");
    
    if( useSanityChecks )
    {
        if( !edgeVerify(edge) )
        {
            print("[augmentedEdgeForEdge] Passed invalid edge, exiting.");
            return();
        }
    }
            
    // Recall that an edge has the format (S,T,L,W) and an augmented edge has the format
    // (S,T,L,W,B,C) where B is an intvec and C a matrix.
    
    // First change ourselves into the ring with only the variables present on the given edge    
    def RRR = basering;
    def nR = ringWithOnlyVars(edgeVariables(edge));
    setring nR;
    
    // Note it's important that we respect the order of the variables in the edge, NOT
    // the natural ordering in the ring.
    list edge = imap(RRR,edge);
    list augEdge = edge;
    
    list L = edgeVariables(edge);
    poly W = edgePotential(edge);
    
    // Build the Jacobi ideal
    ideal J;
    int i;
    for(i=1; i <= size(L); i++)
    {
        J[i] = diff(W,L[i]);
    }
    ideal Jstd = std(J);
    
    list B;
    matrix C;
    for(i=1; i <= size(L); i++)
    {
        // Find a power of the variable va belonging to the ideal generated by the partials of W
        poly va = L[i];
        
        int p;
        int power = 0;
        
        for(p=1; p <= powerCutoff; p++)
        {
            if( reduce(va^p, Jstd) == 0 )
            {
                power = p;
                B[i] = p;
                matrix T = lift(J, va^p);
                
                if( ncols(T) != 1 || nrows(T) != size(J) )
                {
                    print("[augmentedEdgeForEdge] Something wrong with the lift, exiting.");
                    return(0);
                }
                
                if( i == 1 )
                {
                    C = T;
                }
                else
                {
                    C = concat(C,T);
                }
                
                break;
            }
        }
        
        if( power == 0 )
        {
            print("[augmentedEdgeForEdge] Failed to find power of variable in Jacobi ideal, exiting.");
            return(0);
        }
    }
    augEdge[5] = B;
    augEdge[6] = C;
    setring RRR;
    list augEdge = imap(nR, augEdge);
    ideal J = imap(nR,J);
    
    // Perform sanity check
    if( useSanityChecks )
    {
        matrix Q = matrix(J) * augEdge[6];
        list edgeVar = edgeVariables(edge);
        
        int i;
        for(i=1; i <= size(edgeVar); i++)
        {
            if( Q[1,i] != edgeVar[i]^(augEdge[5][i]) )
            {
                print("[augmentedEdgeForEdge] Failed to correct compute augmented edge correctly, exiting.");
                return(0);
            }
        }
    }
    
    return(augEdge);
}

/////////////////////////////////////////////////
// isEdgeInternal

proc isEdgeInternal(list edge)
{
    if( edgeSource(edge) != 0 && edgeTarget(edge) != 0 )
    {
        return(1); 
    }
    
    return(0);
}

/////////////////////////////////////////////////
// webVariables
//
// Returns the list consisting of all variables lying
// on all edges in the web.

proc webVariables(list web)
{
    list edgeList = webEdges(web);
    list vars;
    
    int i;
    for(i=1; i<=size(edgeList); i++)
    {
        list edge = edgeList[i];
        vars = vars + edgeVariables(edge);
        kill edge;
    }
    
    kill edgeList;
    
    return(vars);
}

/////////////////////////////////////////////////
// webInternalVariables
//
// Returns the list consisting of all internal variables

proc webInternalVariables(list web)
{
    list edgeList = webEdges(web);
    list vars;
    
    int i;
    for(i=1; i<=size(edgeList); i++)
    {
        list edge = edgeList[i];
        if( isEdgeInternal(edge) )
        {
            vars = vars + edgeVariables(edge);
        }
        kill edge;
    }
    
    kill edgeList;
    return(vars);
}

/////////////////////////////////////////////////
// webExternalVariables
//
// Returns the list consisting of all external variables

proc webExternalVariables(list web)
{
    list edgeList = webEdges(web);
    list vars;
    
    int i;
    for(i=1; i<=size(edgeList); i++)
    {
        list edge = edgeList[i];
        if( !isEdgeInternal(edge) )
        {
            vars = vars + edgeVariables(edge);
        }
        kill edge;
    }
    
    kill edgeList;
    return(vars);
}

/////////////////////////////////////////////////
// subwebWithVertices
//
// Returns the subweb consisting of the given vertices in the original web
// together with all edges adjacent to those vertices, EVEN THOSE ENDING ON
// VERTICES NOT IN THE LIST. If an edge begins or ends on a vertex in our
// list, but does not end (resp. begin) on a vertex in our list (or 0, which
// is implicitly in every subweb) then we consider it to have ended or begun
// on the boundary; so that in forming the subweb we "shrink" the boundary until
// it cuts across all edges going outside our set of vertices.
//
// NOTE: You should not pass 0 as a vertex to any routines.

proc subwebWithVertices(list web, list vertices)
{
    int useSanityChecks = blowFlags("sanity_checks");
    
    if( size(vertices) == 0 )
    {
        print("[subwebWithVertices] Given no vertices in list, exiting.");
        return();
    }
    
    int useGrading = webIsGraded(web);
    
    // The easy part is building the list of MFs and gradings, so let us do that first
    list MFlist = webMFs(web);
    if( useGrading ){ list MFgrlist = webMFgr(web); }
    list subMFlist;
    if( useGrading ){ list subMFgrlist; }
    
    int i;
    for(i=1; i<=size(vertices); i++)
    {
        subMFlist = subMFlist + list(MFlist[vertices[i]]);
        if( useGrading ){ subMFgrlist = subMFgrlist + list(MFgrlist[vertices[i]]); }
    }
    
    // Now we need to build up the edges. First of all we discard any edges not
    // beginning or ending at one of the vertices in our list.
    // Then in each of these edges replace the integer vertices[i] by i.
    
    list edgeList = webEdges(web);
    list subEdgeList;
    
    for(i=1; i<=size(edgeList); i++)
    {
        list edge = edgeList[i];
        if( checklist(vertices,edgeSource(edge)) || checklist(vertices,edgeTarget(edge)) )
        {
            // Add this edge to our list. The only modification we make is to set 
            // the source or target to 0 if it was previously a vertex not in our list
            if( !checklist(vertices,edgeSource(edge)) )
            {
                edge = setEdgeSource(edge, 0);
            }
            
            if( !checklist(vertices,edgeTarget(edge)) )
            {
                edge = setEdgeTarget(edge, 0);
            }
            
            subEdgeList = subEdgeList + list(edge);
        }
        kill edge;
    }
    
    list finalEdgeList;
    
    // Now go through and modify all the edges in subEdgeList.
    int j;
    for(j=1; j<=size(subEdgeList); j++)
    {
        list edge = subEdgeList[j];
        int fixedSource;
        int fixedTarget;
        for(i=1; i<=size(vertices); i++)
        {
            
            if( edgeSource(edge) == vertices[i] && !fixedSource )
            {
                edge = setEdgeSource(edge, i);
                
                // After we change the source, in a later i iteration a
                // vertex could recognise us as having them as the source,
                // so to avoid this we use a flag.
                fixedSource = 1;
            }
            
            if( edgeTarget(edge) == vertices[i] && !fixedTarget )
            {
                edge = setEdgeTarget(edge, i);
                fixedTarget = 1;
            }
            
        }
        finalEdgeList = finalEdgeList + list(edge);
        kill edge;    
    }
    kill subEdgeList;
    
    list subweb = size(vertices), finalEdgeList, subMFlist;
    kill finalEdgeList, subMFlist;    
    if( useGrading ){ subweb = subweb + list(subMFgrlist); kill subMFgrlist; }
    
    if( useSanityChecks )
    {
        if( !webVerify(subweb) )
        {
            print("[subwebWithVertices] Produced invalid subweb, exiting.");
            return();
        }
    }
    
    return(subweb);
}

/////////////////////////////////////////////////
// compilationStratVerify
//
// Verify that the given list is a valid compilation strategy

proc compilationStratVerify(list web, list compStrat)
{
    int numVertices = webSize(web);
    list edges = webEdges(web);

    list vertexOrder = compStrat[1];

    if( size(vertexOrder) != numVertices )
    {
        dbprint(printlevel, "[compilationStratVerify] Wrong number of vertices in compilation strategy, failed.");
        return(0);
    }
    
    // For each 1 <= i <= n-1 check that one of the first i vertices is connected
    // to the (i+1)st vertex by an edge (vertices indexed according to the comp strat)
    int i;
    for(i=1;i<numVertices;i++)
    {
        int k;
        int foundEdge = 0;
        
        for(k=1;k<=i;k++)
        {
            int j;
            for(j=1; j<=size(edges); j++)
            {
                if( isEdgeBetweenVertices(edges[j],vertexOrder[k],vertexOrder[i+1]) )
                {
                    // Found an edge between one of the first i vertices and the (i+1)st
                    foundEdge = 1;
                    break;
                }
            }
            
            if( foundEdge ){ break; }
        }
        
        if( !foundEdge )
        {
            //dbprint(printlevel, "[compilationStratVerify] Found no edges between the first " + string(i) + " vertices and vertex " + string(i+1) + ", failing.");
            return(0);
        }
    }
    
    return(1);
}

////////////////////////////////////////////////////////////////////
// permuteCompStrat
//
// Tries permuting the order of the vertices in a given compilation
// strategy to find a valid compilation strategy. Since we check
// all possible permutations, this could be very slow!

proc permuteCompStrat(list web, list compStrat)
{
    int numVertices = webSize(web);
    list vertexOrder = compStrat[1];
    list varOrder = compStrat[2];

    list perms = SGroupintvecs(numVertices);
    
    int i;
    for(i=1;i<=size(perms);i++)
    {
        intvec p = perms[i];
        
        list newV;
        int j;
        for(j=1;j<=numVertices;j++)
        {
            newV[j] = vertexOrder[p[j]];
        }
        
        list testStrat = newV, varOrder;
        
        if( compilationStratVerify(web, testStrat) )
        {
            // We found a valid permutation, return it
            dbprint(printlevel, "[permuteCompStrat] Found a valid permutation of the vertices: " + string(p));
            return(testStrat);
        }
        
        kill newV, p, testStrat;
    }
    
    // We failed to find a valid permutation
    dbprint(printlevel, "[permuteCompStrat] Failed to find a valid permutation.");
    return();
}

/////////////////////////////////////////////////
// webDual
//
// Returns the dual web, obtained by reversing the direction
// of all edges, multiplying by -1 the potential at each vertex
// and replacing every MF by its dual. We take care of the gradings
// as well, if they are present.

proc webDual(list web)
{
    int useSanityChecks = blowFlags("sanity_checks");
    
    // Extract the content of the web
    int numVertices = webSize(web);
    list edgeList = webEdges(web);
    list mfList = webMFs(web);
    
    int useGrading = webIsGraded(web);
    if( useGrading ){ list mfGradings = webMFgr(web); }
    
    // First reverse all edges
    list dualEdgeList;
    int i;
    for(i=1;i<=size(edgeList);i++)
    {
        list e = edgeList[i];
        int s = e[1];
        int t = e[2];
        
        e[1] = t;
        e[2] = s;
        
        dualEdgeList = dualEdgeList + list(e);
    }

    // Now dualise all matrices
    list dualMfList;
    list dualMfGradings;
    for(i=1;i<=size(mfList);i++)
    {
        matrix D = MFdual(mfList[i]);
        dualMfList = dualMfList + list(D);
        kill D;
        
        if( useGrading )
        {
            intvec v = MFdualGrading(mfGradings[i]);
            dualMfGradings = dualMfGradings + list(v);
            kill v;
        }
    }
    
    list retWeb = numVertices, dualEdgeList, dualMfList;
    if( useGrading ){ retWeb = retWeb + list(dualMfGradings); }
    
    if( useSanityChecks )
    {
        if( !webVerify(retWeb) )
        {
            dbprint(printlevel, "[webDual] Failed to create valid dual web.");
            return();
        }
    }
    
    return(retWeb);
}

/////////////////////////////////////////////////
// webVerify
//
// Verify that the given list satisfies (some of) the conditions to be a web.

proc webVerify(list web)
{
    if( typeof(web[1]) != "int" || typeof(web[2]) != "list" || typeof(web[3]) != "list" )
    {
        print("[webVerify] Entries in web are not of type int, list, list, failed.");
        return(0);
    }
    
    int useGrading = webIsGraded(web);
    int numVertices = webSize(web);
    
    if( numVertices != size(webMFs(web)) )
    {
        print("[webVerify] There is size mismatch between vertices and MFs, failed.");
        return(0);
    }
    
    int i;
    for( i=1; i<=numVertices; i++ )
    {
        matrix A = webMFatVertex(web, i);
        poly W = vertexPotential(web, i);
        if( useGrading ){ intvec Agr = webMFgratVertex(web, i); }
        
        //if( size(vertexIncomingEdges(web,i)) + size(vertexOutgoingEdges(web,i)) == 0 )
        //{
        //    print("[webVerify] Vertex " + string(i) + " is isolated, failed.");
        //    return(0);
        //}
        
        // Is A a MF of this potential?
        if( A * A != W * unitmat(ncols(A)) )
        {
            print("[webVerify] At vertex " + string(i) + " the MF does not factorise the local potential W=" + string(W) + ", failed.");
            return(0);
        }
        
        // Is A a graded MF?
        if( useGrading )
        {
            if( !isGradingValid(A, Agr) )
            {
                print("[webVerify] At vertex " + string(i) + " the grading is not valid, failed.");
                return(0);
            }
        }
    }
    
    return(1);
}

/////////////////////////////////////////////////
// webCompilePair
//
// The backbone of our compilation routines. We are given a web with only
// two vertices and a compilation strategy. Next we describe the output in
// the standard case, but be aware there is an important special case to be
// explained at the end.
//
// We first return a matrix factorisation in the set of external variables
// of the total potential of the given web. Let X denote the tensor product
// over all vertices of the associated MFs, defined in the ordering determined
// by our compilation strategy. Let pX denote the pushforward/restriction of
// scalars to the ring with all internal variables in the web deleted. Our
// first return value is the differential of a MF Y homotopy equivalent to pX.
//
// The second return value is the grading vector for Y, if we are using gradings.
//
// The third (or second, if there is no grading) return value is a list
// R = (F1,G1,x1,exp1),(F1,G2,x2,exp2),... whose length is equal to the number of
// internal variables in the web. These are the maps returned by mfPushforwardInductive
// applied to the set of internal variables ordered according to our compilation strategy.
//
// In more detail: let p1X denote the pushforward of X just to the ring without the
// "first" internal variable, where "first" means first in the compilation strategy.
// This is a direct summand of Xblow[1], and (F1,G1) is a pair splitting, in the
// homotopy category of finite rank MFs, the corresponding idempotent. Say the
// splitting is Y1. Then F1: Xblow[1] -> Y1 and G1: Y1 -> Xblow[1] with F1 * G1 = 1
// and G1 * Y1 the idempotent corresponding to p1X.
//
// Next, we push forward Y1 to the ring without the first and second internal variables,
// and the same yoga produces maps (F2,G2) splitting an idempotent on Y1blow[1], etc.
//
// This list R is referred to as a "splitting record".
//
// NOTE: Special case
//
// If the total potential of the web is zero, and if in the last step of the compilation
// process we end up with non-square matrices (may not happen) then it is not safe to
// return the differential on Y (it will not be square) so our return list looks like
//
// "nonsquare", [r1,r2,] R
//
// where r1,r2 are grading vectors on the even and odd parts of Y resp. (if we are
// using the grading) and R is as described above, but where the last (F,G,x,exp) pair is
// a tuple f0,f1,g0,g1,x,exp with f0,g0 splitting the even part of the very last idempotent
// and similarly for f1,g1.

proc webCompilePair(list web, list compStrat)
{
    int useSanityChecks = blowFlags("sanity_checks");
    int useGrading = webIsGraded(web);
    
    if( useSanityChecks )
    {
        if( !webVerify(web) )
        {
            print("[webCompilePair] Given web is not valid, exiting.");
            return();
        }
        
        if( webSize(web) != 2 )
        {
            print("[webCompilePair] We only compile pairs. Duh. Exiting.");
            return();
        }
    }
    
    // Extract the content of the web
    int numVertices = webSize(web);
    list edgeList = webEdges(web);
    list mfList = webMFs(web);
    if( useGrading ){ list mfGradings = webMFgr(web); }
    
    // And the compilation strategy
    list vertexOrder = compStrat[1];
    list varOrder = compStrat[2];
    
    if( size(vertexOrder) != 2 )
    {
        print("[webCompilePair] We only compile pairs, exiting.");
        return();
    }
    
    // Order the vertices according to the compilation strategy
    int i;
    for(i=1;i<=2;i++)
    {
        matrix mf(i) = mfList[vertexOrder[i]];
        if( useGrading ){ intvec mfgr(i) = mfGradings[vertexOrder[i]]; }
    }
    
    // Our task is to build up the input to mfPushforwardInductive and
    // then interpret its output.
    matrix A = MFtensor( mf(1), mf(2) );
    if( useGrading ){ intvec Agr = MFtensorGradings( mfgr(1), mfgr(2) ); }
    
    // Build up a list of the internal variables, ordered according to the
    // compilation strategy, a list of powers and a list of homotopies. The
    // powers we collect are the powers of the variables which act null-homotopically
    // on A, with associated homotopies stored in H.
    
    // NOTE: As opposed to the original webCompile, we have to work a little harder
    // to get our homotopies, since a priori we only know that certain partial
    // derivatives of the potentials on edges annihilate, via partial derivatives of
    // suitable differentials, and we have to convert these to null-homotopies of 
    // powers of the internal variables using augmentedEdgeForEdge, which tells us 
    // which powers of the variables on an edge can be written in terms of the partial
    // derivatives of the potential on the edge.
    
    list augEdgeList;
    for(i=1; i<=size(edgeList); i++)
    {
        list edge = edgeList[i];
        list augEdge = augmentedEdgeForEdge(edge);
        augEdgeList = augEdgeList + list(augEdge);
                
        kill augEdge, edge;
    }
    
    list intvars;
    list varPowers;
    list Hlist;
    
    for(i=1;i<=size(varOrder);i++)
    {
        poly u = varOrder[i];
        
        // Find out which edge this variable lies on
        int j;
        for(j=1; j<=size(edgeList); j++)
        {
            list edge = edgeList[j];
            list augEdge = augEdgeList[j];
            
            if( checklist(edgeVariables(edge), u) )
            {
                // Yes, u does lie on this edge
                if( !isEdgeInternal(edge) )
                {
                    print("[webCompilePair] Internal variable lying on external edge, exiting.");
                    print(edge);
                    return();
                }
                
                intvars = intvars + list(u);
                
                poly W = edgePotential(edge);
                int src = edgeSource(edge);
                list varList = edgeVariables(edge);
                list B = augedgePowers(augEdge);
                matrix C = augedgeTransform(augEdge);
                matrix D = mfList[src];
                
                // Find out where in varList u occurs
                int k;
                int uPos;
                for(k=1; k<=size(varList); k++)
                {
                    if( varList[k] == u )
                    {
                        uPos = k;
                        break;
                    }
                }
                
                // Note that D factorises some potential W' with W as a summand, and that
                // for any variable x in the list varList, i.e. a variable on edge, the
                // variable x only occurs in W' as part of W. Hence if we take the
                // partial derivative of D * D = W' we find that
                //
                //   diff(D,x) * D + D * diff(D,x) = diff(W,x).
                //
                // If we sum over the different variables and multiply by the appropriate
                // factors in C, then on the right hand side we have a power of u.
                
                // Recall that by construction
                //
                // u^B[uPos] = C[1,uPos] * diff(W,varList[1]) + ... + C[n,uPos] * diff(W,varList[n]).
                                
                matrix H[nrows(D)][ncols(D)];
                for(k=1; k<=size(varList); k++)
                {
                    matrix h = diff(D,varList[k]);
                    
                    if( h * D + D * h != diff(W,varList[k]) * unitmat(ncols(D)) )
                    {
                        print("[webCompilePair] Inner homotopy construction failure, exiting.");
                        return();
                    }
                    
                    H = H + C[k,uPos] * h;
                    
                    kill h;
                }
                                
                // Convert this from a homotopy on D to a homotopy on A
                if( vertexOrder[1] == src )
                {
                    H = ZZtensor( H, unitmat(ncols(mf(2))) );
                }
                else
                {
                    H = ZZtensor( unitmat(ncols(mf(1))), H );
                }
                
                if( useSanityChecks )
                {
                    if( H * A + A * H != u^(B[uPos]) * unitmat(ncols(A)) )
                    {
                        print("[webCompilePair] Homotopy construction failure, exiting.");
                        return();
                    }
                }
                
                varPowers = varPowers + list(B[uPos]);
                Hlist = Hlist + list(H);
                kill W, src, B, H, varList;
            }
            
            kill edge, augEdge;
        }
    }
    
    // Now we feed all this to mfPushforwardInductive
    if( !useGrading )
    {
        list L = mfPushforwardInductive(A, intvars, varPowers, Hlist);
    }
    else
    {
        list L = mfPushforwardInductive(A, intvars, varPowers, Hlist, Agr);
    }
    kill A, intvars, varPowers, Hlist;
    
    return(L);
}

/////////////////////////////////////////////////
// webCompile
//
// Given a web and compilation strategy we crawl the web, fusing
// vertices in the order dictated by the strategy using webCompilePair.
//
// Our first output is the differential on the compiled MF Y, h.e. to the pushforward
// of the tensor product of all the MFs in the web to the ring without any of the internal
// variables.
//
// If grading is in use, then the second return is the grading vector on this Y.
//
// The third (or second, if grading is not in use) return value is a list Q, which we
// refer to as the "splitting record". Suppose that we write X1, .., Xv for the MFs lying
// at the vertices in our web, ordered according to the compilation strategy. Let
// X1 # X2 denote the result of tensoring X1,X2 and replacing this by a finite rank MF
// homotopy equivalent to the pushforward which eliminates the variables lying between
// vertices 1,2. This reduction generates a splitting record R1 (as defined above) and the
// pair (X1 # X2, R1) is the first entry in the list Q (or (X1 # X2, gr, R1) if we are 
// in the graded setting, gr being the grading vector on X1 # X2).
//
// The second entry in Q is ((X1 # X2) # X3, R2) where R2 is the splitting record of the
// second fusion. All fusions being conducted according to the given compilation strategy.
// The last entry in Q is (Y,Ygr,R?), i.e. we repeat the differential, grading of Y but
// in addition there is the splitting record of the reduction which produced Y.
//
// NOTE: Special case
//
// If the total potential of the web is zero, and if in the last step of the compilation
// process we end up with non-square matrices (may not happen) then it is not safe to
// return the differential on Y (it will not be square) so our return list looks like
//
// "nonsquare", [r1,r2,]  Q
//
// where r1,r2 are grading vectors on the even and odd parts of Y resp. (if we are
// using the grading) and Q is the global splitting record. The last entry in Q is
// a tuple ("nonsquare",[r1,r2,],R) with the last entry of R a tuple (f0,f1,g0,g1,u,NN)
// with f0,g0 splitting the even part of the very last idempotent and similarly for f1,g1.

proc webCompile(list web, list compStrat)
{
    int useSanityChecks = blowFlags("sanity_checks");
    int useGrading = webIsGraded(web);
    
    system("--ticks-per-sec",1000);
    int timeElapsed = timer;
    
    if( useSanityChecks )
    {
        if( !webVerify(web) )
        {
            print("[webCompile] Given web is not valid, exiting.");
            return();
        }
        
        if( !compilationStratVerify(web, compStrat) )
        {
            print("[webCompile] Given compilation strategy is not valid, exiting.");
            return();
        }
    }
    
    def RRR = basering;
    int count = 1;
    def nR(1) = basering;
    list splittingRecord;
    int numInternalVar = size(compStrat[2]);
    
    // Continue fusing until there are no internal variables left. Each time through the
    // loop we will have deleted some internal variables from our current ring and have
    // changed the web and strategy to reflect the fusion of two vertices in the old web.
    // Recall that compStrat[2] is an ordered list of the internal variables.
    int numVertices;
    
    while( numInternalVar > 0 )
    {
        // Our current ring is nR(count)
        dbprint(printlevel, "[webCompile] Entering round " + string(count) + " of compilation.");
        
        // web is being modified each time through the loop, so update our "pointers"
        numVertices = webSize(web);
        list edges = webEdges(web);
        list mfs = webMFs(web);
        if( useGrading ){ list webGr = webMFgr(web); }
                
        // compStrat is also modified at each step
        list vertexOrder = compStrat[1];
        list varOrder = compStrat[2];
        
        if( size(vertexOrder) < 2 )
        {
            print("[webCompile] vertexOrder has impossible size, exiting.");
            return();
        }
        
        // Construct a subweb and substrategy consisting of the vertices to be fused.
        // We always fuse the first two vertices in the compilation strategy
        list subweb, subwebVar, substrat_var;
        subweb = subwebWithVertices(web, list(vertexOrder[1],vertexOrder[2]));
        subwebVar = webVariables(subweb);

        // Intersect the variables in varOrder with the list of internal vars in the subweb
        int i;
        for(i=1; i<=size(varOrder); i++) 
        {
            poly u = varOrder[i];
            
            int j;
            for(j=1; j<=size(edges); j++)
            {
                list e = edges[j];
                if( checklist(edgeVariables(e),u) && isEdgeBetweenVertices(e,vertexOrder[1],vertexOrder[2]) )
                {
                    substrat_var = substrat_var + list(varOrder[i]);
                    break;
                }
                kill e;
            }
            kill u;
        }
            
        if( size(substrat_var) == 0 )
        {
            // In the current intermediate compilation strategy, the first and second vertex
            // are not linked by any edges: that is, the compilation strategy is not valid.
            print("[webCompile] In the intermediate compilation strategy the first and second vertex are not connected by edges, failing.");
            return();
        }
        
        list substrat = list(list(1,2), substrat_var);
        
        // Change to the ring containing only those variables involved in this fusion
        def tempring = ringWithOnlyVars(subwebVar);
        
        kill subwebVar;

        ///////// BEGIN tempring //////////
        setring tempring;
        list subweb = imap(nR(count),subweb);
        list substrat = imap(nR(count),substrat);
            ////// IMPORTANT
            list l = webCompilePair(subweb, substrat);
            ////////////////
        kill subweb, substrat;
        setring nR(count);
        kill subweb, substrat;
        ///////// END tempring ///////////

        list l = imap(tempring, l);
        kill tempring;
        
        if( typeof(l[1]) == "string" )
        {
            // This is the case where webCompilePair produced non-square matrices. This
            // must be our last step, and l looks like
            //
            // "nonsquare", [r1,r2,] R
            //
            // where r1,r2 are grading vectors on the even and odd parts of Y resp. (if we are
            // using the grading) and R is the splitting record, but where the last (F,G) pair in R is
            // a tuple (f0,f1,g0,g1) with f0,g0 splitting the even part of the very last idempotent
            // and similarly for f1,g1.
        
            kill numVertices, edges, mfs, vertexOrder, varOrder;
            if( useGrading ){ kill webGr; }
        
            setring RRR;
            
            // The first time through nR(1) = RRR so there is no need for imap
            if( count > 1 )
            {
                list l = imap(nR(count),l);
            }
            
            kill nR(count);
            list retlist;
            retlist[1] = "nonsquare";
            if( useGrading ){ retlist = retlist + list(l[2],l[3]); }
            splittingRecord = splittingRecord + list(l);
            retlist = retlist + list(splittingRecord);

            dbprint(printlevel, "[webCompile] total elapsed time " + string(timer-timeElapsed) + "ms.");

            if( useGrading )
            {
                dbprint(printlevel, "[webCompile] even grading: " + string(l[2]));
                dbprint(printlevel, "[webCompile]  odd grading: " + string(l[3]));
            }
            kill l, splittingRecord;

            return(retlist);
        }
        
        // The return from webCompilePair is (newvertex [,g], R) where newvertex is the
        // differential on a MF and g is a possible grading vector, R is the splitting record.
        // First we return to RRR and adjoin R to the global splitting record
        setring RRR;
        
            // The first time through nR(1) = RRR so there is no need for imap
            if( count > 1 )
            {
                list l = imap(nR(count),l);
            }
        
            splittingRecord = splittingRecord + list(l);
            
        setring nR(count);
        
        // Now do some engineering on our web. We need to replace the first two vertices
        // by their fusion, and rearrange variables and the compilation strategy accordingly
        // Making the new list of MFs and gradings is easy, so we do that first
        
        // 1. List of matrices
        list fusedWebMfs = mfs;
        fusedWebMfs[vertexOrder[1]] = l[1];
        fusedWebMfs = delete(fusedWebMfs,vertexOrder[2]);
        
        // 2. Gradings
        if( useGrading )
        {
            list fusedWebGr = webGr;
            fusedWebGr[vertexOrder[1]] = l[2];
            fusedWebGr = delete(fusedWebGr,vertexOrder[2]);
        }
        
        // 3. Edges. Remove all squeezed out internal edges, and redirect
        // other edges which end or begin at vertexOrder[2] so that they
        // now end or begin resp. at vertexOrder[1].

        int i;
        list fusedEdges;
        for(i=1; i<=size(edges); i++)
        {
            list e = edges[i];
            
            // Exclude edges which where squeezed out
            if( !isEdgeBetweenVertices(e,vertexOrder[1],vertexOrder[2]) )
            {
                list newEdge = e;
                
                if( edgeSource(newEdge) == vertexOrder[2] )
                {
                    newEdge = setEdgeSource(newEdge, vertexOrder[1]);
                }
                
                if( edgeTarget(newEdge) == vertexOrder[2] )
                {
                    newEdge = setEdgeTarget(newEdge, vertexOrder[1]);
                }
            
                // Finally, if the edge begins or ends on a vertex appearing after
                // vertexOrder[2] we need to decrement the index, as the fused web
                // has one less vertex and all indices shifted
                if( edgeSource(newEdge) > vertexOrder[2] )
                {
                    newEdge = setEdgeSource(newEdge, edgeSource(newEdge) - 1);
                }
                
                if( edgeTarget(newEdge) > vertexOrder[2] )
                {
                    newEdge = setEdgeTarget(newEdge, edgeTarget(newEdge) - 1);
                }
                
                fusedEdges = fusedEdges + list(newEdge);
                kill newEdge;
            }
        }
        
        list fusedWeb = numVertices - 1, fusedEdges, fusedWebMfs;
        kill fusedEdges, fusedWebMfs;
        if( useGrading ){ fusedWeb = fusedWeb + list(fusedWebGr); kill fusedWebGr; }
        
        if( useSanityChecks )
        {
            if( !webVerify( fusedWeb ) )
            {
                print("[webCompile] Fused web is not valid, exiting.");
                return();
            }
        }
        
        // Next we rearrange the compilation strategy. The vertices is simple: we have
        // replaced the first vertex in the strategy by the fusion, so we just drop the
        // entry of vertexOrder in position 2
        list fusedVertexOrder = delete(vertexOrder,2);
        
        for(i=1;i<=size(fusedVertexOrder);i++)
        {
            if( fusedVertexOrder[i] > vertexOrder[2] )
            {
                fusedVertexOrder[i] = fusedVertexOrder[i] - 1;
            }
        }
        
        list fusedVarOrder;
        
        // Our new set of internal variables is the old one with all the variables we
        // squeezed out removed. Note that here is where we update numInternalVar
        int i;
        for(i=1;i<=size(varOrder);i++)
        {
            poly u = varOrder[i];

            if( !checklist(substrat_var,u) )
            {
                fusedVarOrder = fusedVarOrder + list(u);
            }
            else
            {
                // u is an internal variable which was squeezed out, so
                // reduce our internal variable count appropriately
                numInternalVar--;                
            }
            
            kill u;
        }
        
        list fusedStrat = list(fusedVertexOrder, fusedVarOrder);
        kill fusedVertexOrder, fusedVarOrder;
        
        // Replace web by fusedWeb and compStrat by fusedStrat
        web = fusedWeb;
        compStrat = fusedStrat;
        
        // Do some kills before switching rings.
        kill fusedWeb, fusedStrat;
        kill edges, mfs, vertexOrder, varOrder;
        if( useGrading ){ kill webGr; }
        
        // Now we go to the ring nR(count+1), which is nR(count) with the internal
        // variables that we just squeezed out removed. Recall that substrat_var
        // contains a list of these variables. Note we only change rings if there
        // are internal variables left, otherwise ringWithoutVars might try to create
        // a ring with no variables

        if( numInternalVar > 0 )
        {
            def nR(count+1) = ringWithoutVars(substrat_var);
            setring nR(count+1);

            list web = imap(nR(count),web);
            list compStrat = imap(nR(count),compStrat);
            kill nR(count);

            // Increment the counter, get ready to go again!
            count++;
        }
    }
    
    // We exited after exhausting our internal variables, with valid web and compStrat.
    // It remains to tensor together the MFs now (possibly) lying at isolated vertices
    // TODO
    
    // So the first entry in the mfs list of web is our final answer
    setring RRR;
    
    // The first time through nR(1) = RRR so there is no need for imap
    if( count > 1 )
    {
        list web = imap(nR(count),web);
    }
    
    matrix A = webMFatVertex(web,1);
    list retlist = A;
    kill A;
    
    if( useGrading )
    {
        intvec Agr = webMFgratVertex(web,1);
        retlist = retlist + list(Agr);
        kill Agr;
    }
    
    retlist = retlist + list(splittingRecord);
    kill splittingRecord;
    dbprint(printlevel, "[webCompile] total elapsed time " + string(timer-timeElapsed) + "ms.");
    
    return(retlist);
}

/////////////////////////////////////////////////
// webFromLink
//
// Given a link in the format described in linkhom.lib, and an intvec
// giving a "state', i.e. assigning 0 or 1 to every vertex, we return
// the corresponding web (in the general format of this file).

proc webFromLink(list linkData, intvec state, int N, list #)
{
    // If we are passed a fourth parameter we are being asked for a web to
    // compute reduced KR homology, with the given polynomial as the mark
    // to be eliminated. In this case we omit this edge from the web, and
    // substitute reductionMark = 0 in all MFs.
    int useReduced;
    poly reductionMark = 1;
    if( size(#) > 0 )
    {
        useReduced = 1;
        reductionMark = #[1];
    }
    
    intvec vertexParity = linkData[1];
    list edgeList = linkData[2];
    
    int numVertices = size(edgeList);
    
    // Each edge in the link is a list of four integers. The first two are incoming,
    // and the second two are outgoing. Every edge in a web coming from a link has
    // a single variable with exponent N+1, so we can already populate some of the
    // information on each edge
    int i,j;
    list edges;
    
    // Loop through vertices (in the link data)
    for(j=1; j<=numVertices; j++)
    {
        // Loop through entries of vertex
        intvec v = edgeList[j];
        for(i=1; i<=4; i++)
        {
            list edge;
            if( typeof(edges[v[i]]) == "none" )
            {
                // This edge hasn't been created yet
                edge = 0,0,list(var(v[i])),var(v[i])^(N+1);
            }
            else
            {
                edge = edges[v[i]];
            }
            
            if( i <= 2 )
            {
                // This edge begins at this vertex, set the source
                edge = setEdgeSource(edge, j);
            }
            else
            {
                // This edge ends at this vertex, set the target
                edge = setEdgeTarget(edge, j);
            }
            
            edges[v[i]] = edge;
            kill edge;
        }
        kill v;
    }
    
    // Now build the list of MFs and gradings
    list webMfs;
    list webMfsGrading;

    for(j=1; j<=numVertices; j++)
    {
        intvec e = edgeList[j];
    
        matrix X = mfVertex(state[j], list(var(e[1]),var(e[2]),var(e[3]),var(e[4])), N);
        intvec Xgr = mfVertexGrading(state[j], N);
    
        webMfs = webMfs + list(X);
        webMfsGrading = webMfsGrading + list(Xgr);
        kill X,e,Xgr;
    }

    if( useReduced )
    {
        // Delete the edge with the reduced variable
        int i;
        for(i=1; i<=size(edges); i++)
        {
            if( edgeVariables(edges[i])[1] == reductionMark )
            {
                edges = delete(edges,i);
                break;
            }
        }
        
        // Set reductionMark = 0 in all Mfs
        for(i=1; i<=size(webMfs); i++)
        {
            webMfs[i] = subst(webMfs[i],reductionMark,0);
        }
    }
    
    // Define the web
    list web = numVertices, edges, webMfs, webMfsGrading;
    
    if( !webVerify(web) )
    {
        print("[webFromLink] Failed to construct valid web, exiting.");
        return();
    }
    
    return(web);
}


/////////////////////////////////////////////////
// defectTangleFromLink
//
// Given a link in the format described in linkhom.lib, and an intvec
// giving a "state", i.e. assigning 0 or 1 to every vertex, we return a list whose first entry is
// the web that is produced by cutting the cutpos-th edge in the state web. The second entry of the
// list is a compilation strategy that does not contain the cutpos-th variable.

proc defectTangleFromLink(list linkData, intvec state, int N, int cutpos, list #)
{
    // If we are passed a fifth parameter we are being asked for a web to
    // compute reduced KR homology, with the given polynomial as the mark
    // to be eliminated. In this case we omit this edge from the web, and
    // substitute reductionMark = 0 in all MFs.
    int useReduced;
    poly reductionMark = 1;
    if( size(#) > 0 )
    {
        useReduced = 1;
        reductionMark = #[1];
    }
    
    intvec vertexParity = linkData[1];
    list edgeList = linkData[2];    
    int numVertices = size(edgeList);
    
    list compStratVert;
    list compStratVar;
    int i,j;
    for( j=1; j<=size(edgeList); j++ )
    {
        compStratVert[j] = j;
    }
    
    for( j=1; j<=nvars(basering)-1; j++ )
    {
        if( j != cutpos )
        {   
            compStratVar = compStratVar + list(var(j));
        }
    }    
    list compStrat = compStratVert,compStratVar;
    
    // Each edge in the link is a list of four integers. The first two are incoming,
    // and the second two are outgoing. Every edge in a web coming from a link has
    // a single variable with exponent N+1, so we can already populate some of the
    // information on each edge
    list edges;
    
    // Loop through vertices
    for(j=1; j<=numVertices; j++)
    {
        // Loop through entries of vertex
        intvec v = edgeList[j];
        for(i=1; i<=4; i++)
        {
            list edge;
            if( typeof(edges[v[i]]) == "none" )
            {
                // This edge hasn't been created yet
                edge = 0,0,list(var(v[i])),var(v[i])^(N+1);
            }
            else
            {
                edge = edges[v[i]];
            }
            
            if( i <= 2 )
            {
                // This edge begins at this vertex, set the source
                edge = setEdgeSource(edge, j);
            }
            else
            {
                // This edge ends at this vertex, set the target
                edge = setEdgeTarget(edge, j);
            }
            
            edges[v[i]] = edge;
            kill edge;
        }
        kill v;
    }
    
    // The list edges, as defined above, describes the uncut web. Now we cut the cutpos-th edge in
    // two pieces, leaving two loose ends which we associate with the boundary vertex 0. The new edge 
    // whose source is 0 will be labelled by n+1, where n is the number of edges in the original web.
    
    int n = size(edges);
    int target = edges[cutpos][2];
    edges[n+1] = edges[cutpos];

    // Set the target of the cutpos-th edge to 0, set the source of the (n+1)-th edge to 0, and
    // label this last edge with the last variable:
    edges[cutpos][2] = 0;
    edges[n+1][1] = 0;
    edges[n+1][3] = list(var(n+1));
    edges[n+1][4] = var(n+1)^(N+1);
        
    // Now build the list of MFs and gradings
    list webMfs;
    list webMfsGrading;

    for(j=1; j<=numVertices; j++)
    {
        intvec e = edgeList[j];
        
        // If we are looking at the target of the edge cutpos, we need to replace the associated
        // variable by the "last" variable in the MF of vertex j:
        if( j == target )
        {
            for( i=1; i<=4; i++ )
            {
                if( e[i] == cutpos )
                {
                    e[i] = n+1;
                }
            }
        }
    
        matrix X = mfVertex(state[j], list(var(e[1]),var(e[2]),var(e[3]),var(e[4])), N);
        intvec Xgr = mfVertexGrading(state[j], N);
    
        webMfs = webMfs + list(X);
        webMfsGrading = webMfsGrading + list(Xgr);
        kill X,e,Xgr;
    }

    if( useReduced )
    {
        // Delete the edge with the reduced variable
        int i;
        for(i=1; i<=size(edges); i++)
        {
            if( edgeVariables(edges[i])[1] == reductionMark )
            {
                edges = delete(edges,i);
                break;
            }
        }
        
        // Set reductionMark = 0 in all Mfs
        for(i=1; i<=size(webMfs); i++)
        {
            webMfs[i] = subst(webMfs[i],reductionMark,0);
        }
    }
    
    // Define the web
    list web = numVertices + 1, edges, webMfs, webMfsGrading;
    
    // Don't check this because we have one vertex more (namely the boundary vertex 0) than there 
    // are MFs:
    //    if( !webVerify(web) )
    //    {
    //        print("[defectTangleFromLink] Failed to construct valid web, exiting.");
    //        return();
    //    }
    
    list returnlist = web, compStrat;
    return(returnlist);
}



/////////////////////////////////////////////////
// webCompileMorphism
//
// Given a pair of webs, which we assume only differ with regards to the matrix factorisations
// (and possibly their gradings) assigned to the vertices (i.e. we assume that the vertices are
// the same and in the same order), a compilation strategy, and a morphism F from the matrix
// factorisation at vertex mapPos in the first web to the matrix factorisation at the same
// vertex in the second web, we return the induced morphism between the compilations of the
// two webs. We assume that the matrix factorisations in the web agree except possibly at
// the specified vertex mapPos.
//
// NOTE: mapPos is an index into the absolute ordering the vertices, as defined in the
// link encoding, NOT an index into the compilation strategy.
//
// Recall that the first output of webCompile is the differential on the "compilation" MF,
// which we earlier referred to as Y. Let Y1 denote the compilation of web1 and Y2 the
// compilation of web2. Our first return value here is a even block matrix C representing
// the morphism Y1 --> Y2 compiling F.
//
// In the special case where the potential is zero and we end up with nonsquare matrices,
// we return not a block matrix but a pair (c0,c1) which give the compiled map in each
// degree.
//
// You can distinguish the two cases by the fact that typeof( ) in the first case is matrix
// and in the second case is list.

proc webCompileMorphism(list web1, list web2, list compStrat, matrix M, int mapPos)
{
    list compile(1) = webCompile(web1, compStrat);
    list compile(2) = webCompile(web2, compStrat);
    
    def A = webCompileMorphismFromData(web1, compile(1), web2, compile(2), compStrat, M, mapPos);
    
    kill compile(1), compile(2);
    
    return(A);
}

proc webCompileMorphismFromData(list web1, list compile(1), list web2, list compile(2), list compStrat, matrix M, int mapPos)
{
    int useSanityChecks = blowFlags("sanity_checks");
    int linkHomologyChecks = blowFlags("link_homology_checks");
    
    int useGrading = webIsGraded(web1);
    
    list stratVertices = compStrat[1];
    
    int numVertices = webSize(web1);
    
    list mfs(1) = webMFs(web1);
    list mfs(2) = webMFs(web2);
    if( useGrading )
    {
        list mfsgr(1) = webMFgr(web1);
        list mfsgr(2) = webMFgr(web2);
    }
    
    if( useSanityChecks )
    {
        if( mapPos < 1 or mapPos > numVertices )
        {
            print("[webCompileMorphismFromData] mapPos invalid range, exiting.");
            return();
        }
    }
    
    // Figure out which index in the compilation strategy corresponds to mapPos
    int mapPosStrat;
    int i;
    for(i=1;i<=numVertices;i++)
    {
        if( stratVertices[i] == mapPos )
        {
            mapPosStrat = i;
            break;
        }
    }
    
    // Now we walk through the compilation strategy. Here's roughly what's going to happen.
    // Suppose our mfs are mf(1), ..., mf(v), ordered according to the compilation strategy.
    // This means that, somewhat confusingly, mf(i) is mfs(1)[stratVertices[i]].
    // What the compilation process does is tensor mf(1) with mf(2), tensor this with a 
    // finite-dimensional algebra to cut down to the ring without any of the variables that
    // lie between the vertices at which mf(1), mf(2) sit (we call this "inflation") reduce this 
    // inflation, and split an idempotent on this reduction. Let us call the result of this
    // splitting mf(1,2). Note the three steps
    //
    //  1. Inflate, using mablow (let us call this B for blow)
    //  2. Reduce (giving a homotopy equivalence between B and its reduction)
    //  3. Split an idempotent (gives a pair of maps splitting the idempotent).
    //
    // In the end we have a pair F: B[1] --> mf(1,2) and G: mf(1,2) --> B[1] such that
    // F * G = 1 and G * F is the idempotent splitting to give the pushforward of mf(1) * mf(2)
    // (in the infinite category) (these equalities are up to homotopy).
    //
    // It is this pair (F,G) which is recorded in the so-called "splitting records" which
    // are returned by webCompile (see mfPushforward, mfPushforwardInductive and webCompilePair
    // for how these things are put together).
    //
    // These steps are then repeated, to obtain mf(1,2,3) from mf(1,2) * mf(3), and so on.
    //
    // Suppose for the sake of simplicity that mapPosStrat = 3, and write mf'(i) for the MFs occurring
    // in web2. Hence M is a map mf(3) --> mf'(3). Then we have an induced map
    //
    // mf(1,2) * mf(3) --> mf(1,2) * mf'(3).
    //
    // Inflating this using mablow, suspending, and precomposing with the appropriate G and
    // postcomposing with the appropriate F we get a map
    //
    // mf(1,2,3) --> mf'(1,2,3)
    //
    // Tensoring with mf(4) = mf'(4) and inflating, pre and post-composing, etc. we continue
    // to keep track of this map until at the end we have the desired map
    //
    // C: mf(1,..,v) --> mf'(1,..,v).
        
    if( useSanityChecks )
    {
        // First check that apart from mapPos the mfs agree.
        int i;
        for(i=1; i<=numVertices; i++)
        {
            if( i != mapPos && mfs(1)[i] != mfs(2)[i] )
            {
                print("[webCompileMorphism] Matrix factorisations in position " + string(stratVertices[i]) + " differ, exiting.");
                return();
            }
        }
        
        // Check that M is actually a morphism
        if( M * mfs(1)[mapPos] != mfs(2)[mapPos] * M )
        {
            print("[webCompileMorphism] Given matrix is not a morphism, exiting.");
            return();
        }
    }
    
    // The first task is to fastforward to the mapPosStrat - 1 position in the compilation
    // by looking up the differential on mf(1,2,...,mapPosStrat - 1).
    int nonsquare = (typeof(compile(1)[1]) == "string");
    
    // For nonsquare matrices compile looks like string, [intvec,intvec,] list and for
    // square matrices it looks like matrix, [intvec,] list, with the list being the
    // global splitting record.
    
    list splitRec(1) = compile(1)[2 + (nonsquare + 1)*useGrading];
    list splitRec(2) = compile(2)[2 + (nonsquare + 1)*useGrading];
    
    // See webCompile's introduction for the full details on the format of the things
    // in the splitting record. We just recall here that entries are tuples (Z, [gr,] R)
    // where Z is the differential on some mf(1,..,i) and R is the splitting record for
    // the production of Z from mf(1,..,i-1) * mf(i). Note the white lie above; actually
    // we have the inflate-reduce-split triumvirate for _each_ internal variable, and hence
    // a list R of pairs of F,G maps.
    //
    // Note that the ith position in the splitting record speaks about the fusion of
    // mf(1,...,i) with mf(i+1).
    
    int startPoint;
    if( mapPosStrat > 2 )
    {
        list record = splitRec(1)[mapPosStrat - 2];
        
        // The compilations of web1 and web2 should agree up to this point
        if( useSanityChecks )
        {
            if( record[1] != splitRec(2)[mapPosStrat - 2][1] )
            {
                print("[webCompileMorphism] Expected equal compilations, but this was not true. Exiting.");
                return();
            }
        }
        
        // Replace M by mf(1,..,mapPosStrat-1) * mf(mapPosStrat) --> mf(1,..,mapPosStrat-1) * mf'(mapPosStrat)
        M = ZZtensor( unitmat(ncols(record[1])), M );
        startPoint = mapPosStrat - 1;
    }
    
    if( mapPosStrat == 2 )
    {
        // Replace M by mf(1) * mf(2) --> mf(1) * mf'(2).
        M = ZZtensor( unitmat(ncols(mfs(1)[stratVertices[1]])), M );
        startPoint = 1;
    }
    
    if( mapPosStrat == 1 )
    {
        // Replace M by mf(1) * mf(2) --> mf'(1) * mf(2)
        M = ZZtensor( M, unitmat(ncols(mfs(1)[stratVertices[2]])) );
        startPoint = 1;
    }
    
    // Now we're ready to step through the splitting record
    int i;
    for(i=startPoint;i<=size(splitRec(1));i++)
    {
        // Currently M is a map mf(1,...i) * mf(i+1) --> mf'(1,...,i) * mf'(i+1).
        
        // Check whether we are in the step with nonsquare matrices
        if( typeof(splitRec(1)[i][1]) == "string" )
        {
            // This should be the last splitting in the splitting record
            if( i != size(splitRec(1)) )
            {
                print("[webCompileMorphism] Found nonsquare matrices before the end of the splitting record, exiting.");
                return();
            }
            
            // If the second web is not also in this state, something is wrong
            if( typeof(splitRec(2)[i][1]) != "string" )
            {
                print("[webCompileMorphism] Found nonsquare matrices in only one splitting record, exiting.");
                return();
            }
            
            // In this case the local splitting record lives in a different position, as splitRec(1)[i]
            // is a tuple ("nonsquare",[r1,r2,],R) with the last entry of R a tuple (f0,f1,g0,g1,u,N)
            // with f0,g0 splitting the even part of the very last idempotent and similarly for f1,g1.
            
            list R(1) = splitRec(1)[i][2 + 2*useGrading];
            list R(2) = splitRec(2)[i][2 + 2*useGrading];
        }
        else
        {
            // This is the standard case of square matrices
            
            // Extract the R lists for web1 and web2 mentioned above
            list R(1) = splitRec(1)[i][2 + useGrading];
            list R(2) = splitRec(2)[i][2 + useGrading];    
        }
        
        //displaySplittingRecord(R(1));
        
        // A splitting record R looks like a list of tuples (F,G,u,N) where
        // F is a map Ablow[1] --> final, G is a map final --> Ablow[1] and
        // u is the internal variable being squeezed out, NN is the power such that
        // u^NN acts null-homotopically.    
            
        // Unless it is the last step where we get nonsquare matrices, in which case
        // there is a tuple (f0,f1,g0,g1,u,NN).
    
        int j;
        for(j=1;j<=size(R(1));j++)
        {
            // There are two cases: we are either in the main body of the routine,
            // in which the splittings are square, or we are at the last step and
            // we may have nonsquare splittings
            
            if( typeof(R(1)[j][1]) != "string" )
            {
                // Main case. The format of R(1)[j] is (F,G,u,N,D[,r]).
                matrix F = R(2)[j][1];
                matrix G = R(1)[j][2];
                poly u = R(1)[j][3];
                int expo = R(1)[j][4];

                matrix dl = R(1)[j][5]; // differential on the source MF
                matrix dr = R(2)[j][5]; // and the target MF
                
                if( useGrading )
                {
                    intvec grl = R(1)[j][6]; // grading on the source
                    intvec grr = R(2)[j][6]; // and on the target
                }

                // First we need to inflate and suspend.
                matrix Minf = mfInflate(M,u,expo);
                M = mfSuspendMorph( Minf );
                kill Minf;

                // Now premultiply by G and postmultiply by F
                M = F * M * G;

                // Check that M remains a morphism
                if( useSanityChecks )
                {
                    if( dr * M != M * dl )
                    {
                        print("[webCompileMorphismFromData] M is not a morphism in step j = " + string(j) + " and i = " + string(i));
                        return();
                    }
                }

                // Check that M is degree + 1 (this is link homology specific)
                if( linkHomologyChecks )
                {
                    if( !checkMorphismDegree( dl, grl, dr, grr, M, 1 ) )
                    {
                        print("[webCompileMorphismFromData] M is not degree +1 in i = " + string(i) + " j = " + string(j));
                        return();
                    }
                }

                kill u, expo, F, G, dl, dr;  
            }
            else
            {
                // Special case of nonsquare matrices. The format of R(1)[j] is
                // ("nonsquare",f0,f1,g0,g1,intvars[1],N[1][,r1,r2])
                
                // First check this is the last step
                if( j != size(R(1)) || i != size(splitRec(1)) || typeof(splitRec(1)[i][1]) != "string" )
                {
                    print("[webCompileMorphism] Found nonsquare matrices in local splitting record before the end, exiting.");
                    return();
                }
                                        
                matrix f0 = R(2)[j][2]; // Note the F must come from the second web
                matrix f1 = R(2)[j][3];
                matrix g0 = R(1)[j][4]; // and the G from the first
                matrix g1 = R(1)[j][5];
                poly u = R(1)[j][6];
                int expo = R(1)[j][7];
                
                // Inflate and suspend
                matrix Minf = mfInflate(M,u,expo);
                M = mfSuspendMorph( Minf );
                kill Minf;
                    
                list mblock = extractblockmat(M);
                matrix m(0) = mblock[1];
                matrix m(1) = mblock[4];
                
                // Premultiply by G and postmultiply by F
                m(0) = f0 * m(0) * g0;
                m(1) = f1 * m(1) * g1;
                    
                // Check that M has the right dimensions
                if( useSanityChecks && useGrading )
                {
                    // The fact that we are here means that compile(1) = "nonsquare", r1,r2, Q
                    // and compile(2) = "nonsquare", s1,s2, Q' and our maps m0, m1 should therefore
                    // be of size s1xr1 and s2xr2 respectively
                    intvec leftwebgr(0) = compile(1)[2];
                    intvec leftwebgr(1) = compile(1)[3];
                    intvec rightwebgr(0) = compile(2)[2];
                    intvec rightwebgr(1) = compile(2)[3];
                    
                    int z;
                    for(z=0; z<=1; z++)
                    {
                        // Check that this grading agrees with the last step of the global splitting record
                        if( leftwebgr(z) != splitRec(1)[size(splitRec(1))][2 + z]
                            || rightwebgr(z) != splitRec(2)[size(splitRec(2))][2 + z] )
                        {
                            print("[webCompileMorphism] Unexpected grading disagreement, exiting.");
                            return();
                        }
                                                
                        if( nrows(m(z)) != size(rightwebgr(z)) || ncols(m(z)) != size(leftwebgr(z)) )
                        {
                            print("[webCompileMorphism] Size mismatch in nonsquare step, exiting.");
                            print("[webCompileMorphism] m" + string(z) + ": " + string(nrows(m(z))) + "x" + string(ncols(m(z))));
                            print("[webCompileMorphism] Source grading, degree " + string(z) + ": " + string(leftwebgr(z)));
                            print("[webCompileMorphism] Target grading, degree " + string(z) + ": " + string(rightwebgr(z)));
                            return();
                        }
                    }
                }
                    
                list retlist = list(m(0),m(1));
                kill mblock, m(0), m(1), f0, f1, g0, g1, u, expo;
                    
                return(retlist);
            }
        }
        
        // Check that M remains a morphism, after all these splittings
        if( useSanityChecks && typeof(splitRec(1)[i][1]) == "matrix" )
        {
            matrix mfleft = splitRec(1)[i][1];
            matrix mfright = splitRec(2)[i][1];
            
            if( nrows(mfleft) != ncols(M) || nrows(M) != ncols(mfright) )
            {
                print("[webCompileMorphismFromData] Size mismatch in stage " + string(i) + " of the compilation, exiting.");
                print("[webCompileMorphismFromData] mfleft: " + string(nrows(mfleft)) + "x" + string(ncols(mfleft)));
                print("[webCompileMorphismFromData] mfright: " + string(nrows(mfright)) + "x" + string(ncols(mfright)));
                print("[webCompileMorphismFromData] M: " + string(nrows(M)) + "x" + string(ncols(M)));
                return();
            }
            
            if( M * mfleft != mfright * M )
            {
                print("[webCompileMorphismFromData] In stage " + string(i) + " of the compilation M fails to be a morphism, exiting.");
                return();
            }
            
            // Check that M is degree +1 (this is link specific)
            if( linkHomologyChecks )
            {
                intvec mfleftgr = splitRec(1)[i][2];
                intvec mfrightgr = splitRec(2)[i][2];
                
                if( !checkMorphismDegree(mfleft,mfleftgr,mfright,mfrightgr,M,1) )
                {
                    print("[webCompileMorphismFromData] M is not degree +1 in global check for i = " + string(i) + ", exiting.");
                    return();
                }
                
                kill mfleftgr, mfrightgr;
            } 
            
            kill mfleft, mfright;
        }
        
        // Now M is a map mf(1,...,i+1) --> mf'(1,...,i+1). It remains to tensor
        // with the next MF and then we can start again.
        if( i < size(splitRec(1)) )
        {
            M = ZZtensor( M, unitmat(ncols(mfs(1)[stratVertices[i+2]])) );
        }
    }
    
    return(M);
}

////////////////////////////////////////////////////////////////////
// defaultCompStratForWeb
//
// Returns a default compilation strategy for a given web, we order the
// vertices 1,...,V and order the internal variables as given by
// webInternalVariables.

proc defaultCompStratForWeb(list web)
{
    int numVertices = webSize(web);
    list vertexOrder;
    
    int i;
    for(i=1; i<=numVertices; i++)
    {
        vertexOrder[i] = i;
    }
    
    list varOrder = webInternalVariables(web);
    return(list(vertexOrder,varOrder));
}