#Get fat point ideals in various degrees

#Take the derivative of a monomial x^a with respect to y^b and evaluate at a point p
def monomial_derivative_evaluation(a,b,p):
    a = vector(a);
    b = vector(b);
    d = a-b;
    neg_c = false;
    ind = 0;
    while not neg_c and ind<len(d):
        if d[ind]<0:
            neg_c = true
        ind = ind+1
    if neg_c:
        return(0)
    else:
        return(product([factorial(a[i])/factorial(a[i]-b[i])*p[i]^(d[i]) for i in range(0,len(a))]))
        
#Matrix whose kernel consists of all polynomials of degree deg which vanish to order ord at points in ptlist
def fat_matrix(deg,ord,ptlist):
    num_vars = len(ptlist[0]);
    monexp = IntegerVectors(deg,num_vars);
    mat_rows = [];
    for i in range(0,len(ptlist)):
        p = ptlist[i];
        deriv_mons = IntegerVectors(ord[i],num_vars);
        for m in deriv_mons:
            mat_rows=mat_rows+[[monomial_derivative_evaluation(n,m,p) for n in monexp]];
    return(matrix(mat_rows))
                    
#Take matrix with columns consisting of coefficients and turn it into polynomials
def vec_to_poly(coeff,deg,vrs,row=true):
    BS = matrix([[product([(vrs[i])^(a[i]) for i in range(0,len(a))]) for a in IntegerVectors(deg,len(vrs))]]);
    if coeff.nrows()==0 or coeff.ncols()==0:
        return(0)
    if row:
        return(coeff*(BS.transpose()))
    else:
        return(BS*coeff)

#Basis for polynomials of degree deg that vanish to order ord at points in ptlist
def fat_ideal(deg,ord,ptlist,vrs=None,poly_form=true):
    FM = fat_matrix(deg,ord,ptlist);
    K = matrix((fat_matrix(deg,ord,ptlist)).right_kernel().basis());
    if poly_form:
        return(vec_to_poly(K,deg,vrs))
    else:
        return(K)

###The following is a number of commands that allow the computation of (numerical) syzygies.
###The computations rely on passing completely into the realm of linear algebra and using SVD to approximate rank

#monomial exponentiation
def monexp(vrs,expn):
    return(product([v[0]^(v[1]) for v in zip(vrs,expn)]))
    
#Get basis for polynomials homogeneous of degree d
def poly_basis(deg,vrs):
    num=len(vrs);
    basisvecs = [list(m) for m in IntegerVectors(deg,num)];
    return([monexp(vrs,expn) for expn in basisvecs])

#get coefficients of polynomial in preferred order
def coeffs(F,basis):
    return([F.coefficient(m) for m in basis])

#Given a matrix of polynomials, convert each polynomial entry into a coefficient vector (column)
def coeffs_mat(M,Degs,vrs):
    return(block_matrix([[matrix(zip(*[coeffs(M[i,j],poly_basis(Degs[i,j],vrs)) for j in range(0,M.ncols())]))] for i in range(0,M.nrows())]))

#Given a list with polynomial entries, return a list of those polynomials multiplied by all monomials of a given degree (degree may vary for each polynomial)
def mon_mult(PolyList,DegreeDeficits,vrs):
    return(flatten([[expand(m*PolyList[i]) for m in poly_basis(DegreeDeficits[i],vrs)] for i in range(0,len(PolyList))]))

#Given a matrix with polynomial entries, return a matrix whose rows are multiplied by all monomials of a given degree (output of previous function)
def mon_mult_matrix(PolyMat,DegreeDeficits,vrs):
    return(matrix([mon_mult(rw,DegreeDeficits,vrs) for rw in PolyMat.rows()]))

#Index selection functions
def ind_select_less(L,tol):
    iList = [i for i,x in enumerate(L) if x<tol]
    if len(iList)==0:
        return(len(L))
    else:
        return(iList[0])

#def ind_select_greater(L,tol):
#    return([i for i,x in enumerate(L) if x>tol])

#Convert a matrix of coefficients (assumed column vectors) back into a matrix of polynomials
def coeff_to_poly(cfMat,Degs,vrs):
    pMat = [];
    for c in range(0,cfMat.ncols()):
        indx = 0;
        pCol = [];
        for cnt in range(0,Degs.nrows()):
            mns = poly_basis(Degs[cnt,c],vrs);
            cntln = len(mns);
            poly=0;
            for incnt in range(0,cntln):
                poly = poly+vec[indx]*mns[incnt];
                indx = indx+1;
            pCol=pCol+[poly];
        pMat = pMat+[pCol];
    return(transpose(matrix(pMat)))

#Convert a column relation back into a polynomial relation, given DegreeDeficits
def rel_to_poly(rMat,DegreeDeficits,vrs):
        relMat = [];
        for vec in rMat.columns():
            indx = 0;
            relCol = [];
            for cnt in range(0,len(DegreeDeficits)):
                mns = poly_basis(DegreeDeficits[cnt],vrs);
                cntln = len(mns);
                poly=0;
                for incnt in range(0,cntln):
                    poly = poly+vec[indx]*mns[incnt];
                    indx = indx+1;
                relCol=relCol+[poly];
            relMat=relMat+[relCol];
        return(transpose(matrix(relMat)))

#NUMERICAL syzygies of a list of homogeneous polynomials in a specific degree
def num_syz(PolyList,DegreeDeficits,deg,vrs,as_poly=false,tol=10^(-12)):
    bs = poly_basis(deg,vrs);
    exp_mons = mon_mult(PolyList,DegreeDeficits,vrs);
    cfMat = matrix([coeffs(F,bs) for F in exp_mons]);
    [U,S,V] = cfMat.change_ring(CDF).SVD();
    mn = min(S.ncols(),S.nrows());
    singVals = [S[i,i] for i in range(0,mn)];
    IND = ind_select_less(singVals,tol);
    num_ker = conjugate(U[range(0,U.nrows()),range(IND,U.ncols())])
    if as_poly:
        return(rel_to_poly(num_ker,DegreeDeficits,vrs))
    else:
        return(num_ker)

#Numerical intersection of column spaces (same ambient!).  Assumes that matrices have independent columns.
def num_intersect(L,tol=10^(-10)):
    #L is a list of matrices; we will return the numerical intersection of the column spaces
    ocmps = [orth_comp(M) for M in L];
    return(orth_comp(block_matrix([ocmps]),indep_cols=false,tol=tol))

#Numerical orthogonal complement of a column space with respect to another
def orth_comp_subspace(MAT1,MAT2,tol=10^(-12)):
    #Take orthogonal complement of col(MAT2) with respect to col(MAT1)
    [U1,S1,V1] = MAT1.change_ring(CDF).SVD();
    mn1 = min(S1.ncols(),S1.nrows());
    singVals1 = [S1[i,i] for i in range(0,mn1)];
    IND1 = ind_select_less(singVals1,tol);
    O1 = U1[range(0,U1.nrows()),range(0,IND1)];
    ProjM = O1*(O1.conjugate_transpose());
    ID = matrix.identity(ProjM.ncols());
    Orth1 = (ID-ProjM)*MAT2;
    [U,S,V] = Orth1.change_ring(CDF).SVD();
    mn = min(S.ncols(),S.nrows());
    singVals = [S[i,i] for i in range(0,mn)];
    IND = ind_select_less(singVals,tol);
    return(U[range(0,U.nrows()),range(0,IND)])

#Numerical orthogonal complement of column space
def orth_comp(MAT,indep_cols=true,tol=10^(-12)):
    [U,S,V] = MAT.change_ring(CDF).SVD();
    if indep_cols:
        r = MAT.ncols();
        return( (U[range(0,U.nrows()),range(r,U.ncols())]).conjugate())
    else:
        singVals = [S[i,i] for i in range(0,min(S.ncols(),S.nrows()))];
        IND = ind_select_less(singVals,tol);
        return((U[range(0,U.nrows()),range(IND,U.ncols())]).conjugate())            
        
#NUMERICAL syzygies of a matrix of polynomials with homogeneous entries
def num_syz_matrix(PolyMat,DegreeDeficits,row_degs,vrs,tol=10^(-12),as_poly=true):
    row_kers = [num_syz(list(PolyMat.row(i)),DegreeDeficits,row_degs[i],vrs,tol=tol) for i in range(0,PolyMat.nrows())];
    num_ker=num_intersect(row_kers,tol=tol)
    if as_poly:
        return(rel_to_poly(num_ker,DegreeDeficits,vrs))
    else:
        return(num_ker)