-- Goal: given the action of a Lie algebra g on a vector space W, 
-- compute the g action on Sym^d W or Wedge^k W

-- Also a few functions for invariant theory: 
-- casimir
-- DKPcoeffs

-- Load the following packages in order to access a function 
-- for computing the sign of a permutation
loadPackage("SpechtModule");


-- Note: the natural bases of Sym^d W and Wedge^k W are handled differently
-- In Sym^d W, we compute the action on the variable w_j (which of course may show up with degree > 1 in a monomial)
-- In Wedge^k W, we compute the action in each position of a wedge w_i_1 ^ ... ^ w_i_k
 


----------------------
----------------------
-- Sym^d W
----------------------
----------------------



-- Represent a degree d monomial as an exponent vector
-- Compute two hash tables: 
-- Z to E: the exponent vector with integer label i
-- E to Z: the integer label of an exponent vector


-- Let X record the action of a basis element, represented as a list of options (M2 sparse input)
-- Suppose X acts in the variable j
-- Get the coefficient and the monomials of the answer

-- Example: first three nonzero entries from X12 on Wedge^7 S+
-- Xhash = new HashTable from {(0, 4) => 1, (1, 5) => 1, (2, 6) => 1}



XactionOnVariablej = (Xhash,expv,j) -> (
    if expv_j == 0 then return null;
    -- Get the action in coord j
    Xwj:=select(Xhash, p -> p#0#1==j);
    if Xwj=={} then return null;
    newexpv:={};
    answer:= for t in Xwj list (
        if j==t#0#0 then newexpv = expv else newexpv=apply(#expv, i-> if i==j then expv_i-1 else if i==t#0#0 then expv_i+1 else expv_i);	
     	{newexpv, expv_j*(t#1)}
    );
    return answer
);



-- Simplify the expressions obtained
addExpressions = (H) -> (
    K:=unique apply(H, i -> first i);
    P:= for k in K list (
        {k,sum apply(select(H, h -> h_0==k), p -> p_1)}
    );
    select(P, p -> p_1 !=0)
);



-- Let X be an element of the spin basis represented as a hashtable
-- Let the monomial be represented by an exponent vector
-- Get the coefficients and the monomials of the answer
XactionOnMonomial = (Xhash,expv) -> (
    A:=delete(null,apply(#expv, j -> if expv_j > 0 then XactionOnVariablej(Xhash,expv,j)));
    addExpressions(flatten A)
);


-- Now do the above for a basis of Symd
XactionOnSymd = (Xhash,d,ZtoE,EtoZ) -> (
    N:=#keys(ZtoE);
    -- Build the matrix column-by-column
    M:={};
    Mj:={};
    i:=-1;
    for j from 0 to N-1 do (
      Mj=XactionOnMonomial(Xhash,ZtoE#j);
      for p in Mj do (
	i = EtoZ#(p_0);
        M = append(M, (i,j)=> p_1)  
      )
    );
    return M
);


----------------------
----------------------
-- Wedge^k W
----------------------
----------------------


-- Suppose X acts in the variable j
-- Get the coefficient and the monomials of the answer

-- Example: first three nonzero entries from X12 on Wedge^7 S+
-- Xhash = new HashTable from {(0, 4) => 1, (1, 5) => 1, (2, 6) => 1}


-- Input: a list of integer labels, not necessarily increasing
straightenSign = (a) -> (
    b:=sort apply(#a, i -> {a_i,i});
    permutationSign(apply(b, i -> last i))
);



-- Let X be an element of g, represented in M2 sparse form 
-- That is, a list of options, e.g. {(0, 4) => 1, (1, 5) => 1, (2, 6) => 1}
-- Let wedge be a k-subset of {0,...,n-1}
-- Get the answer when X acts in the j^th position of a wedge (where counting starts at 0)
XactionInPositionj = (Xhash,wedge,j) -> (
    -- Get the action in coord j
    Xwj:=select(Xhash, p -> p#0#1==wedge_j);
    if Xwj=={} then return null;
    --if not Xhash#?(ZtoW#(wedge_j)) then return null;
    newwedge:={};
    dj:=0;
    answer:= for t in Xwj list (
	if not member(t#0#0,drop(wedge,{j,j})) then (
            if wedge_j==t#0#0 then newwedge = wedge else newwedge = apply(#wedge, k -> if k==j then t#0#0 else wedge_k);
	    dj = straightenSign(newwedge);
     	    {sort newwedge, dj*(t#1)}
	)
    );
    return delete(null,answer)
); 
    
  

-- Let X be an element of the spin basis represented as a hashtable
-- Let wedge be a k-subset of {0,...,n-1}
-- Get the coefficients and the wedges of the answer
XactionOnWedge = (Xhash,wedge) -> (
    A:=delete(null,apply(#wedge, j -> XactionInPositionj(Xhash,wedge,j)));
    addExpressions(flatten A)
);


-- Now do the above for a basis of wedges
-- Note: ZtoW and WtoZ have a very different role from what they had in "s010 wedge IT.3.m2"
XactionOnWedgek = (Xhash,k,ZtoW,WtoZ) -> (
    N:=#keys(ZtoW);
    -- Build the matrix column-by-column
    M:={};
    Mj:={};
    i:=-1;
    for j from 0 to N-1 do (
      Mj=XactionOnWedge(Xhash,ZtoW#j);
      for p in Mj do (
	i = WtoZ#(p_0);
        M = append(M, (i,j)=> p_1)  
      )
    );
    return M
);



sparseRep = M -> (
    delete(null,flatten apply(numRows(M), i -> apply(numColumns(M), j -> if M_(i,j)!=0 then (i,j)=>M_(i,j))))    
)




checkLieAlgHomOnPair = (B, rhoB, i, j, writeInBasisFunction) -> (
    gbracket := br(B_i,B_j);    
    Wbracket := br(rhoB_i,rhoB_j);
    c := writeInBasisFunction(gbracket);
    Wbracket == sum apply(#rhoB, i -> c_i*rhoB_i)
);


checkLieAlgHom = (B, rhoB, writeInBasisFunction) -> (
    for i from 0 to #B-2 do (
        for j from i+1 to #B-1 do (
	    if not checkLieAlgHomOnPair(B,rhoB,i,j,writeInBasisFunction) then (
	        print concatenate("Brackets not compatible on basis elements ",toString({i,j})) << endl;
		return false
	    )
        )
    );
    return true
);	



isDiagonalMatrix = (M) -> (
    for i from 0 to numRows(M)-1 do (
        for j from 0 to numColumns(M)-1 do (
	    if j!=i and M_(i,j)!=0 then return false
	)
    );
    return true
);



isDiagonalSparseMatrix = (M) -> (
    all(M, o -> (o#0)#0==(o#0)#1)
);



sparseEquals = (A,B) -> (if #A != #B then return false;
    all(#A, i -> A#i#0==B#i#0 and A#i#1==B#i#1)
)



sparseHasZeroColumn = (j,M) -> (
    all(M, o -> (o#0)#1 !=j)
);


exponentVectors = (n,d) -> (
    if d==2 then (
        B2a:=apply(n, i -> apply(n, j -> if i==j then 2 else 0));
        B2b:=apply(subsets(apply(n, i -> i),2), p -> apply(n, j -> if j==p_0 then 1 else if j==p_1 then 1 else 0));
        return join(B2a,B2b)
    );
    SS:=QQ[tt_1..tt_n];
    B:=flatten entries basis(d,SS);
    apply(B, i -> flatten exponents(i))
);



casimir = (B,Bstar) -> (
    sum apply(#B, i -> (B_i)*(Bstar_i))  
);



DKPcoeffs = (f,c) -> (
    E:=f;
    k:=1;
    g:=c*f;
    E = E | g;
    while rank(E)==numColumns(E) do (
        g = c*g;
        E = E | g;
    );
    Pcoeffs:=gens ker E;
    reverse flatten entries Pcoeffs
);


-- Let P be the DKPcoeffs above
reynoldsProjection = (f,c,P) -> (
    Q:=drop(P,{#P-1,#P-1});
    A:={f};
    for i from 1 to #Q-1 do (
        A=prepend(c*(A_0),A)	
    );
    return sum apply(#A, i -> Q_i*A_i)
)



degDInvariantsAsVectors = (invBP) -> (
    apply(invBP, x -> reynoldsProjection(Id_(x_0),c,x_1))
);    


evToMonomial = (ev,R) -> (
    product apply(#ev, i -> if ev_i==0 then 1 else (R_i)^(ev_i))
);


vectorToPolynomial = (w,ZtoE,R) -> (
    v:=flatten entries w;
    sum apply(#v, i -> if v_i==0 then 0_R else (v_i)*evToMonomial(ZtoE#i,R))
);


--R=QQ[w_0..w_(#snk-1)];
degDInvariants = (degDInvariantsAsVectors) -> (
    apply(degDInvariantsAsVectors, v -> vectorToPolynomial(v,ZtoE,R))
)


IdColumn = (n,i) -> (
    transpose matrix {apply(n, j -> if j==i then 1/1 else 0/1)}
);





-*

checkInvariance = (n,k,f,S) -> (
    
R=QQ[a_1..d_n];
I = ideal({sum(apply(n, i -> a_i*b_i))+1,sum(apply(n, i -> c_i*d_i))+1});
S=R/I;
WB= wedgeBasis(n);
E = join(apply(n, i -> contractByemk(i,WB,S)),apply(n, i -> wedgeByek(i,WBS)));
x1 = sum apply(2*n, i -> S_i*E_i);
x2 = sum apply(2*n, i -> S_(2*n+i)*E_i);
s = apply(n, i -> i);
M = (x2*x1)_s^s;
wkM = exteriorPower(k,M);
T=QQ[z_0..z_(binomial(n,k)-1)];
F = map(T,ring(f),gens T];
Ff = F(f);
g = map(T,T,wkM);
g(Ff)==Ff  
);

*-







end

restart

load "LieAlgebraRepresentations.v2.3.m2";


--------------------------------------------------
--------------------------------------------------
-- Example 1: sl2, Sym^2 Std 
--------------------------------------------------
--------------------------------------------------

Wbasis = {matrix {{1,0},{0,-1/1}}, matrix {{0,1},{0,0/1}}, matrix {{0,0},{1,0/1}}};
Wbasis = apply(Wbasis, M -> sparseRep M);

B2 = {{2, 0}, {1, 1}, {0, 2}};
ZtoE = new HashTable from apply(#B2, i -> {i,B2_i});
EtoZ = new HashTable from apply(#B2, i -> {B2_i,i});

map(QQ^3,QQ^3,XactionOnSymd(Wbasis_0,2,ZtoE,EtoZ))
map(QQ^3,QQ^3,XactionOnSymd(Wbasis_1,2,ZtoE,EtoZ))
map(QQ^3,QQ^3,XactionOnSymd(Wbasis_2,2,ZtoE,EtoZ))

-- Check
sparseEquals(XactionOnSymd(Wbasis_0,2,ZtoE,EtoZ),{(0,0) => 2, (2,2) => -2})
sparseEquals(XactionOnSymd(Wbasis_1,2,ZtoE,EtoZ),{(0,1) => 1, (1,2) => 2})
sparseEquals(XactionOnSymd(Wbasis_2,2,ZtoE,EtoZ),{(1,0) => 2, (2,1) => 1})



--------------------------------------------------
--------------------------------------------------
-- Example 2: sl2 acting on degree 2 binary forms
--------------------------------------------------
--------------------------------------------------

-- See Derksen and Kemper, Computational Invariant Theory, Example 4.5.21
-- Note: their basis of sl2 is not the same as Fulton and Harris's
-- so this doesn't match Example 1

gbasis = {matrix {{1,0},{0,-1/1}}, matrix {{0,1},{0,0/1}}, matrix {{0,0},{1,0/1}}};
-- Note: changed implementation from version 0
-- Now the user needs to input the transpose matrix
-- (the dual action is not built into the code)
Wbasis = {{(0,0)=>2, (2,2)=>-2}, {(0,1)=>2,(1,2)=>1}, {(1,0)=>1,(2,1)=>2}};



B2 = {{2, 0, 0}, {1, 1, 0}, {1, 0, 1}, {0, 2, 0}, {0, 1, 1}, {0, 0, 2}};
ZtoE = new HashTable from apply(#B2, i -> {i,B2_i});
EtoZ = new HashTable from apply(#B2, i -> {B2_i,i});
-- Also load Derksen and Kemper's answers for checking
R=QQ[a0,a1,a2];
DKX = (f) -> ( -2*a0*diff(a1,f) - a1*diff(a2,f) );
DKY = (f) -> ( -a1*diff(a0,f) - 2*a2*diff(a1,f) );
DKH = (f) -> ( -2*a0*diff(a0,f) + 2*a2*diff(a2,f));



map(QQ^6,QQ^6,XactionOnSymd(Wbasis_0,2,ZtoE,EtoZ))==transpose matrix apply(flatten entries basis(2,R), f -> apply(flatten entries basis(2,R), b -> -1*coefficient(b,DKH(f))))
map(QQ^6,QQ^6,XactionOnSymd(Wbasis_1,2,ZtoE,EtoZ))==transpose matrix apply(flatten entries basis(2,R), f -> apply(flatten entries basis(2,R), b -> -1*coefficient(b,DKX(f))))
map(QQ^6,QQ^6,XactionOnSymd(Wbasis_2,2,ZtoE,EtoZ))==transpose matrix apply(flatten entries basis(2,R), f -> apply(flatten entries basis(2,R), b -> -1*coefficient(b,DKY(f))))


B3 = {{3, 0, 0}, {2, 1, 0}, {2, 0, 1}, {1, 2, 0}, {1, 1, 1}, {1, 0, 2}, {0, 3, 0}, {0, 2, 1}, {0, 1, 2}, {0, 0, 3}};
ZtoE = new HashTable from apply(#B3, i -> {i,B3_i});
EtoZ = new HashTable from apply(#B3, i -> {B3_i,i});
map(QQ^10,QQ^10,XactionOnSymd(Wbasis_0,3,ZtoE,EtoZ))==transpose matrix apply(flatten entries basis(3,R), f -> apply(flatten entries basis(3,R), b -> -1*coefficient(b,DKH(f))))
map(QQ^10,QQ^10,XactionOnSymd(Wbasis_1,3,ZtoE,EtoZ))==transpose matrix apply(flatten entries basis(3,R), f -> apply(flatten entries basis(3,R), b -> -1*coefficient(b,DKX(f))))
map(QQ^10,QQ^10,XactionOnSymd(Wbasis_2,3,ZtoE,EtoZ))==transpose matrix apply(flatten entries basis(3,R), f -> apply(flatten entries basis(3,R), b -> -1*coefficient(b,DKY(f))))





--------------------------------------------------
--------------------------------------------------
-- Example 3: sl2 acting on degree 4 binary forms
--------------------------------------------------
--------------------------------------------------

-- See Newstead
-- The answers were computed using the code in "Newstead example.m2"

gbasis = {matrix {{1,0},{0,-1/1}}, matrix {{0,1},{0,0/1}}, matrix {{0,0},{1,0/1}}};
Wbasis = {{(0,0) => 4, (1,1) => 2, (3,3) => -2, (4,4) => -4}, {(0,1) => 4, (1,2) => 3, (2,3) => 2, (3,4) => 1}, {(1,0) => 1, (2,1) => 2, (3,2) => 3, (4,3) => 4}};

B2={{2, 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, 2, 0, 0, 0}, {0, 1, 1, 0, 0}, {0, 1, 0, 1, 0}, {0, 1, 0, 0, 1}, {0, 0, 2, 0, 0}, {0, 0, 1, 1, 0}, {0, 0, 1, 0, 1}, {0, 0, 0, 2, 0}, {0, 0, 0, 1, 1}, {0, 0, 0, 0, 2}};
ZtoE = new HashTable from apply(#B2, i -> {i,B2_i});
EtoZ = new HashTable from apply(#B2, i -> {B2_i,i});


map(QQ^15,QQ^15,XactionOnSymd(Wbasis_0,2,ZtoE,EtoZ))==matrix {{8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -8/1}}
map(QQ^15,QQ^15,XactionOnSymd(Wbasis_1,2,ZtoE,EtoZ))==matrix {{0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 3, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 2, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 2, 0, 6, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 3, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/1}}
map(QQ^15,QQ^15,XactionOnSymd(Wbasis_2,2,ZtoE,EtoZ))==matrix {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 2, 0, 6, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 4, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 8, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0/1}}


------------------------------------------------------------
------------------------------------------------------------
-- Example 4: Wedge^2 Std for sl4
------------------------------------------------------------
------------------------------------------------------------

Eijm = (i0,j0,m) -> ( matrix apply(m, i -> apply(m, j -> if i==i0 and j==j0 then 1/1 else 0/1)) );
Him = (i0,m) -> ( Eijm(i0,i0,m) - Eijm(i0+1,i0+1,m));
-- Standard representation of sl4
Wbasis = {Him(0,4),Him(1,4),Him(2,4),Eijm(0,1,4),Eijm(0,2,4),Eijm(0,3,4),Eijm(1,2,4),Eijm(1,3,4),Eijm(2,3,4),Eijm(1,0,4),Eijm(2,0,4),Eijm(3,0,4),Eijm(2,1,4),Eijm(3,1,4),Eijm(3,2,4)};
Wbasis = apply(Wbasis, M -> sparseRep M);


B2 = subsets(apply(15, i -> i),2);
ZtoW = new HashTable from apply(#B2, i -> {i,B2_i});
WtoZ = new HashTable from apply(#B2, i -> {B2_i,i});

time WedgeBasis = apply(Wbasis, X -> XactionOnWedgek(X,2,ZtoW,WtoZ));

-- Check that it's a Lie algebra homomorphism

writeInBasis = (M) -> (
    {M_(0,0),M_(0,0)+M_(1,1),M_(0,0)+M_(1,1)+M_(2,2),M_(0,1),M_(0,2),M_(0,3),M_(1,2),M_(1,3),M_(2,3),M_(1,0),M_(2,0),M_(3,0),M_(2,1),M_(3,1),M_(3,2)}    
);

br = (A,B) -> (A*B-B*A)
WdenseBasis = apply(Wbasis, M -> map(QQ^4,QQ^4,M));
WedgeDenseBasis = apply(WedgeBasis, M -> map(QQ^105,QQ^105,M));
checkLieAlgHom(WdenseBasis,WedgeDenseBasis,writeInBasis)

-- Check that the wedges are eigenvectors
apply(3, i -> isDiagonalSparseMatrix(Wbasis_i))


-- Check that w_0 ^ w_1 is a highest weight vector with weight L_0+L_1
PositiveRootIndices = {3,4,5,6,7,8};
all(PositiveRootIndices, i -> (WedgeDenseBasis_i)_{0} == 0)


--------------------------------------------------------------------------
--------------------------------------------------------------------------
-- Example 5: Wedge^2 of Example 3 (sl2 acting on degree 4 binary forms)
--------------------------------------------------------------------------
--------------------------------------------------------------------------


gbasis = {matrix {{1,0},{0,-1/1}}, matrix {{0,1},{0,0/1}}, matrix {{0,0},{1,0/1}}};
Wbasis = {{(0,0) => 4, (1,1) => 2, (3,3) => -2, (4,4) => -4}, {(1,0) => 4, (2,1) => 3, (3,2) => 2, (4,3) => 1}, {(0,1) => 1, (1,2) => 2, (2,3) => 3, (3,4) => 4}};

s52 = sort subsets(apply(5, i -> i),2);
-- Gives {{0, 1}, {0, 2}, {0, 3}, {0, 4}, {1, 2}, {1, 3}, {1, 4}, {2, 3}, {2, 4}, {3, 4}}

ZtoW = new HashTable from apply(#s52, i -> {i,s52_i});
WtoZ = new HashTable from apply(#s52, i -> {s52_i,i});


XactionOnWedgek(Wbasis_0,2,ZtoW,WtoZ)
XactionOnWedgek(Wbasis_1,2,ZtoW,WtoZ)
XactionOnWedgek(Wbasis_2,2,ZtoW,WtoZ)


-- Check these

-- Action of H
d = {4,2,0,-2,-4};
--Haction = diagonalMatrix apply(s52, i -> (1/1)*(d_(i_0)+d_(i_1)))
--map(QQ^10,QQ^10,XactionOnWedgek(Wbasis_0,2,ZtoW,WtoZ))==Haction

HactionOnPair = p -> (
    if d_(p_0)+d_(p_1)==0 then {} else {{p,1/1*(d_(p_0)+d_(p_1))}}
);
all(s52, p -> XactionOnWedge(Wbasis_0,p)==HactionOnPair(p))


-- Action of X
XactionOnPair = p -> (
    i:=p_0;
    j:=p_1;
    answer:={};
    if i+1!=j then answer = append(answer,{{i+1,j},(4-i)});
    if j<4 then answer = append(answer,{{i,j+1},(4-j)});
    return answer
);
all(s52, p -> XactionOnWedge(Wbasis_1,p)==XactionOnPair(p))

-- Action of Y
YactionOnPair = p -> (
    i:=p_0;
    j:=p_1;
    answer:={};
    if i>0 then answer = append(answer,{{i-1,j},i});
    if j-1!=i then answer = append(answer,{{i,j-1},j});
    return answer
);
all(s52, p -> XactionOnWedge(Wbasis_2,p)==YactionOnPair(p))




-- Check that we get a Lie algebra representation with the expected highest weights


time WedgeBasis = apply(Wbasis, X -> XactionOnWedgek(X,2,ZtoW,WtoZ));

writeInBasis = (M) -> (
    {M_(0,0),M_(0,1),M_(1,0)}
);

--WdenseBasis = apply(Wbasis, M -> map(QQ^4,QQ^4,M));
WedgeDenseBasis = apply(WedgeBasis, M -> (-1)*map(QQ^10,QQ^10,M));
checkLieAlgHom(gbasis,WedgeDenseBasis,writeInBasis)

------------------------------------------------------------
------------------------------------------------------------
-- POTENTIAL BUG ALERT There is an extra minus sign here. Why?
-- (Is it the same reason that we pick up a minus sign in Sym?)
------------------------------------------------------------
------------------------------------------------------------



-- Check that it has the expected highest weight

isDiagonalSparseMatrix(WedgeBasis_0)
(WedgeDenseBasis_1)_{9}==0





------------------------------------------------------------
------------------------------------------------------------
-- Example 6: Wedge^7 S+ for so10
------------------------------------------------------------
------------------------------------------------------------


load "so10.m2";


n=5;
B = so2nBasis(n);
sB = spinBasis(n,QQ);
--Half-spin representation
splusB = apply(sB, M -> M_(apply(2^(n-1), i -> i))^(apply(2^(n-1), i -> i)));


-*
-- First, check that splusB is correct:

-- 1. Check that I have a Lie algebra homomorphism from B to sB
checkLieAlgHom(B,sB,M ->writeInso2nBasis(M,n))

-- 2. Check that {1,2,3,4} is a highest weight vector with weight {1/2,1/2,1/2,1/2,-1/2}:

---- 2a. Check that the e_I's are weight vectors
all(n, i -> isDiagonalMatrix(splusB_i))

---- 2b. Check that {1,2,3,4} has the correct weight
-- In evenWedgeBasis it has position 2
(evenWedgeBasis(5))_2 == {1,2,3,4}
apply(5, i -> (splusB_i)_(2,2))=={1/2,1/2,1/2,1/2,-1/2}

---- 2c. Check that the highest weight vector is killed by all the positive roots
-- n=5:
Phiplus = so2nPositiveRoots(n);
PositiveRootIndices = select(#B, i -> member(B_i,Phiplus));
#PositiveRootIndices==20
all(PositiveRootIndices, i -> (splusB_i)_{2} == 0)

*-


-- Now construct Wedge^7 S+

n=5;
gbasis = so2nBasis(n);
sB = spinBasis(n,QQ);
splusBDense = apply(sB, M -> M_(apply(2^(n-1), i -> i))^(apply(2^(n-1), i -> i)));
splusB = apply(splusBDense, M -> sparseRep(M));

s167 = sort subsets(apply(16, i -> i),7);
ZtoW = new HashTable from apply(#s167, i -> {i,s167_i});
WtoZ = new HashTable from apply(#s167, i -> {s167_i,i});
time WedgeBasis = apply(splusB, X -> XactionOnWedgek(X,7,ZtoW,WtoZ));


-- Now check that Wedge^7 S+ is correct
-- Adapted from "so10 wedge IT.4.m2"

-- Not performed yet: Check that we have a Lie algebra homomorphism






-- Are the wedges eigenvectors?
apply(n, i -> isDiagonalSparseMatrix(WedgeBasis_i))
-- Yes



-- Check that v1 ^ v2 ^ v3 ^ v4 ^ v5 ^ v6 ^ v10 and v1^ v2 ^ v3 ^ v4 ^ v5 ^ v6 ^ v7 are highest weight vectors with weights {5/2,3/2,3/2,1/2,1/2} and {7/2,1/2,1/2,1/2,-1/2}

-- From Sage: Wedge^7 S+ decomposes as 
-- D5(5/2,3/2,3/2,1/2,1/2) + D5(7/2,1/2,1/2,1/2,-1/2)

-- What are the eigenvectors with these weights?
Splusweights = apply(16, i -> apply(5, j -> (splusBDense_j)_(i,i)));
select(s167, w -> sum(apply(w, i -> Splusweights_i))=={5/2,3/2,3/2,1/2,1/2})
select(s167, w -> sum(apply(w, i -> Splusweights_i))=={7/2,1/2,1/2,1/2,-1/2})
-- What position are these on the list?
select(#s167, i -> s167_i=={1, 2, 3, 4, 5, 6, 10})
-- Answer: 5008
select(#s167, i -> s167_i=={1, 2, 3, 4, 5, 6, 7})
-- Answer: 5005



Phiplus = so2nPositiveRoots(n);
PositiveRootIndices = select(#gbasis, i -> member(gbasis_i,Phiplus));
#PositiveRootIndices
all(PositiveRootIndices, i -> sparseHasZeroColumn(5008,WedgeBasis_i))
all(PositiveRootIndices, i -> sparseHasZeroColumn(5005,WedgeBasis_i))






------------------------------------------------------------
------------------------------------------------------------
-- Now do invariant theory
------------------------------------------------------------
------------------------------------------------------------


-- Look for an invariant quadric
-- Compute Sym^2 Wedge^7 S+
-- Output these matrices to some software that supports sparse linear algebra (Magma, Matlab, etc)

exponentVectors = (n,d) -> (
    if d!=2 then error "Only implemented for d=2";
    B2a:=apply(n, i -> apply(n, j -> if i==j then 2 else 0));
    B2b:=apply(subsets(apply(n, i -> i),2), p -> apply(n, j -> if j==p_0 then 1 else if j==p_1 then 1 else 0));
    join(B2a,B2b)
);

B2 = exponentVectors(#s167,2);
-- Can't complete this in a reasonable amount of time

-*
Here's what I would have done next:
ZtoE = new HashTable from apply(#B2, i -> {i,B2_i});
EtoZ = new HashTable from apply(#B2, i -> {B2_i,i});

-- How long does it take to compute one basis element?
time M=XactionOnSymd(WedgeBasis_0,2,ZtoE,EtoZ)
time Sym2Wedge7Basis = apply(WedgeBasis, X -> XactionOnSymd(X,2,ZtoE,EtoZ));
*-


-- Try some smaller examples
-- What about so(6)?
-- Or so(10) but on a smaller wedge?