Merge pull request #3129 from mahrud/VirtualResolutions

VirtualResolutions minor update

M2/Macaulay2/m2/matrix.m2

@@ -634,12 +634,13 @@ inducedMap(Module,Module,Matrix) := Matrix => opts -> (N',M',f) -> (
      N := target f;
      M := source f;
      if ring N' =!= ring M' or ring N' =!= ring f then error "inducedMap: expected modules and map over the same ring";
-     if isFreeModule N and isFreeModule M and (N =!= ambient N' and rank N === rank ambient N' or M =!= ambient M' and rank M === rank ambient M')
+    if isFreeModule N and isFreeModule M and (
+	N != ambient N' and rank N === rank ambient N' or
+	M != ambient M' and rank M === rank ambient M')
      then f = map(N = ambient N', M = ambient M', f)
      else (
-     	  if ambient N' =!= ambient N then error "inducedMap: expected new target and target of map provided to be subquotients of same free module";
-     	  if ambient M' =!= ambient M then error "inducedMap: expected new source and source of map provided to be subquotients of same free module";
-	  );
+	if ambient N' != ambient N then error "inducedMap: expected new target and target of map provided to be subquotients of same free module";
+	if ambient M' != ambient M then error "inducedMap: expected new source and source of map provided to be subquotients of same free module");
      c := runHooks((inducedMap, Module, Module, Matrix), (opts, N', M', f));
      (f', g, gbN', gbM) := if c =!= null then c else error "inducedMap: no method implemented for this type of input";
      if opts.Verify then (

M2/Macaulay2/m2/matrix2.m2

@@ -389,8 +389,7 @@ remainder'(Matrix,Matrix) := Matrix => (f,g) -> (
      or not isFreeModule source g or not isFreeModule source g then error "expected maps between free modules";
      dual remainder(dual f, dual g))
 remainder(Matrix,Matrix) := Matrix % Matrix := Matrix => (n,m) -> (
-     R := ring n;
-     if target m =!= target n then error "expected matrices with the same target";
+     if target m != target n then error "expected matrices with the same target";
      if not isFreeModule source n or not isFreeModule source m then error "expected maps from free modules";
      if not isQuotientModule target m then error "expected maps to a quotient module";
      c := runHooks((remainder, Matrix, Matrix), (n, m));

M2/Macaulay2/m2/ringmap.m2

@@ -213,7 +213,9 @@ RingMap Module := Module => (f, M) -> (
 -- misc
 tensor(RingMap, Module) := Module => {} >> opts -> (f, M) -> (
     if source f =!= ring M then error "expected module over source ring";
-    cokernel f presentation M);
+    subquotient(f ambient M,
+	if M.?generators then f M.generators,
+	if M.?relations  then f M.relations))
 RingMap ** Module := Module => (f, M) -> tensor(f, M)
 
 tensor(RingMap, Matrix) := Matrix => {} >> opts -> (f, m) -> (

M2/Macaulay2/packages/FourTiTwo.m2

@@ -49,6 +49,7 @@ if not programPaths#?"4ti2" and FourTiTwo#Options#Configuration#"path" != ""
     then programPaths#"4ti2" = FourTiTwo#Options#Configuration#"path"
 
 fourTiTwo = null
+debugLimit = 5
 
 run4ti2 = (exe, args) -> (
     if fourTiTwo === null then
@@ -110,7 +111,7 @@ toBinomial(Matrix,Ring) := (M,S) -> (
 toricMarkov = method(Options=> {InputType => null})
 toricMarkov Matrix := Matrix => o -> (A) -> (
      filename := getFilename();
-     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
+     if debugLevel >= debugLimit then << "using temporary file name " << filename << endl;
      if o.InputType === "lattice" then
      	  F := openOut(filename|".lat")
      else 
@@ -125,7 +126,7 @@ toricMarkov(Matrix,Ring) := o -> (A,S) -> toBinomial(toricMarkov(A,o), S)
 toricGroebner = method(Options=>{Weights=>null})
 toricGroebner Matrix := o -> (A) -> (
      filename := getFilename();
-     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
+     if debugLevel >= debugLimit then << "using temporary file name " << filename << endl;
      F := openOut(filename|".mat");
      putMatrix(F,A);
      close F;
@@ -141,7 +142,7 @@ toricGroebner(Matrix,Ring) := o -> (A,S) -> toBinomial(toricGroebner(A,o), S)
 toricCircuits = method()
 toricCircuits Matrix := Matrix => (A ->(
      filename := getFilename();
-     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
+     if debugLevel >= debugLimit then << "using temporary file name " << filename << endl;
      F := openOut(filename|".mat");
      putMatrix(F,A);
      close F;
@@ -152,7 +153,7 @@ toricCircuits Matrix := Matrix => (A ->(
 toricGraver = method()
 toricGraver Matrix := Matrix => (A ->(
      filename := getFilename();
-     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
+     if debugLevel >= debugLimit then << "using temporary file name " << filename << endl;
      F := openOut(filename|".mat");
      putMatrix(F,A);
      close F;
@@ -164,7 +165,7 @@ toricGraver (Matrix,Ring) := Ideal => ((A,S)->toBinomial(toricGraver(A),S))
 hilbertBasis = method(Options=> {InputType => null})
 hilbertBasis Matrix := Matrix => o -> (A ->(
      filename := getFilename();
-     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
+     if debugLevel >= debugLimit then << "using temporary file name " << filename << endl;
      if o.InputType === "lattice" then
      	  F := openOut(filename|".lat")
      else 
@@ -177,7 +178,7 @@ hilbertBasis Matrix := Matrix => o -> (A ->(
 
 rays Matrix := Matrix => (A ->(
      filename := getFilename();
-     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
+     if debugLevel >= debugLimit then << "using temporary file name " << filename << endl;
      F := openOut(filename|".mat");
      putMatrix(F,A);
      close F;
@@ -193,7 +194,7 @@ rays Matrix := Matrix => (A ->(
 toricGraverDegrees = method()
 toricGraverDegrees Matrix := Matrix => (A ->(
      filename := getFilename();
-     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
+     if debugLevel >= debugLimit then << "using temporary file name " << filename << endl;
      F := openOut(filename|".mat");
      putMatrix(F,A);
      close F;

M2/Macaulay2/packages/LinearTruncations.m2

@@ -303,6 +303,9 @@ isQuasiLinear(List,Module) := opts -> (d,M) -> (
     irr := if opts.IrrelevantIdeal =!= null
     then opts.IrrelevantIdeal else irrelevantIdeal ring M;
     N := truncate(d,M);
+    -- first check whether N is generated only by degree d elements
+    if {d} =!= unique degrees source mingens N then return false;
+    -- then check if the rest of the resolution is quasilinear
     degs := supportOfTor N;
     allowed := supportOfTor comodule irr;
     if (len := #degs) > #allowed then return false;

M2/Macaulay2/packages/VirtualResolutions.m2

@@ -489,8 +489,8 @@ isComputationDone = method(TypicalValue => Boolean, Options => true)
 isComputationDone MultigradedRegularityComputation := Boolean => options multigradedRegularity >> opts -> container -> (
     -- this function determines whether we can use the cached result, or further computation is necessary
     instance(container.Result, List)
-    and (opts.LowerLimit === null or container.LowerLimit <= opts.LowerLimit)
-    and (opts.UpperLimit === null or container.UpperLimit >= opts.UpperLimit))
+    and container.LowerLimit <= min \ transpose{container.LowerLimit, opts.LowerLimit}
+    and container.UpperLimit >= max \ transpose{container.UpperLimit, opts.UpperLimit})
 
 cacheHit := type -> if debugLevel > 0 then printerr("Cache hit on a ", synonym type, "! 🎉");
 
@@ -528,8 +528,16 @@ multigradedRegularityHelper = (X, S, M, opts) -> (
     else if S.?TateData then X = normalToricVarietyFromTateData S
     -- start from module over multigraded polynomial ring and get module over Cox ring of a product of toricProjectiveSpaces
     else X = normalToricVarietyFromTateData imbueRingWithTateData S;
+    r := regularity M;
+    n := degreeLength ring M;
     -- the multigraded regularity of the zero module is -infinity in every component
-    if M == 0 then return {toList(degreeLength ring X : -infinity)};
+    -- TODO: use Hilbert polynomial to detect irrelevant modules quickly
+    if M == 0 then return {toList(n : -infinity)};
+    opts = opts ++ {
+	-- from Proposition 3.7 of [BCHS22] we know reg M \subset mindegs + Eff X
+	LowerLimit => if opts.LowerLimit =!= null then opts.LowerLimit else compMin degrees M,
+	-- Note: an upper limit that works for all examples isn't known
+	UpperLimit => if opts.UpperLimit =!= null then opts.UpperLimit else compMax join(degrees M, {toList(n : r)})};
     -- store a cached computation object in M
     --   MultigradedRegularityOptions{} => MultigradedRegularityComputation{ LowerLimit, UpperLimit, Result }
     container := new MultigradedRegularityComputation from M;
@@ -543,13 +551,15 @@ multigradedRegularityHelper = (X, S, M, opts) -> (
     then error("no applicable strategy for ", toString key)
     else error("assumptions for computing multigraded regularity with strategy ", toString strategy, " are not met"))
 
+-- This is the old strategy for products of projective spaces.
+-- It is based on a direct sheaf cohomology calculation.
+-- See [ABLS20]: https://msp.org/jsag/2020/10-1/p06.xhtml
 multigradedRegularityCohomologySearchStrategy = (X, M, opts) -> (
     S := ring X;
     -- TODO: also check that X and S are indeed a product of
     -- projective spaces and its Cox ring, otherwise return null
     if instance(M, Ideal) then M = comodule M;
     if ring M =!= S then M = map(S, ring M, gens S) ** M;
-    -- This is the default strategy, outlined in https://msp.org/jsag/2020/10-1/p06.xhtml
     debugInfo := if debugLevel < 1 then identity else printerr;
     -- For products of projective space, the dimension is the
     -- number of variables minus the rank of the Picard group
@@ -566,21 +576,17 @@ multigradedRegularityCohomologySearchStrategy = (X, M, opts) -> (
         "reg M = " | toString r,
         "mindegs = " | toString mindegs};
     H := hilbertPolynomial(X, M);
-    HP := x -> (map(QQ, ring H, x))(H);
     debugInfo \ {
 	"HP M = " | toString H,
 	"degs = " | toString degs};
-    -- TODO: why is this the right upper bound?
-    high := if opts.UpperLimit =!= null then opts.UpperLimit else apply(n, i -> max({r} | degs / (deg -> deg_i)));
-    -- TODO: why is mindegs - toList(n:d) the right lower bound?
-    low  := if opts.LowerLimit =!= null then opts.LowerLimit else mindegs - toList(n:d);
-    --
+    (low, high) := (opts.LowerLimit, opts.UpperLimit);
     debugInfo("Computing cohomologyHashTable from ", toString low, " to ", toString high);
     L := pairs cohomologyHashTable(M, low, high);
     --
     gt := new MutableHashTable;
     debugInfo("Beginning search in Picard group");
     -- TODO: rewrite this loop
+    HP := x -> (map(QQ, ring H, x))(H);
     apply(L, ell -> (
             -- Check that Hilbert function and Hilbert polynomial match
             -- (this imposes a condition on the alternating sum of local cohomology dimensions)

M2/Macaulay2/packages/VirtualResolutions/development.m2

@@ -31,6 +31,10 @@ isChiH0(List, Module) := opts -> (d, M) -> (
     H := hilbertPolynomial(variety ring opts.IrrelevantIdeal, M);
     hilbertFunction(d, M) == sub(H, sub(matrix{d}, QQ)))
 
+
+-- This is the new strategy for products of projective spaces.
+-- It is based on quasilinearity of truncations of modules.
+-- See [BCHS21]: https://arxiv.org/abs/2110.10705
 multigradedRegularityTruncationSearchStrategy = (X, M, opts) -> (
     S := ring X;
     -- TODO: also check that X and S are indeed a product of
@@ -56,19 +60,16 @@ multigradedRegularityTruncationSearchStrategy = (X, M, opts) -> (
     debugInfo \ {
         "HP M = " | toString H,
         "degs = " | toString degs};
-    -- TODO: why is this the right upper bound?
-    high := if opts.UpperLimit =!= null then opts.UpperLimit else apply(n, i -> max({r} | degs / (deg -> deg_i)));
-    -- this is just used for shifting the degrees
-    low  := mindegs - toList(n:d);
+    (low, high) := (opts.LowerLimit, opts.UpperLimit);
     debugInfo("Searching for the quasi-linear region from " | toString low | " to " | toString high);
     -- the combinatorial upperbound on regularity from betti numbers
     U0 := regularityBound M;
     if debugLevel > 2 then plotRegion(U0, low, high);
+    -- lowerbound from degrees where H_B^1 would need to vanish if M was regular
+    L0 := findRegion({low, high}, M, isChiH0,       Inner => U0, IrrelevantIdeal => B);
+    if debugLevel > 2 then plotRegion(L0, low, high);
     -- extend U0 to degrees where the truncation is quasi-linear (see Theorem 2.9 of BES)
-    U0  = findRegion({mindegs, high}, M, isQuasiLinear, Inner => U0, IrrelevantIdeal => B);
-    if debugLevel > 2 then plotRegion(U0, low, high);
-    -- limit U0 to degrees where H_B^1 vanishes
-    U0  = findRegion({mindegs, high}, M, isChiH0,       Outer => U0, IrrelevantIdeal => B);
+    U0  = findRegion({low, high}, M, isQuasiLinear, Inner => U0, IrrelevantIdeal => B, Outer => L0);
     if debugLevel > 2 then plotRegion(U0, low, high);
     debugInfo("Upper bound from LinearTruncations: " | toString U0);
     -* old code before the conjecture was proved

M2/Macaulay2/packages/VirtualResolutions/helpers.m2

@@ -1,7 +1,7 @@
 plotRegion = method()
 plotRegion(Function, List, List) := (func, low, high) -> printerr netList(Boxes => false,
-    table(min(high - low) + 1, max(high - low) + 1,
-	(i, j) -> if func(j + first low, first high - i) then "." else "x"))
+    table(last(high - low) + 1, first(high - low) + 1,
+	(i, j) -> if func(j + first low, last high - i) then "." else "x"))
 plotRegion(List, List, List) := (L, low, high) -> plotRegion(
     (i, j) -> any(L, ell -> i >= ell_0 and j >= ell_1), low, high)
 
@@ -42,7 +42,7 @@ normalToricVarietyWithTateData = X -> (
 -- input: multigraded polynomial ring with Tate Data
 -- output: NormalToricVariety, with the given ring cached in it
 normalToricVarietyFromTateData = S -> (
-    if S.?variety and S.?variety.?ring and S.?variety.cache.ring.?TateData then return S.variety;
+    if S.?variety and S.variety.cache.?ring and S.variety.cache.ring.?TateData then return S.variety;
     if not S.?TateData then error "expected a ring with TateData";
     X := cartesianProduct apply(toSequence dimVector S, n -> toricProjectiveSpace(n, CoefficientRing => coefficientRing S));
     X.cache.ring = S; S.variety = X)

M2/Macaulay2/packages/VirtualResolutions/tests.m2

@@ -54,6 +54,13 @@ TEST ///
     S = ZZ/32003[x_0,x_1,x_2,x_3,x_4, Degrees=>{2:{1,0},3:{0,1}}];
     irr = intersect(ideal(x_0,x_1),ideal(x_2,x_3,x_4));
     I = ideal(x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2, x_0^3*x_4+x_1^3*(x_2+x_3));
+    --
+    d1 = matrix{{x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2,
+        x_0*x_1*x_2^3+x_0*x_1*x_2^2*x_3-x_0^2*x_3^2*x_4+x_1^2*x_2*x_4^2+x_1^2*x_3*x_4^2,
+        x_1^2*x_2^3+x_1^2*x_2^2*x_3-x_0*x_1*x_3^2*x_4-x_0^2*x_4^3}};
+    C = chainComplex({d1});
+    elapsedTime assert(isVirtual(irr,C) == false) -- 0.022
+    --
     d1 = matrix{{x_1^3*x_2+x_1^3*x_3+x_0^3*x_4,
             x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2,
             x_0*x_1*x_2^3+x_0*x_1*x_2^2*x_3-x_0^2*x_3^2*x_4+x_1^2*x_2*x_4^2+x_1^2*x_3*x_4^2,
@@ -66,28 +73,6 @@ TEST ///
     elapsedTime assert(isVirtual(irr,C) == true) -- 0.028
 ///
 
-TEST ///
-    S = ZZ/32003[x_0,x_1,x_2,x_3,x_4, Degrees=>{2:{1,0},3:{0,1}}];
-    irr = intersect(ideal(x_0,x_1),ideal(x_2,x_3,x_4));
-    I = ideal(x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2, x_0^3*x_4+x_1^3*(x_2+x_3));
-    d1 = matrix{{x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2,
-        x_0*x_1*x_2^3+x_0*x_1*x_2^2*x_3-x_0^2*x_3^2*x_4+x_1^2*x_2*x_4^2+x_1^2*x_3*x_4^2,
-        x_1^2*x_2^3+x_1^2*x_2^2*x_3-x_0*x_1*x_3^2*x_4-x_0^2*x_4^3}};
-    C = chainComplex({d1});
-    elapsedTime assert(isVirtual(irr,C) == false) -- 0.022
-///
-
-TEST ///
-    S = ZZ/32003[x_0,x_1,x_2,x_3,x_4, Degrees=>{2:{1,0},3:{0,1}}];
-    irr = intersect(ideal(x_0,x_1),ideal(x_2,x_3,x_4));
-    I = ideal(x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2, x_0^3*x_4+x_1^3*(x_2+x_3));
-    d1 = matrix{{x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2,
-        x_0*x_1*x_2^3+x_0*x_1*x_2^2*x_3-x_0^2*x_3^2*x_4+x_1^2*x_2*x_4^2+x_1^2*x_3*x_4^2,
-        x_1^2*x_2^3+x_1^2*x_2^2*x_3-x_0*x_1*x_3^2*x_4-x_0^2*x_4^3}};
-    C = chainComplex({d1});
-    elapsedTime assert(isVirtual(irr,C) == false) -- 0.027
-///
-
 TEST ///
     S = ZZ/101[x_0,x_1,x_2,x_3,x_4, Degrees=>{2:{1,0},3:{0,1}}];
     irr = intersect(ideal(x_0,x_1),ideal(x_2,x_3,x_4));
@@ -97,15 +82,6 @@ TEST ///
 ///
 
 ----- Tests for idealSheafGens
-TEST ///
-    debug needsPackage "VirtualResolutions"
-    S = ZZ/32003[x_0,x_1,x_2,x_3,x_4, Degrees=>{2:{1,0},3:{0,1}}];
-    irr = intersect(ideal(x_0,x_1),ideal(x_2,x_3,x_4));
-    I = ideal(x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2, x_0^3*x_4+x_1^3*(x_2+x_3));
-    J = ourSaturation(I,irr);
-    elapsedTime assert(idealSheafGens(2,J,irr) == {I}) -- 4.2
-///
-
 TEST ///
     debug needsPackage "VirtualResolutions"
     S = ZZ/32003[x_0,x_1,y_0,y_1, Degrees=>{2:{1,0},2:{0,1}}];
@@ -122,6 +98,7 @@ TEST ///
     irr = intersect(ideal(x_0,x_1),ideal(x_2,x_3,x_4));
     I = ideal(x_0^2*x_2^2+x_1^2*x_3^2+x_0*x_1*x_4^2, x_0^3*x_4+x_1^3*(x_2+x_3));
     J = ourSaturation(I,irr);
+    elapsedTime assert(idealSheafGens(2,J,irr) == {I}) -- 4.2
     elapsedTime output = idealSheafGens(2,J,irr,GeneralElements=>true); -- 6.8
     elapsedTime assert(J == ourSaturation(output_0, irr)) -- 0.06
 ///
@@ -173,6 +150,11 @@ TEST ///
     elapsedTime assert(multigradedRegularity(S, I) == {{2,2},{4,1},{1,5}}) -- woohoo cache hit!!
     -- test for weird ring problems
     assert(ring x_0 === value getSymbol "S")
+    -- test for slow resolutions
+    elapsedTime M = module saturate(I^2, B); -- ~0.3s
+    elapsedTime assert(multigradedRegularity(S, M) == {{6, 4}, {5, 6}, {4, 8}}) -- ~19s
+    -- elapsedTime M = module saturate(I^3, B); -- ~0.8s
+    -- elapsedTime assert(multigradedRegularity(S, M) == {{8, 7}, {9, 6}, {7, 9}, {6, 11}}) -- ~26min
 ///
 
 TEST ///
@@ -316,7 +298,7 @@ TEST ///
   assert(multigradedRegularity(S, comodule I) == {{0,4}})
 ///
 
-TEST /// -- test of returning -infinity for irrelevant ideals
+/// -- test of returning -infinity for irrelevant ideals
   debug needsPackage "VirtualResolutions"
   X = toricProjectiveSpace(1)**toricProjectiveSpace(2);
   --X = normalToricVarietyWithTateData X