Code and implementation

The file makingWaves.m2 implements Algorithm 1 in our paper via the command wavePairs(A, r). The input A is an \(\ell \times k\) matrix in a polynomial ring in \(n\) variables over \(\mathbb{Q}\), and r is an integer between \(0\) and \(n-1\). The output is the saturated ideal of the wave pair variety \(\mathcal{P}_r^A\). This ideal lies in the coordinate ring of \(\operatorname{Gr}(n-r, n) \times \mathbb{P}^{k-1}\). This is represented in Macaulay2 by the ring \(\mathbb{Q}[p_I, z_1,\dotsc, z_k]/G\), where the \(p_I\) represent Plücker coordinates, (\(I\) ranges over all subsets of \(\{1,2,\dotsc,n\}\) of size \(n-r\)) and \(G\) is the ideal of Plücker relations.

By default, the computation is performed globally, and thus can be rather slow. This can be sped up by performing the computation locally on an affine patch of the Grassmannian, which can be achieved by specifying the optional argument Patch => I. Valid types for I are List, ZZ, or Boolean. If \(I\) is list representing a subset of \(\{1,2,\dotsc,n\}\) of size \(n-r\), setting Patch => I will restrict to the affine patch where the coordinate \(p_I = 1\). Alternatively, \(I\) can be an integer between \(0\) and \(\binom{n}{n-r} - 1\). The option Patch => true is equivalent to Patch => 0 and Patch => toList(0..<n-r).

For convenience, we also provide the option ModPlucker => false, which returns the lift of the wave pair ideal in \(\mathbb{Q}[p_I, z_1,\dotsc,z_k]\).

For an example showcasing options and usage, see the page Example computation.

Source code

 1wavePairs = method(Options => { ModPlucker => true, Patch => null })
 2wavePairs(Matrix, ZZ) := o -> (A, r) -> (
 3  if r == 0 then return waveIncidenceRZero(A);
 4  R := ring A;
 5  KK := coefficientRing R;
 6  n := numgens R;
 7  if r < 0 or r >= n then error("expected argument between 0 and " | toString(n-1) | ".");
 8
 9  s := symbol s;
10  w := symbol w;
11  R' := local R';
12  S := local S;
13  if o.Patch === null then (
14    R' = KK[s_(1,1)..s_(n-r, n), w_1..w_(n-r)];
15    S = transpose genericMatrix(R',n,n-r);
16  )
17  else (
18    R' = KK[s_(1,1)..s_(n-r, r), w_1..w_(n-r)];
19    S = transpose genericMatrix(R',r,n-r);
20    Id := map(R'^(n-r), R'^(n-r), 1);
21    patch := null;
22    if instance(o.Patch, List) then (
23      if #o.Patch != n-r then error("option Patch expects a list with " | toString(n-r) | " entries");
24      patch = o.Patch;
25    ) else if instance(o.Patch, ZZ) and o.Patch != 0 then (
26      patch = (subsets(n,n-r))#(o.Patch);
27    );
28
29    if patch =!= null then (
30      S = transpose entries S;
31      Id = transpose entries Id;
32      ids := apply(Id, sort patch, (i,j) -> {i,j});
33      scan(ids, i -> S = insert(last i, first i, S));
34      S = matrix transpose S;
35    ) else S = Id | S;
36  );
37
38  phi := map(R', R, matrix{{w_1..w_(n-r)}} * S);
39
40  k := numColumns A;
41  z := symbol z;
42  T := R'[z_1..z_k];
43  M := sub(phi(A), T) * transpose vars T;
44
45  wVars := take(gens R', -(n-r)) / (wi -> sub(wi, T));
46  J' := ideal last coefficients(M, Variables => wVars);
47
48  -- remove wVars
49  SS := KK[s_(1,1)..s_(n-r,n), z_1..z_k, MonomialOrder => {(n-r)*n, k}];
50  J := sub(J', SS);
51
52  SSbar := SS/J;
53
54  G := Grassmannian(n-r-1, n-1, CoefficientRing => KK);
55  RG := ring G; -- ring of Grassmannian
56
57  T = RG ** KK[z_1..z_k];
58  idxs := gens RG / last @@ baseName / myToList;
59  pMappings := gens RG /
60                (p -> det submatrix(S, myToList last baseName p)) /
61                (f -> sub(f, SSbar));
62  zMappings := take(gens SSbar, -k);
63
64  Tbar := if o.ModPlucker then T/sub(G,T) else T;
65
66  psi := map(SSbar, Tbar, pMappings | zMappings);
67  I := ker psi;
68
69  zIrrelevant := ideal take(gens Tbar, -k);
70  pIrrelevant := ideal take(gens Tbar, #gens Tbar - k);
71
72  saturate(saturate(I, pIrrelevant), zIrrelevant)
73)
74
75waveIncidenceRZero = A -> (
76  R := ring A;
77  KK := coefficientRing R;
78  k := numColumns A;
79
80  z := symbol z;
81  S := R[z_1..z_k];
82
83  M := vars S * transpose sub(A, S);
84  MM := last coefficients(M, Variables => gens R / (x -> sub(x, S)));
85
86  T := KK[z_1..z_k];
87  saturate(ideal sub(MM, T), ideal gens T)
88)
89
90myToList = L -> (if instance(L,ZZ) then {L} else toList L)