needsPackage("Polyhedra")

-- basis(4,I) does not necessarily return a basis where the leading monomials of each term are distinct
-- So we need to start by upper triangularizing the basis
utBasis = (B) -> (
    C:={B_0};
    f:=0;
    c:=0;
    g:=0;
    for i from 1 to #B-1 do (
     f=B_i;	
     for j from 0 to i-1 do ( 
      g = C_j;	 
      c=coefficient(leadMonomial(g),f)/leadCoefficient(g);
      f = f - c*g
     );
     C=append(C,f)
    );
    return C    
);

-- Now it is easy to write any element of I4 in this basis
writeInUTBasis = (f, utI4) -> (
    coeffs:={};
    c:=0;
    g:=0;
    for i from 0 to #utI4-1 do (
	g = utI4_i;
        c=coefficient(leadMonomial(g),f)/leadCoefficient(g);
	coeffs = append(coeffs, c);
	f = f - c*g	
    );    
    return coeffs
);


quadraticFormOnI2 = (I2) -> (
    I:=ideal(I2);
    I4:=reverse sort flatten entries super basis(4,I);
    utI4 := utBasis(I4);
    Sym2I2 := flatten apply(#I2, i -> apply(i+1, j -> I2_i*I2_j));
    M := transpose matrix apply(Sym2I2, f -> writeInUTBasis(f,utI4)); 
    q := getSymbol "q";
    QFR:=(coefficientRing(ring(I2_0)))[q_0..q_9];
    Q2 := matrix {flatten apply(numgens QFR, i -> apply(i+1, j -> QFR_i*QFR_j))};
    print concatenate("rank ker(Sym2I2 -> I4) = ",toString(rank ker M)) << endl;
    return (Q2*( gens ker M ))_(0,0)
);


intersectionDimension = (A,B) -> (
    rank(A)+rank(B)-rank(A||B)    
);


eval = (p,f) -> (
    R:=ring(f);
    CR:=coefficientRing(R);
    lift(substitute(f, apply(numgens R, i -> R_i=>p_i)),CR)
);


-- Claim: the rowspace of the returned matrix O is W_p^perp
-- w = sum c_j*g_j is in W_p if all seven of the equations 
-- encoded by O evaluate to 0
-- So the rows of O are orthogonal to w for all w in W_p
-- So the rowspace of O is W_p^perp
Wpperp = (p,I2) -> (
    R:=ring(I2_0);
    O:=matrix apply(gens R, i -> apply(I2, j -> eval(p,diff(i,j))));
    S5:=subsets({0,1,2,3,4,5,6},5);
    i:=0;
    while(rank(O^(S5_i))<5) do i=i+1;
    return O^(S5_i)
);

-*
Replace this version of spinorOfPoint with the version that works for non-generic U below

pfaffian = (M) -> (
    M_(0,3)*M_(1,2)-M_(0,2)*M_(1,3)+M_(0,1)*M_(2,3)
);


-- The -1 in the definition of M below is because
-- Mukai writes things in the basis e_{-i} - sum a_{ij} e_j
spinorOfPoint = (y,I2,U0,Uinfty,Q) -> (
    B := (x,y) -> ( Q(x+y) - Q(x) - Q(y));
    BPM:=matrix apply(5, i -> apply(5, j -> B(flatten entries(U0^{i}),flatten entries(Uinfty^{j}))));
    E := transpose(U0 || (BPM^-1*Uinfty));
    U:=Wpperp(y,I2);
    if intersectionDimension(U,Uinfty) != 0 then error "dim(U cap Uinfty) != 0";
    Gammaf := transpose(E^-1*(transpose(U)));
    Gammaf = (Gammaf_{0,1,2,3,4})^-1*Gammaf;
    M := (-1)*Gammaf_{5,6,7,8,9};
    return {1,M_(0,1),M_(0,2),M_(0,3),M_(0,4),M_(1,2),M_(1,3),M_(1,4),M_(2,3),M_(2,4),M_(3,4),pfaffian((M_{1,2,3,4}^{1,2,3,4})),pfaffian((M_{0,2,3,4}^{0,2,3,4})),pfaffian((M_{0,1,3,4}^{0,1,3,4})),pfaffian((M_{0,1,2,4}^{0,1,2,4})),pfaffian((M_{0,1,2,3}^{0,1,2,3}))}
);
*-


-- U is a nu x 2*nu matrix w.r.t. the basis of V 
-- Q is the quadratic form on V
-- Check whether the rowspace of U is a Lagrangian of (V,Q)
isLagrangian = (U,Q) -> (
    K:=ring(U_(0,0));
    E:=entries U;
    nc:=numColumns(U);
    nr:=numRows(U);
    R:=K[cc_0..cc_(nc-1)];
    u:=sum apply(nr, i -> cc_i*E_i);
    Q(u) == 0
)



S = subsets({1,2,3,4,5});
evenWedgeBasis = {{}, {1, 2}, {1, 3}, {1, 4}, {1, 5}, {2, 3}, {2, 4}, {2, 5}, {3, 4}, {3, 5}, {4, 5}, {1, 2, 3, 4}, {1, 2, 3, 5}, {1, 2, 4, 5}, {1, 3, 4, 5}, {2, 3, 4, 5}};
oddWedgeBasis = {{1}, {2}, {3}, {4}, {5}, {1, 2, 3}, {1, 2, 4}, {1, 2, 5}, {1, 3, 4}, {1, 3, 5}, {1, 4, 5}, {2, 3, 4}, {2, 3, 5}, {2, 4, 5}, {3, 4, 5}, {1, 2, 3, 4, 5}};
wedgeBasis = join(evenWedgeBasis,oddWedgeBasis);



contracteIByemk = (k,I,oddWedgeBasis,R) -> (    
    if not member(k,I) then return apply(#oddWedgeBasis, i -> 0_R);
    Imk := select(I, j -> j!=k);
    p:=position(I, j -> j==k);
    apply(oddWedgeBasis, J -> if J==Imk then (-1_R)^p else 0_R)
);



contractByemk = (k,evenWedgeBasis,oddWedgeBasis,R) ->( 
    transpose matrix apply(evenWedgeBasis, I -> contracteIByemk(k,I,oddWedgeBasis,R))    
)


wedgeeIByek = (k,I,oddWedgeBasis,R) -> (    
    if member(k,I) then return apply(#oddWedgeBasis, i -> 0_R);
    Icupk := sort prepend(k,I);
    p:=position(Icupk, j -> j==k);
    apply(oddWedgeBasis, J -> if J==Icupk then (-1_R)^p else 0_R)
);



wedgeByek = (k,evenWedgeBasis,oddWedgeBasis,R) ->( 
    transpose matrix apply(evenWedgeBasis, I -> wedgeeIByek(k,I,oddWedgeBasis,R))    
)


Gammaf = (y,I2,U0,Uinfty,Q) -> (
    B := (x,y) -> ( Q(x+y) - Q(x) - Q(y));
    BPM:=matrix apply(5, i -> apply(5, j -> B(flatten entries(U0^{i}),flatten entries(Uinfty^{j}))));
    E := transpose(U0 || (BPM^-1*Uinfty));
    U:=Wpperp(y,I2);
    return transpose(E^-1*(transpose(U)))
);



phiSubv = (v) -> (
    R:=ring(v_0);
    M:=matrix apply(16, i -> apply(16, j -> 0_R));  
    for i from 0 to 4 do (
        M = M + (v_i)*contractByemk(i+1,evenWedgeBasis,oddWedgeBasis,R)
    );
    for i from 0 to 4 do (
        M = M + (v_(i+5))*wedgeByek(i+1,evenWedgeBasis,oddWedgeBasis,R)
    );
    return M
);


pfaffian = (M) -> (
    M_(0,3)*M_(1,2)-M_(0,2)*M_(1,3)+M_(0,1)*M_(2,3)
);



-- Assume U is the rowspace of a matrix of vectors written in the basis e_(-1)...e_(5)
spinorOfPoint = (y,I2,U0,Uinfty,Q) -> (
    U:=Gammaf(y,I2,U0,Uinfty,Q);
    if intersectionDimension(U,Uinfty) == 0 then (
         U = (U_{0,1,2,3,4})^-1*U;	
	 M := (-1)*U_{5,6,7,8,9};
	 return {1_(ring M),M_(0,1),M_(0,2),M_(0,3),M_(0,4),M_(1,2),M_(1,3),M_(1,4),M_(2,3),M_(2,4),M_(3,4),pfaffian((M_{0,1,2,3}^{0,1,2,3})),pfaffian((M_{0,1,2,4}^{0,1,2,4})),pfaffian((M_{0,1,3,4}^{0,1,3,4})),pfaffian((M_{0,2,3,4}^{0,2,3,4})),pfaffian((M_{1,2,3,4}^{1,2,3,4}))}	
    );
    if intersectionDimension(U,Uinfty) > 0 then (
        v:=intersect apply(entries U, u -> ker(phiSubv(u)));
        return flatten entries gens v
    )
);



-- Given the general point on a line in P6 with parameters [a:b], give the ideal of the image of this line under the spinor embedding 
spinorRepOfLine = (p,It,U0,Uinfty,Q) -> (
    s:= spinorOfPoint(p,It,U0,Uinfty,Q);
    d:=lcm(apply(s, j -> denominator(j)));
    s = d*s;  
    CR:=baseRing(ring(s_0));
    s = apply(s, i -> lift(i,CR));
    ElimR:=QQ[a,b,t_0,t_1,x_0,x_12,x_13,x_14,x_15,x_23,x_24,x_25,x_34,x_35,x_45,x_1234,x_1235,x_1245,x_1345,x_2345, MonomialOrder=>Eliminate(2)];
    spinvars := {x_0,x_12,x_13,x_14,x_15,x_23,x_24,x_25,x_34,x_35,x_45,x_1234,x_1235,x_1245,x_1345,x_2345};
    F:=map(ElimR,CR,{a,b,t_0,t_1});
    ElimI:=ideal(apply(16, i -> spinvars_i - F(s_i)));
    I:=selectInSubring(1,gens gb ElimI);
    return flatten entries I
);






writeInI2Basis = (f,g,F) -> (
    M = gens kernel(matrix {prepend(F(f),g)});
    v = flatten entries(M_{0});
    v = apply(v, i -> lift(i,K));
    -- Write the relation as F(f) - sum c_i*g_i and then return the c_i's
    return drop((-1)*apply(v, i -> i/v_0),{0,0})
);


alphaxv = (x,v) -> ( ((Q(v+x)-Q(v)-Q(x))/Q(x))*x - v)


alpha = (x) -> ( 
    n:=lift((#x)/2,ZZ);
    I:= apply(2*n, i -> apply(2*n, j -> if i==j then 1 else 0));
    M:=matrix apply(I, j -> -alphaxv(x,j));
    transpose M
);


factorInReflections = (T) -> (
    n:=lift(numRows(T)/2,ZZ);
    K:=ring(T);
    WL:={};
    RL:={};
    Twi:={};
    wi:={};
    ci:={};
    indexList:=flatten apply(n, i -> {i,i+n});
    for i in indexList do (
        Twi=flatten entries((product reverse RL)*T_{i});
	wi=apply(2*n, j -> if j==i then 1_K else 0_K);
	ci=Twi-wi;
	if ci==apply(2*n, j -> 0_K) then continue;
	if Q(ci)==0 then error concatenate("Q(ci)=0 at step ",toString(i)) << endl;
	WL=append(WL,ci);
	print toString({i,ci}) << endl;
	RL=append(RL,alpha(ci))
    );
    return WL
);




actionOnSplus = (L,E,s) -> (
    (product (apply(#L, i -> sum apply(#(L_i), j -> (L_i_j)*(E_j)))))_s^s
);



isInFultonHarrisSO2n = (A) -> (
    n:=lift(numRows(A)/2,ZZ);
    K:=ring(A);
    In:=matrix apply(n, i -> apply(n, j -> if i==j then 1_K else 0));
    Zn:=matrix apply(n, i -> apply(n, j -> 0_K));
    Ma:=Zn | In;
    Mb:=In | Zn;
    M:=Ma || Mb;
    I2n:=matrix apply(2*n, i -> apply(2*n, j -> if i==j then 1_K else 0));
    return (entries(transpose(A)*M*A) == entries(M))
);


------------------------------------------------------------------
------------------------------------------------------------------
-- Equations of special curves and surfaces
------------------------------------------------------------------
------------------------------------------------------------------


---------------------------------------------------------------------
-- Code for equations of the tangent developable/g cuspidal curve
---------------------------------------------------------------------




-- Eisenbud's equations from "Green's Conjecture: An Orientation" page 13
eisenbudTangentDevelopableQuadrics = (v) -> (  
    M := matrix {drop(v,{#v-1,#v-1}),drop(v,{0,0})};
    g:=#v-1;
    I:={};
    for a from 0 to g-4 do (
     for b from a+1 to g-3 do (
        I=append(I,det(M_{a+2,b})-2*det(M_{a+1,b+1})+det(M_{a,b+2}))
     )
    );
    return I
)


-- I found this in Bayer and Stillman, "Some matrices related to Green's Conjecture", Prop. 2.9
schreyerTangentDevelopableQuadrics = (v) -> (
  L:={};
  g=#v-1;
  v=append(v,v_0);
  for i from 2 to g-2 do (
    for j from i to g-2 do (
      if i+j<=g-1 then (
        L=append(L,(i+j-1)*v_i*v_j-i*j*v_(i+j-1)*v_1-(i+j-i*j-1)*v_(i+j)*v_0)
      ) else (      
        L=append(L,(2*g-i-j-1)*v_i*v_j-(g-i)*(g-j)*v_(i+j-g+1)*v_(g-1)+(g-i-1)*(g-j-1)*v_(i+j-g)*v_g)
      )
    )  
  );	  
  return L
);


-- Below we are taking a hyperplane section of the tangent developable by setting v_g = v_0
eisenbudCuspidalCurveQuadrics = (v) -> (
    w:=append(drop(v,{0,0}),v_0);    
    M := matrix {v,w};
    g:=#v;
    I:={};
    for a from 0 to g-4 do (
     for b from a+1 to g-3 do (
        I=append(I,det(M_{a+2,b})-2*det(M_{a+1,b+1})+det(M_{a,b+2}))
     )
    );
    return I
);



-- Below we are taking a hyperplane section of the tangent developable by setting v_g = v_0
schreyerCuspidalCurveQuadrics = (v) -> (
  L:={};
  g=#v;
  v=append(v,v_0);
  for i from 2 to g-2 do (
    for j from i to g-2 do (
      if i+j<=g-1 then (
        L=append(L,{(i,j),(i+j-1)*v_i*v_j-i*j*v_(i+j-1)*v_1-(i+j-i*j-1)*v_(i+j)*v_0})
      ) else (    
        L=append(L,{(i,j),(2*g-i-j-1)*v_i*v_j-(g-i)*(g-j)*v_(i+j-g+1)*v_(g-1)+(g-i-1)*(g-j-1)*v_(i+j-g)*v_g})
      )
    )  
  );	  
  return L
);




---------------------------------------------------------------------
-- Code for equations of the balanced ribbon/K3 carpet
---------------------------------------------------------------------



RNC = (k) -> (
    R:=QQ[y_0..y_(k)];
    --get RNC equations from catalecticant matrices
    RNC1Row1:=apply(k, i -> x_i);
    RNC1Row2:=apply(k, i -> x_(i+1));
    RNC1:=flatten entries gens minors(2,matrix {RNC1Row1,RNC1Row2});
    return RNC1
);


balancedRibbon = (k) -> (
    if not instance(k,ZZ) then error "the input must be an integer";
    if k < 2 then error "k must be at least 2"; 
    x:=getSymbol "x";
    R:=QQ(monoid [y_0..y_(2*k)]);
    --get RNC equations from catalecticant matrices
    RNC1Row1:=apply(k, i -> R_i);
    RNC1Row2:=apply(k, i -> R_(i+1));
    RNC1:=flatten entries gens minors(2,matrix {RNC1Row1,RNC1Row2});
    RNC2Row1:=apply(k, i -> R_(k+i));
    RNC2Row2:=apply(k, i -> R_(k+i+1));
    RNC2:=flatten entries gens minors(2,matrix {RNC2Row1,RNC2Row2});
    --get trinomials
    T:=flatten apply(k-1, i -> apply(k-1, j -> R_(i+2)*R_(k+j) - 2*R_(i+1)*R_(k+j+1) + R_i*R_(k+j+2) ) );
    return flatten {T,RNC1,RNC2}
);


K3carpet = (k) -> (
    if not instance(k,ZZ) then error "the input k must be an integer";
    if k < 2 then error "the input k must be at least 2"; 
    x:=getSymbol "x";
    y:=getSymbol "y";
    L:=flatten apply(k+1, i -> {x_i,y_i});
    R:=QQ(monoid [L]);
    x = apply(k+1, i -> R_(2*i));
    y = apply(k+1, i -> R_(2*i+1));
    --get RNC equations from catalecticant matrices
    RNC1Row1:=apply(k, i -> x_i);
    RNC1Row2:=apply(k, i -> x_(i+1));
    RNC1:=flatten entries gens minors(2,matrix {RNC1Row1,RNC1Row2});
    RNC2Row1:=apply(k, i -> y_i);
    RNC2Row2:=apply(k, i -> y_(i+1));
    RNC2:=flatten entries gens minors(2,matrix {RNC2Row1,RNC2Row2});
    --get trinomials
    T:=flatten apply(k-1, i -> apply(k-1, j -> x_i*y_j-2*x_(i+1)*y_(j+1)+x_(i+2)*y_(j+2)) );
    return flatten {T,RNC1,RNC2}
);





---------------------------------------------------------------------
-- Code for equations of trivalent graph curves
---------------------------------------------------------------------

-*

Goal: given the adjacency matrix of a 3-edge-connected trivalent graph, compute the canonical ideal of the associated graph curve (in the sense of Bayer and Eisenbud)

*-


edgeList = M -> (
    -- Make a list of edges
    n:=numRows(M);
    vG:=apply(n, i -> i);
    eG:=sort select(subsets(vG,2), p -> M_(p_0,p_1)==1);
    return eG
);


coboundaryHomomorphism = M -> (
    -- Make a list of edges
    n:=numRows(M);
    vG:=apply(n, i -> i);
    eG:=edgeList(M);
    return matrix apply(vG, v -> apply(eG, e -> if v==min(e) then -1 else if v==max(e) then 1 else 0))
);


disjoint = (p) -> (
    #(unique flatten p)==4
);

triangle = p -> (
    #(unique flatten p)==3
);


BEGraphCurveIdeal = (M) -> (
    Z1G := gens ker coboundaryHomomorphism M;
    gG := numColumns Z1G;
    eG := edgeList(M);    
    e:=getSymbol "e";
    yy:=getSymbol "yy";
    y:=getSymbol "y";
    edgeVarList := apply(eG, i -> e_i);
    cycleVarList := apply(gG, i -> yy_i);
    fullVarList := join(edgeVarList, cycleVarList);
    R := QQ[fullVarList, MonomialOrder=>Eliminate(#edgeVarList)];
    mIdG := matrix apply(gG, i -> apply(gG, j -> if i==j then -1 else 0));
    A := Z1G || mIdG;
    -- Get monomials for pairs of disjoint edges
    T := new HashTable from apply(#eG, i -> {eG_i,R_i});
    J1 := ideal(vars(R)*A);
    J2 := ideal delete(null,apply( subsets(eG,2), p -> if disjoint(p) then (T#(p_0))*(T#(p_1))));
    J3 := ideal delete(null,apply( subsets(eG,3), p -> if triangle(p) then (T#(p_0))*(T#(p_1))*(T#(p_2))));
    J := J1+J2+J3;
    Igens := flatten entries selectInSubring(1,gens gb J);
    cycleVarList2 := apply(gG, i -> y_i);
    S := QQ[cycleVarList2];
    F := map(S,R,apply(numgens R, i -> if i < #eG then 0_S else S_(i-#eG)));
    I := ideal apply(Igens, f -> F(f));
    return I
);


BEGraphCurveIdealGivenZ1G = (M,Z1G) -> (
    gG := numColumns Z1G;
    eG := edgeList(M);
    e:=getSymbol "e";
    yy:=getSymbol "yy";
    y:=getSymbol "y";
    edgeVarList := apply(eG, i -> e_i);
    cycleVarList := apply(gG, i -> yy_i);
    fullVarList := join(edgeVarList, cycleVarList);
    R := QQ[fullVarList, MonomialOrder=>Eliminate(#edgeVarList)];
    mIdG := matrix apply(gG, i -> apply(gG, j -> if i==j then -1 else 0));
    A := Z1G || mIdG;
    -- Get monomials for pairs of disjoint edges
    T := new HashTable from apply(#eG, i -> {eG_i,R_i});
    J1 := ideal(vars(R)*A);
    J2 := ideal delete(null,apply( subsets(eG,2), p -> if disjoint(p) then (T#(p_0))*(T#(p_1))));
    J3 := ideal delete(null,apply( subsets(eG,3), p -> if triangle(p) then (T#(p_0))*(T#(p_1))*(T#(p_2))));
    J := J1+J2+J3;
    Igens := flatten entries selectInSubring(1,gens gb J);
    cycleVarList2 := apply(gG, i -> y_i);
    S := QQ[cycleVarList2];
    F := map(S,R,apply(numgens R, i -> if i < #eG then 0_S else S_(i-#eG)));
    I := ideal apply(Igens, f -> F(f));
    return I
);







---------------------------------------------------------------------
-- Code related to T state polytopes
---------------------------------------------------------------------

charactersOfT = (G) -> (
    if G=="SL16" then return apply(16, i -> apply(16, j -> if i==j then 1 else 0));   
    if G=="Spin10" then return {{0,0,0,0,0},{1,1,0,0,0},{1,0,1,0,0},{1,0,0,1,0},{1,0,0,0,1},{0,1,1,0,0},{0,1,0,1,0},{0,1,0,0,1},{0,0,1,1,0},{0,0,1,0,1},{0,0,0,1,1},{1,1,1,1,0},{1,1,1,0,1},{1,1,0,1,1},{1,0,1,1,1},{0,1,1,1,1}};  
);


-- M should be a 9x16 matrix
stateForTinG = (M,G) -> (
    characterList:=charactersOfT(G);   
    k:=numRows(M);
    N:=numColumns(M);
    sk := subsets(apply(N, i -> i),k);
    St := {};
    for i from 0 to #sk-1 do (
        if rank(M_(sk_i)) == k then (
	     St = append(St, sum apply(sk_i, j -> characterList_j))
        )
    );
    return unique St
);
  

isTsemistable = (M,G) -> (
    St:=stateForTinG(M,G);
    characterList:=charactersOfT(G);
    chi0:=(numRows M)*(1/#characterList)*(sum characterList);
    print concatenate("chi0 = ",toString(chi0)) << endl;
    StPolytope:=convexHull(transpose matrix St);
    contains(StPolytope, transpose matrix {chi0})
);