brick_functions.sage

#### Functions for working with reflections, hyperplanes, cutting, and shards


# INPUT: A quiver.
# OUTPUT: The Cartan matrix associated to the undirected graph obtained by
# forgetting the orientation of the quiver.
def cartan(Q):
    G = Q.to_undirected()
    return 2*identity_matrix(len(G.vertices(sort=True))) - G.adjacency_matrix()

# INPUT: A Cartan matrix
# OUTPUT: A list of matrices representing the reflections associated to the Cartan matrix
def reflections(cartan):
    rank = cartan.nrows()
    reflection_list = []
    for i in range(rank):
        reflection = identity_matrix(QQ, rank) - (cartan * diagonal_matrix(QQ, [j == i for j in range(rank)]))
        reflection_list.append(reflection)
    return reflection_list

# INPUT: A Cartan matrix, an index (representing a simple root alpha_0),
# and a sequence of indices (representing a sequence of reflections s_1, s_2, ..., s_l)
# OUTPUT: If the expression which applies the reflections to the simple root is not positive, throws a ValueError.
# Otherwise, returns a pair consisting of the root s_l s_{l-1} ... s_1(alpha_0) (the "plane", as a 1 x n matrix)
# and a list of truncations s_l s_{l-1} ... s_{t+1}(alpha_t), which represent the roots cutting the hyperplane.
# (the "fractures", as 1 x n matrices)
def plane_and_fractures(cartan, starting_root, reflect_seq):
    rank = cartan.nrows()
    # Base case: no reflections
    if len(reflect_seq) == 0:
        return (Matrix([i == starting_root - 1 for i in range(rank)]), [])
    # Recursive step
    else:
        reflects = reflections(cartan)
        old_fractures = plane_and_fractures(cartan, starting_root, reflect_seq[:-1])
        next_reflection = reflects[reflect_seq[-1] - 1]
        new_plane = old_fractures[0] * next_reflection
        difference = new_plane - old_fractures[0]
        if not difference[0, reflect_seq[-1] - 1] > 0:
            raise ValueError('This is not a positive expression, witnessed by the {0}th reflection.'.format(len(reflect_seq)))
        else:
            new_fractures = [s * next_reflection for s in old_fractures[1]] + [Matrix([i == reflect_seq[-1] - 1 for i in range(rank)])]
            return (new_plane, new_fractures)

# INPUT: A pair (hyperplane H, a list of hyperplanes) where each hyperplane is
# specified by a 1 x n matrix.
# OUTPUT: A list of vectors defining a hyperplane arrangement in R^{n-1}
# isomorphic to the arrangement within H defined by its intersections
# with the other hyperplanes
def restricted_fractures(planes):
    ambient_plane = list(list(planes[0])[0])
    fractures = [list(list(p)[0]) for p in planes[1]]
    scale = max(ambient_plane)
    # find index at which coefficient of H is nonzero
    normal_pivot = ambient_plane.index(scale)
    ambient_plane.pop(normal_pivot)
    ambient_plane_vector_trunc = vector(ambient_plane)
    restricted_fractures = []
    # for each other hyperplane, subtract appropriate multiple of H to zero out
    # previously selected index
    for plane in fractures:
        rescale = plane[normal_pivot] / scale
        plane.pop(normal_pivot)
        restricted_fractures.append(vector(plane) - rescale*ambient_plane_vector_trunc)
    return restricted_fractures

# INPUT: A list of vectors representing equations of hyperplanes
# and a point in the space those hyperplanes occupy
# OUTPUT: If the point lies on a hyperplane, throws a ValueError.
# Otherwise, returns a string of '+' and '-', indicating for each hyperplane
# (in the order they were provided) which side the point lies on
def manual_sign_vector(planes, point):
    vector = ''
    for p in planes:
        eval = (Matrix(p) * point)[0]
        if eval > 0:
            vector += '+'
        elif eval < 0:
            vector += '-'
        else:
            raise ValueError('This point lies on a hyperplane.')
    return vector

# INPUT: A Cartan matrix, an index (representing a simple root alpha_0),
# and a sequence of indices (representing a sequence of reflections s_1, s_2, ..., s_l)
# OUTPUT: A list of sign vectors describing the regions of the arrangement of shards within
# the hyperplane defined by s_l s_{l-1} ... s_1(alpha_0). The jth sign describes
# which side of the jth fracture (as output by plane_and_fractures) the region lies on.
def region_signs(cartan, starting_root, reflect_seq):
    fractures = restricted_fractures(plane_and_fractures(cartan, starting_root, reflect_seq))
    H = HyperplaneArrangements(QQ, tuple('x{0}'.format(i+1) for i in range(cartan.nrows() - 1)))
    arr = H([[tuple(f), 0] for f in fractures])
    sign_vectors = [manual_sign_vector(fractures, r.representative_point()) for r in arr.regions()]
    return sign_vectors

# INPUT: A Cartan matrix, an index (representing a simple root alpha_0),
# a sequence of indices (representing a sequence of reflections s_1, s_2, ..., s_l),
# and a sequence of signs specifying a region of the arrangement of shards
# within the hyperplane defined by s_l s_{l-1} ... s_1(alpha_0), as output by region_signs.
# OUTPUT: A Polyhedron, the dual cone to the specified shard.
def dual_shard(cartan, start_root, reflect_seq, signs):
    p_and_f = plane_and_fractures(cartan, start_root, reflect_seq)
    root = p_and_f[0]
    fracs = p_and_f[1]
    ray_list = []
    # specifying that the shard lies on one side of a hyperplane corresponds to
    # adding a ray dual to that half-plane to the dual cone of the shard
    for i in range(len(fracs)):
        if signs[i] == '+':
            ray_list.append(list(-fracs[i][0]))
        elif signs[i] == '-':
            ray_list.append(list(fracs[i][0]))
    return Polyhedron(rays=ray_list, lines=[list(root[0])])

#### Functions for working with representations of preprojective algebras and reflection functors

# An enriched class for directed graphs in which each edge has a counterpart pointing the other direction.
class DoubleQuiver(sage.graphs.digraph.DiGraph):
    
    # Constructor can take an existing DoubleQuiver (to prevent repeated edge doubling) or a general DiGraph.
    # In the latter case it adds a reverse for every edge. In the new quiver, the original edge is renamed
    # by adding a '0' to the end of its name, while its reverse counterpart is named by adding a '1' to the
    # end of the original edge's name.
    def __init__(self, digraph, multiedges=True, **kwargs):
        if isinstance(digraph, DoubleQuiver):
            super(DoubleQuiver, self).__init__(digraph, multiedges=True, **kwargs)
            self._edge_pairs = digraph.edge_pairs()
        else:
            _edges = []
            self._edge_pairs = {}
            for e in digraph.edges(sort=True):
                self._edge_pairs[(e[0], e[1], e[2] + '0')] = (e[1], e[0], e[2] + '1')
                self._edge_pairs[(e[1], e[0], e[2] + '1')] = (e[0], e[1], e[2] + '0')
                _edges.append((e[0], e[1], e[2] + '0'))
                _edges.append((e[1], e[0], e[2] + '1'))
            super(DoubleQuiver, self).__init__(_edges, multiedges=True, **kwargs)

    # Take an edge as input and return its counterpart.
    def flip(self, edge):
        return self._edge_pairs[edge]

    def edge_pairs(self):
        return self._edge_pairs


# INPUT: A representation of a DoubleQuiver.
# OUTPUT: Whether the representation satisfies the relation necessary to be a representation of the preprojective algebra.
def satisfies_preproj_relation(rep):
    q = rep.quiver()
    for v in q.vertices(sort=True):
        dim = rep.dimension(v)
        # check preprojective algebra relation at vertex v
        if not sum((rep.get_map(q.flip(e)) * rep.get_map(e) * (1 if e[2][-1] == '0' else -1)).matrix() for e in q.outgoing_edge_iterator(v)) == 0:
            return False
    return True


# INPUT: A list of vector spaces and an optional base ring.
# OUTPUT: A tuple (direct sum of spaces, list of inclusion maps as matrices, list of projection maps as matrices)
def direct_sum_with_maps(spaces, ring=QQ):
    # empty sum
    if len(spaces) == 0:
        return (VectorSpace(ring, 0), [], [])

    # proceed recursively, using sum of two vector spaces
    else:
        #ring = spaces[0].base_ring()

        # take sum of all spaces but the last one
        prelim_sum = direct_sum_with_maps(spaces[:-1])

        # take sum of all spaces including the last one
        full_sum = prelim_sum[0].direct_sum(spaces[-1])

        # two useful block matrices: inclusion of prelim_sum into full_sum, and inclusion of last space into full_sum.
        binary_sum_inclusion_maps = (End(prelim_sum[0]).identity().matrix().augment(Hom(prelim_sum[0], spaces[-1]).zero().matrix()), Hom(spaces[-1], prelim_sum[0]).zero().matrix().augment(End(spaces[-1]).identity().matrix()))

        # compose the first of the above with inclusion maps into prelim_sum
        list_of_inclusion_maps = [M * binary_sum_inclusion_maps[0] for M in prelim_sum[1]] + [binary_sum_inclusion_maps[1]]

        # same as above, but for projection maps
        binary_sum_projection_maps = ((End(prelim_sum[0]).identity().matrix().augment(Hom(prelim_sum[0], spaces[-1]).zero().matrix())).transpose(), (Hom(spaces[-1], prelim_sum[0]).zero().matrix().augment(End(spaces[-1]).identity().matrix())).transpose())
        list_of_projection_maps = [binary_sum_projection_maps[0] * M for M in prelim_sum[2]] + [binary_sum_projection_maps[1]]

        return (full_sum, list_of_inclusion_maps, list_of_projection_maps)


# INPUT: A representation of a DoubleQuiver, (the label of) a vertex of the quiver, 
# and whether to require that the representation is in NoQuot (assumed by default)
# OUTPUT: If no_quot_check == True and the representation has S_i as a quotient,
# throws a ValueError. Otherwise, returns the representation obtained by applying 
# the positive reflection functor at that vertex.
def reflection_functor_plus(rep, vertex, no_quot_check=True):
    # check for preprojectivity
    if not satisfies_preproj_relation(rep):
        raise ValueError('This representation does not satisfy the preprojective relation.')

    # extract underlying parameters of representation
    base_ring = rep.base_ring()
    Q = rep.quiver()
    PQ = Q.path_semigroup()
    spaces = {v: rep.get_space(v) for v in Q.vertices(sort=True)}
    maps = {e: rep.get_map(e) for e in Q.edges(sort=True)}


    # pull out the given vertex's space and the arrows / maps incident to the given vertex
    Vi = spaces[vertex]
    in_edges = Q.incoming_edges([vertex])
    maps_in = [maps[e] for e in in_edges]
    # make sure paired edges have the same index
    out_edges = [Q.flip(e) for e in in_edges]
    maps_out = [maps[e] for e in out_edges]

    # sum together all incident spaces
    Vdi = direct_sum_with_maps([spaces[e[0]] for e in in_edges], ring=base_ring)
    
    
    # if the space at the given vertex is 0, there's nothing to do.
    if Vdi[0].dimension() == 0 and rep.get_space(vertex).dimension() == 0:
        return rep

    # concatenate inward maps to get a map Vdi -> Vi
    concat_matrix_in = sum(Vdi[2][i] * maps[in_edges[i]].matrix() * (1 if in_edges[i][2][-1] == '0' else -1) for i in range(len(in_edges)))
    map_in = linear_transformation(Vdi[0], Vi, concat_matrix_in)

    # concatenate outward maps to get a map Vi -> Vdi
    concat_matrix_out = sum(maps[out_edges[i]].matrix() * Vdi[1][i] for i in range(len(out_edges)))
    map_out = linear_transformation(Vi, Vdi[0], concat_matrix_out)

    # if the sum of inward maps is not surjective, then rep has the simple at vertex as a quotient
    if (no_quot_check and not map_in.is_surjective()):
        raise ValueError('This representation has S_{0} as a quotient.'.format(vertex))

    # compute kernel of map Vdi -> Vi
    kernel = map_in.kernel()

    # define new maps out of and into kernel
    # UPDATE: the ".matrix()" at the end of each expression may be unnecessary for Sage's quiver representation
    # constructor. I don't think this messes things up either way, though. Leaving it in for now.
    new_maps_out = {out_edges[i]: linear_transformation(Vdi[0], spaces[out_edges[i][1]], Vdi[2][i]).restrict_domain(kernel).matrix() for i in range(len(out_edges))}
    new_maps_in = {in_edges[i]: linear_transformation(spaces[in_edges[i][0]], Vdi[0], Vdi[1][i] * concat_matrix_in * concat_matrix_out * (1 if in_edges[i][2][-1] == '0' else -1)).restrict_codomain(kernel).matrix() for i in range(len(in_edges))}

    # update space at given vertex
    spaces.update({vertex: kernel})

    # update maps incident to given vertex
    maps.update(new_maps_out)
    maps.update(new_maps_in)

    # return reflected representation
    new_rep = PQ.representation(base_ring, spaces, maps)
    return new_rep


# INPUT: A representation of a DoubleQuiver, (the label of) a vertex of the quiver, 
# and whether to require that the representation is in NoSub (assumed by default)
# OUTPUT: If no_sub_check == True and the representation has S_i as a submodule,
# throws a ValueError. Otherwise, returns the representation obtained by applying 
# the negative reflection functor at that vertex.
def reflection_functor_minus(rep, vertex, no_sub_check=True):
    # check for preprojectivity
    if not satisfies_preproj_relation(rep):
        raise ValueError('This representation does not satisfy the preprojective relation.')

    # extract underlying parameters of representation
    base_ring = rep.base_ring()
    Q = rep.quiver()
    PQ = Q.path_semigroup()
    spaces = {v: rep.get_space(v) for v in Q.vertices(sort=True)}
    maps = {e: rep.get_map(e) for e in Q.edges(sort=True)}


    # pull out the given vertex's space and the arrows incident to the given vertex
    Vi = spaces[vertex]
    out_edges = Q.outgoing_edges([vertex])
    maps_out = [maps[e] for e in out_edges]
    # make sure paired edges have the same index
    in_edges = [Q.flip(e) for e in out_edges]
    maps_in = [maps[e] for e in in_edges]

    # sum together all incident spaces
    Vdi = direct_sum_with_maps([spaces[e[1]] for e in out_edges], ring=base_ring)
    
    # if the space at the given vertex is 0, there's nothing to do.
    if Vdi[0].dimension() == 0 and rep.get_space(vertex).dimension() == 0:
        return rep

    # concatenate inward maps (twisted by sign) to get a map Vdi -> Vi
    concat_matrix_in = sum(Vdi[2][i] * maps[in_edges[i]].matrix() * (1 if in_edges[i][2][-1] == '0' else -1) for i in range(len(in_edges)))
    map_in = linear_transformation(Vdi[0], Vi, concat_matrix_in)

    # concatenate outward maps to get a map Vi -> Vdi
    concat_matrix_out = sum(maps[out_edges[i]].matrix() * Vdi[1][i] for i in range(len(out_edges)))
    map_out = linear_transformation(Vi, Vdi[0], concat_matrix_out)

    # if the sum of outward maps is not injective, then rep has the simple at vertex as a submodule
    if (no_sub_check and not map_out.is_injective()):
        raise ValueError('This representation has S_{0} as a subrepresentation.'.format(vertex))

    # compute cokernel of map Vi -> Vdi
    cokernel, proj, lift = Vdi[0].quotient_abstract(map_out.image())

    # define new maps out of and into kernel
    # UPDATE: the ".matrix()" at the end of each expression may be unnecessary for Sage's quiver representation
    # constructor. I don't think this messes things up either way, though. Leaving it in for now.
    new_maps_out = {out_edges[i]: (linear_transformation(Vdi[0], spaces[out_edges[i][1]], Vdi[2][i]) * map_out * map_in * lift).matrix() for i in range(len(out_edges))}
    new_maps_in = {in_edges[i]: (proj * linear_transformation(spaces[in_edges[i][0]], Vdi[0], Vdi[1][i]) * (1 if in_edges[i][2][-1] == '0' else -1)).matrix() for i in range(len(in_edges))}

    # update space at given vertex
    spaces.update({vertex: cokernel})

    # update maps incident to given vertex
    maps.update(new_maps_out)
    maps.update(new_maps_in)

    # return reflected representation
    new_rep = PQ.representation(base_ring, spaces, maps)
    return new_rep


# INPUT: A representation of a DoubleQuiver and a list of pairs (vertex, sign), where sign is either '+' or '-'.
# OUTPUT: The result of applying the sequence of reflection functors in left-to-right order.
def reflection_functor_seq(rep, vertices):
    if len(vertices) == 0:
        return rep
    elif vertices[-1][1] == '+':
        return reflection_functor_plus(reflection_functor_seq(rep, vertices[:-1]), vertices[-1][0])
    elif vertices[-1][1] == '-':
        return reflection_functor_minus(reflection_functor_seq(rep, vertices[:-1]), vertices[-1][0])


# INPUT: A (not necessarily doubled) quiver, a vertex, a list of vertices, and a string of '+' and '-'.
# OUTPUT: The representation obtained by starting with the simple at the first vertex and applying
# the sequence of reflection functors given by the other vertices and signs.
# NOTE: The output of this function may not actually be a shard module!
def shard_module(Q, start_vertex, reflect_vertices, sign_seq):
    DQ = DoubleQuiver(Q)
    PDQ = DQ.path_semigroup()
    reflect_seq = list(zip(reflect_vertices, list(sign_seq)))
    return reflection_functor_seq(PDQ.S(QQ, start_vertex), reflect_seq)



# INPUT: A pair (hyperplane H, a list of hyperplanes) where each hyperplane is
# given by a 1 x n matrix.
# OUTPUT: A list of lists, each of which consists of indices into the list in the input
# corresponding to hyperplanes which give the same result when intersected with H.
def walls(p_and_f):
    wall_grouping = []
    for i in range(len(p_and_f[1])):
        independent = True
        for j in range(len(wall_grouping)):
            w = wall_grouping[j]
            subsystem_matrix = Matrix([p_and_f[0][0], p_and_f[1][w[0]][0], p_and_f[1][i][0]])
            if rank(subsystem_matrix) == 2:
                wall_grouping[j] = w + [i]
                independent = False
        if independent:
            wall_grouping = wall_grouping + [[i]]
    return wall_grouping


# INPUT: A representation of a DoubleQuiver and a sequence of vertices of the quiver
# OUTPUT: A string of '+' and '-' such that the given representation can be obtained
# by another representation by applying reflection functors at the given indices with
# these signs.
# (In other words, we can apply reflection functors at the given sequence of indices
# in reverse order with the opposite signs, and the intermediate representations thus
# obtained will be NoQuot or NoSub, as appropriate.)
def undo_reflections(rep, reflect_seq):
    if len(reflect_seq) == 0:
        return ''
    try: 
        new_rep = reflection_functor_minus(rep, reflect_seq[-1])
        return undo_reflections(new_rep, reflect_seq[:-1]) + '+'
    except ValueError:
        new_rep = reflection_functor_plus(rep, reflect_seq[-1])
        return undo_reflections(new_rep, reflect_seq[:-1]) + '-'


# Given a shard module:
# - finds reduced expression for root
# - finds module's sign sequence with respect to this reduced expression
# - for each wall associated to the reduced expression, flips those signs and tries computing
# the ensuing shard module (NOTE! This may be super inefficient, and for larger situations it may be 
# better to implement this with convex geometry jiggery-pokery!)
# - also it depends on the truth of David & Hugh's conjecture, but oh well
# - return the reduced expression, the sign sequence, and the list of walls for which this works

# This method determines, given an expression for a representation in terms of reflection functors,
# the minimal flips we can make in the signs of those reflection functors
# such that the results are also legitimate sequences of reflection functors.
#
# INPUT: A representation of a DoubleQuiver
# and a pair (vertex, sequence of vertices) representing a positive expression 
# (simple root, sequence of reflections) for beta, the dimension vector of the representation
# OUTPUT: 
# - the provided expression (for legacy reasons)
# - a sequence of signs such that, if we start with the simple module associated
#    to the simple root in our expression and apply the reflection functors
#    associated to the reflections in our expression with these signs, we obtain
#    the given representation
# - a collection of lists of indices representing signs we can flip in this sequence
#    such that we produce another valid sequence of reflection functors.
#    Each index represents an entry in the output of walls(), representing one or more hyperplanes
#    which cut the hyperplane defined by our dimension vector in the same way.
# - for each list of indices in the previous output, the hyperplanes produced by plane_and_fractures
#    which correspond to those indices
def shard_module_walls(module, chosen_expression):
    root = module.dimension_vector()
    G = module.quiver().to_undirected()
    cartan_matrix = cartan(G)/2 + identity_matrix(len(G.vertices(sort=True)))
    #if chosen_expression == None:
    #    expression = get_reduced_expression(cartan_matrix, root)
    #else:
    expression = chosen_expression
    signs = undo_reflections(module, expression[1])
    p_and_f = plane_and_fractures(cartan_matrix, expression[0], expression[1])
    all_walls = walls(p_and_f)

    border_walls = []
    for w in all_walls:
        new_signs = ''
        for i in range(len(signs)):
            if i in w:
                if signs[i] == '+':
                    new_signs = new_signs + '-'
                elif signs[i] == '-':
                    new_signs = new_signs + '+'
            else:
                new_signs = new_signs + signs[i]
        try:
            flipped_rep = shard_module(module.quiver(), expression[0], expression[1], new_signs)
            border_walls.append(w)
        except:
            pass
    return (expression, signs, border_walls, [p_and_f[1][w[0]] for w in border_walls])