needsPackage("LieAlgebraRepresentations");
debug LieAlgebraRepresentations;


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.7.m2";


-- Example 1
g = simpleLieAlgebra("A",3);
lambda = {2,0,0};
Vlambda = irreducibleLieAlgebraModule(lambda,g);
LAB = lieAlgebraBasis(g);
time Llambda = GTrepresentationMatrices(Vlambda);
rhoVlambda = lieAlgebraRepresentation(Vlambda,LAB,Llambda);
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
N = dim Vlambda;
MC = apply(Blambda, p -> molevCoefficient(gtp(p)));
Q1 = diagonalMatrix apply(#MC, i-> sqrt(1/(MC_i)));
LOI = LAB#"LoweringOperatorIndices";
ROI = LAB#"RaisingOperatorIndices";

-- Check the X_(a_k)
all(#lambda, k -> matApprox(inverse(Q1)*(Llambda_(ROI_k))*Q1,AXkrepresentationMatrix(Vlambda,k+1,Blambda),10^-12))
-- Check the Y_(a_k)
all(#lambda, k -> matApprox(inverse(Q1)*(Llambda_(LOI_k))*Q1,AYkrepresentationMatrix(Vlambda,k+1,Blambda),10^-12))





-- Example 2
g = simpleLieAlgebra("A",3);
lambda = {3,2,1};
Vlambda = irreducibleLieAlgebraModule(lambda,g);
LAB = lieAlgebraBasis(g);
time Llambda = GTrepresentationMatrices(Vlambda);
rhoVlambda = lieAlgebraRepresentation(Vlambda,LAB,Llambda);
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
N = dim Vlambda;
MC = apply(Blambda, p -> molevCoefficient(gtp(p)));
Q1 = diagonalMatrix apply(#MC, i-> sqrt(1/(MC_i)));
LOI = LAB#"LoweringOperatorIndices";
ROI = LAB#"RaisingOperatorIndices";

-- Check the X_(a_k)
all(#lambda, k -> matApprox(inverse(Q1)*(Llambda_(ROI_k))*Q1,AXkrepresentationMatrix(Vlambda,k+1,Blambda),10^-12))
-- Check the Y_(a_k)
all(#lambda, k -> matApprox(inverse(Q1)*(Llambda_(LOI_k))*Q1,AYkrepresentationMatrix(Vlambda,k+1,Blambda),10^-12))





-- Example 3
g = simpleLieAlgebra("A",5);
lambda = {1,0,0,0,1};
Vlambda = irreducibleLieAlgebraModule(lambda,g);
LAB = lieAlgebraBasis(g);
time Llambda = GTrepresentationMatrices(Vlambda);
rhoVlambda = lieAlgebraRepresentation(Vlambda,LAB,Llambda);
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
N = dim Vlambda;
MC = apply(Blambda, p -> molevCoefficient(gtp(p)));
Q1 = diagonalMatrix apply(#MC, i-> sqrt(1/(MC_i)));
LOI = LAB#"LoweringOperatorIndices";
ROI = LAB#"RaisingOperatorIndices";

-- Check the X_(a_k)
all(#lambda, k -> matApprox(inverse(Q1)*(Llambda_(ROI_k))*Q1,AXkrepresentationMatrix(Vlambda,k+1,Blambda),10^-12))
-- Check the Y_(a_k)
all(#lambda, k -> matApprox(inverse(Q1)*(Llambda_(LOI_k))*Q1,AYkrepresentationMatrix(Vlambda,k+1,Blambda),10^-12))





-- Example 4
g = simpleLieAlgebra("A",6);
lambda = {2,0,0,0,1,0};
Vlambda = irreducibleLieAlgebraModule(lambda,g);
LAB = lieAlgebraBasis(g);
time Llambda = GTrepresentationMatrices(Vlambda);
rhoVlambda = lieAlgebraRepresentation(Vlambda,LAB,Llambda);
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
N = dim Vlambda;
MC = apply(Blambda, p -> molevCoefficient(gtp(p)));
Q1 = diagonalMatrix apply(#MC, i-> sqrt(1/(MC_i)));
LOI = LAB#"LoweringOperatorIndices";
ROI = LAB#"RaisingOperatorIndices";

-- Check the X_(a_k)
all(#lambda, k -> matApprox(inverse(Q1)*(Llambda_(ROI_k))*Q1,AXkrepresentationMatrix(Vlambda,k+1,Blambda),10^-12))
-- Check the Y_(a_k)
all(#lambda, k -> matApprox(inverse(Q1)*(Llambda_(LOI_k))*Q1,AYkrepresentationMatrix(Vlambda,k+1,Blambda),10^-12))