ExternalActivityComplexes.m2

newPackage(
    "ExternalActivityComplexes",
    AuxiliaryFiles => false,
    Version => "0.1", 
    Date => "",
    Authors => {
        {Name => "Ayah Almousa", 
            Email => "aalmousa@sc.edu", 
            HomePage => "http://sites.google.com/view/ayah-almousa"}
    },
    Headline => "functions for investigating affine Schubert varieties of a pair of linear spaces",
    PackageExports => {
        "Posets",
	"Matroids",
	"Graphs"
    },
    DebuggingMode => true
)

export{
    "diagonalDilworthTruncation",	  --documented, has tests
    "kempfCollapsing",			  --no doc, no tests
    "affineSchubertVariety",		  --no doc, no tests
    "kempfCollapsingBeta",		  --no doc, no tests
    "kempfCollapsingSQvee",		  --no doc, no tests
    "multigradedKCBeta"			  --no doc, no tests
}
    
------------------------------------------------------------------------------
------------------------------------------------------------------------------
-- **CODE** --
------------------------------------------------------------------------------
------------------------------------------------------------------------------

-----------------------------------------------
--**Part I: Dilworth truncation** --
-----------------------------------------------


diagonalDilworthTruncation = method()

--maybe cache I1 and I2 in a hash table for each C?
diagonalDilworthTruncation(Matroid, Matroid) := Matroid => (M1,M2) -> (
    G := groundSet M1; --assuming they have the same G for now
    possibleCircuits := drop(subsets G,1); --drop empty set
    circuitList := new MutableList;
    for C in possibleCircuits do(
	if (rank(M1,C) + rank(M2,C) == #C) then circuitList#(#circuitList) = C;
	);
    circuitsFinal := apply(toList circuitList, i-> toList i);
    matroid(toList G, circuitsFinal, EntryMode => "circuits")
    )

-----------------------------------------------
--**Part II: Products of Chern Classes** --
-----------------------------------------------

--eventually do this with 2 linear spaces

--inputs:
    --M1: matroid on n elements
    --M2: matroid on n elements
    --Wx: vector in R^n
    --Wy: vector in R^n
    --p : point in R^n/R

-*
--THIS IS WRONG RIGHT NOW
--E currently doesn't distinguish between E1 and E2
--actually I don't really feel like this needs to be a graph
intersectionGraph(List, List, List) := Graph => (p, Wx, Wy) -> (
    --TODO: add check that length p == length Wx == length Wy
    n := #p;
    V1 := unique (p-Wx);
    V2 := unique (p-Wy);
    E := apply(n, i-> (p-Wx)_i,(p-Wy)_i);
    graph E
    )
*-

-----------------------------------------------
--**Part III: Kempf Collapsing functions** --
-----------------------------------------------

kempfCollapsing = method(
    Options => {
	CoefficientRing => QQ
	}
    )

kempfCollapsing(Matrix, Matrix) := o -> (A1,A2) -> (
    if (numcols A1 != numcols A2) then error("A1 and A2 should have the same number of columns.");
    kk := o.CoefficientRing;
    n := numcols A1;
    m1 := numrows A1;
    m2 := numrows A2;
    (t,u,v,x,y) := (getSymbol "t", getSymbol "u", getSymbol "v",getSymbol "x", getSymbol "y");
    R := kk[u_1..u_m1,v_1..v_m2,t_1..t_n];
    S := kk[x_1..x_n,y_1..y_n];
    use R;
    T := diagonalMatrix((gens R)_{numgens R-n..numgens R-1});
    tInvVars := apply(1..n, i-> product select((gens R)_{numgens R-n..numgens R-1}, j-> last baseName j != i));
    Tinv := diagonalMatrix(toList tInvVars);
    L1 := matrix{(gens R)_{0..m1-1}}*A1*T;
    L2 := matrix{(gens R)_{m1..m1+m2-1}}*A2*Tinv;
    phi := map(R,S,L1|L2);
    kernel phi
    );


affineSchubertVariety = method(
    Options => {
	CoefficientRing => QQ
	}
    )

affineSchubertVariety(Matrix, Matrix) := o -> (A1,A2) -> (
    if (numcols A1 != numcols A2) then error("A1 and A2 should have the same number of columns.");
    kk := o.CoefficientRing;
    n := numcols A1;
    m1 := numrows A1;
    m2 := numrows A2;
    (t,u,v,x,y) := (getSymbol "t", getSymbol "u", getSymbol "v", getSymbol "x", getSymbol "y");
    R := kk[u_1..u_m1,v_1..v_m2,t_1..t_n];
    S := kk[x_1..x_n,y_1..y_n];
    use R;
    T := diagonalMatrix((gens R)_{numgens R-n..numgens R-1});
    L1 := matrix{(gens R)_{0..m1-1}}*A1*T;
    L2 := matrix{(gens R)_{m1..m1+m2-1}}*A2*T;
    phi := map(R,S,L1|L2);
    kernel phi
    );


kempfCollapsingBeta = method(
    Options => {
	CoefficientRing => QQ
	}
    )

kempfCollapsingBeta(Matrix, Matrix) := o -> (A1,A2) -> (
    if (numcols A1 != numcols A2) then error("A1 and A2 should have the same number of columns.");
    kk := o.CoefficientRing;
    n := numcols A1;
    m1 := numrows A1;
    m2 := numrows A2;
    (t,u,v,x,y,z) := (getSymbol "t", getSymbol "u", getSymbol "v",getSymbol "x", getSymbol "y", getSymbol "z");
    R := kk[u_1..u_m1,v_1..v_m2,t_1..t_n];
    S := kk[x_1..x_n,y_1..y_n,z_1..z_n];
    use R;
    T := diagonalMatrix((gens R)_{numgens R-n..numgens R-1});
    tInvVars := apply(1..n, i-> product select((gens R)_{numgens R-n..numgens R-1}, j-> last baseName j != i));
    Tinv := diagonalMatrix(toList tInvVars);
    L1 := matrix{(gens R)_{0..m1-1}}*A1*T;
    L2 := matrix{(gens R)_{m1..m1+m2-1}}*A2*Tinv;
    C := matrix {(gens R)_{numgens R-n..numgens R-1}};
    phi := map(R,S,L1|L2|C);
    kernel phi
    );

--kempf collapsing of S ++ Q*

kempfCollapsingSQvee = method (
    Options => {
	CoefficientRing => QQ
	}
    )

kempfCollapsingSQvee(Matrix, Matrix) := o -> (A1,A2) -> (
    if (numcols A1 != numcols A2) then error("A1 and A2 should have the same number of columns.");
    kk := o.CoefficientRing;
    n := numcols A1;
    m1 := numrows A1;
    m2 := numrows A2;
    (t,u,v,x,y) := (getSymbol "t", getSymbol "u", getSymbol "v", getSymbol "x", getSymbol "y");
    R := kk[u_1..u_m1, v_1..v_m2, t_1..t_n];
    S := kk[x_1..x_n, y_1..y_n];
    T := diagonalMatrix((gens R)_{numgens R-n..numgens R-1});
    tInvVars := apply(1..n, i-> product select((gens R)_{numgens R-n..numgens R-1}, j-> last baseName j != i));
    Tinv := diagonalMatrix(toList tInvVars);
    L1 := matrix{(gens R)_{0..(m1-1)}}*A1*T;
    L2 := matrix{(gens R)_{m1..(m1+m2-1)}}*A2*Tinv;
    phi := map(R,S,L1|L2);
    kernel phi
);



------------------------
--test: a ZZ^3-multigraded version of kempf collapsing
--goal: try to use torus grading to guess representations
multigradedKCBeta = method(
    Options => {
	CoefficientRing => QQ
	}
    )

multigradedKCBeta(Matrix, Matrix) := o -> (A1,A2) -> (
    if (numcols A1 != numcols A2) then error("A1 and A2 should have the same number of columns.");
    kk := o.CoefficientRing;
    n := numcols A1;
    m1 := numrows A1;
    m2 := numrows A2;
    (t,u,v,x,y,z) := (getSymbol "t", getSymbol "u", getSymbol "v",getSymbol "x", getSymbol "y", getSymbol "z");
    Dx := n:{1,0,0};
    Dy := n:{0,1,0};
    Dz := n:{0,0,1};
    DegList = toList(Dx)|toList(Dy)|toList(Dz);
    R := kk[u_1..u_m1,v_1..v_m2,t_1..t_n];
    S := kk[x_1..x_n,y_1..y_n,z_1..z_n, Degrees => DegList];
    use R;
    T := diagonalMatrix((gens R)_{numgens R-n..numgens R-1});
    tInvVars := apply(1..n, i-> product select((gens R)_{numgens R-n..numgens R-1}, j-> last baseName j != i));
    Tinv := diagonalMatrix(toList tInvVars);
    L1 := matrix{(gens R)_{0..m1-1}}*A1*T;
    L2 := matrix{(gens R)_{m1..m1+m2-1}}*A2*Tinv;
    C := matrix {(gens R)_{numgens R-n..numgens R-1}};
    phi := map(R,S,L1|L2|C);
    kernel phi
    );

------------------------------------------------------------------------------
------------------------------------------------------------------------------
-- **DOCUMENTATION** --
------------------------------------------------------------------------------
------------------------------------------------------------------------------
beginDocumentation ()    

doc ///
    Key
        ExternalActivityComplexes
    Headline
        a package for investigating affine Schubert varieties of a pair of linear spaces.
    Description
        Text
            This package provides functions for constructing and investigating coordinate rings of
	    affine Schubert varieties associated to a pair of linear spaces and several combinatorial objects
	    associated to them, including their initial ideals and their external activity complexes.
        Text
            @UL {
	    {"[BF24] Andy Berget and Alex Fink, ",
	    HREF("https://arxiv.org/abs/2412.11759", EM "The external activity complex of a pair of matroids"),
	    ", arXiv 2412.11759"},
            }@  
///

doc ///
    Key
        diagonalDilworthTruncation
        (diagonalDilworthTruncation, Matroid, Matroid)
    Headline
        compute diagonal Dilworth truncation of a pair of matroids
    Usage
        diagonalDilworthTruncation(M1,M2)
    Inputs
        M1:Matroid
	M2:Matroid
    Outputs
	:Matroid
    Description
        Text
            Given a pair of matroids, computes their diagonal Dilworth truncation as
	    defined in [BF24]. By example 4.8 in [BF24], the diagonal Dilworth truncation
	    of the fano matroid and the nonfano matroid is the uniform matroid of rank 5
	    on 7 elements.
        Example
	    F = specificMatroid "fano"
	    F' = specificMatroid "nonfano"
	    diagonalDilworthTruncation(F,F') == uniformMatroid(5,7)	    
///

------------------------------------------------------------------------------
------------------------------------------------------------------------------
-- **TESTS** --
------------------------------------------------------------------------------
------------------------------------------------------------------------------


TEST ///
--diagonal Dilworth of any M with U(1,n) should yield M
U35 = uniformMatroid(3,5);
U15 = uniformMatroid(1,5);
D = diagonalDilworthTruncation(U35,U15)
assert(diagonalDilworthTruncation(U35,U15) == U35)
///

TEST ///
--BF24 Example 4.8: Fano matroid and non-Fano matroid give U(5,7)
F7 = specificMatroid "fano"
F7' = specificMatroid "nonfano"
assert(diagonalDilworthTruncation(F7,F7') == uniformMatroid(5,7))
///

end---------------------------------------------------------------------------     

------------------------------------------------------------------------------
------------------------------------------------------------------------------
-- **SCRATCH SPACE** --
------------------------------------------------------------------------------
------------------------------------------------------------------------------

    G = groundSet M1; --assuming they have the same G for now
    possibleCircuits = drop(subsets G,1); --drop empty set
    circuitList = new MutableList;
    circuitDecomps = new MutableHashTable;
    for C in possibleCircuits do(
	if (rank(M1,C) + rank(M2,C) == #C) then (
	    circuitList#(#circuitList) = C
	    (M1C, M2C) = (M1|C, M2|C)
	    basesM1C = apply(bases M1C, i-> M1C_(toList i))
	    basesM2C = apply(bases M2C, i-> M2C_(toList i))
	    for i in basesM1C do(
	        Cdiff = sort toList((set C)-(set i))
		if any(basesM2C, i-> i == Cdiff) then (
		    circuitDecomps#C = (i,Cdiff);
		    break;
		    );
		);
	    );
	);
    circuitsFinal = apply(toList circuitList, i-> toList i);
    D = matroid(toList G, circuitsFinal, EntryMode => "circuits");
    D.cache.circuitDecomposition = circuitDecomps;
    D


    diagonalDilworthTruncation = method()

--maybe cache I1 and I2 in a hash table for each C?
diagonalDilworthTruncation(Matroid, Matroid) := Matroid => (M1,M2) -> (
    G := groundSet M1; --assuming they have the same G for now
    possibleCircuits := drop(subsets G,1); --drop empty set
    circuitList := new MutableList;
    circuitDecomps := new MutableHashTable;
    for C in possibleCircuits do(
	if (rank(M1,C) + rank(M2,C) == #C) then (
	    circuitList#(#circuitList) = C;
	    (M1C, M2C) := (M1|C, M2|C);
	    basesM1C := apply(bases M1C, i-> M1C_(toList i));
	    basesM2C := apply(bases M2C, i-> M2C_(toList i));
	    for i in basesM1C do(
	        Cdiff := sort toList((set C)-(set i));
		if any(basesM2C, i-> i == Cdiff) then (
		    circuitDecomps#C = (i,Cdiff);
		    break;
		    );
		);
	    );
	);
    circuitsFinal := apply(toList circuitList, i-> toList i);
    D := matroid(toList G, circuitsFinal, EntryMode => "circuits");
    D.cache.circuitDecomposition = circuitDecomps;
    D
    )
*-
---------------------
--Ayah's sandbox
---------------------

--------------------------
--smaller example that has the sextic relation
--note that this one is NOT uniform
L' = matrix{
    {1,0,0,1,1,1},
    {0,1,0,1,2,2},
    {0,0,1,1,2,3}
    }
I' = kempfCollapsingBeta(L',L')
transpose mingens I'
g = last flatten entries mingens I'
netList terms g

S = QQ[x_0..x_20];
A = genericSymmetricMatrix(S,6)
B = L'*A*(transpose L')
phi = map(S,S,{x_0,0,0,0,0,0,x_6,0,0,0,0,x_11,0,0,0,x_15,0,0,x_18,0,x_20})
det phi (B)

A = genericSymmetricMatrix(S,6)
B = L'*A*(transpose L')
phi = map(S,S,{x_0,0,0,0,0,0,x_6,0,0,0,0,x_11,0,0,0,x_15,0,0,x_18,0,x_20})
Y = phi (B)



------------
R = QQ[u_1..u_3,v_1..v_3]
T = QQ[l_(1,1)..l_(1,3),l_(2,2),l_(2,3),l_(3,3)]
A = genericSymmetricMatrix(T,3)
u = matrix{{u_1,u_2,u_3}}
v = matrix{{v_1,v_2,v_3}}

X = (transpose u)*v + (transpose v)*u


mapVec = (transpose sym2L')*(vector gens T)
omega = map(T,T,transpose matrix mapVec)
omega (ker psi)

psi = map(R,T,{X_(0,0),X_(0,1),X_(0,2),X_(1,1),X_(1,2),X_(2,2)})
--psi' = map(S,T,{Y_(0,0),Y_(0,1),Y_(0,2),Y_(1,1),Y_(1,2),Y_(2,2)})

kernel psi'

kernel psi

------------------
--finding a uniform matroid closer to the last one
L' = matrix{
    {1,0,0,1,1,1},
    {0,1,0,1,2,3},
    {0,0,1,1,3,4}
    }
I' = kempfCollapsingBeta(L',L')
transpose mingens I'
g = last flatten entries mingens I'
netList terms g
gensI' = (ideal mingens I')_*
g = last gensI'
#terms g
binomial(6+3-1,3)
-----------
-----------
L0 = matrix{
    {1,1,1,1,1,1},
    {1,2,3,4,5,6},
    {1,4,9,16,25,36} --36 vs 37
    };

I0 = kempfCollapsingBeta(L0,L0);
mingens I0

--second veronese

-------



(reduceHilbert hilbertSeries I0, reduceHilbert hilbertSeries I)
M0 = matroid L0; M1 = matroid L;


------------------------------------
--Development Section
------------------------------------

restart
uninstallPackage "ExternalActivityComplexes"
restart
installPackage "ExternalActivityComplexes"
restart
needsPackage "ExternalActivityComplexes"
elapsedTime check "ExternalActivityComplexes"
viewHelp "ExternalActivityComplexes"