source: git/Singular/LIB/gitfan.lib @ 05a147

jengelh-datetimespielwiese
Last change on this file since 05a147 was 05a147, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: version string (MUST be $Id$ or $Id: <name> <ver> <date> ..)
  • Property mode set to 100644
File size: 12.0 KB
Line 
1////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Algebraic Geometry";
4info="
5LIBRARY:  gitfan.lib       Compute GIT-fans.
6
7AUTHORS:  Janko Boehm      boehm@mathematik.uni-kl.de
8@*        Simon Keicher    keicher@mail.mathematik.uni-tuebingen.de
9@*        Yue Ren          ren@mathematik.uni-kl.de
10
11OVERVIEW:
12This library computes GIT-fans, torus orbits and GKZ-fans.
13It uses the package 'gfanlib' by Anders N. Jensen
14and some algorithms have been outsourced to C++ to improve the performance.
15Check https://github.com/skeicher/gitfan_singular for updates.
16
17KEYWORDS: library; gitfan; GIT; geometric invariant theory; quotients
18
19PROCEDURES:
20afaces(ideal);                      Returns a list of intvecs that correspond to all a-faces
21gitCone(ideal,bigintmat,bigintmat); Returns the GIT-cone around the given weight vector w
22gitFan(ideal,bigintmat);            Returns the GIT-fan of the H-action defined by Q on V(a)
23gkzFan(bigintmat);                  Returns the GKZ-fan of the matrix Q
24isAface(ideal,intvec);              Checks whether intvec corresponds to an ideal-face
25orbitCones(ideal,bigintmat);        Returns the list of all projected a-faces
26";
27
28LIB "parallel.lib"; // for parallelWaitAll
29
30////////////////////////////////////////////////////
31
32static proc int2face(int n, int r)
33{
34  int k = r-1;
35  intvec v;
36  int n0 = n;
37
38  while(n0 > 0){
39    while(2^k > n0){
40      k--;
41      //v[size(v)+1] = 0;
42    }
43
44    v = k+1,v;
45    n0 = n0 - 2^k;
46    k--;
47  }
48  v = v[1..size(v)-1];
49  return(v);
50}
51
52/////////////////////////////////
53
54proc isAface(ideal a, intvec gam0)
55"USAGE:  isAface(a,gam0); a: ideal, gam0:intvec
56PURPOSE: Checks whether the face of the positive orthant,
57         that is spanned by all i-th unit vectors,
58         where i ranges amongst the entries of gam0,
59         is an a-face.
60RETURN:  int
61EXAMPLE: example isaface; shows an example
62"
63{
64  poly pz;
65
66  // special case: gam0 is the zero-cone:
67  if (size(gam0) == 1 and gam0[1] == 0){
68    ideal G;
69
70    // is an a-face if and only if RL0(0,...,0) = const
71    // set all entries to 0:
72    int i;
73    for (int k = 1; k <= size(a); k++) {
74      pz = subst(a[k], var(1), 0);
75      for (i = 2; i <= nvars(basering); i++) {
76        pz = subst(pz, var(i), 0);
77      }
78      G = G, pz;
79    }
80
81    G = std(G);
82
83    // monomial inside?:
84    if(1 == G){
85      return(0);
86    }
87
88    return(1);
89  }
90
91
92  // ring is too big: Switch to KK[T_i; e_i\in gam0] instead:
93  string initNewRing = "ring Rgam0 = 0,(";
94
95  for (int i=1; i<size(gam0); i++){
96    initNewRing = initNewRing + string(var(gam0[i])) + ",";
97  }
98
99  initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),dp;";
100  def R = basering;
101  execute(initNewRing);
102
103  ideal agam0 = imap(R,a);
104
105  poly p = var(1); // first entry of g; p = prod T_i with i element of g
106  for ( i = 2; i <= nvars(basering); i++ ) {
107      p = p * var(i);
108    }
109  // p is now the product over all T_i, with e_i in gam0
110
111  agam0 = agam0, p - 1; // rad-membership
112  ideal G = std(agam0);
113
114  // does G contain 1?, i.e. is G = 1?
115  if(G <> 1) {
116    return(1); // true
117  }
118
119  return(0); // false
120}
121example
122{
123  echo = 2;
124
125  ring R = 0,(T(1..4)),dp;
126  ideal I = T(1)*T(2)-T(4);
127
128  intvec w = 1,4;
129  intvec v = 1,2,4;
130
131  isAface(I,w); // should be 0
132  "-----------";
133  isAface(I,v); // should be 1
134}
135
136////////////////////////////////////////////////////
137
138proc afacesPart(ideal a, int d, int start, int end, int r){
139  intvec gam0;
140  int i;
141  list AF;
142
143  for(i = start; i <= end; i++){
144    if(i < 2^r){
145      gam0 = int2face(i,r);
146
147      // take gam0 only if it has
148      // at least d rays:
149      if(size(gam0) >= d){
150        if (isAface(a,gam0)){
151          AF[size(AF) + 1] = gam0;
152        }
153      }
154    }
155  }
156
157  return(AF);
158}
159
160////////////////////////////////////////////////////
161
162proc afaces(ideal a, list #)
163"USAGE:  afaces(a, b, c); a: ideal, d: int, c: int
164PURPOSE: Returns a list of all a-faces (represented by intvecs).
165         Moreover, it is possible to specify a dimensional bound b,
166         upon which only cones of that dimension and above are returned.
167         Lastly, as the computation is parallizable, one can specify c,
168         the number of cores to be used by the computation.
169RETURN:  a list of intvecs
170EXAMPLE: example afaces; shows an example
171"
172{
173  int d = 1;
174  int ncores = 1;
175
176  if ((size(#) > 0) and (typeof(#[1]) == "int")){
177    d = #[1];
178  }
179
180  if ((size(#) > 1) and (typeof(#[2]) == "int")){
181    ncores = #[2];
182  }
183
184  list AF;
185  intvec gam0;
186  int r = nvars(basering);
187
188  // check if 0 is an a-face:
189  gam0 = 0;
190  if (isAface(a,gam0)){
191      AF[size(AF) + 1] = gam0;
192  }
193
194  // check for other a-faces:
195  // make ncores processes:
196  int step = 2^r div ncores;
197  int i;
198
199  list args;
200  for(int k = 0; k < ncores; k++){
201    args[size(args) + 1] = list(a, d, k * step + 1, (k+1) * step, r);
202  }
203
204  string command = "afacesPart";
205  list out = parallelWaitAll(command, args);
206
207  // do remaining ones:
208  for(i = ncores * step +1; i < 2^r; i++){
209    "another one needed";
210    gam0 = int2face(i,r);
211
212    // take gam0 only if it has
213    // at least d rays:
214    if(size(gam0) >= d){
215      if (isAface(a,gam0)){
216        AF[size(AF) + 1] = gam0;
217      }
218    }
219  }
220
221  // read out afaces of out into AF:
222  for(i = 1; i <= size(out); i++){
223    AF = AF + out[i];
224  }
225
226  return(AF);
227}
228example
229{
230
231  echo = 2;
232  ring R = 0,T(1..3),dp;
233  ideal a = T(1)+T(2)+T(3);
234
235  list F = afaces(a,3,4);
236  print(F);
237  print(size(F));
238
239  // 2nd ex //
240  ring R2 = 0,T(1..3),dp;
241  ideal a2 = T(2)^2*T(3)^2+T(1)*T(3);
242
243  list F2 = afaces(a2,3,4);
244  print(F2);
245  print(size(F2));
246
247  // 3rd ex //
248  ring R3 = 0,T(1..3),dp;
249  ideal a3 = 0;
250
251  list F3 = afaces(a3,3,4);
252  print(F3);
253  print(size(F3));
254
255  // bigger example //
256  ring R = 0,T(1..15),dp;
257  ideal a =
258    T(1)*T(10)-T(2)*T(7)+T(3)*T(6),
259    T(1)*T(11)-T(2)*T(8)+T(4)*T(6),
260    T(1)*T(12)-T(2)*T(9)+T(5)*T(6),
261    T(1)*T(13)-T(3)*T(8)+T(4)*T(7),
262    T(1)*T(14)-T(3)*T(9)+T(5)*T(7),
263    T(1)*T(15)-T(4)*T(9)+T(5)*T(8),
264    T(2)*T(13)-T(3)*T(11)+T(4)*T(10),
265    T(2)*T(14)-T(3)*T(12)+T(5)*T(10),
266    T(2)*T(15)-T(4)*T(12)+T(5)*T(11),
267    T(3)*T(15)-T(4)*T(14)+T(5)*T(13),
268    T(6)*T(13)-T(7)*T(11)+T(8)*T(10),
269    T(6)*T(14)-T(7)*T(12)+T(9)*T(10),
270    T(6)*T(15)-T(8)*T(12)+T(9)*T(11),
271    T(7)*T(15)-T(8)*T(14)+T(9)*T(13),
272    T(10)*T(15)-T(11)*T(14)+T(12)*T(13);
273
274  int t = timer;
275  list F4 = afaces(a,0,2);
276  print(size(F4));
277  timer - t;
278
279  int t = timer;
280  list F4 = afaces(a,0);
281  print(size(F4));
282  timer - t;
283
284}
285
286///////////////////////////////////////
287
288proc orbitCones(ideal a, bigintmat Q, list #)
289"USAGE:  orbitCones(a, Q, b, c); a: ideal, Q: bigintmat, b: int, c: int
290PURPOSE: Returns a list consisting of all cones Q(gam0) where gam0 is an a-face.
291         Moreover, it is possible to specify a dimensional bound b,
292         upon which only cones of that dimension and above are returned.
293         Lastly, as the computation is parallizable, one can specify c,
294         the number of cores to be used by the computation.
295RETURN:  a list of cones
296EXAMPLE: example orbitCones; shows an example
297"
298{
299  list AF;
300
301  if((size(#) > 1) and (typeof(#[2]) == "int")){
302    AF = afaces(a, nrows(Q), #[2]);
303  } else
304  {
305    AF = afaces(a, nrows(Q));
306  }
307
308  int dimensionBound = 0;
309  if((size(#) > 0) and (typeof(#[1]) == "int")){
310    dimensionBound = #[1];
311  }
312
313  list OC;
314  intvec gam0;
315  int j;
316
317  for(int i = 1; i <= size(AF); i++){
318    gam0 = AF[i];
319
320    if(gam0 == 0){
321      bigintmat M[1][nrows(Q)];
322    } else {
323      bigintmat M[size(gam0)][nrows(Q)];
324      for (j = 1; j <= size(gam0); j++){
325        M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]];
326      }
327    }
328    cone c = coneViaPoints(M);
329
330    if((dimension(c) >= dimensionBound) and (!(listContainsCone(OC, c)))){
331      OC[size(OC)+1] = c;
332    }
333
334    kill M, c;
335  }
336
337  return(OC);
338}
339example
340{
341  echo=2;
342  intmat Q[3][4] =
343    1,0,1,0,
344    0,1,0,1,
345    0,0,1,1;
346
347  ring R = 0,T(1..4),dp;
348  ideal a = 0;
349
350  orbitCones(a, Q);
351}
352
353///////////////////////////////////////
354
355proc gitCone(ideal a, bigintmat Q, bigintmat w)
356"USAGE: gitCone(a, Q, w); a: ideal, Q:bigintmat, w:bigintmat
357PURPOSE: Returns the GIT-cone lambda(w), i.e. the intersection of all
358orbit cones containing the vector w.
359NOTE: call this only if you are interested in a single GIT-cone.
360RETURN: a cone.
361EXAMPLE: example gitCone; shows an example
362"
363{
364  list OC =  orbitCones(a, Q);
365  cone lambda = nrows(Q);
366
367  for(int i = 1; i <= size(OC); i++){
368    cone c = OC[i];
369
370    if(containsInSupport(c, w)){
371      lambda = convexIntersection(lambda, c);
372    }
373
374    kill c;
375  }
376
377  return(lambda);
378}
379example
380{
381  echo=2;
382  intmat Q[3][4] =
383    1,0,1,0,
384    0,1,0,1,
385    0,0,1,1;
386
387  ring R = 0,T(1..4),dp;
388  ideal a = 0;
389
390  bigintmat w[1][3] = 3,3,1;
391  cone lambda = gitCone(a, Q, w);
392  rays(lambda);
393
394  bigintmat w2[1][3] = 1,1,1;
395  cone lambda2 = gitCone(a, Q, w2);
396  rays(lambda2);
397}
398
399/////////////////////////////////////
400
401proc gitFan(ideal a, bigintmat Q, list #)
402"USAGE: gitFan(a, Q); a: ideal, Q:bigintmat
403PURPOSE: Returns the GIT-fan of the H-action defined by Q on V(a).
404An optional third parameter of type 'int' is interpreted as the
405number of CPU-cores you would like to use.
406Note that 'system("cpu");' returns the number of cpu available
407in your system.
408RETURN: a fan.
409EXAMPLE: example gitFan; shows an example
410"
411{
412  list OC = orbitCones(a, Q, #);
413
414  fan f = refineCones(OC, Q);
415  return(f);
416}
417example
418{
419  echo=2;
420  intmat Q[3][4] =
421    1,0,1,0,
422    0,1,0,1,
423    0,0,1,1;
424
425  ring R = 0,T(1..4),dp;
426  ideal a = 0;
427
428  gitFan(a, Q);
429
430  // 2nd example //
431  kill Q;
432  intmat Q[3][6] =
433    1,1,0,0,-1,-1,
434    0,1,1,-1,-1,0,
435    1,1,1,1,1,1;
436
437  ring R = 0,T(1..6),dp;
438  ideal a = T(1)*T(6) + T(2)*T(5) + T(3)*T(4);
439
440  int t = rtimer;
441  fan F = gitFan(a, Q);
442  t = rtimer - t;
443
444  int tt = rtimer;
445  fan F = gitFan(a, Q, 4);
446  tt = rtimer - tt;
447
448  t;
449  tt;
450  "--------";
451  kill R, Q, t, tt;
452  // next example //
453  ring R = 0,T(1..10),dp;
454  ideal a = T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
455    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
456    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
457    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
458    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
459
460  bigintmat Q[4][10] =
461    1,0,0,0,1,1,1,0,0,0,
462    0,1,0,0,1,0,0,1,1,0,
463    0,0,1,0,0,1,0,1,0,1,
464    0,0,0,1,0,0,1,0,1,1;
465
466  int t = rtimer;
467  fan F = gitFan(a, Q);
468  t = rtimer - t;
469
470  int tt = rtimer;
471  fan F = gitFan(a, Q, 4);
472  tt = rtimer - tt;
473
474  t;
475  tt;
476
477  "--------";
478  kill R, Q, t, tt;
479  // next example //
480  ring R = 0,T(1..15),dp;
481  ideal a =
482    T(1)*T(10)-T(2)*T(7)+T(3)*T(6),
483    T(1)*T(11)-T(2)*T(8)+T(4)*T(6),
484    T(1)*T(12)-T(2)*T(9)+T(5)*T(6),
485    T(1)*T(13)-T(3)*T(8)+T(4)*T(7),
486    T(1)*T(14)-T(3)*T(9)+T(5)*T(7),
487    T(1)*T(15)-T(4)*T(9)+T(5)*T(8),
488    T(2)*T(13)-T(3)*T(11)+T(4)*T(10),
489    T(2)*T(14)-T(3)*T(12)+T(5)*T(10);
490
491  bigintmat Q[5][15] =
492    1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,
493    0,1,0,0,0,1,0,0,0,1,1,1,0,0,0,
494    0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,
495    0,0,0,1,0,0,0,1,0,0,1,0,1,0,1,
496    0,0,0,0,1,0,0,0,1,0,0,1,0,1,1;
497
498  int t = rtimer;
499  fan F = gitFan(a, Q);
500  t = rtimer - t;
501
502  int tt = rtimer;
503  fan F = gitFan(a, Q, 4);
504  tt = rtimer - tt;
505
506  t;
507  tt;
508
509}
510
511/////////////////////////////////////
512// Computes all simplicial orbit cones
513// w.r.t. the 0-ideal:
514
515static proc simplicialToricOrbitCones(bigintmat Q){
516  intvec gam0;
517  list OC;
518  cone c;
519  int r = ncols(Q);
520  int j;
521
522  for(int i = 1; i < 2^r; i++ ){
523    gam0 = int2face(i,r);
524
525    // each simplicial cone is generated by
526    // exactly nrows(Q) many columns of Q:
527    if(size(gam0) == nrows(Q)){
528      bigintmat M[size(gam0)][nrows(Q)];
529
530      for(j = 1; j <= size(gam0); j++){
531        M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]];
532      }
533
534      c = coneViaPoints(M);
535
536      if((dimension(c) == nrows(Q)) and (!(listContainsCone(OC, c)))){
537        OC[size(OC)+1] = c;
538      }
539
540      kill M;
541    }
542  }
543
544  return(OC);
545}
546
547/////////////////////////////////////
548
549proc gkzFan(bigintmat Q)
550"USAGE: gkzFan(Q); a: ideal, Q:bigintmat
551PURPOSE: Returns the GKZ-fan of the matrix Q.
552RETURN: a fan.
553EXAMPLE: example gkzFan; shows an example
554"
555{
556  // only difference to gitFan:
557  // it suffices to consider all faces
558  // that are simplicial:
559  list OC = simplicialToricOrbitCones(Q);
560
561  fan f = refineCones(OC, Q);
562  return(f);
563}
564example
565{
566  echo=2;
567  intmat Q[3][4] =
568    1,0,1,0,
569    0,1,0,1,
570    0,0,1,1;
571
572  gkzFan(Q);
573}
Note: See TracBrowser for help on using the repository browser.