-- Macaulay2
-- Purpose of this talk:
--  a. Get a feeling for what Macaulay2 can do in this area
--  b. See some of the new features in the latest version 1.7
--  c. Give you some tools to investigate interesting problems

-- Let's investigate some 3-tensor spaces A ** B ** C (A,B,C, vector spaces)
-- 1. Equations of the Segre embedding
kk = ZZ/32003 -- Could take QQ as well.
A = kk^2
B = kk^2
C = kk^2
V = A ** B ** C

R = kk[a_1,a_2,b_1,b_2,c_1,c_2]
S = kk[x_(1,1,1)..x_(2,2,2)]
aa = matrix{{a_1,a_2}}
bb = matrix{{b_1,b_2}}
cc = matrix{{c_1,c_2}}
phi = map(R,S,aa**bb**cc)
I = ker phi
dim I
codim I
degree I
hilbertSeries I
reduceHilbert oo
res I
betti oo

-- Here is a naive way to compute the secant locus, but in general, it doesn't work very well!
-- In this case, sec_2(P1 x P1 x P2) = codim 2 in P11, complete 
R2 = kk[a_1,a_2,b_1,b_2,c_1,c_2,c_3, d_1,d_2,e_1,e_2,f_1,f_2,f_3,s,t]
S = kk[x_(1,1,1)..x_(2,2,3)]
aa = matrix{{a_1,a_2}}
bb = matrix{{b_1,b_2}}
cc = matrix{{c_1,c_2,c_3}}
a2 = matrix{{d_1,d_2}}
b2 = matrix{{e_1,e_2}}
c2 = matrix{{f_1,f_2,f_3}}
phi = map(R2,S, s**(aa**bb**cc) + t**(a2**b2**c2))
Isec2 = ker phi;
res Isec2
oo.dd_2

-----------------------------
-- Deficient triple Segres --
-----------------------------
-- 2. If dim A = nA, dim B = nB, dim C = nC, then
--    the expected dimension of sec_r(PA x PB x PC) is
--    r * (nA-1 + nB-1 + nC-1 + 1) - 1
-- A variety is called r-deficient, if Sec_r has lower 
-- than the expected dimension (can't have higher).
--
-- Our goal: find all of the deficient triple 
--   Segre's up to some size, e.g. dim A,B,C <= 5, or 10, or 15...

-- Our tool: Terracini's lemma.  
--   If X \subset PV is a projective variety, 
--   and p in X is a smooth point,
--   let T_(X,p) \subset V denote the embedded tangent space of the cone of X, at p, 
--   this is a vector subspace of V
-- Theorem: If p1, ...pr in X are smooth points of X, and P is a smooth point 
--   of Sec_r(X), in the linear space span(p1, ..., pr), then
--   T_(Sec_r(X), P) = T_(X,p1) + ... + T_(X,pr).

-- Our plan: 
--    a. choose random points on X, 
--    b. compute the (embedded) tangent space T_(X,p)
--    c. sum these.
-- The dimension will be heuristically the dimension of Sec_r(X),
-- although it is possible that the dimension is higher than we compute.
-- In practice, this is very rare, especially if we choose a larger finite field, or QQ.

-- Step (a): Sampling from a variety.
-- In this case, it is easy to sample: we choose random points in 
-- PA, PB, PC, and map them into PV.
kk = ZZ/32003
A = kk^2
B = kk^2
C = kk^2
V = A ** B ** C

sampleSegre = (A,B,C) -> (
    kk := ring A;
    a := random(A,kk^1);
    b := random(B,kk^1);
    c := random(C,kk^1);
    (a,b,c, a**b**c)
    )

pt1 = sampleSegre(A,B,C)

-- If the variety is NOT rational, this is much harder to do.
-- Numerical Algebraic Geometry can help out here!

-- Step (b): The embedded tangent space to a point a**b**c of the Segre.
-- This is just the image of the following vector space map
jacobianMap = (pt,A,B,C) -> (
    (pt_0 ** pt_1 ** C) | (pt_0 ** B ** pt_2) | (A ** pt_1 ** pt_2)
    )

jacobianMap(pt1,A,B,C)
-- it should have rank 4 (dim of Segre is 3, and this is the cone over that, so 4):
rank oo

-- Step (c): Using Terracini, find the tangent space to the secant locus:
-- first, do it by hand:
pt2 = sampleSegre(A,B,C)
J = jacobianMap(pt1,A,B,C) | jacobianMap(pt2,A,B,C)
rank oo -- 8, the whole space!
-- so, as expected, dim sec_2(P1 x P1 x P1) = 7 (8-1).
matrix{{jacobianMap(pt1,A,B,C), jacobianMap(pt2,A,B,C)}}
rank oo
-- Now: let's make a function that takes r, and the dimensions of
-- A,B,C, and returns the dimension of the secant locus (of course,
-- this is a random algorithm).  We also pass in the field to use:

dimSecantToSegre = (r,nA,nB,nC,kk) -> (
    A := kk^nA;
    B := kk^nB;
    C := kk^nC;
    pts := for i from 1 to r list sampleSegre(A,B,C);
    tanspaces := for pt in pts list jacobianMap(pt,A,B,C);
    sumtan := matrix{tanspaces};
    rank sumtan - 1
    )

for r from 1 to 10 list dimSecantToSegre(r,3,3,3,kk)

expectedDim = (r,nA,nB,nC) -> (
    dimX := nA-1 + nB-1 + nC-1;
    min(r * (dimX + 1) - 1, nA*nB*nC-1)
    )
deficient = (nA,nB,nC,kk) -> (
    r := 2;
    d := dimSecantToSegre(r,nA,nB,nC,kk);
    expected := expectedDim(r,nA,nB,nC);
    )
for r from 1 to 10 list expectedDim(r,3,3,3)

triples = (N) -> (
    -- returns all triples (i,j,k), with 2 <= i <= j <= k <= N
    flatten for i from 2 to N list 
    flatten for j from i to N list
    flatten for k from j to N list (i,j,k)
    )

triples 5

N = 5
flatten for ijk in triples N list (
    (i,j,k) := ijk;
    topr := ceiling((i*j*k-1)/(i+j+k-2));
    for r from 2 to topr list (
        expected := expectedDim(r,i,j,k);
        d := dimSecantToSegre(r,i,j,k,kk); 
        if d == expected then continue else {i,j,k,r}=>(d, expected)
        )
    )
netList oo

-- now redo:
N = 9
-- <retype in code above>
#oo
netList ooo

-------------------------------------------
-- Equations of secand loci, flattenings --
-------------------------------------------
-- secant locus sec_2(P2 x P2 x P2 inside P26)
-- notice: this one has the expected dim 2*7-1 = 13.  (but Sec_4 does not).
kk = ZZ/32003
R = kk[x_(1,1,1)..x_(3,3,3)]
X3 = matrix for i in (1,1)..(3,3) list for j in 1..3 list x_(i#0,i#1,j)
X2 = matrix for i in (1,1)..(3,3) list for j in 1..3 list x_(i#0,j,i#1)
X1 = matrix for i in (1,1)..(3,3) list for j in 1..3 list x_(j,i#0,i#1)
sec2 = minors(3,X1) + minors(3,X2) + minors(3,X3);
codim sec2 -- codim 13, therefore dim 13, which is the expected dim: 2*7-1
betti sec2
sec2 = trim sec2;
betti sec2
  hilbertSeries sec2
  reduceHilbert oo
  degree sec2

  needsPackage "MinimalPrimes"
  elapsedTime minprimes(sec2);
  #oo == 1
  ooo_0 == sec2
  -- One minimal prime, which matches sec2, therefore sec2 is prime.
  installMinprimes()
  isPrime sec2

-- Sec_3, Sec_4, Sec_5: there is knowledge of the ideals of these 
-- (see notes of Luke Oeding).  BUT: some of the equations are 
-- interesting and slightly obscure.  Enter representation theory!

--------------------------------------------
-- SchurRings and plethysm -----------------
--------------------------------------------
needsPackage "SchurRings"
A = schurRing(QQ,s,4)
s_{3}
dim s_{3} -- 20 (3rd symmetric power of 4 diml vector space)
dim s_1
s_1 * s_1 
s_2 * s_2
(dim s_2)^2
(dim s_4, dim s_{3,1}, dim s_{2,2})
symmetricPower(2,s_{2})
symmetricPower(3,s_{2})
symmetricPower(4,s_{2})
toE s_{2}
toH s_{2}
toP s_{2}
toE p_3
toH p_3

-------------------------------------------------------------
-- In our case, we often have GL(A) x GL(B) x GL(C) acting.--
-- How to deal with that?                                  --
-------------------------------------------------------------
B = schurRing(A,t,3)
-- Let's verify Cauchy's formula for how S^d(V ** W), 
-- exterior^d(V ** W) break up into irreducibles:
symmetricPower(3, s_1 * t_1)
symmetricPower(4, s_1 * t_1)
exteriorPower(3, s_1 * t_1)

-------------------------------------------------------------
-- Permanental ideals ---------------------------------------
-------------------------------------------------------------
restart
needsPackage "MinimalPrimes"
R = QQ[vars(0..11)]
M = genericMatrix(R,a,3,4)
I = permanents(3,M)
C = minprimes I -- this is using the new minprimes code
#C
C = C/trim;
netList C
J = C_7
netList J_*
codim J -- 4
degree J -- 66
res J
betti res J
  -- Question: understand this component!
viewHelp "Permanents"
--------------------------------------------
-- Annihilators: Macaulay dual, apolarity---
--------------------------------------------
restart
R = QQ[x,y]
F = 193*x^3-1002*x^2*y+1776*x*y^2-1021*y^3 -- not so random cubic

-- we will compute the vector space of polynomials p(x,y)
-- of degree 2 s.t. 
--   p(del_x, del_y)(F) = 0
B = basis(2,R)
D = diff(B, F)
coefficients D
coeffs = last coefficients(D)
-- linear relations on the columns gives a poly of degree 2 which annihilates F
syz coeffs
duals = B * (syz coeffs)
J = ideal duals
diff(gens J, F)

-- Let's make a function out of this
annd = method()
annd(ZZ, RingElement) := (d, F) -> (
    B := basis(d,R);
    D := diff(B, F);
    coeffs := last coefficients(D);
    duals := B * (syz coeffs);
    ideal duals
    )

annd(2,F)
annd(1,F)
annd(4,F)

-- back to the example in hand:
J = annd(2,F)
factor J_0
-- 2x-y annihilates the linear form x+2y, 
-- and 7x+4y annihilates 4x-74
diff(2*x-y, x+2*y)
diff(7*x+4*y, 4*x-7*y)
diff(2*x-y, (x+2*y)^10)
-- Sylvester's method to solve a univariate cubic proceeds as follows:
-- F = u L1^3 + v L2^3, for some u,v, where L1,L2 are above linear forms
-- Solve this:
F // matrix{{(x+2*y)^3,(4*x-7*y)^3}}
-- so 
F == (x+2*y)^3 + 3 * (4*x-7*y)^3

-- Let's apply this to cubics in more variables.
-- Here is an interesting special case (5 variables).
-- see Alexander-Hirshowitz, Ranestad-Schreyer.
R = QQ[a,b,c,d,e]
N = 4
lins = for i from 1 to N list random(1,R)
F = sum for f in lins list f^3
annd(1,F)
J = annd(2,F)
netList J_*
res J
betti res J

-- now do: N=5, 6, 7, 8
N = 5
N = 6 -- 9 linear 2nd syzygies
N = 7 -- 3 linear 2nd syzygies
N = 8 -- general!
  -- in fact sec_7(veronese(3,P^4) in P^34) has expected dimension 7*5-1 = 34
  -- BUT: it is not! it has dim 33, and it is on the small list of AH exceptions.

-------------------------------------
-- Numerical computations -----------
-------------------------------------
-- A. lapack: eigenvalues and eigenvectors
restart
kk = CC_53
M = random(kk^4, kk^4)
eigenvalues M
eigenvectors M
SVD M
det M
M^-1 * M
clean(1.0e-10, oo)

M = mutableMatrix(kk, 500, 500);
fillMatrix M;
elapsedTime SVD M;
elapsedTime eigenvalues M;
elapsedTime eigenvectors M;
B = mutableMatrix(kk,500, 3)
fillMatrix B;
elapsedTime solve(M,B);
oo

needsPackage "NumericalAlgebraicGeometry"
RC = CC_53[x]
F = x^10-3*x^9-4*x^7-x^6+13*x^5-1
solveSystem{F}
oo/coordinates//sort//matrix
-- let's analyze this using numerical algebraic geometry methods

restart
needsPackage "NumericalAlgebraicGeometry"
needsPackage "MinimalPrimes"
ndecomp = numericalIrreducibleDecomposition
NAGtrace 1
R = QQ[vars(0..11)]
M = genericMatrix(R,a,3,4)
I = permanents(3,M)
C = minprimes I
RC = CC[gens R]
IC = sub(I,RC)
elapsedTime W = ndecomp IC
peek W
components W
W7 = last components W
sample W7