HsurfHomology.m2

newPackage("HsurfHomology")
export{
  "countCritPoints",
  "prepareIdeals",
  "findImageEquation",
  "partials",

  "compareSaturation",
  "hasRadicalCritPoints"
}

needsPackage "Depth"

countCritPoints = method()
    ---- desc: given an unfolding, count the no. of its divergent critical points
    ----
    ---- ii01 : countCritPoints(exampleB1())
    ---- oo01 = 2
countCritPoints(HashTable) := (unfolding) -> (
  -- find the ideal of Sigma (trajectories of divergent critical points)
  Sigma := saturate prepareIdeals(unfolding, "default");
  print "found the ideal of Sigma";

  -- verify the Cohen-Macaulay property
  targetRing := unfolding#"R" / ideal(unfolding#"sourceCoords");
  Sigma = substitute(Sigma, targetRing);
  if(not isCM(targetRing / Sigma)) then ( print "not Cohen-Macaulay"; return -7; );
  print "checked for Cohen-Macaulay";

  -- calculate and return the number of divergent critical points
  return degree(targetRing / (Sigma + ideal(unfolding#"targetTime")))
)


prepareIdeals = method()
    ---- desc: prepare the ideals critPoints and zeroLocus which give Sigma (trajectories of DCP)
    ---- 
    ---- 3 strategies are available:
    ----   default: zeroLocus is generated by the equation G
    ----   partials:  ..by the full Jacobian of G
    ----   radical partials: .. by the radical of the ideal generated by the full Jacobian of G
    ----
    ---- ii01 : prepareIdeals(exampleA1())
    ---- ii02 : prepareIdeals(exampleA1(), "radical_partials")
prepareIdeals(HashTable, String) := (unfolding, strategy) -> (
  equation := findImageEquation(unfolding);
  critPoints := ideal partials({equation}, unfolding#"targetSpace"); 
  
  -- default strategy
  zeroLocus := ideal equation;

  if strategy == "partials" then (
    targetCoords := join(unfolding#"targetSpace", unfolding#"targetTime");
    zeroLocus = ideal partials({equation}, targetCoords); 
  );

  if strategy == "radical_partials" then (
    targetCoords = join(unfolding#"targetSpace", unfolding#"targetTime");
    zeroLocus = radical ideal partials({equation}, targetCoords); 
  );

  print "successfully prepared ideals";
  return (critPoints, zeroLocus)
)


findImageEquation = method()
    ---- desc: given an unfolding, return the equation of its image
    ----
    ---- ii01 : findImageEquation(exampleA1())
    ---- oo01 = t^2*x + 2t*x^2 + x^3 - y^2
findImageEquation(HashTable) := (unfolding) -> (
  -- eliminate source variables from the parametrisation
  param := unfolding#"parametrisation";
  apply(unfolding#"sourceCoords", var -> (
    param = eliminate(var, param)
  ));
  imageIdeal := param;
  imageEquation := first first entries gens imageIdeal;
  return imageEquation
)


partials = method()
    ---- desc: return partial derivatives of a function wrt given variables.
    ---- 
    ---- ii01 : QQ[x,y]; 
    ---- ii02 : partials({x^3 + x*y + y^2}, {x,y})
    ---- oo02 = {3x^2 + y, x + 2y}
partials(List, List) := (funcData, vars) -> (
  func := first funcData;
  apply(vars, var -> diff(var, func))
)

compareSaturation = method()
    ---- find what power of the transporter (I : J^k) is required for the saturation
    ---- 
    ---- ii01 : compareSaturation(exampleA1())
    ---- oo01 = 1
compareSaturation(HashTable) := (unfolding) -> (
  return compareSaturation(unfolding, "default");
)
compareSaturation(HashTable, String) := (unfolding, strategy) -> (
  print concatenate("loaded example: ", toString(unfolding#"parametrisation"), "; method: ", strategy);
  (critPoints, zeroLocus) := prepareIdeals(unfolding, strategy);
  saturation := saturate(critPoints, zeroLocus);
  
  k := 1;
  for k from 1 to 5 do (
    transK := ideal quotient(critPoints : zeroLocus^k);
    print concatenate("(k=", toString(k), ") (I:J^", toString(k), ") = ", toString(transK));
    print concatenate("saturation: ", toString(saturation));

    if(saturation == transK) then break;
  );

  return k;
)


hasRadicalCritPoints = method()
    ---- desc: given an unfolding, test if the critPoints ideal is radical
    ----
    ---- ii01 : hasRadicalCritPoints(exampleA1())
hasRadicalCritPoints(HashTable) := (unfolding) -> (
  (critPoints, zeroLocus) := prepareIdeals(unfolding, "default");
  critPointsRadical := radical critPoints;
  print concatenate("relative Jacobian: ", toString(critPoints) , "\n", "its radical: ", toString(critPointsRadical));
  return (critPoints == critPointsRadical);
)


TEST /// -- countCritPoints
  load("unfolding_examples.m2");
  assert (1 == countCritPoints(exampleA1()));
  assert (4 == countCritPoints(exampleC1(4)));
///

TEST /// -- prepareIdeals
  QQ[s,t,x,y]
  unfolding = new HashTable from {
    "parametrisation" => ideal(s^2 - x, s^3 + s*t - y),
    "sourceCoords" => {s},
    "targetSpace" => {x,y},
    "targetCoords" => {t,x,y}
  };
  critPointsExp = ideal(t^2 + 4*t*x + 3*x^2, -2*y)
  zeroLocusExp = ideal(t^2*x + 2*t*x^2 + x^3 - y^2)
  (critPoints, zeroLocus) = prepareIdeals(unfolding, "default")
  assert (critPoints == critPointsExp)
  assert (zeroLocus  == zeroLocusExp )
///

TEST /// -- partials
  QQ[x,y]
  expect = {3*x^2 + 4*y, 4*x + 2*y}
  assert (expect == partials({x^3 + 4*x*y+y^2}, {x,y}))
///

TEST /// -- findImageEquation
  QQ[s,t,x,y]
  unfolding = new HashTable from {
    "parametrisation" => ideal(s^2 - x, s^3 + s*t - y),
    "sourceCoords" => {s}
  };
  expect = t^2*x + 2*t*x^2 + x^3 - y^2;
  assert (expect == findImageEquation(unfolding)) 
///

TEST /// -- compareSaturation
  load("unfolding_examples.m2")
  assert (1 == compareSaturation(exampleA2()))
///

TEST /// -- hasRadicalCritPoints
  load("unfolding_examples.m2");
  assert (true == hasRadicalCritPoints(exampleA1()))
  assert (false == hasRadicalCritPoints(exampleA2()))
///