Changeset a1ef3a2 in git for kernel/combinatorics

Timestamp:

Mar 27, 2021, 7:21:22 PM (3 years ago)

Author:

Karim Abou Zeid <karim23697@…>

Branches:

(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'c5facdfddea2addfd91babd8b9019161dea4b695')

Children:

80dee2b8505d045ce3d5bf4fa6bee7b4383af86c

Parents:

2daaaec2c36328eb8638fe3a6b61d9971c661dc4

git-author:

Karim Abou Zeid <karim23697@gmail.com>2021-03-27 19:21:22+01:00

git-committer:

Hans Schoenemann <hannes@mathematik.uni-kl.de>2021-03-30 00:45:51+02:00

Message:

implement kdim in kernel as vdim

Location:

kernel/combinatorics

Files:

• hdegree.cc (modified) (8 diffs)
• stairc.h (modified) (1 diff)

kernel/combinatorics/hdegree.cc

@@ -19 +19 @@
 #include "kernel/combinatorics/stairc.h"
 #include "reporter/reporter.h"
+
+#ifdef HAVE_SHIFTBBA
+#include <vector>
+#include "Singular/libsingular.h"
+#endif
 
 VAR int  hCo, hMu, hMu2;
@@ -1668 +1673 @@
 }
 
-static void _lp_computeStandardWords(ideal words, int n, ideal M, int& last)
-{
-  // assume <M> != <1>
-  if (n <= 0){
-    words->m[0] = pOne();
-    last = 0;
+// ATTENTION:
+//  - `words` contains the words normal modulo M of length n
+//  - `numberOfNormalWords` contains the number of words normal modulo M of length 0 ... n
+static void _lp_computeNormalWords(ideal words, int& numberOfNormalWords, int length, ideal M, int minDeg, int& last)
+{
+  if (length <= 0){
+    poly one = pOne();
+    if (p_LPDivisibleBy(M, one, currRing)) // 1 \in M => no normal words at all
+    {
+      pDelete(&one);
+      last = -1;
+      numberOfNormalWords = 0;
+    }
+    else
+    {
+      words->m[0] = one;
+      last = 0;
+      numberOfNormalWords = 1;
+    }
     return;
   }
 
-  _lp_computeStandardWords(words, n - 1, M, last);
+  _lp_computeNormalWords(words, numberOfNormalWords, length - 1, M, minDeg, last);
 
   int nVars = currRing->isLPring - currRing->LPncGenCount;
+  int numberOfNewNormalWords = 0;
 
   for (int j = nVars - 1; j >= 0; j--)
@@ -1693 +1712 @@
         }
 
-        int varOffset = ((n - 1) * currRing->isLPring) + 1;
+        int varOffset = ((length - 1) * currRing->isLPring) + 1;
         pSetExp(words->m[index], varOffset + j, 1);
         pSetm(words->m[index]);
         pTest(words->m[index]);
 
-        if (p_LPDivisibleBy(M, words->m[index], currRing))
+        if (length >= minDeg && p_LPDivisibleBy(M, words->m[index], currRing))
         {
           pDelete(&words->m[index]);
           words->m[index] = NULL;
         }
+        else
+        {
+          numberOfNewNormalWords++;
+        }
       }
     }
@@ -1708 +1731 @@
 
   last = nVars * last + nVars - 1;
-}
-
-static ideal lp_computeStandardWords(int n, ideal M)
-{
+
+  numberOfNormalWords += numberOfNewNormalWords;
+}
+
+static ideal lp_computeNormalWords(int length, ideal M)
+{
+  long minDeg = IDELEMS(M) > 0 ? pTotaldegree(M->m[0]) : 0;
+  for (int i = 1; i < IDELEMS(M); i++)
+  {
+    minDeg = si_min(minDeg, pTotaldegree(M->m[i]));
+  }
+
   int nVars = currRing->isLPring - currRing->LPncGenCount;
 
   int maxElems = 1;
-  for (int i = 0; i < n; i++) // maxElems = nVars^n
+  for (int i = 0; i < length; i++) // maxElems = nVars^n
     maxElems *= nVars;
   ideal words = idInit(maxElems);
-  int last;
-  _lp_computeStandardWords(words, n, M, last);
+  int last, numberOfNormalWords;
+  _lp_computeNormalWords(words, numberOfNormalWords, length, M, minDeg, last);
   idSkipZeroes(words);
   return words;
+}
+
+static int lp_countNormalWords(int upToLength, ideal M)
+{
+  long minDeg = IDELEMS(M) > 0 ? pTotaldegree(M->m[0]) : 0;
+  for (int i = 1; i < IDELEMS(M); i++)
+  {
+    minDeg = si_min(minDeg, pTotaldegree(M->m[i]));
+  }
+
+  int nVars = currRing->isLPring - currRing->LPncGenCount;
+
+  int maxElems = 1;
+  for (int i = 0; i < upToLength; i++) // maxElems = nVars^n
+    maxElems *= nVars;
+  ideal words = idInit(maxElems);
+  int last, numberOfNormalWords;
+  _lp_computeNormalWords(words, numberOfNormalWords, upToLength, M, minDeg, last);
+  idDelete(&words);
+  return numberOfNormalWords;
 }
 
@@ -1738 +1789 @@
   int lV = currRing->isLPring;
 
-  standardWords = lp_computeStandardWords(l, G);
+  standardWords = lp_computeNormalWords(l, G);
 
   int n = IDELEMS(standardWords);
@@ -1823 +1874 @@
     int ncGenCount = currRing->LPncGenCount;
     if (lV - ncGenCount == 0)
+    {
+      idDelete(&G);
       return 0;
+    }
     if (lV - ncGenCount == 1)
+    {
+      idDelete(&G);
       return 1;
+    }
     if (lV - ncGenCount >= 2)
+    {
+      idDelete(&G);
       return -1;
+    }
   }
 
@@ -1840 +1900 @@
     {
       WerrorS("GK-Dim not defined for 0-ring");
+      idDelete(&G);
       return -2;
     }
@@ -1850 +1911 @@
     int ncGenCount = currRing->LPncGenCount;
     if (IDELEMS(G) == lV - ncGenCount) // V = {1} no edges
+    {
+      idDelete(&G);
       return 0;
+    }
     if (IDELEMS(G) == lV - ncGenCount - 1) // V = {1} with loop
+    {
+      idDelete(&G);
       return 1;
+    }
     if (IDELEMS(G) <= lV - ncGenCount - 2) // V = {1} with more than one loop
+    {
+      idDelete(&G);
       return -1;
+    }
   }
 
   ideal standardWords;
   intvec* UG = lp_ufnarovskiGraph(G, standardWords);
-  if (errorreported || UG == NULL) return -2;
-  return graphGrowth(UG);
+  if (UG == NULL)
+  {
+    idDelete(&G);
+    return -2;
+  }
+  if (errorreported)
+  {
+    delete UG;
+    idDelete(&G);
+    return -2;
+  }
+  int gkDim = graphGrowth(UG);
+  delete UG;
+  idDelete(&G);
+  return gkDim;
+}
+
+// converts an intvec matrix to a vector<vector<int> >
+static std::vector<std::vector<int> > iv2vv(intvec* M)
+{
+  int rows = M->rows();
+  int cols = M->cols();
+
+  std::vector<std::vector<int> > mat(rows, std::vector<int>(cols));
+
+  for (int i = 0; i < rows; i++)
+  {
+    for (int j = 0; j < cols; j++)
+    {
+      mat[i][j] = IMATELEM(*M, i + 1, j + 1);
+    }
+  }
+
+  return mat;
+}
+
+static void vvPrint(const std::vector<std::vector<int> >& mat)
+{
+  for (int i = 0; i < mat.size(); i++)
+  {
+    for (int j = 0; j < mat[i].size(); j++)
+    {
+      Print("%d ", mat[i][j]);
+    }
+    PrintLn();
+  }
+}
+
+static void vvTest(const std::vector<std::vector<int> >& mat)
+{
+  if (mat.size() > 0)
+  {
+    int cols = mat[0].size();
+    for (int i = 1; i < mat.size(); i++)
+    {
+      if (cols != mat[i].size())
+        WerrorS("number of cols in matrix inconsistent");
+    }
+  }
+}
+
+static void vvDeleteRow(std::vector<std::vector<int> >& mat, int row)
+{
+  mat.erase(mat.begin() + row);
+}
+
+static void vvDeleteColumn(std::vector<std::vector<int> >& mat, int col)
+{
+  for (int i = 0; i < mat.size(); i++)
+  {
+    mat[i].erase(mat[i].begin() + col);
+  }
+}
+
+static BOOLEAN vvIsRowZero(const std::vector<std::vector<int> >& mat, int row)
+{
+  for (int i = 0; i < mat[row].size(); i++)
+  {
+    if (mat[row][i] != 0)
+      return FALSE;
+  }
+  return TRUE;
+}
+
+static BOOLEAN vvIsColumnZero(const std::vector<std::vector<int> >& mat, int col)
+{
+  for (int i = 0; i < mat.size(); i++)
+  {
+    if (mat[i][col] != 0)
+      return FALSE;
+  }
+  return TRUE;
+}
+
+static BOOLEAN vvIsZero(const std::vector<std::vector<int> >& mat)
+{
+  for (int i = 0; i < mat.size(); i++)
+  {
+    if (!vvIsRowZero(mat, i))
+      return FALSE;
+  }
+  return TRUE;
+}
+
+static std::vector<std::vector<int> > vvMult(const std::vector<std::vector<int> >& a, const std::vector<std::vector<int> >& b)
+{
+  int ra = a.size();
+  int rb = b.size();
+  int ca = a.size() > 0 ? a[0].size() : 0;
+  int cb = b.size() > 0 ? b[0].size() : 0;
+
+  if (ca != rb)
+  {
+    WerrorS("matrix dimensions do not match");
+    return std::vector<std::vector<int> >();
+  }
+
+  std::vector<std::vector<int> > res(ra, std::vector<int>(cb));
+  for (int i = 0; i < ra; i++)
+  {
+    for (int j = 0; j < cb; j++)
+    {
+      int sum = 0;
+      for (int k = 0; k < ca; k++)
+        sum += a[i][k] * b[k][j];
+      res[i][j] = sum;
+    }
+  }
+  return res;
+}
+
+static BOOLEAN isAcyclic(const intvec* G)
+{
+  // init
+  int n = G->cols();
+  std::vector<int> path;
+  std::vector<BOOLEAN> visited;
+  std::vector<BOOLEAN> cyclic;
+  std::vector<int> cache;
+  visited.resize(n, FALSE);
+  cyclic.resize(n, FALSE);
+  cache.resize(n, -2);
+
+  for (int v = 0; v < n; v++)
+  {
+    cache = countCycles(G, v, path, visited, cyclic, cache);
+    // check that there are 0 cycles from v
+    if (cache[v] != 0)
+      return FALSE;
+  }
+  return TRUE;
+}
+
+/*
+ * Computation of the K-Dimension
+ */
+
+// -1 is infinity, -2 is error
+int lp_kDim(const ideal _G)
+{
+  if (rField_is_Ring(currRing)) {
+      WerrorS("K-Dim not implemented for rings");
+      return -2;
+  }
+
+  for (int i=IDELEMS(_G)-1;i>=0; i--)
+  {
+    if (_G->m[i] != NULL)
+    {
+      if (pGetComp(_G->m[i]) != 0)
+      {
+        WerrorS("K-Dim not implemented for modules");
+        return -2;
+      }
+      if (pGetNCGen(_G->m[i]) != 0)
+      {
+        WerrorS("K-Dim not implemented for bi-modules");
+        return -2;
+      }
+    }
+  }
+
+  ideal G = id_Head(_G, currRing); // G = LM(G) (and copy)
+  if (TEST_OPT_PROT)
+    Print("%d original generators\n", IDELEMS(G));
+  idSkipZeroes(G); // remove zeros
+  id_DelLmEquals(G, currRing); // remove duplicates
+  if (TEST_OPT_PROT)
+    Print("%d non-zero unique generators\n", IDELEMS(G));
+
+  // check if G is the zero ideal
+  if (IDELEMS(G) == 1 && G->m[0] == NULL)
+  {
+    // NOTE: this is needed because if the ideal is <0>, then idSkipZeroes keeps this element, and IDELEMS is still 1!
+    int lV = currRing->isLPring;
+    int ncGenCount = currRing->LPncGenCount;
+    if (lV - ncGenCount == 0)
+    {
+      idDelete(&G);
+      return 1;
+    }
+    if (lV - ncGenCount == 1)
+    {
+      idDelete(&G);
+      return -1;
+    }
+    if (lV - ncGenCount >= 2)
+    {
+      idDelete(&G);
+      return -1;
+    }
+  }
+
+  // get the max deg
+  long maxDeg = 0;
+  for (int i = 0; i < IDELEMS(G); i++)
+  {
+    maxDeg = si_max(maxDeg, pTotaldegree(G->m[i]));
+
+    // also check whether G = <1>
+    if (pIsConstantComp(G->m[i]))
+    {
+      WerrorS("K-Dim not defined for 0-ring"); // TODO is it minus infinity ?
+      idDelete(&G);
+      return -2;
+    }
+  }
+  if (TEST_OPT_PROT)
+    Print("max deg: %ld\n", maxDeg);
+
+
+  // for normal words of length minDeg ... maxDeg-1
+  // brute-force the normal words
+  if (TEST_OPT_PROT)
+    PrintS("Computing normal words normally...\n");
+  long numberOfNormalWords = lp_countNormalWords(maxDeg - 1, G);
+
+  if (TEST_OPT_PROT)
+    Print("%ld normal words up to length %ld\n", numberOfNormalWords, maxDeg - 1);
+
+  // early termination if G \subset X
+  if (maxDeg <= 1)
+  {
+    int lV = currRing->isLPring;
+    int ncGenCount = currRing->LPncGenCount;
+    if (IDELEMS(G) == lV - ncGenCount) // V = {1} no edges
+    {
+      idDelete(&G);
+      return numberOfNormalWords;
+    }
+    if (IDELEMS(G) == lV - ncGenCount - 1) // V = {1} with loop
+    {
+      idDelete(&G);
+      return -1;
+    }
+    if (IDELEMS(G) <= lV - ncGenCount - 2) // V = {1} with more than one loop
+    {
+      idDelete(&G);
+      return -1;
+    }
+  }
+
+  if (TEST_OPT_PROT)
+    PrintS("Computing Ufnarovski graph...\n");
+
+  ideal standardWords;
+  intvec* UG = lp_ufnarovskiGraph(G, standardWords);
+  if (UG == NULL)
+  {
+    idDelete(&G);
+    return -2;
+  }
+  if (errorreported)
+  {
+    delete UG;
+    idDelete(&G);
+    return -2;
+  }
+
+  if (TEST_OPT_PROT)
+    Print("Ufnarovski graph is %dx%d.\n", UG->rows(), UG->cols());
+
+  if (TEST_OPT_PROT)
+    PrintS("Checking whether Ufnarovski graph is acyclic...\n");
+
+  if (!isAcyclic(UG))
+  {
+    // in this case we have infinitely many normal words
+    return -1;
+  }
+
+  std::vector<std::vector<int> > vvUG = iv2vv(UG);
+  for (int i = 0; i < vvUG.size(); i++)
+  {
+    if (vvIsRowZero(vvUG, i) && vvIsColumnZero(vvUG, i)) // i is isolated vertex
+    {
+      vvDeleteRow(vvUG, i);
+      vvDeleteColumn(vvUG, i);
+      i--;
+    }
+  }
+  if (TEST_OPT_PROT)
+    Print("Simplified Ufnarovski graph to %ldx%ld.\n", vvUG.size(), vvUG.size());
+
+  // for normal words of length >= maxDeg
+  // use Ufnarovski graph
+  if (TEST_OPT_PROT)
+    PrintS("Computing normal words via Ufnarovski graph...\n");
+  std::vector<std::vector<int> > UGpower = vvUG;
+  long nUGpower = 1;
+  while (!vvIsZero(UGpower))
+  {
+    if (TEST_OPT_PROT)
+      PrintS("Start count graph entries.\n");
+    for (int i = 0; i < UGpower.size(); i++)
+    {
+      for (int j = 0; j < UGpower[i].size(); j++)
+      {
+        numberOfNormalWords += UGpower[i][j];
+      }
+    }
+
+    if (TEST_OPT_PROT)
+    {
+      PrintS("Done count graph entries.\n");
+      Print("%ld normal words up to length %ld\n", numberOfNormalWords, maxDeg - 1 + nUGpower);
+    }
+
+    if (TEST_OPT_PROT)
+      PrintS("Start mat mult.\n");
+    UGpower = vvMult(UGpower, vvUG); // TODO: avoid creation of new intvec
+    if (TEST_OPT_PROT)
+      PrintS("Done mat mult.\n");
+    nUGpower++;
+  }
+
+  delete UG;
+  idDelete(&G);
+  return numberOfNormalWords;
 }
 #endif

kernel/combinatorics/stairc.h

@@ -34 +34 @@
 #if HAVE_SHIFTBBA
 int lp_gkDim(const ideal G);
+int lp_kDim(const ideal G);
 intvec * lp_ufnarovskiGraph(ideal G, ideal &standardWords);
 #endif