source: git/Singular/LIB/gitfan.lib @ dc8eb2

///////////////////////////////////////////////////////////////////
version="version gitfan.lib 4.0.0.0 Jun_2013 "; // $Id$
category="Algebraic Geometry";
info="
LIBRARY:  gitfan.lib       Compute GIT-fans.

AUTHORS:  Janko Boehm      boehm@mathematik.uni-kl.de
@*        Simon Keicher    keicher@mail.mathematik.uni-tuebingen.de
@*        Yue Ren          ren@mathematik.uni-kl.de

OVERVIEW:
This library computes GIT-fans, torus orbits and GKZ-fans.
It uses the package 'gfanlib' by Anders N. Jensen
and some algorithms have been outsourced to C++ to improve the performance.
Check https://github.com/skeicher/gitfan_singular for updates.

KEYWORDS: library; gitfan; GIT; geometric invariant theory; quotients

PROCEDURES:
afaces(ideal);                      Returns a list of intvecs that correspond to all a-faces
gitCone(ideal,bigintmat,bigintmat); Returns the GIT-cone around the given weight vector w
gitFan(ideal,bigintmat);            Returns the GIT-fan of the H-action defined by Q on V(a)
gkzFan(bigintmat);                  Returns the GKZ-fan of the matrix Q
isAface(ideal,intvec);              Checks whether intvec corresponds to an ideal-face
orbitCones(ideal,bigintmat);        Returns the list of all projected a-faces
";

LIB "parallel.lib"; // for parallelWaitAll
LIB "general.lib";

////////////////////////////////////////////////////

proc int2face(int n, int r)
{
  int k = r-1;
  intvec v;
  int n0 = n;

  while(n0 > 0){
    while(2^k > n0){
      k--;
    }

    v = k+1,v;
    n0 = n0 - 2^k;
    k--;
  }
  v = v[1..size(v)-1];
  return(v);
}

/////////////////////////////////

proc isAface(ideal a, intvec gam0, int n)
"USAGE:  isAface(a,gam0); a: ideal, gam0:intvec
PURPOSE: Checks whether the face of the positive orthant,
         that is spanned by all i-th unit vectors,
         where i ranges amongst the entries of gam0,
         is an a-face.
RETURN:  int
EXAMPLE: example isaface; shows an example
"
{
  // special case: gam0 is the zero-cone:
  if (size(gam0) == 1 and gam0[1] == 0){
    poly pz; ideal G;

    // is an a-face if and only if RL0(0,...,0) = const
    // set all entries to 0:
    int i;
    for (int k = 1; k <= size(a); k++) {
      pz = subst(a[k], var(1), 0);
      for (i = 2; i <= nvars(basering); i++) {
        pz = subst(pz, var(i), 0);
      }
      G = G, pz;
    }

    G = std(G);

    // monomial inside?:
    if(1 == G){
      return(0);
    }

    return(1);
  }

  // ring is too big: Switch to KK[T_i | e_i\in gam0] instead:
  intvec w=ringlist(basering)[3][1][2];
  intvec w0;
  string initNewRing = "ring Rgam0 = 0,(";

  for (int i=1; i<size(gam0); i++){
    initNewRing = initNewRing + string(var(gam0[i])) + ",";
    w0[i]=w[gam0[i]];
  }

  w0 = w0,w[gam0[i]],1;
  initNewRing = initNewRing + string(var(gam0[size(gam0)])) + ",U),Wp("+string(w0)+");";
  def R = basering;
  execute(initNewRing);

  ideal agam0 = imap(R,a);

  for (i = 1; i<=size(agam0); i=i+1)
  {
    if (size(agam0[i]) == 1)
    { return(0,0); }
  }

  poly p = var(1); // first entry of g; p = prod T_i with i element of g
  for ( i = 2; i <= nvars(basering); i++ ) {
      p = p * var(i);
    }
  // p is now the product over all T_i, with e_i in gam0

  agam0 = agam0, p - 1; // rad-membership
  agam0 = timeStd(agam0,5);
  // "timestd finished after: "+string(timer-t);
  // int timeout = 0;
  if (attrib(agam0,"isSB") < 1)
  {
    kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
    return(0,1);
    // // "timestd failed in "+string(gam0)+", falling back to saturation!";
    // setring R; kill Rgam0;

    // initNewRing = "ring Rgam0 = 0,(";

    // w0 = 0;
    // for (i=1; i<size(gam0); i++){
    //   initNewRing = initNewRing + string(var(gam0[i])) + ",";
    //   w0[i]=w[gam0[i]];
    // }

    // w0 = w0,w[gam0[i]];
    // initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),Wp("+string(w0)+");";
    // execute(initNewRing);

    // ideal G = imap(R,a);

    // timeout = 1;
    // int t=rtimer;
    // for(int k=nvars(basering); k>0; k--) { G=sat(G,var(k))[1]; }
    // t = (rtimer - t) div 1000;
    // string(n)+": saturation successful after "+string(t)+" with: "+string(G<>1);
  }

  // does G contain 1?, i.e. is G = 1?
  if(agam0 <> 1) {
    kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
    return(1,0); // true
  }

  kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
  return(0,0); // false
}
example
{
  echo = 2;

  ring R = 0,(T(1..4)),dp;
  ideal I = T(1)*T(2)-T(4);

  intvec w = 1,4;
  intvec v = 1,2,4;

  isAface(I,w); // should be 0
  "-----------";
  isAface(I,v); // should be 1
}


proc isAfaceNonZero(ideal a, intvec gam0)
{
  string initNewRing = "ring Rgam0 = 0,(";
  for (int i=1; i<size(gam0); i++)
    { initNewRing = initNewRing + string(var(gam0[i])) + ","; }
  initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),dp;";
  def R = basering;
  execute(initNewRing);

  ideal agam0 = imap(R,a);

  for ( i = 1; i<=size(agam0); i=i+1)
    { if (size(agam0[i]) == 1) { return(0); } }

  poly p = var(1);
  for ( i = 2; i <= nvars(basering); i++ )
    { p = p * var(i); }

  agam0 = agam0, p - 1;
  ideal G = std(agam0);

  if(G <> 1)
    { return(1); }

  return(0);
}
////////////////////////////////////////////////////

proc afacesPart(ideal a, int d, int start, int end, int r){
  intvec gam0;
  int i;
  list AF;

  for(i = start; i <= end; i=i+1){
    if(i < 2^r){
      string(i)+"    "+string(size(AF));
      gam0 = int2face(i,r);

      // take gam0 only if it has
      // at least d rays:
      if(size(gam0) >= d){
        if (isAface(a,gam0)){
          AF[size(AF) + 1] = gam0;
        }
      }
    }
  }

  return(AF);
}

////////////////////////////////////////////////////

proc afaces(ideal a, list #)
"USAGE:  afaces(a, b, c); a: ideal, b: int, c: int
PURPOSE: Returns a list of all a-faces (represented by intvecs).
         Moreover, it is possible to specify a dimensional bound b,
         upon which only cones of that dimension and above are returned.
         Lastly, as the computation is parallizable, one can specify c,
         the number of cores to be used by the computation.
RETURN:  a list of intvecs
EXAMPLE: example afaces; shows an example
"
{
  int d = 1;
  int ncores = 1;

  if ((size(#) > 0) and (typeof(#[1]) == "int")){
    d = #[1];
  }

  if ((size(#) > 1) and (typeof(#[2]) == "int")){
    ncores = #[2];
  }

  list AF;
  intvec gam0;
  int r = nvars(basering);

  // check if 0 is an a-face:
  gam0 = 0;
  if (isAface(a,gam0)){
      AF[size(AF) + 1] = gam0;
  }

  // check for other a-faces:
  // make ncores processes:
  int step = 2^r div ncores;
  int i;

  list args;
  for(int k = 0; k < ncores; k++){
    args[size(args) + 1] = list(a, d, k * step + 1, (k+1) * step, r);
  }

  string command = "afacesPart";
  list out = parallelWaitAll(command, args);

  // do remaining ones:
  for(i = ncores * step +1; i < 2^r; i=i+1){
    "another one needed";
    gam0 = int2face(i,r);

    // take gam0 only if it has
    // at least d rays:
    if(size(gam0) >= d){
      if (isAface(a,gam0)){
        AF[size(AF) + 1] = gam0;
      }
    }
  }

  int l;
  // read out afaces of out into AF:
  for(l = 1; l <= size(out); l++){
    AF = AF + out[l];
  }

  return(AF);
}
example
{

  echo = 2;
  ring R = 0,T(1..3),dp;
  ideal a = T(1)+T(2)+T(3);

  list F = afaces(a,3,4);
  print(F);
  print(size(F));

  // 2nd ex //
  ring R2 = 0,T(1..3),dp;
  ideal a2 = T(2)^2*T(3)^2+T(1)*T(3);

  list F2 = afaces(a2,3,4);
  print(F2);
  print(size(F2));

  // 3rd ex //
  ring R3 = 0,T(1..3),dp;
  ideal a3 = 0;

  list F3 = afaces(a3,3,4);
  print(F3);
  print(size(F3));

  // bigger example //
  ring R = 0,T(1..15),dp;
  ideal a =
    T(1)*T(10)-T(2)*T(7)+T(3)*T(6),
    T(1)*T(11)-T(2)*T(8)+T(4)*T(6),
    T(1)*T(12)-T(2)*T(9)+T(5)*T(6),
    T(1)*T(13)-T(3)*T(8)+T(4)*T(7),
    T(1)*T(14)-T(3)*T(9)+T(5)*T(7),
    T(1)*T(15)-T(4)*T(9)+T(5)*T(8),
    T(2)*T(13)-T(3)*T(11)+T(4)*T(10),
    T(2)*T(14)-T(3)*T(12)+T(5)*T(10),
    T(2)*T(15)-T(4)*T(12)+T(5)*T(11),
    T(3)*T(15)-T(4)*T(14)+T(5)*T(13),
    T(6)*T(13)-T(7)*T(11)+T(8)*T(10),
    T(6)*T(14)-T(7)*T(12)+T(9)*T(10),
    T(6)*T(15)-T(8)*T(12)+T(9)*T(11),
    T(7)*T(15)-T(8)*T(14)+T(9)*T(13),
    T(10)*T(15)-T(11)*T(14)+T(12)*T(13);

  int t = timer;
  list F4 = afaces(a,0,2);
  print(size(F4));
  timer - t;

  int t = timer;
  list F4 = afaces(a,0);
  print(size(F4));
  timer - t;

}

///////////////////////////////////////

proc orbitCones(ideal a, bigintmat Q, list #)
"USAGE:  orbitCones(a, Q, b, c); a: ideal, Q: bigintmat, b: int, c: int
PURPOSE: Returns a list consisting of all cones Q(gam0) where gam0 is an a-face.
         Moreover, it is possible to specify a dimensional bound b,
         upon which only cones of that dimension and above are returned.
         Lastly, as the computation is parallizable, one can specify c,
         the number of cores to be used by the computation.
RETURN:  a list of cones
EXAMPLE: example orbitCones; shows an example
"
{
  list AF;

  if((size(#) > 1) and (typeof(#[2]) == "int")){
    AF = afaces(a, nrows(Q), #[2]);
  } else
  {
    AF = afaces(a, nrows(Q));
  }

  int dimensionBound = 0;
  if((size(#) > 0) and (typeof(#[1]) == "int")){
    dimensionBound = #[1];
  }

  list OC;
  intvec gam0;
  int j;

  for(int i = 1; i <= size(AF); i++){
    gam0 = AF[i];

    if(gam0 == 0){
      bigintmat M[1][nrows(Q)];
    } else {
      bigintmat M[size(gam0)][nrows(Q)];
      for (j = 1; j <= size(gam0); j++){
        M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]];
      }
    }
    cone c = coneViaPoints(M);

    if((dimension(c) >= dimensionBound) and (!(listContainsCone(OC, c)))){
      OC[size(OC)+1] = c;
    }

    kill M, c;
  }

  return(OC);
}
example
{
  echo=2;
  intmat Q[3][4] =
    1,0,1,0,
    0,1,0,1,
    0,0,1,1;

  ring R = 0,T(1..4),dp;
  ideal a = 0;

  orbitCones(a, Q);
}

///////////////////////////////////////

proc gitCone(ideal a, bigintmat Q, bigintmat w)
"USAGE: gitCone(a, Q, w); a: ideal, Q:bigintmat, w:bigintmat
PURPOSE: Returns the GIT-cone lambda(w), i.e. the intersection of all
orbit cones containing the vector w.
NOTE: call this only if you are interested in a single GIT-cone.
RETURN: a cone.
EXAMPLE: example gitCone; shows an example
"
{
  list OC =  orbitCones(a, Q);
  cone lambda = nrows(Q);

  for(int i = 1; i <= size(OC); i++){
    cone c = OC[i];

    if(containsInSupport(c, w)){
      lambda = convexIntersection(lambda, c);
    }

    kill c;
  }

  return(lambda);
}
example
{
  echo=2;
  intmat Q[3][4] =
    1,0,1,0,
    0,1,0,1,
    0,0,1,1;

  ring R = 0,T(1..4),dp;
  ideal a = 0;

  bigintmat w[1][3] = 3,3,1;
  cone lambda = gitCone(a, Q, w);
  rays(lambda);

  bigintmat w2[1][3] = 1,1,1;
  cone lambda2 = gitCone(a, Q, w2);
  rays(lambda2);
}

/////////////////////////////////////

proc gitFan(ideal a, bigintmat Q, list #)
"USAGE: gitFan(a, Q); a: ideal, Q:bigintmat
PURPOSE: Returns the GIT-fan of the H-action defined by Q on V(a).
An optional third parameter of type 'int' is interpreted as the
number of CPU-cores you would like to use.
Note that 'system("cpu");' returns the number of cpu available
in your system.
RETURN: a fan.
EXAMPLE: example gitFan; shows an example
"
{
  list OC = orbitCones(a, Q, #);

  fan f = refineCones(OC, Q);
  return(f);
}
example
{
  echo=2;
  intmat Q[3][4] =
    1,0,1,0,
    0,1,0,1,
    0,0,1,1;

  ring R = 0,T(1..4),dp;
  ideal a = 0;

  gitFan(a, Q);

  // 2nd example //
  kill Q;
  intmat Q[3][6] =
    1,1,0,0,-1,-1,
    0,1,1,-1,-1,0,
    1,1,1,1,1,1;

  ring R = 0,T(1..6),dp;
  ideal a = T(1)*T(6) + T(2)*T(5) + T(3)*T(4);

  int t = rtimer;
  fan F = gitFan(a, Q);
  t = rtimer - t;

  int tt = rtimer;
  fan F = gitFan(a, Q, 4);
  tt = rtimer - tt;

  t;
  tt;
  "--------";
  kill R, Q, t, tt;
  // next example //
  ring R = 0,T(1..10),dp;
  ideal a = T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);

  bigintmat Q[4][10] =
    1,0,0,0,1,1,1,0,0,0,
    0,1,0,0,1,0,0,1,1,0,
    0,0,1,0,0,1,0,1,0,1,
    0,0,0,1,0,0,1,0,1,1;

  int t = rtimer;
  fan F = gitFan(a, Q);
  t = rtimer - t;

  int tt = rtimer;
  fan F = gitFan(a, Q, 4);
  tt = rtimer - tt;

  t;
  tt;

  "--------";
  kill R, Q, t, tt;
  // next example //
  ring R = 0,T(1..15),dp;
  ideal a =
    T(1)*T(10)-T(2)*T(7)+T(3)*T(6),
    T(1)*T(11)-T(2)*T(8)+T(4)*T(6),
    T(1)*T(12)-T(2)*T(9)+T(5)*T(6),
    T(1)*T(13)-T(3)*T(8)+T(4)*T(7),
    T(1)*T(14)-T(3)*T(9)+T(5)*T(7),
    T(1)*T(15)-T(4)*T(9)+T(5)*T(8),
    T(2)*T(13)-T(3)*T(11)+T(4)*T(10),
    T(2)*T(14)-T(3)*T(12)+T(5)*T(10);

  bigintmat Q[5][15] =
    1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,
    0,1,0,0,0,1,0,0,0,1,1,1,0,0,0,
    0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,
    0,0,0,1,0,0,0,1,0,0,1,0,1,0,1,
    0,0,0,0,1,0,0,0,1,0,0,1,0,1,1;

  int t = rtimer;
  fan F = gitFan(a, Q);
  t = rtimer - t;

  int tt = rtimer;
  fan F = gitFan(a, Q, 4);
  tt = rtimer - tt;

  t;
  tt;

}

/////////////////////////////////////
// Computes all simplicial orbit cones
// w.r.t. the 0-ideal:

static proc simplicialToricOrbitCones(bigintmat Q){
  intvec gam0;
  list OC;
  cone c;
  int r = ncols(Q);
  int j;

  for(int i = 1; i < 2^r; i++ ){
    gam0 = int2face(i,r);

    // each simplicial cone is generated by
    // exactly nrows(Q) many columns of Q:
    if(size(gam0) == nrows(Q)){
      bigintmat M[size(gam0)][nrows(Q)];

      for(j = 1; j <= size(gam0); j++){
        M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]];
      }

      c = coneViaPoints(M);

      if((dimension(c) == nrows(Q)) and (!(listContainsCone(OC, c)))){
        OC[size(OC)+1] = c;
      }

      kill M;
    }
  }

  return(OC);
}

/////////////////////////////////////

proc gkzFan(bigintmat Q)
"USAGE: gkzFan(Q); a: ideal, Q:bigintmat
PURPOSE: Returns the GKZ-fan of the matrix Q.
RETURN: a fan.
EXAMPLE: example gkzFan; shows an example
"
{
  // only difference to gitFan:
  // it suffices to consider all faces
  // that are simplicial:
  list OC = simplicialToricOrbitCones(Q);

  fan f = refineCones(OC, Q);
  return(f);
}
example
{
  echo=2;
  intmat Q[3][4] =
    1,0,1,0,
    0,1,0,1,
    0,0,1,1;

  gkzFan(Q);
}