"""

Gives the Alexander matrix of the universal abelian cover of a link using Fox Calculus.  That is, let X_L be the universal abelian cover of S^3 \setminus L and x_0 a point in S^3 \setminus L.  
Then this computes a presentation of the module H_1(X_L,\tilde{x_0}) where \tilde{x_0} is the preimage in the cover. 
Thanks to Nathan Dunfield and Sherry Gong for sharing their Fox Milnor algorithm with me on which much of this is based.

"""
import snappy
from snappy.snap.nsagetools import (join_lists,
                                    minimum_exponents,
                                    uniform_poly_exponents,
                                    fox_derivative,
                                    convert_laurent_to_poly,
                                    MapToFreeAbelianization,
                                    MapToGroupRingOfFreeAbelianization)
from sage.all import matrix, QQ, ZZ, PolynomialRing, Integers, vector, LaurentPolynomialRing, prod


def standard_map_to_group_ring_of_abelianization_for_link(manifold_group):
    """

    From Dunfield and Gong
    >>> G = snappy.Link('L5a1').exterior().fundamental_group()
    >>> phi = standard_map_to_group_ring_of_abelianization_for_link(G)
    >>> phi(''.join(m for m, l in G.peripheral_curves()))
    a*b*c
    """
    G = manifold_group
    phi = MapToGroupRingOfFreeAbelianization(G)
    R = phi.range()
    betti = R.ngens()
    meridians = [m for m, l in G.peripheral_curves()]
    assert len(meridians) == phi.U.nrows() == betti
    A = matrix([MapToFreeAbelianization.__call__(phi, m) for m in meridians])
    A = A.transpose().inverse()
    phi = MapToGroupRingOfFreeAbelianization(G)
    phi.U = A * phi.U
    assert all(phi(m) == g for m, g in zip(meridians, R.gens()))
    return phi


def alexander_matrix_link_group(G, simplify):
    """
    Compute a presentation of the pointed Alexander module from the Fox derivatives from the group.  Here the columns represent the generators of the module and the rows the relations. 
    Input is a finitely presented group for G. Simplify can be True or False.  If True, we multiply each row by a monomial to make each entry a polynomial rather than a Laurent polynomial.  
    If False, it returns matrix from the Fox derivatives. 
    >>> L = snappy.Link('L5a1')
    >>> Man = L.exterior()
    >>> G = Man.fundamental_group()
    >>> M = alexander_matrix_link_group(G, True)
    """
    
    phi = standard_map_to_group_ring_of_abelianization_for_link(G)
    R = phi.range()

    if G.num_relators() == 0:
        return matrix(R, 1, G.num_generators())
    M = [[fox_derivative(rel, phi, var) for var in G.generators()] for rel in G.relators()]
    if simplify == False:
        return matrix(R,M)
    for i in range(G.num_relators()):
        entries = sum([M[i]], [])
        minexp = minimum_exponents(entries)
        for j in range(G.num_generators()):
            M[i][j] = convert_laurent_to_poly(M[i][j], minexp, R)
    return matrix(R,M)


def alexander_matrix_link(L,simplify):
    """
    Compute a presentation of the pointed Alexander module from the Fox derivatives from a snappy link.  Here the columns represent the generators of the module and the rows the relations. 
    Simplify can be True or False.  If True, we multiply each row by a monomial to make each entry a polynomial rather than a Laurent polynomial.  
    If False, it returns matrix from the Fox derivatives. 
    >>> L = snappy.Link('L5a1')
    >>> Man = L.exterior()
    >>> M = alexander_matrix_link(Man, True)
    """
    manifold = L.exterior()
    G = manifold.fundamental_group()
    return alexander_matrix_link_from_group(G, simplify)