needsPackage("LieAlgebraRepresentations");
debug LieAlgebraRepresentations;


matApprox = (A,B,e) -> (
    max apply(flatten entries(A-B), x->abs(x) ) < e
);

matApprox = (A,B,e) -> (
    max apply(flatten entries(A-B), x->abs(x) ) < e
);


aIJXk = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num1 := product apply(toList(1..k-1), q -> ((GTP#(k,i)-i+1)-(GTP#(k-1,q)-q+1)+1));
    num2 := product apply(toList(1..k+1), q -> ((GTP#(k,i)-i+1)-(GTP#(k+1,q)-q+1)));
    denom1 := product apply(toList(1..k), q -> ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)+1));
    denom2 := product apply(toList(1..k), q -> if q==i then 1 else ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)));
    (-1)*(num1*num2)/(denom1*denom2)
);



AXk = (V, GTP, k) -> (
    num:=1;
    denom:=1;
    GTPPlusDeltakiEntries:={};
    GTPPlusDeltaki:={};
    L:=for i from 1 to k list (
	GTPPlusDeltakiEntries:=gtpPMDeltakiEntries(GTP,1,k,i);
	if not isValidEntryList("A",GTPPlusDeltakiEntries) then continue;
	GTPPlusDeltaki=gtpA(GTPPlusDeltakiEntries);
	c = aIJXk(GTP,i,k);
	{GTPPlusDeltaki,sqrt(c)}
    );
    lieAlgebraModuleElement(V,L)
);



AXkrepresentationMatrix = (V,k,BGT) -> (
    (1/1)*(transpose matrix apply(BGT, p -> writeInGTBasisA(AXk(V,gtpA(p),k),BGT)))
);



aIJYk = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num1 := product apply(toList(1..k-1), q -> ((GTP#(k,i)-i+1)-(GTP#(k-1,q)-q+1)));
    num2 := product apply(toList(1..k+1), q -> ((GTP#(k,i)-i+1)-(GTP#(k+1,q)-q+1)-1));
    denom1 := product apply(toList(1..k), q -> ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)-1));
    denom2 :=product apply(toList(1..k), q -> if q==i then 1 else ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)));
    (num1*num2)/(denom1*denom2)
);


AYk = (V, GTP, k) -> (
    num:=1;
    denom:=1;
    GTPMinusDeltakiEntries:={};
    GTPMinusDeltaki:={};
    L:=for i from 1 to k list (
	GTPMinusDeltakiEntries:=gtpPMDeltakiEntries(GTP,-1,k,i);
	if not isValidEntryList("A",GTPMinusDeltakiEntries) then continue;
	GTPPlusDeltaki=gtpA(GTPMinusDeltakiEntries);
	c = aIJYk(GTP,i,k);
	{GTPPlusDeltaki,sqrt(c)}
    );
    lieAlgebraModuleElement(V,L)
);


AYkrepresentationMatrix = (V,k,BGT) -> (
    (1/1)*(transpose matrix apply(BGT, p -> writeInGTBasisA(AYk(V,gtpA(p),k),BGT)))
);


end


load "Proposition_A.4_statement_for_Bstarplus.m2";




-- Example 1
g = simpleLieAlgebra("A",3);
lambda = {2,0,0};
Vlambda = irreducibleLieAlgebraModule(lambda,g);
N = dim Vlambda;
LAB = lieAlgebraBasis(g);
lambdastar = starInvolution(lambda,g);
Vlambdastar = irreducibleLieAlgebraModule(lambdastar,g);
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
Blambdastar = gtPatterns("A",dynkinToPartition("A",lambdastar));
sigma = apply(N, i -> first select(N, k -> Blambdastar_k==dualGTPattern(Blambda_i)));
Blambdastar= apply(sigma, j -> Blambdastar_j);


all(for k from 1 to #lambda list (
   N1 = AXkrepresentationMatrix(Vlambda,k,Blambda);
   N2 = AXkrepresentationMatrix(Vlambdastar,k,Blambdastar);
   matApprox(N2,transpose(N1),10^-12)
))

all(for k from 1 to #lambda list (
   N1 = AYkrepresentationMatrix(Vlambda,k,Blambda);
   N2 = AYkrepresentationMatrix(Vlambdastar,k,Blambdastar);
   matApprox(N2,transpose(N1),10^-12)
))





-- Example 2
g = simpleLieAlgebra("A",3);
lambda = {3,2,1};
Vlambda = irreducibleLieAlgebraModule(lambda,g);
N = dim Vlambda;
LAB = lieAlgebraBasis(g);
lambdastar = starInvolution(lambda,g);
Vlambdastar = irreducibleLieAlgebraModule(lambdastar,g);
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
Blambdastar = gtPatterns("A",dynkinToPartition("A",lambdastar));
sigma = apply(N, i -> first select(N, k -> Blambdastar_k==dualGTPattern(Blambda_i)));
Blambdastar= apply(sigma, j -> Blambdastar_j);



all(for k from 1 to #lambda list (
   N1 = AXkrepresentationMatrix(Vlambda,k,Blambda);
   N2 = AXkrepresentationMatrix(Vlambdastar,k,Blambdastar);
   matApprox(N2,transpose(N1),10^-12)
))

all(for k from 1 to #lambda list (
   N1 = AYkrepresentationMatrix(Vlambda,k,Blambda);
   N2 = AYkrepresentationMatrix(Vlambdastar,k,Blambdastar);
   matApprox(N2,transpose(N1),10^-12)
))






-- Example 3
g = simpleLieAlgebra("A",5);
lambda = {1,0,0,0,1};
Vlambda = irreducibleLieAlgebraModule(lambda,g);
N = dim Vlambda;
LAB = lieAlgebraBasis(g);
lambdastar = starInvolution(lambda,g);
Vlambdastar = irreducibleLieAlgebraModule(lambdastar,g);
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
Blambdastar = gtPatterns("A",dynkinToPartition("A",lambdastar));
sigma = apply(N, i -> first select(N, k -> Blambdastar_k==dualGTPattern(Blambda_i)));
Blambdastar= apply(sigma, j -> Blambdastar_j);



all(for k from 1 to #lambda list (
   N1 = AXkrepresentationMatrix(Vlambda,k,Blambda);
   N2 = AXkrepresentationMatrix(Vlambdastar,k,Blambdastar);
   matApprox(N2,transpose(N1),10^-12)
))

all(for k from 1 to #lambda list (
   N1 = AYkrepresentationMatrix(Vlambda,k,Blambda);
   N2 = AYkrepresentationMatrix(Vlambdastar,k,Blambdastar);
   matApprox(N2,transpose(N1),10^-12)
))






-- Example 4
g = simpleLieAlgebra("A",6);
lambda = {2,0,0,0,1,0};
Vlambda = irreducibleLieAlgebraModule(lambda,g);
N = dim Vlambda;
LAB = lieAlgebraBasis(g);
lambdastar = starInvolution(lambda,g);
Vlambdastar = irreducibleLieAlgebraModule(lambdastar,g);
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
Blambdastar = gtPatterns("A",dynkinToPartition("A",lambdastar));
sigma = apply(N, i -> first select(N, k -> Blambdastar_k==dualGTPattern(Blambda_i)));
Blambdastar= apply(sigma, j -> Blambdastar_j);



all(for k from 1 to #lambda list (
   N1 = AXkrepresentationMatrix(Vlambda,k,Blambda);
   N2 = AXkrepresentationMatrix(Vlambdastar,k,Blambdastar);
   matApprox(N2,transpose(N1),10^-12)
))

all(for k from 1 to #lambda list (
   N1 = AYkrepresentationMatrix(Vlambda,k,Blambda);
   N2 = AYkrepresentationMatrix(Vlambdastar,k,Blambdastar);
   matApprox(N2,transpose(N1),10^-12)
))