Changeset e653c9 in git
- Timestamp:
- Jul 7, 2017, 2:46:14 PM (7 years ago)
- Branches:
- (u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
- Children:
- fb98c8dd692771969a7b893425029e42cd04a51d
- Parents:
- 67f9e8e1fc9856b995aae2d27c9a3848419e3d13
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/gitfan.lib
r67f9e8 re653c9 1 /////////////////////////////////////////////////////////////////// /2 version="version gitfan.lib 4 -1-0-3 Jun_2017 ";1 /////////////////////////////////////////////////////////////////// 2 version="version gitfan.lib 4.0.2.0 Apr_2015 "; // $Id$ 3 3 category="Algebraic Geometry"; 4 4 info=" 5 LIBRARY: 6 7 AUTHORS: Janko Boehm, boehm at mathematik.uni-kl.de @*8 Simon Keicher, keicher at mail.mathematik.uni-tuebingen.de @* 9 Yue Ren, ren at mathematik.uni-kl.de @* 5 LIBRARY: gitfan.lib Compute GIT-fans. 6 7 AUTHORS: 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 10 11 11 OVERVIEW: 12 This library allows you to calculate GIT-fans, torus orbits and GKZ-fans. 13 14 In provides features to make use of symmetries of the torus action under consideration. 15 16 The main procedure is GITfan which can be directly applied to an ideal and a grading matrix 17 encoding the torus action, and returns a fan, the associated GIT-fan. 18 We also provide various procedures implementing substeps of the algorithm 19 to deal with large computations. 20 21 The library uses the package 'gfanlib' by Anders N. Jensen. 22 23 For notation, background, and algorithms see [BKR16]. 24 25 Functions produce debug output if printlevel is positive. 26 27 Elements of the symmetric group Sn of type permutation can be created by 28 the function permutationFromIntvec. 29 The images of 1,...,n can be obtained by permutationToIntvec. 30 Composition of permutations can be done by the *-Operator, 31 also powers can be computed in the usual way. 32 33 34 REFERENCES: 35 36 [BKR16] J. Boehm, S. Keicher, Y. Ren: Computing GIT-Fans with Symmetry and the Mori Chamber Decomposition of M06bar, https://arxiv.org/abs/1603.09241 37 38 TYPES: 39 permutation; Permutation in map representation. 12 This library computes GIT-fans, torus orbits and GKZ-fans. 13 It uses the package 'gfanlib' by Anders N. Jensen 14 and some algorithms have been outsourced to C++ to improve the performance. 15 Check https://github.com/skeicher/gitfan_singular for updates. 16 17 KEYWORDS: library; gitfan; GIT; geometric invariant theory; quotients 40 18 41 19 PROCEDURES: 42 isAface(ideal,intvec); Checks whether the given face is an a-face. 43 afaces(ideal); Returns a list of intvecs that correspond to the set of all a-faces, optionally for given list of simplex faces. 44 fullDimImages(list,intmat); Finds the afaces which have a full-dimensional projection. 45 minimalAfaces(list); compute the minimal a-faces among the a-faces with full dimensional projection. 46 orbitCones(list,intmat); Returns the list of all orbit cones. 47 GITcone(list,bigintmat); Returns the GIT-cone containing the given weight vector. 48 GITfan(ideal,intmat); Compute GIT-fan. 49 GITfanFromOrbitCones(list,intmat,cone); Compute GIT-fan from orbit cones. 50 GITfanParallel(list,intmat,cone); Compute GIT-fan in parallel from orbit cones. 51 GKZfan(intmat); Returns the GKZ-fan of the matrix Q. 52 movingCone(intmat); Compute the moving cone. 53 computeAfaceOrbits(list,list); Compute orbits of a-faces under a permutation group action. 54 minimalAfaceOrbits(list); Find the minimal a-face orbits. 55 orbitConeOrbits(list,intmat); Project the a-face orbits to orbit cone orbits. 56 minimalOrbitConeOrbits(list); Find the minimal orbit cone orbits. 57 intersectOrbitsWithMovingCone(list,cone); Intersect orbit cone orbits with moving cone. 58 groupActionOnQImage(list,intmat); Determine the induced group action in the target of the grading matrix. 59 groupActionOnHashes(list,list); Determine the induced group action on the set of orbit cones. 60 storeActionOnOrbitConeIndices(list, string); Write the group action on the set of orbit cones to a file. 61 permutationFromIntvec(intvec); Create a permutation from an intvec of images. 62 permutationToIntvec(permutation); Return the intvec of images. 63 evaluateProduct(list,list); Evaluate a list of products of group elements in terms of a given representation of the elements as permutations. 64 GITfanSymmetric(list,intmat,cone,list); Compute GIT-fan from orbit cones by determining a minimal representing set for the orbits of maximal dimensional GIT-cones. 65 GITfanParallelSymmetric(list, intmat, cone, list); Compute GIT-fan in parallel from orbit cones by determining a minimal representing set for the orbits of maximal dimensional GIT-cones. 66 bigintToBinary(bigint,int); Convert bigint into a sparse binary represenation specifying the indices of the one-entries 67 binaryToBigint(intvec); Convert sparse binary represenation specifying the indices of the one-entries to bigint 68 applyPermutationToIntvec(intvec,permutation); Apply permutation to a set of integers represented as an intvec 69 hashToCone(bigint,list); Convert a bigint hash to a GIT-cone 70 hashesToFan(list hashes,list OC) 71 72 KEYWORDS: library, gitfan, GIT, geometric invariant theory, quotients 20 afaces(ideal); Returns a list of intvecs that correspond to all a-faces 21 gitCone(ideal,bigintmat,bigintmat); Returns the GIT-cone around the given weight vector w 22 gitFan(ideal,bigintmat); Returns the GIT-fan of the H-action defined by Q on V(a) 23 gkzFan(bigintmat); Returns the GKZ-fan of the matrix Q 24 isAface(ideal,intvec); Checks whether intvec corresponds to an ideal-face 25 orbitCones(ideal,bigintmat); Returns the list of all projected a-faces 73 26 "; 74 27 75 76 LIB "customstd.lib"; 77 LIB "linalg.lib"; 78 LIB "multigrading.lib"; 79 LIB "parallel.lib"; 80 81 static proc mod_init() 82 { 83 LIB "gfanlib.so"; 84 LIB "gitfan.so"; 85 newstruct("permutation","intvec image"); 86 system("install","permutation","*",composePermutations,2); 87 system("install","permutation","^",permutationPower,2); 88 system("install","permutation","print",printPermutation,1); 89 } 90 91 static proc emptyString(int n) 92 { 93 string st; 94 for (int i = 1; i<=n;i++){ 95 st=st+" "; 96 } 97 return(st);} 98 99 proc printPermutation(permutation sigma) 100 { 101 intvec v = permutationToIntvec(sigma); 102 string vsrc,vimg; 103 for (int i = 1; i<=size(v);i++){ 104 vsrc=vsrc+emptyString(1+size(string(size(v)))-size(string(i)))+string(i); 105 } 106 for (i = 1; i<=size(v);i++){ 107 vimg=vimg+emptyString(1+size(string(size(v)))-size(string(v[i])))+string(v[i]); 108 } 109 print("|"+vsrc+"|"); 110 print("|"+vimg+"|"); 111 } 112 28 LIB "parallel.lib"; // for parallelWaitAll 113 29 114 30 //////////////////////////////////////////////////// 115 // converts e.g. n=5 to its binary representation, i.e. 0,1,0,1 116 // if r = 3. 117 // and stores it in an intvec. 118 // r gives the bound for n <= 2^r: 31 proc mod_init() 32 { 33 LIB"customstd.so"; 34 LIB"gfanlib.so"; 35 } 36 119 37 static proc int2face(int n, int r) 120 38 { … … 123 41 int n0 = n; 124 42 125 while(n0 > 0){ 126 while(2^k > n0){ 43 while(n0 > 0) 44 { 45 while(2^k > n0) 46 { 127 47 k--; 128 48 //v[size(v)+1] = 0; … … 136 56 return(v); 137 57 } 58 59 ///////////////////////////////// 60 61 proc isAface(ideal a, intvec gam0) 62 "USAGE: isAface(a,gam0); a: ideal, gam0:intvec 63 PURPOSE: Checks whether the face of the positive orthant, 64 that is spanned by all i-th unit vectors, 65 where i ranges amongst the entries of gam0, 66 is an a-face. 67 RETURN: int 68 EXAMPLE: example isaface; shows an example 69 " 70 { 71 poly pz; 72 73 // special case: gam0 is the zero-cone: 74 if (size(gam0) == 1 and gam0[1] == 0) 75 { 76 ideal G; 77 78 // is an a-face if and only if RL0(0,...,0) = const 79 // set all entries to 0: 80 int i; 81 for (int k = 1; k <= ncols(a); k++) 82 { 83 pz = subst(a[k], var(1), 0); 84 for (i = 2; i <= nvars(basering); i++) 85 { 86 pz = subst(pz, var(i), 0); 87 } 88 G = G, pz; 89 } 90 91 G = std(G); 92 93 // monomial inside?: 94 if(1 == G) 95 { 96 return(0); 97 } 98 99 return(1); 100 } 101 102 103 // ring is too big: Switch to KK[T_i; e_i\in gam0] instead: 104 def R = basering; 105 string initNewRing = "ring Rgam0 = 0,("; 106 for (int i=1; i<size(gam0); i++) 107 { 108 initNewRing = initNewRing + string(var(gam0[i])) + ","; 109 } 110 initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),dp;"; 111 execute(initNewRing); 112 kill i; 113 114 ideal agam0 = imap(R,a); 115 116 poly p = var(1); // first entry of g; p = prod T_i with i element of g 117 for (int i = 2; i <= nvars(basering); i++ ) 118 { 119 p = p * var(i); 120 } 121 // p is now the product over all T_i, with e_i in gam0 122 123 agam0 = agam0, p - 1; // rad-membership 124 ideal G = std(agam0); 125 126 // does G contain 1?, i.e. is G = 1? 127 if(G <> 1) 128 { 129 return(1); // true 130 } 131 132 return(0); // false 133 } 138 134 example 139 135 { 140 136 echo = 2; 141 int n = 5; 142 int r = 4; 143 int2face(n, r); 144 145 n = 1; 146 r = 1; 147 int2face(n,r); 148 } 149 150 //////// 151 152 137 138 ring R = 0,(T(1..4)),dp; 139 ideal I = T(1)*T(2)-T(4); 140 141 intvec w = 1,4; 142 intvec v = 1,2,4; 143 144 isAface(I,w); // should be 0 145 "-----------"; 146 isAface(I,v); // should be 1 147 } 153 148 154 149 //////////////////////////////////////////////////// 155 150 151 proc afacesPart(ideal a, int d, int start, int end, int r) 152 { 153 intvec gam0; 154 int i; 155 list AF; 156 157 for(i = start; i <= end; i++) 158 { 159 if(i < 2^r) 160 { 161 gam0 = int2face(i,r); 162 163 // take gam0 only if it has 164 // at least d rays: 165 if(size(gam0) >= d) 166 { 167 if (isAface(a,gam0)) 168 { 169 AF[size(AF) + 1] = gam0; 170 } 171 } 172 } 173 } 174 return(AF); 175 } 176 177 //////////////////////////////////////////////////// 178 156 179 proc afaces(ideal a, list #) 157 "USAGE: afaces(a [,L]); a: ideal, L: list of intvecs 158 PURPOSE: Returns a list of all a-faces (considered as intvecs of 0 and 1, where the i-th entry is 1 if the cone has the i-th unit basis vector as a generator), if L is specified only the faces of the simplex listed in L are considered (e.g. representatives with respect to a group action). 159 RETURN: a list of intvecs 180 "USAGE: afaces(a, b, c); a: ideal, d: int, c: int 181 PURPOSE: Returns a list of all a-faces (represented by intvecs). 182 Moreover, it is possible to specify a dimensional bound b, 183 upon which only cones of that dimension and above are returned. 184 Lastly, as the computation is parallizable, one can specify c, 185 the number of cores to be used by the computation. 186 RETURN: a list of intvecs 160 187 EXAMPLE: example afaces; shows an example 161 188 " 162 189 { 190 int d = 1; 191 int ncores = 1; 192 193 if ((size(#) > 0) and (typeof(#[1]) == "int")) 194 { 195 d = #[1]; 196 } 197 198 if ((size(#) > 1) and (typeof(#[2]) == "int")) 199 { 200 ncores = #[2]; 201 } 202 163 203 list AF; 164 if ((size(#)>0) and (typeof(#[1])=="intvec")){165 list L = #;166 for(int i = 1; i<=size(L); i++ ){167 dbprint("representative "+string(i)+" of "+string(size(L)));168 if (isAface(a,L[i])){169 AF[size(AF) + 1] = L[i];170 }171 }172 return(AF);173 }174 175 204 intvec gam0; 176 205 int r = nvars(basering); 177 206 178 207 // check if 0 is an a-face: 179 int bd;180 if (size(#)>0){bd=#[1];}181 208 gam0 = 0; 182 if ( size(gam0)>=bd){183 if (isAface(a,gam0)){209 if (isAface(a,gam0)) 210 { 184 211 AF[size(AF) + 1] = gam0; 185 186 } 212 } 213 187 214 // check for other a-faces: 188 for(int i = 1; i < 2^r; i++ ){ 215 // make ncores processes: 216 int step = 2^r div ncores; 217 int i; 218 219 list args; 220 for(int k = 0; k < ncores; k++) 221 { 222 args[size(args) + 1] = list(a, d, k * step + 1, (k+1) * step, r); 223 } 224 225 string command = "afacesPart"; 226 list out = parallelWaitAll(command, args); 227 228 // do remaining ones: 229 for(i = ncores * step +1; i < 2^r; i++) 230 { 231 "another one needed"; 189 232 gam0 = int2face(i,r); 190 if (size(gam0)>=bd){ 191 if (isAface(a,gam0)){ 233 234 // take gam0 only if it has 235 // at least d rays: 236 if(size(gam0) >= d) 237 { 238 if (isAface(a,gam0)) 239 { 192 240 AF[size(AF) + 1] = gam0; 193 241 } … … 195 243 } 196 244 197 //"done checking a-faces!"; 245 // read out afaces of out into AF: 246 for(i = 1; i <= size(out); i++) 247 { 248 AF = AF + out[i]; 249 } 250 198 251 return(AF); 199 252 } … … 205 258 ideal a = T(1)+T(2)+T(3); 206 259 207 list F = afaces(a );260 list F = afaces(a,3,4); 208 261 print(F); 209 262 print(size(F)); … … 213 266 ideal a2 = T(2)^2*T(3)^2+T(1)*T(3); 214 267 215 list F2 = afaces(a2 );268 list F2 = afaces(a2,3,4); 216 269 print(F2); 217 270 print(size(F2)); … … 221 274 ideal a3 = 0; 222 275 223 list F3 = afaces(a3 );276 list F3 = afaces(a3,3,4); 224 277 print(F3); 225 278 print(size(F3)); 226 279 227 228 // 4th ex // 229 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 230 ideal J = 231 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 232 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 233 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 234 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 235 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 236 list F4 = afaces(J); 280 // bigger example // 281 ring R = 0,T(1..15),dp; 282 ideal a = 283 T(1)*T(10)-T(2)*T(7)+T(3)*T(6), 284 T(1)*T(11)-T(2)*T(8)+T(4)*T(6), 285 T(1)*T(12)-T(2)*T(9)+T(5)*T(6), 286 T(1)*T(13)-T(3)*T(8)+T(4)*T(7), 287 T(1)*T(14)-T(3)*T(9)+T(5)*T(7), 288 T(1)*T(15)-T(4)*T(9)+T(5)*T(8), 289 T(2)*T(13)-T(3)*T(11)+T(4)*T(10), 290 T(2)*T(14)-T(3)*T(12)+T(5)*T(10), 291 T(2)*T(15)-T(4)*T(12)+T(5)*T(11), 292 T(3)*T(15)-T(4)*T(14)+T(5)*T(13), 293 T(6)*T(13)-T(7)*T(11)+T(8)*T(10), 294 T(6)*T(14)-T(7)*T(12)+T(9)*T(10), 295 T(6)*T(15)-T(8)*T(12)+T(9)*T(11), 296 T(7)*T(15)-T(8)*T(14)+T(9)*T(13), 297 T(10)*T(15)-T(11)*T(14)+T(12)*T(13); 298 299 int t = timer; 300 list F4 = afaces(a,0,2); 237 301 print(size(F4)); 238 } 239 240 241 242 static proc saturateWithRespectToVariable(ideal I, int k) 243 { 244 // "saturating with respect to variable "+string(k); 245 ASSUME(1,k>=1); 246 ASSUME(1,k<=nvars(basering)); 247 248 def origin = basering; 249 int n = nvars(basering); 250 intvec weightVector = ringlist(origin)[3][1][2]; 251 252 string newVars; 253 intvec newWeightVector; 254 for (int i=1; i<k; i++) 255 { 256 newVars = newVars+string(var(i))+","; 257 newWeightVector[i]=weightVector[i]; 258 } 259 for (i=k+1; i<=n; i++) 260 { 261 newVars = newVars+string(var(i))+","; 262 newWeightVector[i-1]=weightVector[i]; 263 } 264 newVars = newVars+string(var(k)); 265 newWeightVector[n]=weightVector[k]; 266 execute("ring ringForSaturation = ("+charstr(origin)+"),("+newVars+"),(wp(newWeightVector));"); 267 268 ideal I = satstd(imap(origin,I)); 269 I = simplify(I,2+4+32); 270 271 // "finished saturating with respect to variable "+string(k); 272 setring origin; 273 return (imap(ringForSaturation,I)); 274 } 275 276 static proc stepwiseSaturation(ideal I) 277 { 278 if (I!=1) 279 { 280 list variablesToBeSaturated; 281 int l = nvars(basering); 282 for (int i=1; i<=l; i++) 283 { variablesToBeSaturated[i]=l-i+1; } 284 285 while (size(variablesToBeSaturated)>0) 286 { 287 I = saturateWithRespectToVariable(I,variablesToBeSaturated[1]); 288 variablesToBeSaturated = delete(variablesToBeSaturated,1); 289 if ((I==1) || (I==-1)) 290 { 291 break; 292 } 293 } 294 } 295 296 return (I); 297 } 298 299 300 301 302 303 304 proc isAface(ideal a, intvec gam0) 305 "USAGE: isAface(a,gam0); a: ideal, gam0:intvec 306 PURPOSE: Checks whether gam0 is an a-face w.r.t. the ideal a. 307 RETURN: int 308 EXAMPLE: example isaface; shows an example 309 " 310 { 311 // special case: gam0 is the zero-cone: 312 if (size(gam0) == 1 and gam0[1] == 0){ 313 poly pz; 314 ideal G; 315 int i; 316 for (int k = 1; k <= size(a); k++) { 317 pz = subst(a[k], var(1), 0); 318 for (i = 2; i <= nvars(basering); i++) { 319 pz = subst(pz, var(i), 0); 320 } 321 G = G, pz; 322 } 323 G = std(G); 324 // monomial inside?: 325 if(G == 1){ 326 return(0); 327 } 328 return(1); 329 } 330 331 string initNewRing = "ring Rgam0 = 0,("; 332 333 intvec w = ringlist(basering)[3][1][2]; 334 intvec w0; 335 for (int i=1; i<size(gam0); i++) 336 { 337 // take only entries i of w0 with i in gam0 338 initNewRing = initNewRing + string(var(gam0[i])) + ","; 339 w0[i] = w[gam0[i]]; 340 } 341 w0[size(gam0)] = w[gam0[size(gam0)]]; 342 343 initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),wp(w0);"; 344 def R = basering; 345 execute(initNewRing); 346 347 ideal agam0 = imap(R,a); 348 349 ideal G = stepwiseSaturation(agam0); 350 351 if(G == 1) 352 { 353 return(0); // true 354 } 355 return(1); // false 356 } 357 example 358 { 359 echo = 2; 360 361 ring R = 0,(T(1..4)),dp; 362 ideal I = T(1)*T(2)-T(4); 363 364 intvec w = 1,4; 365 intvec v = 1,2,4; 366 367 isAface(I,w); // should be 0 368 isAface(I,v); // should be 1 369 } 370 371 372 proc applyPermutationToIntvec(intvec v, permutation g) 373 "USAGE: applyPermutationToIntvec(v,g); v: intvec, g:permutation 374 PURPOSE: Apply g to the set of entries of v. The result is a sorted intvec. We assume that the entries of v are valid arguments of g. We do not require that the input is sorted. 375 RETURN: intvec. 376 EXAMPLE: example applyPermutationToIntvec; shows an example 377 " 378 { 379 intvec gv = composeIntvecs(permutationToIntvec(g),v); 380 //for (int i=1;i<=size(v);i++){ 381 // gv[i]=g[v[i]]; 382 //} 383 return(sort(gv)[1]); 384 } 385 example 386 { 387 permutation g = permutationFromIntvec(intvec(10, 9, 7, 4, 8, 6, 3, 5, 2, 1)); 388 applyPermutationToIntvec(intvec(1, 3, 4, 6, 8),g); 389 } 390 391 static proc isContained(intvec v, list orbit){ 392 for (int i=1;i<=size(orbit);i++){ 393 if (v==orbit[i]){return(1);} 394 } 395 return(0); 396 } 397 398 399 static proc computeSimplexOrbit(intvec v, list G){ 400 list orbit; 401 intvec gv; 402 for (int i=1;i<=size(G);i++){ 403 gv=applyPermutationToIntvec(v,G[i]); 404 if (isContained(gv,orbit)==0){orbit[size(orbit)+1]=gv;} 405 } 406 return(orbit); 407 } 408 409 410 proc computeAfaceOrbits(list AF, list G) 411 "USAGE: computeAfaceOrbits(AF,G); AF list of intvecs, G: list of permutations 412 PURPOSE: Computes the orbits of the afaces in the list AF under the group action in G, where G is a list of permutations. We assume that the elements of G form a group and the first entry corresponds to the neutral element. 413 RETURN: a list of lists of intvecs 414 EXAMPLE: example computeAfaceOrbits; shows an example 415 " 416 { 417 list listorbits; 418 for (int i=1;i<=size(AF);i++){ 419 dbprint(i); 420 listorbits[i]=computeSimplexOrbit(AF[i],G); 421 } 422 return(listorbits);} 423 example 424 { 425 echo = 2; 426 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 427 ideal J = 428 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 429 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 430 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 431 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 432 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 433 intmat Q[5][10] = 434 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 435 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 436 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 437 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 438 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 439 list simplexSymmetryGroup = G25Action(); 440 list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ), 441 intvec( 1, 2, 3, 5, 6 ), 442 intvec( 1, 2, 3, 5, 7 ), 443 intvec( 1, 2, 3, 5, 10 ), 444 intvec( 1, 2, 3, 7, 9 ), 445 intvec( 1, 2, 6, 9, 10 ), 446 intvec( 1, 2, 3, 4, 5, 6 ), 447 intvec( 1, 2, 3, 4, 5, 10 ), 448 intvec( 1, 2, 3, 5, 6, 8 ), 449 intvec( 1, 2, 3, 5, 6, 9 ), 450 intvec( 1, 2, 3, 5, 7, 10 ), 451 intvec( 1, 2, 3, 7, 9, 10 ), 452 intvec( 1, 2, 3, 4, 5, 6, 7 ), 453 intvec( 1, 2, 3, 4, 5, 6, 8 ), 454 intvec( 1, 2, 3, 4, 5, 6, 9 ), 455 intvec( 1, 2, 3, 5, 6, 9, 10 ), 456 intvec( 1, 2, 3, 4, 5, 6, 7, 8 ), 457 intvec( 1, 2, 3, 4, 5, 6, 9, 10 ), 458 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ), 459 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ); 460 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 461 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 462 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 463 apply(afaceOrbits,size);} 464 465 466 467 468 static proc isSubset(intvec v, intvec w) 469 { 470 int i,j; 471 int jStart=1; 472 for (i=1; i<=size(v); i++) 473 { 474 for (j=jStart; j<=size(w); j++) 475 { 476 if (v[i]==w[j]) 477 { 478 break; 479 } 480 } 481 if (j<=size(w)) 482 { 483 jStart = j+1; 484 } 485 else 486 { 487 return (0); 488 } 489 } 490 return (1); 491 } 492 493 494 495 static proc containsOrbitd(list A, list B) 496 { 497 intvec a = A[1]; 498 for (int i=1; i<=size(B); i++) 499 { 500 if (isSubset(a,B[i])) 501 { 502 return (1); 503 } 504 } 505 return (0); 506 } 507 508 509 proc minimalAfaceOrbits(list listOfOrbits) 510 "USAGE: minimalAfaceOrbits(afaceOrbits); afaceOrbits: list 511 PURPOSE: Returns a list of all minimal a-face orbits. 512 RETURN: a list of intvecs 513 EXAMPLE: example minimalAfaceOrbits; shows an example 514 " 515 { 516 int i,j; 517 list L = listOfOrbits; 518 for (i=1; i<=size(L); i++) 519 { 520 dbprint(string(i)+" of "+string(size(L))); 521 for (j=1; j<=size(L); j++) 522 { 523 if (i!=j) 524 { 525 if (containsOrbitd(L[i],L[j])) 526 { 527 dbprint("werfe raus (Laenge): " + string(size(L[j]))); 528 L = delete(L,j); 529 if (j<i) 530 { 531 i = i-1; 532 } 533 j = j-1; 534 } 535 } 536 } 537 } 538 return(L); 539 } 540 example 541 { 542 echo = 2; 543 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 544 ideal J = 545 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 546 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 547 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 548 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 549 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 550 intmat Q[5][10] = 551 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 552 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 553 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 554 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 555 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 556 list simplexSymmetryGroup = G25Action(); 557 list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ), 558 intvec( 1, 2, 3, 5, 6 ), 559 intvec( 1, 2, 3, 5, 7 ), 560 intvec( 1, 2, 3, 5, 10 ), 561 intvec( 1, 2, 3, 7, 9 ), 562 intvec( 1, 2, 6, 9, 10 ), 563 intvec( 1, 2, 3, 4, 5, 6 ), 564 intvec( 1, 2, 3, 4, 5, 10 ), 565 intvec( 1, 2, 3, 5, 6, 8 ), 566 intvec( 1, 2, 3, 5, 6, 9 ), 567 intvec( 1, 2, 3, 5, 7, 10 ), 568 intvec( 1, 2, 3, 7, 9, 10 ), 569 intvec( 1, 2, 3, 4, 5, 6, 7 ), 570 intvec( 1, 2, 3, 4, 5, 6, 8 ), 571 intvec( 1, 2, 3, 4, 5, 6, 9 ), 572 intvec( 1, 2, 3, 5, 6, 9, 10 ), 573 intvec( 1, 2, 3, 4, 5, 6, 7, 8 ), 574 intvec( 1, 2, 3, 4, 5, 6, 9, 10 ), 575 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ), 576 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ); 577 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 578 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 579 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 580 apply(afaceOrbits,size); 581 list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits); 582 apply(minAfaceOrbits,size); 583 } 584 585 586 static proc list2string(list L){ 587 string s = ""; 588 589 for(int i = 1; i <=size(L); i++){ 590 s = s + string(size(L[i])) + ", "; 591 } 592 593 return(s); 594 } 595 596 597 598 proc fullDimImages(list afaces,intmat Q) 599 "USAGE: fullDimImages(afaces, Q); afaces: list, Q: intmat 600 PURPOSE: Determines the a-faces (represented as intvecs) from the list afaces which have a full-dimensional projection with respect to Q. 601 RETURN: a list of intvecs 602 EXAMPLE: example fullDimImages; shows an example 603 " 604 { 605 list L; 606 for (int i=1; i<=size(afaces); i++) 607 { 608 intvec gam0 = afaces[i]; 609 intmat Qgam0[nrows(Q)][size(gam0)] = Q[intvec(1..nrows(Q)),gam0]; 610 //cone c = coneViaPoints(Qgam0); 611 if(rank(Qgam0) == nrows(Q)){ 612 L[size(L)+1]=afaces[i]; 613 } 614 kill gam0,Qgam0; 615 } 616 return(L); 617 } 618 example 619 { 620 echo=2; 621 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 622 ideal J = 623 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 624 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 625 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 626 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 627 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 628 intmat Q[5][10] = 629 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 630 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 631 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 632 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 633 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 634 635 list AF= afaces(J,nrows(Q)); 636 size(AF); 637 size(fullDimImages(AF,Q)); 638 } 639 640 641 642 643 proc orbitConeOrbits(list F, intmat Q) 644 "USAGE: orbitConeOrbits(F, Q); F: list, Q: intmat 645 PURPOSE: Projects a list F of a-face orbits to the orbit cones with respect to Q. The function checks whether the projections are of full dimension and returns an error otherwise. 646 RETURN: a list of lists of cones 647 EXAMPLE: example orbitConeOrbits; shows an example 648 " 649 { 650 list OC; 651 int j; 652 intmat Qt = transpose(Q); 653 for(int i = 1; i <= size(F); i++){ 654 list Orbit = F[i]; 655 list U; 656 for(j = 1; j <= size(Orbit); j++){ 657 intvec aface = Orbit[j]; 658 intmat Qgam0[size(aface)][ncols(Qt)] = Qt[aface,intvec(1..ncols(Qt))]; 659 cone c = coneViaPoints(Qgam0); 660 if (dimension(c)!=ncols(Qt)){ERROR("cone should have full dimension");} 661 // insert it even if it is already inside (mincones will take care of this) 662 U[size(U)+1] = c; 663 kill aface; 664 kill Qgam0; 665 kill c; 666 } 667 OC[size(OC)+1] = U; 668 dbprint("current size of OC:"); 669 dbprint(size(OC)); 670 kill U; 671 kill Orbit; 672 } 673 return(OC); 674 } 675 example 676 { 677 echo = 2; 678 // Note that simplexOrbitRepresentatives and simplexSymmetryGroup are subsets of the actual sets for G25. For the full example see the examples in the documentation 679 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 680 ideal J = 681 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 682 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 683 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 684 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 685 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 686 intmat Q[5][10] = 687 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 688 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 689 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 690 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 691 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 692 list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ), 693 intvec( 1, 2, 3, 5, 6 ), 694 intvec( 1, 2, 3, 5, 7 ), 695 intvec( 1, 2, 3, 5, 10 ), 696 intvec( 1, 2, 3, 7, 9 ), 697 intvec( 1, 2, 3, 4, 5, 6, 9, 10 ), 698 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ), 699 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ); 700 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 701 list simplexSymmetryGroup = permutationFromIntvec(intvec( 1 .. 10 )), 702 permutationFromIntvec(intvec( 1, 2, 4, 3, 5, 7, 6, 9, 8, 10 )), 703 permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )), 704 permutationFromIntvec(intvec( 1, 3, 4, 2, 6, 7, 5, 10, 8, 9 )), 705 permutationFromIntvec(intvec( 1, 4, 2, 3, 7, 5, 6, 9, 10, 8 )), 706 permutationFromIntvec(intvec( 1, 4, 3, 2, 7, 6, 5, 10, 9, 8 )); 707 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 708 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 709 apply(afaceOrbits,size); 710 list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits); 711 apply(minAfaceOrbits,size); 712 list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q); 713 } 714 715 716 proc minimalOrbitConeOrbits(list listOfOrbits) 717 "USAGE: minimalOrbitConeOrbits(listOfOrbits); listOfOrbits: list of lists of cones 718 PURPOSE: Minimizes a list of orbit cone orbits. 719 RETURN: a list of lists of cones 720 EXAMPLE: example minimalOrbitConeOrbits; shows an example 721 " 722 { 723 int i,j; 724 list L = listOfOrbits; 725 726 727 for (i=1; i<=size(L); i++) 728 { 729 dbprint("::::: " + string(i)+" of "+string(size(L))); 730 731 for (j=1; j<=size(L); j++) 732 { 733 dbprint("outer: "+string(i)+" inner: " + string(j)); 734 if (i!=j) 735 { 736 if (containsOrbitdCone(L[i],L[j])) 737 { 738 dbprint("werfe raus, i="+string(i)+", j="+string(j)+" (Laenge): " + string(size(L[j]))); 739 L = delete(L,j); 740 if (j<i) 741 { 742 i = i-1; 743 } 744 j = j-1; 745 } 746 } 747 } 748 } 749 return(L); 750 } 751 example 752 { 753 echo = 2; 754 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 755 ideal J = 756 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 757 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 758 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 759 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 760 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 761 intmat Q[5][10] = 762 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 763 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 764 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 765 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 766 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 767 list simplexSymmetryGroup = G25Action(); 768 list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ), 769 intvec( 1, 2, 3, 5, 6 ), 770 intvec( 1, 2, 3, 5, 7 ), 771 intvec( 1, 2, 3, 5, 10 ), 772 intvec( 1, 2, 3, 7, 9 ), 773 intvec( 1, 2, 6, 9, 10 ), 774 intvec( 1, 2, 3, 4, 5, 6 ), 775 intvec( 1, 2, 3, 4, 5, 10 ), 776 intvec( 1, 2, 3, 5, 6, 8 ), 777 intvec( 1, 2, 3, 5, 6, 9 ), 778 intvec( 1, 2, 3, 5, 7, 10 ), 779 intvec( 1, 2, 3, 7, 9, 10 ), 780 intvec( 1, 2, 3, 4, 5, 6, 7 ), 781 intvec( 1, 2, 3, 4, 5, 6, 8 ), 782 intvec( 1, 2, 3, 4, 5, 6, 9 ), 783 intvec( 1, 2, 3, 5, 6, 9, 10 ), 784 intvec( 1, 2, 3, 4, 5, 6, 7, 8 ), 785 intvec( 1, 2, 3, 4, 5, 6, 9, 10 ), 786 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ), 787 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ); 788 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 789 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 790 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 791 apply(afaceOrbits,size); 792 list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits); 793 apply(minAfaceOrbits,size); 794 list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q); 795 apply(listOfOrbitConeOrbits,size); 796 list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits); 797 size(listOfMinimalOrbitConeOrbits); 798 } 799 800 801 static proc containsOrbitdCone(list A, list B) 802 { 803 cone a = A[1]; 804 for (int i=1; i<=size(B); i++) 805 { 806 if (isSubcone(a,B[i])) 807 { 808 dbprint("found subcone i="+string(i)); 809 return (1); 810 } 811 } 812 return (0); 813 } 814 815 static proc isSubcone(cone c1, cone c2) 816 { 817 return(containsInSupport(c2, c1)); 818 } 819 820 821 822 823 824 825 proc intersectOrbitsWithMovingCone(list OCmin,cone mov) 826 "USAGE: intersectOrbitsWithMovingCone(OCmin, mov); OCmin: list of lists of cones, mov: cone 827 PURPOSE: Intersects all cones in the orbits in OCmin with mov discarting all orbits of cones which are not of full dimension. The resulting orbits are duplicate free. 828 RETURN: a list of lists of cones 829 EXAMPLE: example intersectOrbitsWithMovingCone; shows an example 830 " 831 { 832 list OCmov; 833 cone ocmov; 834 int i,j; 835 int fulldim=dimension(mov); 836 for (i=1; i<=size(OCmin); i++) 837 { 838 dbprint("intersecting orbit "+string(i)+" with moving cone"); 839 ocmov = convexIntersection(OCmin[i][1],mov); 840 if (dimension(ocmov)==fulldim) 841 { 842 list currentOrbit; 843 for (j=1; j<=size(OCmin[i]); j++) 844 { 845 dbprint("checking cone "+string(j)+" of "+string(size(OCmin[i]))); 846 ocmov = convexIntersection(OCmin[i][j],mov); 847 ocmov = canonicalizeCone(ocmov); 848 if (!listContainsCone(currentOrbit,ocmov)) 849 { 850 dbprint("inserting cone"); 851 currentOrbit[size(currentOrbit)+1] = ocmov; 852 } 853 } 854 OCmov[size(OCmov)+1] = currentOrbit; 855 kill currentOrbit; 856 } 857 else 858 { 859 dbprint("intersection with moving cone lower-dimensional, abandoning orbit"); 860 } 861 } 862 return(OCmov); 863 } 864 example { 865 echo=2; 866 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 867 ideal J = 868 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 869 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 870 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 871 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 872 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 873 intmat Q[5][10] = 874 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 875 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 876 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 877 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 878 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 879 list simplexSymmetryGroup = G25Action(); 880 list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ), 881 intvec( 1, 2, 3, 5, 6 ), 882 intvec( 1, 2, 3, 5, 7 ), 883 intvec( 1, 2, 3, 5, 10 ), 884 intvec( 1, 2, 3, 7, 9 ), 885 intvec( 1, 2, 6, 9, 10 ), 886 intvec( 1, 2, 3, 4, 5, 6 ), 887 intvec( 1, 2, 3, 4, 5, 10 ), 888 intvec( 1, 2, 3, 5, 6, 8 ), 889 intvec( 1, 2, 3, 5, 6, 9 ), 890 intvec( 1, 2, 3, 5, 7, 10 ), 891 intvec( 1, 2, 3, 7, 9, 10 ), 892 intvec( 1, 2, 3, 4, 5, 6, 7 ), 893 intvec( 1, 2, 3, 4, 5, 6, 8 ), 894 intvec( 1, 2, 3, 4, 5, 6, 9 ), 895 intvec( 1, 2, 3, 5, 6, 9, 10 ), 896 intvec( 1, 2, 3, 4, 5, 6, 7, 8 ), 897 intvec( 1, 2, 3, 4, 5, 6, 9, 10 ), 898 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ), 899 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ); 900 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 901 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 902 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 903 apply(afaceOrbits,size); 904 list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits); 905 apply(minAfaceOrbits,size); 906 list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q); 907 cone mov = movingCone(Q); 908 intersectOrbitsWithMovingCone(listOfOrbitConeOrbits,mov); 909 } 910 911 912 static proc coneToString(cone C){ 913 bigintmat F = facets(C); 914 ring R=0,(x),dp; 915 matrix FF[nrows(F)][ncols(F)]; 916 int i,j; 917 for (i=1; i<=nrows(F);i++){ 918 for (j=1; j<=ncols(F);j++){ 919 FF[i,j]=number(F[i,j]); 920 } 921 } 922 return("bigintmat F["+string(nrows(F))+"]["+string(ncols(F))+"]="+string(FF)+";cone C=canonicalizeCone(coneViaInequalities(F));"); 923 } 924 925 proc storeOrbitConeOrbits(list OC, string fn) 926 "USAGE: storeOrbitConeOrbits(OC, fn); OC: list of lists of cones, fn: string 927 PURPOSE: Writes OC to file fn in Singular readable format. Can be importet using <. 928 RETURN: nothing 929 " 930 { 931 int i,j; 932 write(":w "+fn,"list OC;"); 933 for (i=1;i<=size(OC);i++){ 934 write(":a "+fn,"list OCj;"); 935 for (j=1;j<=size(OC[i]);j++){ 936 write(":a "+fn,coneToString(OC[i][j])); 937 write(":a "+fn,"OCj[size(OCj)+1]=C;"); 938 write(":a "+fn,"kill C;kill F;"); 939 } 940 write(":a "+fn,"OC[size(OC)+1]=OCj;"); 941 write(":a "+fn,"kill OCj;"); 942 } 943 } 944 945 static proc storeCone(cone C,string fn){ 946 intmat H = intmat(inequalities(C)); 947 write(":w "+fn,"bigintmat H["+string(nrows(H))+"]["+string(ncols(H))+"] ="); 948 write(":w "+fn,string(H)+";"); 949 write(":a "+fn,"cone C= coneViaInequalities(H); kill H;"); 950 } 951 952 static proc listToIntvec(list L){ 953 intvec v; 954 for (int i = 1;i<=size(L);i++){ 955 v[i]=L[i];} 956 return(v); 957 } 958 959 static proc intvecToList(intvec L){ 960 list v; 961 for (int i = 1;i<=size(L);i++){ 962 v[i]=L[i];} 963 return(v); 964 } 965 966 static proc stackMatrices(bigintmat M, bigintmat N){ 967 if (ncols(M)!=ncols(N)){ERROR("matrices should have the same number of columns");} 968 bigintmat MN[nrows(M)+nrows(N)][ncols(N)]; 969 MN[1..nrows(M),1..ncols(M)]=M[1..nrows(M),1..ncols(M)]; 970 MN[(nrows(M)+1)..(nrows(M)+nrows(N)),1..ncols(M)]=N[1..nrows(N),1..ncols(N)]; 971 return(MN); 972 } 973 974 975 proc movingCone(intmat Q) 976 "USAGE: movingCone(Q); Q: intmat 977 PURPOSE: Computes the moving cone from the grading matrix, with the degrees in the columns of Q. 978 RETURN: a cone 979 EXAMPLE: example movingCone; shows an example 980 " 981 { 982 cone Qgammadual = coneViaInequalities(transpose(Q)); 983 bigintmat R = facets(Qgammadual); 984 list idx1=intvecToList(intvec(1..nrows(R))); 985 intvec idx2=1..ncols(R); 986 intvec idx1iv; 987 bigintmat Ri[nrows(R)-1][ncols(R)]; 988 list C; 989 for (int i = 1; i<=nrows(R);i++){ 990 idx1iv = listToIntvec(delete(idx1,i)); 991 Ri=R[idx1iv,idx2]; 992 C[i]=coneViaPoints(Ri); 993 dbprint(string(i)+" "+string(nrows(facets(C[i])))); 994 } 995 cone mov = convexIntersection(C); 996 return(mov); 997 } 998 example 999 { 1000 echo=2; 1001 intmat Q[5][10] = 1002 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1003 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1004 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 1005 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 1006 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 1007 cone mov = movingCone(Q); 1008 mov; 1009 rays(mov); 1010 } 1011 1012 1013 static proc pivotIndices(matrix H) 1014 { 1015 intvec p; 1016 int pp; 1017 int i,j; 1018 int l=1; 1019 for (i=1; i<=ncols(H); i++) 1020 { 1021 for (j=nrows(H); j>=0; j--) 1022 { 1023 if (H[j,i]!=0) 1024 { 1025 if (j>pp) 1026 { 1027 p[l] = i; 1028 l++; 1029 pp = j; 1030 } 1031 break; 1032 } 1033 } 1034 } 1035 return (p); 1036 } 1037 1038 1039 1040 proc groupActionOnQImage(list G,intmat Q) 1041 "USAGE: groupActionOnQImage(G,Q); G: list of permutations, Q: intmat 1042 PURPOSE: Given the group G of permutations acting on the simplex on ncols(Q) objects, computes the corresponding group action on the image of Q. We assume that the basering has characteristic 0. 1043 RETURN: list of matrices 1044 EXAMPLE: example groupActionOnQImage; shows an example 1045 " 1046 { 1047 matrix Qmat = transpose(matrix(Q)); 1048 matrix H = gauss_nf(Qmat); 1049 intvec indices = pivotIndices(H); 1050 intmat Qbasis[nrows(Q)][size(indices)]=Q[1..nrows(Q),indices]; 1051 matrix QbasisInv = inverse(Qbasis); 1052 list L; 1053 intmat Qt = transpose(Q); 1054 for(int i = 1; i <= size(G); i++){ 1055 intvec sig = permutationToIntvec(G[i]); 1056 intmat Bsig = perm2mat(sig, indices,Qt); 1057 matrix Asig = Bsig * QbasisInv; 1058 L[size(L)+1] = matrixToIntmat(Asig); 1059 kill sig; 1060 kill Bsig; 1061 kill Asig; 1062 } 1063 return(L); 1064 } 1065 example { 1066 echo=2; 1067 ring R = 0,(x),dp; 1068 intmat Q[5][10] = 1069 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1070 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1071 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 1072 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 1073 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 1074 list generatorsG = permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )), 1075 permutationFromIntvec(intvec( 5, 7, 1, 6, 9, 2, 8, 4, 10, 3 )); 1076 groupActionOnQImage(generatorsG,Q); 1077 } 1078 1079 1080 1081 static proc perm2mat(intvec sig, intvec indices, intmat Q){ 1082 intvec sigind = 1:size(indices); 1083 for(int i = 1; i <= size(indices); i++){ 1084 if(indices[i] > size(sig)){ 1085 sigind[i] = indices[i]; 1086 } else { 1087 sigind[i] = sig[indices[i]]; 1088 } 1089 } 1090 intmat Asig[size(indices)][ncols(Q)]; 1091 for(i = 1; i <= size(sigind); i++){ 1092 Asig[i,1..ncols(Q)] = Q[sigind[i], 1..ncols(Q)]; 1093 } 1094 // the new basis must be in the cols: 1095 return(transpose(Asig)); 1096 } 1097 1098 /////// 1099 1100 1101 1102 1103 static proc imageCone(cone c, intmat A){ 1104 bigintmat ineqs = inequalities(c); 1105 cone cc = coneViaInequalities (ineqs * A); 1106 return(cc); 1107 } 1108 1109 1110 static proc matrixToIntmat(matrix A){ 1111 int i,j; 1112 intmat Aint[nrows(A)][ncols(A)]; 1113 for(i = 1; i<=nrows(A); i++){ 1114 for(j = 1; j<=ncols(A); j++){ 1115 if (deg(A[i,j])>0){ERROR("entries should be constants");} 1116 if (denominator(number(A[i,j]))!=1){ERROR("entries should be integers");} 1117 Aint[i,j]=int(number(A[i,j])); 1118 } 1119 } 1120 return(Aint); 1121 } 1122 1123 1124 proc groupActionOnHashes(list Asigma, list OCmov) 1125 "USAGE: groupActionOnHashes(Asigma,OCmov); Asigma: list, OCmov: list of list of cones 1126 PURPOSE: From the list of orbits of orbitcones, and the symmetry group represenation given by the matrices in Asigma, compute the corresponding permutation representation of the symmetry group on the orbit cones. The permutations are specified in a map representation of length the sum of the size of the orbits of OCmov. 1127 RETURN: list of permutations 1128 EXAMPLE: example groupActionOnHashes; shows an example 1129 " 1130 { 1131 // for each A in S6: 1132 // for each c in OC90: 1133 // store index of A*c 1134 // --> intvec v_A 1135 // store this in the list Ind 1136 list Ind; 1137 int i,j,b,k,sizepreviousorbits; 1138 list remaining; 1139 intmat A; 1140 cone c,Ac; 1141 for(i = 1; i<=size(Asigma); i++){ 1142 intvec vA; 1143 A = intmat(Asigma[i]); 1144 dbprint("element "+string(i)+" of symmetry group"); 1145 sizepreviousorbits=0; 1146 for(b = 1; b <= size(OCmov); b++){ 1147 remaining = intvecToList(intvec(1..size(OCmov[b]))); 1148 for(j = 1; j <= size(OCmov[b]); j++){ 1149 Ac = imageCone(OCmov[b][j], A); 1150 1151 // find out index: 1152 for(k= 1; k <= size(remaining); k++){ 1153 if(OCmov[b][remaining[k]] == Ac){ 1154 vA[j+sizepreviousorbits] = remaining[k]+sizepreviousorbits; 1155 remaining = delete(remaining,k); 1156 break; 1157 } 1158 } 1159 } 1160 1161 sizepreviousorbits=size(vA); 1162 } 1163 Ind[size(Ind)+1] = permutationFromIntvec(vA); 1164 dbprint("vA: "+string(vA)); 1165 kill vA; 1166 } 1167 return(Ind); 1168 } 1169 example 1170 { 1171 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 1172 ideal J = 1173 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 1174 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 1175 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 1176 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 1177 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 1178 intmat Q[5][10] = 1179 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1180 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1181 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 1182 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 1183 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 1184 list AF= afaces(J); 1185 list OC = orbitCones(AF,Q); 1186 list generatorsG = permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )), 1187 permutationFromIntvec(intvec( 5, 7, 1, 6, 9, 2, 8, 4, 10, 3 )); 1188 list Asigmagens = groupActionOnQImage(generatorsG,Q); 1189 groupActionOnHashes(Asigmagens,list(OC)); 1190 1191 1192 list simplexSymmetryGroup = G25Action(); 1193 list orb = findOrbits(simplexSymmetryGroup,nrows(Q)); 1194 list simplexOrbitRepresentatives; 1195 for (int i=1;i<=size(orb);i++){simplexOrbitRepresentatives[i]=orb[i][1];} 1196 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 1197 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 1198 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 1199 list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits); 1200 list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q); 1201 list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits); 1202 list Asigma = groupActionOnQImage(simplexSymmetryGroup,Q); 1203 groupActionOnHashes(Asigma,listOfOrbitConeOrbits); 1204 1205 } 1206 1207 1208 1209 1210 1211 1212 1213 1214 static proc composePermutationsGAP(permutation sigma, permutation tau){ 1215 intvec sigmaTauImage; 1216 1217 for (int i=1;i<=size(sigma.image);i++){ 1218 sigmaTauImage[i]=tau.image[sigma.image[i]]; 1219 } 1220 1221 permutation sigmaTau; 1222 sigmaTau.image = sigmaTauImage; 1223 return(sigmaTau); 1224 } 1225 1226 static proc composePermutations(permutation sigma, permutation tau){ 1227 intvec sigmaTauImage; 1228 1229 for (int i=1;i<=size(sigma.image);i++){ 1230 sigmaTauImage[i]=sigma.image[tau.image[i]]; 1231 } 1232 1233 permutation sigmaTau; 1234 sigmaTau.image = sigmaTauImage; 1235 return(sigmaTau); 1236 } 1237 1238 1239 static proc permutationPower(permutation sigma, int k){ 1240 int i; 1241 if (k==0) 1242 { 1243 permutation identity; 1244 identity.image = intvec(1..size(sigma.image)); 1245 return (identity); 1246 } 1247 if (k<0) 1248 { 1249 // if k<0 invert sigma and replace k with -k 1250 intvec sigmaImage = sigma.image; 1251 for (i=1; i<=size(sigmaImage); i++) 1252 { 1253 sigma.image[sigmaImage[i]] = i; 1254 } 1255 k = -k; 1256 } 1257 permutation sigmaToPowerK = sigma; 1258 for (i=2; i<=k; i++) 1259 { 1260 sigmaToPowerK = composePermutations(sigmaToPowerK, sigma); 1261 } 1262 return (sigmaToPowerK); 1263 } 1264 1265 proc permutationFromIntvec(intvec sigmaImage) 1266 "USAGE: permutationFromIntvec(sigmaImage); sigmaImage: intvec 1267 PURPOSE: Create a permutation from an intvec of images. 1268 RETURN: permutation 1269 EXAMPLE: example permutationFromIntvec; shows an example 1270 " 1271 { 1272 permutation sigma; 1273 sigma.image = sigmaImage; 1274 return (sigma); 1275 } 1276 1277 1278 proc evaluateProduct(list generatorsGperm, string st) 1279 "USAGE: evaluateProduct(generatorsGperm,st); generatorsGperm: list, st: string 1280 PURPOSE: Evaluates a formal product of variables xi in st, where xi corresponds to the permutation generatorsGperm[i]. 1281 RETURN: permutation 1282 EXAMPLE: example evaluateProduct; shows an example 1283 " 1284 { 1285 for (int i=1;i<=size(generatorsGperm);i++){ 1286 execute("permutation x"+string(i)+"= generatorsGperm[i];"); 1287 } 1288 if (st==""){ 1289 st="x1^0"; 1290 } 1291 execute("permutation sigma = "+st); 1292 return(sigma);} 1293 example 1294 { 1295 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 1296 ideal J = 1297 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 1298 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 1299 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 1300 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 1301 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 1302 intmat Q[5][10] = 1303 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1304 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1305 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 1306 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 1307 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 1308 list AF= afaces(J); 1309 list OC = orbitCones(AF,Q); 1310 list generatorsG = permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )), 1311 permutationFromIntvec(intvec( 5, 7, 1, 6, 9, 2, 8, 4, 10, 3 )); 1312 list Asigmagens = groupActionOnQImage(generatorsG,Q); 1313 //list actionOnOrbitconeIndicesForGenerators = groupActionOnHashes(Asigmagens,OC); 1314 string elementInTermsOfGenerators = 1315 "(x2^-1*x1^-1)^3*x1^-1"; 1316 //evaluateProduct(actionOnOrbitconeIndicesForGenerators, elementInTermsOfGenerators); 1317 } 1318 1319 1320 proc permutationToIntvec(permutation sigma) 1321 "USAGE: permutationToIntvec(sigma); sigma: permutation 1322 PURPOSE: Convert a permutation to an intvec of images. 1323 RETURN: intvec 1324 EXAMPLE: example permutationToIntvec; shows an example 1325 " 1326 {return(sigma.image);} 1327 example 1328 { 1329 permutation sigma = permutationFromIntvec(intvec( 1, 2, 4, 3, 5, 7, 6, 9, 8, 10 )); 1330 sigma; 1331 permutationToIntvec(sigma); 1332 } 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 proc storeActionOnOrbitConeIndices(list Ind,string fn) 1346 "USAGE: storeActionOnOrbitConeIndices(generatorsGperm,st); generatorsGperm: list, fn: string 1347 PURPOSE: Write the action on the set of orbit cones to the file fn in Singular readable format. 1348 RETURN: nothing 1349 " 1350 { 1351 string s = "list actionOnOrbitconeIndices;"; 1352 for(int i =1; i <= size(Ind); i++){ 1353 s = s + "intvec v = " + string(Ind[i]) + ";" + newline; 1354 s = s + "actionOnOrbitconeIndices[size(actionOnOrbitconeIndices)+1] = permutationFromIntvec(v);" + newline + "kill v;" + newline; 1355 } 1356 write(":w "+fn, s); 1357 } 1358 1359 1360 1361 1362 1363 static proc getNeighborHash(list OC, bigintmat w, bigintmat v, int mu) 1364 { 1365 int success = 0; 1366 int zz; 1367 intvec Jtmp; 1368 bigintmat wtmp; 1369 1370 while(!success) 1371 { 1372 mu = mu*2; 1373 wtmp = mu*w - v; 1374 success = 1; 1375 for(zz = 1; zz <= size(OC); zz++) 1376 { 1377 if(containsInSupport(OC[zz], wtmp)) 1378 { 1379 if(!containsInSupport(OC[zz], w)) 1380 { 1381 success = 0; 1382 Jtmp = 0; 1383 break; 1384 } 1385 // insert index zz: 1386 if(size(Jtmp) ==1 && Jtmp[1] == 0) 1387 { 1388 Jtmp[1] = zz; 1389 } 1390 else 1391 { 1392 Jtmp[size(Jtmp)+1] = zz; 1393 } 1394 } 1395 } 1396 } 1397 1398 return(Jtmp); 1399 } 1400 1401 1402 302 timer - t; 303 304 int t = timer; 305 list F4 = afaces(a,0); 306 print(size(F4)); 307 timer - t; 308 309 } 1403 310 1404 311 /////////////////////////////////////// 1405 312 1406 proc GITfanSymmetric(list OC, bigintmat Q, cone Qgamma, list actiononorbitcones, list #) 1407 "USAGE: GITfanSymmetric(OC, Q, Qgamma, actiononorbitcones [, file1, file2]); OC:list, Q:bigintmat, Qgamma:cone, actiononorbitcones: list of intvec, file1:string, file2:string 1408 PURPOSE: Returns the common refinement of the cones given in 1409 the list OC which is supposed to contain the orbit cones intersected with Qgamma. The list actiononorbitcones is supposed to contain the symmetry group acting as permutations of on the list of orbit cones in OC. The optional argument can be used to specify one or two strings with file names, where the first file will contain the hashes of the GIT-cones and the second argument the actual cones in their H-representation. 1410 To obtain the whole GIT-fan Qgamma has to be take the cone generated by the columns of Q. 1411 RETURN: a list containing the bigint hashes of the GIT cones. 1412 EXAMPLE: example GITfanSymmetric; shows an example 1413 " 1414 { 1415 actiononorbitcones=apply(actiononorbitcones,permutationToIntvec); 1416 /** 1417 * stores the hashes of all maximal GIT cones computed 1418 */ 1419 list hashesOfCones; 1420 /** 1421 * stores to the corresponding maximal GIT cone in hashesOfCones 1422 * - 0, if guaranteed that all adjacent GIT cones are known 1423 * or are to be computed in the next iteration (*) 1424 * - hash as intvec, otherwise 1425 */ 1426 list workingList; 1427 1428 1429 /** 1430 * compute starting cone 1431 */ 1432 bigintmat w,v; 1433 cone lambda; 1434 intvec lambdaHash; 1435 while(dimension(lambda) <> nrows(Q)){ 1436 w = randConeEl(transpose(Q),100); 1437 dbprint("testing "+string(w)); 1438 if (containsRelatively(Qgamma,w)) { 1439 lambda,lambdaHash = GITcone(OC,w); 1440 dbprint("computed cone of dimension "+string(dimension(lambda))); 1441 } 1442 } 1443 int nCones = 1; // essentially size(hashesOfCones) 1444 int nConesOpen = 1; // essentially size(workingList)+1, see (*) above 1445 1446 1447 /** 1448 * initialize lists 1449 */ 1450 int posToInsert = 1; 1451 bigint hashToInsert = binaryToBigint(lambdaHash); 1452 intvec sigLambdaHash; 1453 for(int i = 2; i <= size(actiononorbitcones); i++){ 1454 sigLambdaHash = composeIntvecs(actiononorbitcones[i],lambdaHash); 1455 hashToInsert = min(hashToInsert,binaryToBigint(sigLambdaHash)); 1456 } 1457 hashesOfCones[1] = hashToInsert; 1458 workingList[1] = int(0); 1459 if (size(#)>0) {write(":w "+#[1],string(hashesOfCones[1]) + ",");} 1460 if (size(#)>1) {write(":w "+#[2],"");writeGitconeToFile(lambda,#[2]);} 1461 1462 1463 /** 1464 * traverse fan 1465 */ 1466 int t,tt; 1467 list FL; 1468 intvec neighbourHash; 1469 int mu = 1024; 1470 while (lambdaHash>0) 1471 { 1472 tt=timer; 1473 1474 /** 1475 * compute all facets of lambda 1476 */ 1477 t = timer; 1478 FL = listOfFacetsAndInteriorVectors(lambda, Qgamma); 1479 dbprint("time for facets: "+string(timer - t)); 1480 1481 /** 1482 * compute all neighbours of lambda 1483 */ 1484 for (i=size(FL[1]); i>0; i--) 1485 { 1486 v = FL[1][i][1]; // interior facet normal 1487 w = FL[1][i][2]; // interior facet point 1488 neighbourHash = getNeighborHash(OC,w,v,mu); 1489 posToInsert,hashToInsert = findPosToInsertSymmetric(hashesOfCones,neighbourHash,actiononorbitcones); 1490 if(posToInsert > 0) 1491 { 1492 if (size(#)>0){write(":a "+#[1],string(binaryToBigint(neighbourHash)) + ",");} 1493 hashesOfCones = insertToList(hashesOfCones,hashToInsert,posToInsert); 1494 workingList = insertToList(workingList,neighbourHash,posToInsert); 1495 nConesOpen++; 1496 nCones++; 1497 } 1498 } 1499 1500 /** 1501 * pick lambdaHash and lambda for next iteration, 1502 * set respective entry in workingList to 0 1503 */ 1504 lambdaHash = 0; 1505 for (i=size(workingList); i>0; i--) 1506 { 1507 if (typeof(workingList[i])=="intvec") 1508 { 1509 lambdaHash = workingList[i]; 1510 lambda = gitConeFromHash(OC,lambdaHash); 1511 if (size(#)>1) {writeGitconeToFile(lambda,#[2]);} 1512 workingList[i] = int(0); 1513 break; 1514 } 1515 } 1516 nConesOpen--; 1517 1518 dbprint("overall: "+string(nCones)+" open: "+string(nConesOpen)+" time for loop: "+string(timer-tt)); 1519 } 1520 return(hashesOfCones); 1521 } 1522 example 1523 { 1524 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 1525 ideal J = 1526 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 1527 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 1528 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 1529 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 1530 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 1531 intmat Q[5][10] = 1532 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1533 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1534 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 1535 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 1536 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 1537 list simplexSymmetryGroup = G25Action(); 1538 list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ), 1539 intvec( 1, 2, 3, 5, 6 ), 1540 intvec( 1, 2, 3, 5, 7 ), 1541 intvec( 1, 2, 3, 5, 10 ), 1542 intvec( 1, 2, 3, 7, 9 ), 1543 intvec( 1, 2, 6, 9, 10 ), 1544 intvec( 1, 2, 3, 4, 5, 6 ), 1545 intvec( 1, 2, 3, 4, 5, 10 ), 1546 intvec( 1, 2, 3, 5, 6, 8 ), 1547 intvec( 1, 2, 3, 5, 6, 9 ), 1548 intvec( 1, 2, 3, 5, 7, 10 ), 1549 intvec( 1, 2, 3, 7, 9, 10 ), 1550 intvec( 1, 2, 3, 4, 5, 6, 7 ), 1551 intvec( 1, 2, 3, 4, 5, 6, 8 ), 1552 intvec( 1, 2, 3, 4, 5, 6, 9 ), 1553 intvec( 1, 2, 3, 5, 6, 9, 10 ), 1554 intvec( 1, 2, 3, 4, 5, 6, 7, 8 ), 1555 intvec( 1, 2, 3, 4, 5, 6, 9, 10 ), 1556 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ), 1557 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ); 1558 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 1559 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 1560 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 1561 apply(afaceOrbits,size); 1562 list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits); 1563 apply(minAfaceOrbits,size); 1564 list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q); 1565 apply(listOfOrbitConeOrbits,size); 1566 list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits); 1567 size(listOfMinimalOrbitConeOrbits); 1568 list Asigma = groupActionOnQImage(simplexSymmetryGroup,Q); 1569 list actionOnOrbitconeIndices = groupActionOnHashes(Asigma,listOfOrbitConeOrbits); 1570 list OClist = listOfOrbitConeOrbits[1]; 1571 for (int i =2;i<=size(listOfOrbitConeOrbits);i++){ 1572 OClist = OClist + listOfOrbitConeOrbits[i]; 1573 } 1574 cone mov = coneViaPoints(transpose(Q)); 1575 mov = canonicalizeCone(mov); 1576 printlevel = 3; 1577 list Sigma = GITfanSymmetric(OClist, Q, mov, actionOnOrbitconeIndices); 1578 Sigma; 1579 } 1580 1581 1582 proc GITfanParallelSymmetric(list OC, bigintmat Q, cone Qgamma, list actiononorbitcones, list #) 1583 "USAGE: GITfanParallelSymmetric(OC, Q, Qgamma, actiononorbitcones [, file1]); OC:list, Q:bigintmat, Qgamma:cone, actiononorbitcones: list of intvec, file1:string 1584 PURPOSE: Returns the common refinement of the cones given in 1585 the list OC which is supposed to contain the orbit cones intersected with Qgamma. The list actiononorbitcones is supposed to contain the symmetry group acting as permutations of on the list of orbit cones in OC. The optional argument can be used to specify a name for a file which will contain the hashes of the GIT-cones. 1586 To obtain the whole GIT-fan Qgamma has to be take the cone generated by the columns of Q. 1587 RETURN: a list containing the bigint hashes of the GIT cones. 1588 NOTE: The proceduce uses parallel computation for the construction of the GIT-cones. 1589 EXAMPLE: example GITfanParallelSymmetric; shows an example 1590 " 1591 { 1592 actiononorbitcones=apply(actiononorbitcones,permutationToIntvec); 1593 /** 1594 * stores the hashes of all maximal GIT cones computed 1595 */ 1596 list hashesOfCones; 1597 /** 1598 * stores to the corresponding maximal GIT cone in hashesOfCones 1599 * - 0, if guaranteed that all adjacent GIT cones are known 1600 * or are to be computed in the next iteration (*) 1601 * - hash as intvec, otherwise 1602 */ 1603 list workingList; 1604 1605 1606 /** 1607 * compute starting cone 1608 */ 1609 bigintmat w,v; 1610 cone lambda; 1611 intvec lambdaHash; 1612 while(dimension(lambda) <> nrows(Q)){ 1613 w = randConeEl(transpose(Q),100); 1614 dbprint("testing "+string(w)); 1615 if (containsRelatively(Qgamma,w)) { 1616 lambda,lambdaHash = GITcone(OC,w); 1617 dbprint("computed cone of dimension "+string(dimension(lambda))); 1618 } 1619 } 1620 int nCones = 1; // essentially size(hashesOfCones) 1621 int nConesOpen = 1; // essentially size(workingList)+1, see (*) above 1622 1623 1624 /** 1625 * initialize lists 1626 */ 1627 bigint hashToInsert = binaryToBigint(lambdaHash); 1628 int posToInsert = 1; 1629 intvec sigLambdaHash; 1630 for(int i = 2; i <= size(actiononorbitcones); i++){ 1631 sigLambdaHash = composeIntvecs(actiononorbitcones[i],lambdaHash); 1632 hashToInsert = min(hashToInsert,binaryToBigint(sigLambdaHash)); 1633 } 1634 hashesOfCones[1] = hashToInsert; 1635 workingList[1] = int(0); 1636 if (size(#)>0) {write(":w "+#[1],string(binaryToBigint(lambdaHash)) + ",");} 1637 //if (size(#)>1) {write(":w "+#[2],"list listOfMaximalCones;");writeGitconeToFile(lambda,#[2]);} 1638 1639 1640 list iterationArgs = list(list(lambdaHash,OC,Qgamma,actiononorbitcones)); 1641 list iterationRes; 1642 1643 /** 1644 * traverse fan 1645 */ 1646 int j,t,tloop; 1647 list FL; 1648 list neighbourHashes; 1649 intvec neighbourHash; 1650 while (size(iterationArgs)>0) 1651 { 1652 tloop=rtimer; 1653 1654 /** 1655 * compute all neighbours of lambda 1656 */ 1657 t = rtimer; 1658 iterationRes = parallelWaitAll("computeNeighbourMinimalHashes",iterationArgs); 1659 dbprint("time neighbours: "+string(rtimer - t)); 1660 1661 /** 1662 * central book keeping 1663 */ 1664 t = rtimer; 1665 for (i=1; i<=size(iterationRes); i++) 1666 { 1667 neighbourHashes = iterationRes[i]; 1668 for (j=1; j<=size(neighbourHashes); j++) 1669 { 1670 neighbourHash = neighbourHashes[j]; 1671 hashToInsert = binaryToBigint(neighbourHash); 1672 posToInsert = findPlaceToInsert(hashesOfCones,hashToInsert); 1673 if(posToInsert > 0) 1674 { 1675 if (size(#)>0){write(":a "+#[1],string(binaryToBigint(neighbourHash)) + ",");} 1676 hashesOfCones = insertToList(hashesOfCones,hashToInsert,posToInsert); 1677 workingList = insertToList(workingList,neighbourHash,posToInsert); 1678 nConesOpen++; 1679 nCones++; 1680 } 1681 } 1682 } 1683 nConesOpen = nConesOpen - size(iterationArgs); 1684 dbprint("time bookkeeping: "+string(rtimer - t)); 1685 1686 /** 1687 * pick arguments for next iteration, 1688 * set respective entry in workingList to 0 1689 */ 1690 t = rtimer; 1691 iterationArgs = list(); 1692 for (i=size(workingList); i>0; i--) 1693 { 1694 if (typeof(workingList[i])=="intvec") 1695 { 1696 iterationArgs[size(iterationArgs)+1] = list(workingList[i],OC,Qgamma,actiononorbitcones); 1697 workingList[i] = int(0); 1698 if (size(iterationArgs) >= getcores()) 1699 { 1700 break; 1701 } 1702 } 1703 } 1704 dbprint("time preparation: "+string(rtimer - t)); 1705 1706 dbprint("overall: "+string(nCones)+" open: "+string(nConesOpen)+" time loop: "+string(rtimer-tloop)); 1707 } 1708 return(hashesOfCones); 1709 } 1710 example 1711 { 1712 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 1713 ideal J = 1714 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 1715 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 1716 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 1717 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 1718 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 1719 intmat Q[5][10] = 1720 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1721 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1722 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 1723 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 1724 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 1725 list simplexSymmetryGroup = G25Action(); 1726 list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ), 1727 intvec( 1, 2, 3, 5, 6 ), 1728 intvec( 1, 2, 3, 5, 7 ), 1729 intvec( 1, 2, 3, 5, 10 ), 1730 intvec( 1, 2, 3, 7, 9 ), 1731 intvec( 1, 2, 6, 9, 10 ), 1732 intvec( 1, 2, 3, 4, 5, 6 ), 1733 intvec( 1, 2, 3, 4, 5, 10 ), 1734 intvec( 1, 2, 3, 5, 6, 8 ), 1735 intvec( 1, 2, 3, 5, 6, 9 ), 1736 intvec( 1, 2, 3, 5, 7, 10 ), 1737 intvec( 1, 2, 3, 7, 9, 10 ), 1738 intvec( 1, 2, 3, 4, 5, 6, 7 ), 1739 intvec( 1, 2, 3, 4, 5, 6, 8 ), 1740 intvec( 1, 2, 3, 4, 5, 6, 9 ), 1741 intvec( 1, 2, 3, 5, 6, 9, 10 ), 1742 intvec( 1, 2, 3, 4, 5, 6, 7, 8 ), 1743 intvec( 1, 2, 3, 4, 5, 6, 9, 10 ), 1744 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ), 1745 intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ); 1746 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 1747 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 1748 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 1749 apply(afaceOrbits,size); 1750 list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits); 1751 apply(minAfaceOrbits,size); 1752 list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q); 1753 apply(listOfOrbitConeOrbits,size); 1754 list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits); 1755 size(listOfMinimalOrbitConeOrbits); 1756 list Asigma = groupActionOnQImage(simplexSymmetryGroup,Q); 1757 list actionOnOrbitconeIndices = groupActionOnHashes(Asigma,listOfOrbitConeOrbits); 1758 list OClist = listOfOrbitConeOrbits[1]; 1759 for (int i =2;i<=size(listOfOrbitConeOrbits);i++){ 1760 OClist = OClist + listOfOrbitConeOrbits[i]; 1761 } 1762 cone mov = coneViaPoints(transpose(Q)); 1763 mov = canonicalizeCone(mov); 1764 list Sigma = GITfanParallelSymmetric(OClist, Q, mov, actionOnOrbitconeIndices); 1765 Sigma; 1766 } 1767 1768 1769 // not static to be used in parallel.lib 1770 proc computeNeighbourMinimalHashes(intvec lambdaHash, list OC, cone Qgamma, list actiononorbitcones) 1771 { 1772 /** 1773 * compute all facets of lambda 1774 */ 1775 cone lambda = gitConeFromHash(OC,lambdaHash); 1776 list FL = listOfFacetsAndInteriorVectors(lambda, Qgamma); 1777 1778 /** 1779 * compute all minimal hashes of neighbours of lambda 1780 */ 1781 int i,j; 1782 bigintmat v,w; 1783 intvec neighbourHash; 1784 intvec neighbourHashPerm; 1785 bigint nPerm; 1786 intvec neighbourHashMin; 1787 bigint nMin; 1788 list neighbourHashes; 1789 for (i=size(FL[1]); i>0; i--) 1790 { 1791 v = FL[1][i][1]; // interior facet normal 1792 w = FL[1][i][2]; // interior facet point 1793 neighbourHash = getNeighborHash(OC,w,v,1024); 1794 neighbourHashMin = neighbourHash; 1795 nMin = binaryToBigint(neighbourHash); 1796 for (j=size(actiononorbitcones); j>1; j--) 1797 { 1798 neighbourHashPerm = composeIntvecs(actiononorbitcones[j],neighbourHash); 1799 nPerm = binaryToBigint(neighbourHashPerm); 1800 if (nPerm < nMin) 1801 { 1802 nMin = nPerm; 1803 neighbourHashMin = neighbourHashPerm; 1804 } 1805 } 1806 neighbourHashes[i] = neighbourHashMin; 1807 } 1808 1809 return (neighbourHashes); 1810 } 1811 1812 1813 1814 static proc writeGitconeToFile(cone lambda,string fn) 1815 { 1816 bigintmat H = facets(lambda); 1817 1818 int rows = nrows(H); 1819 int cols = ncols(H); 1820 string toBeWritten = "bigintmat H["+string(rows)+"]["+string(cols)+"]="; 1821 int i,j; 1822 for (i=1; i<=rows; i++) 1823 { 1824 for (j=1; j<=cols; j++) 1825 { 1826 toBeWritten = toBeWritten + string(H[i,j]) + ","; 1827 } 1828 } 1829 toBeWritten = toBeWritten[1..size(toBeWritten)-1]; 1830 toBeWritten = toBeWritten + ";"; 1831 write(":a "+fn,toBeWritten); 1832 1833 toBeWritten = "listOfMaximalCones[size(listOfMaximalCones)+1] = coneViaInequalities(V);"; 1834 write(":a "+fn,toBeWritten); 1835 1836 toBeWritten = "kill V;"; 1837 write(":a "+fn,toBeWritten); 1838 1839 toBeWritten = ""; 1840 write(":a "+fn,toBeWritten); 1841 } 1842 1843 1844 1845 static proc listOfFacetsAndInteriorVectors(cone lambda, cone Qgamma, list #){ 1846 list FL; 1847 int numboundary; 1848 1849 bigintmat FL0 = facets(lambda); 1850 bigintmat H[1][ncols(FL0)]; 1851 bigintmat FL1[nrows(FL0)-1][ncols(FL0)]; // delete row of H from FL0 1852 1853 bigintmat w; 1854 cone facetCone; 1855 for(int i = 1; i <= nrows(FL0); i++){ 1856 H = FL0[i,1..ncols(FL0)]; 1857 1858 if(i > 1 and i < nrows(FL0)){ 1859 FL1 = FL0[1..i-1, 1..ncols(FL0)], FL0[i+1..nrows(FL0), 1..ncols(FL0)]; 1860 } else { 1861 if(i == nrows(FL0)){ 1862 FL1 = FL0[1..i-1, 1..ncols(FL0)]; 1863 } else { // i = 1: 1864 FL1 = FL0[2..nrows(FL0), 1..ncols(FL0)]; 1865 } 1866 } 1867 1868 facetCone = coneViaInequalities(FL1, H); 1869 1870 if (size(#)>0) 1871 { 1872 if (containsInSupport(facetCone,#[1])) 1873 { 1874 i++; 1875 continue; 1876 } 1877 } 1878 1879 w = relativeInteriorPoint(facetCone); 1880 1881 if(containsRelatively(Qgamma,w)){ 1882 FL[size(FL) + 1] = list(H,w); 1883 } else { 1884 numboundary=numboundary+1; 1885 } 1886 } 1887 1888 return(list(FL,numboundary)); 1889 } 1890 1891 1892 1893 1894 1895 //////////////////////////////// 1896 1897 // CAREFUL: we assume that actiononorbitcones[1] = id. 1898 static proc findPosToInsertSymmetric(list hashesOfCones, intvec J, list actiononorbitcones){ 1899 1900 // 1.: compute minimal hash sigJ 1901 bigint n = binaryToBigint(J); 1902 intvec sigJ; 1903 for(int i = 2; i <= size(actiononorbitcones); i++){ 1904 sigJ = composeIntvecs(actiononorbitcones[i],J); 1905 n = min(n,binaryToBigint(sigJ)); 1906 } 1907 1908 // 2.:check if minimal hash is already in list 1909 return(findPlaceToInsert(hashesOfCones,n),n); 1910 } 1911 1912 ////////////////////// 1913 1914 // insert n at position pos and move bigger elements by one index 1915 proc insertToList(list L, def elementToInsert, int posToInsert){ 1916 if(posToInsert == 1){ 1917 return(list(elementToInsert)+L); 1918 } 1919 if(posToInsert == size(L)+1){ 1920 return(L+list(elementToInsert)); 1921 } 1922 return(list(L[1..(posToInsert-1)],elementToInsert,L[posToInsert..size(L)])); 1923 } 1924 1925 proc GITcone(list OCcones, bigintmat w) 1926 "USAGE: GITcone(OCcones, w); OCcones: list of orbit cones, w: bigintmat with one row 1927 PURPOSE: Returns the intersection of all orbit cones containing w. 1928 RETURN: cone,intvec with the GIT cone containing w, and the hash of this cone (the indices of the orbit cones contributing to the intersection) 1929 EXAMPLE: example GITcone; shows an example 1930 " 1931 { 1932 list OCrelcones; 1933 intvec Hash; 1934 for(int i = 1; i <= size(OCcones); i++){ 1935 if(containsInSupport(OCcones[i], w)){ 1936 OCrelcones[size(OCrelcones)+1]=OCcones[i]; 1937 Hash[size(Hash)+1] = i; 1938 } 1939 } 1940 Hash = intvec(Hash[2..size(Hash)]); 1941 return(convexIntersection(OCrelcones),Hash); 1942 } 1943 example { 1944 echo=2; 1945 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1); 1946 ideal J = 1947 T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 1948 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 1949 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 1950 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 1951 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 1952 intmat Q[5][10] = 1953 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1954 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1955 0, 1, 1, 0, 0, 0, -1, 1, 0, 0, 1956 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 1957 0, 0, 1, 1, -1, 0, 0, 0, 0, 1; 1958 1959 list AF= afaces(J,nrows(Q)); 1960 AF=fullDimImages(AF,Q); 1961 AF = minimalAfaces(AF); 1962 list OC = orbitCones(AF,Q); 1963 bigintmat w[1][nrows(Q)]; 1964 for(int i = 1; i <= nrows(Q); i++){ 1965 w = w + Q[1..nrows(Q),i]; 1966 } 1967 GITcone(OC,w); 1968 } 1969 1970 1971 static proc gitConeFromHash(list OC, intvec Hash) 1972 { 1973 list OCtoIntersect = OC[Hash]; 1974 return (convexIntersection(OCtoIntersect)); 1975 } 1976 1977 1978 1979 1980 1981 static proc randConeEl(bigintmat Q, int bound){ 1982 bigintmat w[1][ncols(Q)]; 1983 1984 for(int i = 1; i <= nrows(Q); i++){ 1985 bigintmat v[1][ncols(Q)] = Q[i,1..ncols(Q)]; 1986 1987 w = w + random(1,bound) * v; 1988 1989 kill v; 1990 } 1991 1992 return(w); 1993 } 1994 1995 1996 1997 static proc subdividelist(list OC,int ncores){ 1998 if (ncores>size(OC)){ncores=size(OC);} 1999 int percore = size(OC) div (ncores); 2000 int lastcore = size(OC) mod (ncores); 2001 list OCLL; 2002 int j=1; 2003 int starti=1; 2004 int endi; 2005 for (int i=1;i<=ncores;i++){ 2006 endi = starti+percore-1; 2007 if (i<=lastcore) {endi=endi+1;} 2008 list OCLLi=OC[starti..endi]; 2009 starti=endi+1; 2010 OCLL[i]=OCLLi; 2011 kill OCLLi; 2012 } 2013 return(OCLL); 2014 } 2015 //list OC = 1,2,3,4,5,6,7,8,9,10,11,12; 2016 //subdividelist(OC,4); 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029 proc orbitCones(list AF, intmat Q, list #) 2030 "USAGE: orbitCones(AF, Q[, d]); AF: list of intvecs, Q: intmat, d: int 2031 PURPOSE: Returns the list consisting of all cones Q(gam0) where gam0 in AF. If the optional argument d is given then the function returns only the orbit cones of dimension at least d 2032 RETURN: a list of cones 313 proc orbitCones(ideal a, bigintmat Q, list #) 314 "USAGE: orbitCones(a, Q, b, c); a: ideal, Q: bigintmat, b: int, c: int 315 PURPOSE: Returns a list consisting of all cones Q(gam0) where gam0 is an a-face. 316 Moreover, it is possible to specify a dimensional bound b, 317 upon which only cones of that dimension and above are returned. 318 Lastly, as the computation is parallizable, one can specify c, 319 the number of cores to be used by the computation. 320 RETURN: a list of cones 2033 321 EXAMPLE: example orbitCones; shows an example 2034 322 " 2035 323 { 324 list AF; 325 326 if((size(#) > 1) and (typeof(#[2]) == "int")) 327 { 328 AF = afaces(a, nrows(Q), #[2]); 329 } 330 else 331 { 332 AF = afaces(a, nrows(Q)); 333 } 334 335 int dimensionBound = 0; 336 if((size(#) > 0) and (typeof(#[1]) == "int")) 337 { 338 dimensionBound = #[1]; 339 } 340 2036 341 list OC; 2037 342 intvec gam0; 2038 343 int j; 2039 cone c; 2040 int d = 0; 2041 if (size(#)>0) { 2042 d = #[1]; 2043 } 2044 2045 for(int i = 1; i <= size(AF); i++){ 344 345 for(int i = 1; i <= size(AF); i++) 346 { 2046 347 gam0 = AF[i]; 2047 348 2048 if(gam0 == 0){ 349 if(gam0 == 0) 350 { 2049 351 bigintmat M[1][nrows(Q)]; 2050 } else { 352 } 353 else 354 { 2051 355 bigintmat M[size(gam0)][nrows(Q)]; 2052 2053 for (j = 1; j <= size(gam0); j++){356 for (j = 1; j <= size(gam0); j++) 357 { 2054 358 M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]]; 2055 359 } 2056 360 } 2057 2058 c = coneViaPoints(M); 2059 2060 if(listContainsCone(OC, c) == 0 && dimension(c)>=d){361 cone c = coneViaPoints(M); 362 363 if((dimension(c) >= dimensionBound) and (!(listContainsCone(OC, c)))) 364 { 2061 365 OC[size(OC)+1] = c; 2062 366 } 2063 367 2064 kill M ;368 kill M, c; 2065 369 } 2066 370 2067 371 return(OC); 2068 }2069 example2070 {2071 echo=2;2072 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);2073 ideal J =2074 T(5)*T(10)-T(6)*T(9)+T(7)*T(8),2075 T(1)*T(9)-T(2)*T(7)+T(4)*T(5),2076 T(1)*T(8)-T(2)*T(6)+T(3)*T(5),2077 T(1)*T(10)-T(3)*T(7)+T(4)*T(6),2078 T(2)*T(10)-T(3)*T(9)+T(4)*T(8);2079 intmat Q[5][10] =2080 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,2081 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,2082 0, 1, 1, 0, 0, 0, -1, 1, 0, 0,2083 0, 1, 0, 1, 0, -1, 0, 0, 1, 0,2084 0, 0, 1, 1, -1, 0, 0, 0, 0, 1;2085 2086 list AF= afaces(J);2087 print(size(AF));2088 list OC = orbitCones(AF,Q);2089 size(OC);2090 }2091 2092 2093 2094 ///////////////////////////////////////2095 2096 proc GITfanFromOrbitCones(list OC, bigintmat Q, cone Qgamma, list #)2097 "USAGE: GITfanFromOrbitCones(OC, Q, Qgamma [, file1, file2]); OC:list, Q:bigintmat, Qgamma:cone, file1:string, file2:string2098 PURPOSE: Returns the common refinement of the cones given in2099 the list OC which is supposed to contain the orbit cones intersected with Qgamma. The optional argument can be used2100 to specify one or two strings with file names, where the first file will contain the hashes of the2101 GIT-cones and the second argument the actual cones in their H-representation.2102 To obtain the whole GIT-fan Qgamma has to be take the cone generated by the columns of Q.2103 RETURN: a list containing the bigint hashes of the GIT cones.2104 EXAMPLE: example GITfanFromOrbitCones; shows an example2105 "2106 {2107 /**2108 * stores the hashes of all maximal GIT cones computed2109 */2110 list hashesOfCones;2111 /**2112 * stores to the corresponding maximal GIT cone in hashesOfCones2113 * - 0, if guaranteed that all adjacent GIT cones are known2114 * or are to be computed in the next iteration (*)2115 * - hash as intvec, otherwise2116 */2117 list workingList;2118 2119 2120 /**2121 * compute starting cone2122 */2123 bigintmat w,v;2124 cone lambda;2125 intvec lambdaHash;2126 while(dimension(lambda) <> nrows(Q)){2127 w = randConeEl(transpose(Q),100);2128 dbprint("testing "+string(w));2129 if (containsRelatively(Qgamma,w)) {2130 lambda,lambdaHash = GITcone(OC,w);2131 dbprint("computed cone of dimension "+string(dimension(lambda)));2132 }2133 }2134 int nCones = 1; // essentially size(hashesOfCones)2135 int nConesOpen = 1; // essentially size(workingList)+1, see (*) above2136 2137 2138 /**2139 * initialize lists2140 */2141 int posToInsert = 1;2142 bigint hashToInsert = binaryToBigint(lambdaHash);2143 hashesOfCones[1] = hashToInsert;2144 workingList[1] = int(0);2145 if (size(#)>0) {write(":w "+#[1],string(hashesOfCones[1]) + ",");}2146 if (size(#)>1) {write(":w "+#[2],"list listOfMaximalCones;");writeGitconeToFile(lambda,#[2]);}2147 2148 2149 /**2150 * traverse fan2151 */2152 int i,t,tt;2153 list FL;2154 intvec neighbourHash;2155 int mu = 1024;2156 while (lambdaHash>0)2157 {2158 tt=timer;2159 2160 /**2161 * compute all facets of lambda2162 */2163 t = timer;2164 FL = listOfFacetsAndInteriorVectors(lambda, Qgamma);2165 dbprint("time for facets: "+string(timer - t));2166 2167 /**2168 * compute all neighbours of lambda2169 */2170 for (i=size(FL[1]); i>0; i--)2171 {2172 v = FL[1][i][1]; // interior facet normal2173 w = FL[1][i][2]; // interior facet point2174 neighbourHash = getNeighborHash(OC,w,v,mu);2175 hashToInsert = binaryToBigint(neighbourHash);2176 posToInsert = findPlaceToInsert(hashesOfCones,hashToInsert);2177 2178 if(posToInsert > 0)2179 {2180 if (size(#)>0){write(":a "+#[1],string(binaryToBigint(neighbourHash)) + ",");}2181 hashesOfCones = insertToList(hashesOfCones,hashToInsert,posToInsert);2182 workingList = insertToList(workingList,neighbourHash,posToInsert);2183 nConesOpen++;2184 nCones++;2185 }2186 }2187 2188 /**2189 * pick lambdaHash and lambda for next iteration,2190 * set respective entry in workingList to 02191 */2192 lambdaHash = 0;2193 for (i=size(workingList); i>0; i--)2194 {2195 if (typeof(workingList[i])=="intvec")2196 {2197 lambdaHash = workingList[i];2198 lambda = gitConeFromHash(OC,lambdaHash);2199 if (size(#)>1) {writeGitconeToFile(lambda,#[2]);}2200 workingList[i] = int(0);2201 break;2202 }2203 }2204 nConesOpen--;2205 2206 dbprint("overall: "+string(nCones)+" open: "+string(nConesOpen)+" time for loop: "+string(timer-tt));2207 }2208 return(hashesOfCones);2209 }2210 example2211 {2212 echo=2;2213 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);2214 ideal J =2215 T(5)*T(10)-T(6)*T(9)+T(7)*T(8),2216 T(1)*T(9)-T(2)*T(7)+T(4)*T(5),2217 T(1)*T(8)-T(2)*T(6)+T(3)*T(5),2218 T(1)*T(10)-T(3)*T(7)+T(4)*T(6),2219 T(2)*T(10)-T(3)*T(9)+T(4)*T(8);2220 intmat Q[5][10] =2221 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,2222 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,2223 0, 1, 1, 0, 0, 0, -1, 1, 0, 0,2224 0, 1, 0, 1, 0, -1, 0, 0, 1, 0,2225 0, 0, 1, 1, -1, 0, 0, 0, 0, 1;2226 2227 list AF= afaces(J,nrows(Q));2228 AF=fullDimImages(AF,Q);2229 AF = minimalAfaces(AF);2230 list OC = orbitCones(AF,Q);2231 cone Qgamma = coneViaPoints(transpose(Q));2232 list GIT = GITfanFromOrbitCones(OC,Q,Qgamma);2233 size(GIT);2234 }2235 2236 2237 2238 2239 2240 proc GITfanParallel(list OC, intmat Q, cone Qgamma, list #)2241 "USAGE: GITfanParallel(OC, Q, Qgamma [, file1]); OC:list, Q:intmat, Qgamma:cone, file1:string2242 PURPOSE: Returns the common refinement of the cones given in2243 the list OC which is supposed to contain the orbit cones intersected with Qgamma. The optional argument can be used to specify a name for a file which will contain the hashes of the GIT-cones.2244 To obtain the whole GIT-fan Qgamma has to be take the cone generated by the columns of Q.2245 RETURN: a list containing the bigint hashes of the GIT cones.2246 NOTE: The proceduce uses parallel computation for the construction of the GIT-cones.2247 EXAMPLE: example GITfanParallel; shows an example2248 "2249 {2250 /**2251 * stores the hashes of all maximal GIT cones computed2252 */2253 list hashesOfCones;2254 /**2255 * stores to the corresponding maximal GIT cone in hashesOfCones2256 * - 0, if guaranteed that all adjacent GIT cones are known2257 * or are to be computed in the next iteration (*)2258 * - hash as intvec, otherwise2259 */2260 list workingList;2261 2262 2263 /**2264 * compute starting cone2265 */2266 bigintmat w,v;2267 cone lambda;2268 intvec lambdaHash;2269 while(dimension(lambda) <> nrows(Q)){2270 w = randConeEl(transpose(Q),100);2271 dbprint("testing "+string(w));2272 if (containsRelatively(Qgamma,w)) {2273 lambda,lambdaHash = GITcone(OC,w);2274 dbprint("computed cone of dimension "+string(dimension(lambda)));2275 }2276 }2277 int nCones = 1; // essentially size(hashesOfCones)2278 int nConesOpen = 1; // essentially size(workingList)+1, see (*) above2279 2280 2281 /**2282 * initialize lists2283 */2284 bigint hashToInsert = binaryToBigint(lambdaHash);2285 int posToInsert = 1;2286 hashesOfCones[1] = hashToInsert;2287 workingList[1] = int(0);2288 if (size(#)>0) {write(":w "+#[1],string(binaryToBigint(lambdaHash)) + ",");}2289 //if (size(#)>1) {write(":w "+#[2],"list listOfMaximalCones;");writeGitconeToFile(lambda,#[2]);}2290 2291 2292 list iterationArgs = list(list(lambdaHash,OC,Qgamma));2293 list iterationRes;2294 2295 /**2296 * traverse fan2297 */2298 int i,j,t,tloop;2299 list FL;2300 list neighbourHashes;2301 intvec neighbourHash;2302 while (size(iterationArgs)>0)2303 {2304 tloop=rtimer;2305 2306 /**2307 * compute all neighbours of lambda2308 */2309 t = rtimer;2310 iterationRes = parallelWaitAll("computeNeighbourHashes",iterationArgs);2311 dbprint("time neighbours: "+string(rtimer - t));2312 2313 /**2314 * central book keeping2315 */2316 t = rtimer;2317 for (i=1; i<=size(iterationRes); i++)2318 {2319 neighbourHashes = iterationRes[i];2320 for (j=1; j<=size(neighbourHashes); j++)2321 {2322 neighbourHash = neighbourHashes[j];2323 hashToInsert = binaryToBigint(neighbourHash);2324 posToInsert = findPlaceToInsert(hashesOfCones,hashToInsert);2325 if(posToInsert > 0)2326 {2327 if (size(#)>0){write(":a "+#[1],string(binaryToBigint(neighbourHash)) + ",");}2328 hashesOfCones = insertToList(hashesOfCones,hashToInsert,posToInsert);2329 workingList = insertToList(workingList,neighbourHash,posToInsert);2330 nConesOpen++;2331 nCones++;2332 }2333 }2334 }2335 nConesOpen = nConesOpen - size(iterationArgs);2336 dbprint("time bookkeeping: "+string(rtimer - t));2337 2338 /**2339 * pick arguments for next iteration,2340 * set respective entry in workingList to 02341 */2342 t = rtimer;2343 iterationArgs = list();2344 for (i=size(workingList); i>0; i--)2345 {2346 if (typeof(workingList[i])=="intvec")2347 {2348 iterationArgs[size(iterationArgs)+1] = list(workingList[i],OC,Qgamma);2349 workingList[i] = int(0);2350 if (size(iterationArgs) >= getcores())2351 {2352 break;2353 }2354 }2355 }2356 dbprint("time preparation: "+string(rtimer - t));2357 2358 dbprint("overall: "+string(nCones)+" open: "+string(nConesOpen)+" time loop: "+string(rtimer-tloop));2359 }2360 2361 if (size(#)==0)2362 {2363 return(hashesOfCones);2364 }2365 }2366 example2367 {2368 echo=2;2369 setcores(4);2370 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);2371 ideal J =2372 T(5)*T(10)-T(6)*T(9)+T(7)*T(8),2373 T(1)*T(9)-T(2)*T(7)+T(4)*T(5),2374 T(1)*T(8)-T(2)*T(6)+T(3)*T(5),2375 T(1)*T(10)-T(3)*T(7)+T(4)*T(6),2376 T(2)*T(10)-T(3)*T(9)+T(4)*T(8);2377 intmat Q[5][10] =2378 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,2379 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,2380 0, 1, 1, 0, 0, 0, -1, 1, 0, 0,2381 0, 1, 0, 1, 0, -1, 0, 0, 1, 0,2382 0, 0, 1, 1, -1, 0, 0, 0, 0, 1;2383 2384 list AF= afaces(J);2385 print(size(AF));2386 list OC = orbitCones(AF,Q);2387 cone Qgamma = coneViaPoints(transpose(Q));2388 list GIT = GITfanParallel(OC,Q,Qgamma);2389 size(GIT);2390 }2391 2392 // not static to be used in parallel.lib2393 proc computeNeighbourHashes(intvec lambdaHash, list OC, cone Qgamma)2394 {2395 /**2396 * compute all facets of lambda2397 */2398 cone lambda = gitConeFromHash(OC,lambdaHash);2399 list FL = listOfFacetsAndInteriorVectors(lambda, Qgamma);2400 2401 /**2402 * compute all minimal hashes of neighbours of lambda2403 */2404 bigintmat v,w;2405 intvec neighbourHash;2406 list neighbourHashes;2407 for (int i=size(FL[1]); i>0; i--)2408 {2409 v = FL[1][i][1]; // interior facet normal2410 w = FL[1][i][2]; // interior facet point2411 neighbourHash = getNeighborHash(OC,w,v,1024);2412 neighbourHashes[i] = neighbourHash;2413 }2414 2415 return (neighbourHashes);2416 }2417 2418 2419 2420 proc minimalAfaces(list listOfAfaces)2421 "USAGE: minimalAfaces(listOfAfaces); listOfAfaces: list2422 PURPOSE: Returns a list of all minimal a-faces. Note that listOfAfaces must only contain2423 afaces which project to full dimension.2424 RETURN: a list of intvecs2425 EXAMPLE: example minimalAfaces; shows an example2426 "2427 {2428 int i,j;2429 for (i=1; i<=size(listOfAfaces); i++)2430 {2431 for (j=1; j<=size(listOfAfaces); j++)2432 {2433 if (i!=j)2434 {2435 if (isSubset(listOfAfaces[i],listOfAfaces[j]))2436 {2437 listOfAfaces = delete(listOfAfaces,j);2438 if (j<i)2439 {2440 i = i-1;2441 }2442 j = j-1;2443 }2444 }2445 }2446 }2447 return(listOfAfaces);2448 }2449 example2450 {2451 echo=2;2452 setcores(4);2453 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);2454 ideal J =2455 T(5)*T(10)-T(6)*T(9)+T(7)*T(8),2456 T(1)*T(9)-T(2)*T(7)+T(4)*T(5),2457 T(1)*T(8)-T(2)*T(6)+T(3)*T(5),2458 T(1)*T(10)-T(3)*T(7)+T(4)*T(6),2459 T(2)*T(10)-T(3)*T(9)+T(4)*T(8);2460 intmat Q[5][10] =2461 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,2462 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,2463 0, 1, 1, 0, 0, 0, -1, 1, 0, 0,2464 0, 1, 0, 1, 0, -1, 0, 0, 1, 0,2465 0, 0, 1, 1, -1, 0, 0, 0, 0, 1;2466 list AF= afaces(J,nrows(Q));2467 size(AF);2468 AF=fullDimImages(AF,Q);2469 size(AF);2470 AF=minimalAfaces(AF);2471 size(AF);2472 }2473 2474 2475 2476 proc GITfan(ideal J, intmat Q, list #)2477 "USAGE: GITfan(J,Q [, G]); J:ideal, Q:intmat, G:list2478 PURPOSE: Computes the GIT fan associated to J and Q. Optionally a symmetry group action on the column space of Q can be specified.2479 RETURN: a fan, the GIT fan.2480 NOTE: The proceduce uses parallel computation for the construction of the GIT-cones. The a-faces are not computed in parallel. This can be done by calling the aface procedure specifying a list of simplex faces. If used with the optional argument G, the orbit decomposition of the simplex of columns of Q is computed. Refer to the Singular documentation on how to do this more efficiently using GAP.2481 EXAMPLE: example GITfan; shows an example2482 "2483 {2484 list GIT;2485 list OC;2486 if (size(#)>0){2487 (GIT,OC) = GITfanWrapperWithSymmetry(J,Q,#);2488 } else {2489 list AF= afaces(J,nrows(Q));2490 AF=fullDimImages(AF,Q);2491 AF = minimalAfaces(AF);2492 OC = orbitCones(AF,Q);2493 cone Qgamma = coneViaPoints(transpose(Q));2494 GIT = GITfanParallel(OC,Q,Qgamma);2495 }2496 fan Sigma = hashesToFan(GIT,OC);2497 return(Sigma);2498 }2499 example2500 {2501 echo=2;2502 setcores(4);2503 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);2504 ideal J =2505 T(5)*T(10)-T(6)*T(9)+T(7)*T(8),2506 T(1)*T(9)-T(2)*T(7)+T(4)*T(5),2507 T(1)*T(8)-T(2)*T(6)+T(3)*T(5),2508 T(1)*T(10)-T(3)*T(7)+T(4)*T(6),2509 T(2)*T(10)-T(3)*T(9)+T(4)*T(8);2510 intmat Q[5][10] =2511 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,2512 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,2513 0, 1, 1, 0, 0, 0, -1, 1, 0, 0,2514 0, 1, 0, 1, 0, -1, 0, 0, 1, 0,2515 0, 0, 1, 1, -1, 0, 0, 0, 0, 1;2516 fan GIT = GITfan(J,Q);2517 2518 intmat Q[5][10] =2519 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,2520 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,2521 0, 1, 1, 0, 0, 0, -1, 1, 0, 0,2522 0, 1, 0, 1, 0, -1, 0, 0, 1, 0,2523 0, 0, 1, 1, -1, 0, 0, 0, 0, 1;2524 list simplexSymmetryGroup = G25Action();2525 fan GIT2 = GITfan(J,Q,simplexSymmetryGroup);2526 2527 }2528 2529 2530 proc hashToCone(bigint v, list OC)2531 "USAGE: hashToCone(v, OC): v bigint, OC list of cones.2532 ASSUME: the elements of OC are the orbit cones used in the hash representation of the GIT cones.2533 RETURN: a cone, the intersection of the cones in OC according to the binary representation of the hash v.2534 EXAMPLE: example hashToCone; shows an example2535 "2536 {2537 intvec J = bigintToBinary(v, size(OC));2538 return(convexIntersection(list(OC[J])));2539 }2540 example2541 {2542 echo = 2;2543 setcores(4);2544 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);2545 ideal J =2546 T(5)*T(10)-T(6)*T(9)+T(7)*T(8),2547 T(1)*T(9)-T(2)*T(7)+T(4)*T(5),2548 T(1)*T(8)-T(2)*T(6)+T(3)*T(5),2549 T(1)*T(10)-T(3)*T(7)+T(4)*T(6),2550 T(2)*T(10)-T(3)*T(9)+T(4)*T(8);2551 intmat Q[5][10] =2552 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,2553 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,2554 0, 1, 1, 0, 0, 0, -1, 1, 0, 0,2555 0, 1, 0, 1, 0, -1, 0, 0, 1, 0,2556 0, 0, 1, 1, -1, 0, 0, 0, 0, 1;2557 list AF= afaces(J,nrows(Q));2558 AF=fullDimImages(AF,Q);2559 AF = minimalAfaces(AF);2560 list OC = orbitCones(AF,Q);2561 bigint v = 21300544;2562 hashToCone(v, OC);2563 }2564 2565 2566 proc bigintToBinary(bigint n, int r)2567 "2568 USAGE: bigintToBinary(n, r): n bigint, r int.2569 ASSUME: n is smaller then 2^r.2570 RETURN: an intvec, with entries the positions of 1 in the binary representation of n with r bits.2571 EXAMPLE: example bigintToBinary; shows an example2572 "2573 {2574 int k = r-1;2575 2576 intvec v;2577 bigint n0 = n;2578 2579 while(n0 > 0){2580 bigint tmp = bigint(2)^k;2581 2582 while(tmp > n0){2583 k--;2584 tmp = bigint(2)^k;2585 }2586 2587 v = k+1,v;2588 n0 = n0 - tmp;2589 k--;2590 2591 kill tmp;2592 }2593 2594 v = v[1..size(v)-1];2595 return(v);2596 }2597 example2598 {2599 echo = 2;2600 bigintToBinary(bigint(2)^90-1, 90);2601 }2602 2603 2604 2605 2606 proc hashesToFan(list hashes, list OC)2607 "USAGE: hashesToFan(hashes, OC): hashes list of bigint, OC list of cones.2608 ASSUME: the elements of OC are the orbit cones used in the hash representation of the GIT cones.2609 RETURN: a fan, with maximal cones the intersections of the cones in OC according to the binary representation of the hashes.2610 EXAMPLE: example hashesToFan; shows an example2611 "2612 {2613 fan Sigma = emptyFan(ambientDimension(OC[1]));2614 for (int i=1;i<=size(hashes);i++){2615 insertCone(Sigma,hashToCone(hashes[i],OC),0);2616 }2617 return(Sigma);2618 }2619 example2620 {2621 echo = 2;2622 setcores(4);2623 ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);2624 ideal J =2625 T(5)*T(10)-T(6)*T(9)+T(7)*T(8),2626 T(1)*T(9)-T(2)*T(7)+T(4)*T(5),2627 T(1)*T(8)-T(2)*T(6)+T(3)*T(5),2628 T(1)*T(10)-T(3)*T(7)+T(4)*T(6),2629 T(2)*T(10)-T(3)*T(9)+T(4)*T(8);2630 intmat Q[5][10] =2631 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,2632 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,2633 0, 1, 1, 0, 0, 0, -1, 1, 0, 0,2634 0, 1, 0, 1, 0, -1, 0, 0, 1, 0,2635 0, 0, 1, 1, -1, 0, 0, 0, 0, 1;2636 list AF= afaces(J,nrows(Q));2637 AF=fullDimImages(AF,Q);2638 AF = minimalAfaces(AF);2639 list OC = orbitCones(AF,Q);2640 cone Qgamma = coneViaPoints(transpose(Q));2641 list GIT = GITfanParallel(OC,Q,Qgamma);2642 fan Sigma = hashesToFan(GIT,OC);2643 }2644 2645 2646 2647 2648 /////////////////////////////////////2649 2650 proc gkzFan(intmat Q)2651 "USAGE: gkzFan(Q); a: ideal, Q:intmat2652 PURPOSE: Returns the GKZ-fan of the matrix Q.2653 RETURN: a fan.2654 EXAMPLE: example gkzFan; shows an example2655 "2656 {2657 // only difference to gitFan:2658 // it suffices to consider all faces2659 // that are simplicial:2660 list OC = simplicialToricOrbitCones(Q);2661 print(size(OC));2662 cone Qgamma = coneViaPoints(transpose(Q));2663 list GIT = GITfanParallel(OC,Q,Qgamma);2664 fan Sigma = hashesToFan(GIT,OC);2665 return(Sigma);2666 372 } 2667 373 example … … 2673 379 0,0,1,1; 2674 380 2675 gkzFan(Q); 2676 } 2677 2678 381 ring R = 0,T(1..4),dp; 382 ideal a = 0; 383 384 orbitCones(a, Q); 385 } 386 387 /////////////////////////////////////// 388 389 proc gitCone(ideal a, bigintmat Q, bigintmat w) 390 "USAGE: gitCone(a, Q, w); a: ideal, Q:bigintmat, w:bigintmat 391 PURPOSE: Returns the GIT-cone lambda(w), i.e. the intersection of all 392 orbit cones containing the vector w. 393 NOTE: call this only if you are interested in a single GIT-cone. 394 RETURN: a cone. 395 EXAMPLE: example gitCone; shows an example 396 " 397 { 398 list OC = orbitCones(a, Q); 399 cone lambda = nrows(Q); 400 401 for(int i = 1; i <= size(OC); i++) 402 { 403 cone c = OC[i]; 404 405 if(containsInSupport(c, w)) 406 { 407 lambda = convexIntersection(lambda, c); 408 } 409 410 kill c; 411 } 412 413 return(lambda); 414 } 415 example 416 { 417 echo=2; 418 intmat Q[3][4] = 419 1,0,1,0, 420 0,1,0,1, 421 0,0,1,1; 422 423 ring R = 0,T(1..4),dp; 424 ideal a = 0; 425 426 bigintmat w[1][3] = 3,3,1; 427 cone lambda = gitCone(a, Q, w); 428 rays(lambda); 429 430 bigintmat w2[1][3] = 1,1,1; 431 cone lambda2 = gitCone(a, Q, w2); 432 rays(lambda2); 433 } 2679 434 2680 435 ///////////////////////////////////// 2681 436 2682 // Computes all simplicial orbit cones 2683 // w.r.t. the 0-ideal: 2684 static proc simplicialToricOrbitCones(bigintmat Q){ 2685 intvec gam0; 2686 list OC; 2687 cone c; 2688 int r = ncols(Q); 2689 int j; 2690 2691 for(int i = 1; i < 2^r; i++ ){ 2692 gam0 = int2face(i,r); 2693 2694 // each simplicial cone is generated by 2695 // exactly nrows(Q) many columns of Q: 2696 if(size(gam0) == nrows(Q)){ 2697 bigintmat M[size(gam0)][nrows(Q)]; 2698 2699 for(j = 1; j <= size(gam0); j++){ 2700 M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]]; 2701 } 2702 2703 c = coneViaPoints(M); 2704 2705 if((dimension(c) == nrows(Q)) and (!(listContainsCone(OC, c)))){ 2706 OC[size(OC)+1] = c; 2707 } 2708 2709 kill M; 2710 } 2711 } 2712 2713 return(OC); 437 proc gitFan(ideal a, bigintmat Q, list #) 438 "USAGE: gitFan(a, Q); a: ideal, Q:bigintmat 439 PURPOSE: Returns the GIT-fan of the H-action defined by Q on V(a). 440 An optional third parameter of type 'int' is interpreted as the 441 number of CPU-cores you would like to use. 442 Note that 'system("--cpus");' returns the number of cpu available 443 in your system. 444 RETURN: a fan. 445 EXAMPLE: example gitFan; shows an example 446 " 447 { 448 list OC = orbitCones(a, Q, #); 449 450 fan f = refineCones(OC, Q); 451 return(f); 2714 452 } 2715 453 example 2716 454 { 2717 echo = 2; 2718 bigintmat Q[3][4] = 2719 1,0,1,0, 2720 0,1,0,1, 2721 0,0,1,1; 2722 2723 list OC = simplicialToricOrbitCones(Q); 2724 print(OC); 455 echo=2; 456 intmat Q[3][4] = 457 1,0,1,0, 458 0,1,0,1, 459 0,0,1,1; 460 461 ring R = 0,T(1..4),dp; 462 ideal a = 0; 463 464 gitFan(a, Q); 465 466 // 2nd example // 467 kill Q; 468 intmat Q[3][6] = 469 1,1,0,0,-1,-1, 470 0,1,1,-1,-1,0, 471 1,1,1,1,1,1; 472 473 ring R = 0,T(1..6),dp; 474 ideal a = T(1)*T(6) + T(2)*T(5) + T(3)*T(4); 475 476 int t = rtimer; 477 fan F = gitFan(a, Q); 478 t = rtimer - t; 479 480 int tt = rtimer; 481 fan F = gitFan(a, Q, 4); 482 tt = rtimer - tt; 483 484 t; 485 tt; 486 "--------"; 487 kill R, Q, t, tt; 488 // next example // 489 ring R = 0,T(1..10),dp; 490 ideal a = T(5)*T(10)-T(6)*T(9)+T(7)*T(8), 491 T(1)*T(9)-T(2)*T(7)+T(4)*T(5), 492 T(1)*T(8)-T(2)*T(6)+T(3)*T(5), 493 T(1)*T(10)-T(3)*T(7)+T(4)*T(6), 494 T(2)*T(10)-T(3)*T(9)+T(4)*T(8); 495 496 bigintmat Q[4][10] = 497 1,0,0,0,1,1,1,0,0,0, 498 0,1,0,0,1,0,0,1,1,0, 499 0,0,1,0,0,1,0,1,0,1, 500 0,0,0,1,0,0,1,0,1,1; 501 502 int t = rtimer; 503 fan F = gitFan(a, Q); 504 t = rtimer - t; 505 506 int tt = rtimer; 507 fan F = gitFan(a, Q, 4); 508 tt = rtimer - tt; 509 510 t; 511 tt; 512 513 "--------"; 514 kill R, Q, t, tt; 515 // next example // 516 ring R = 0,T(1..15),dp; 517 ideal a = 518 T(1)*T(10)-T(2)*T(7)+T(3)*T(6), 519 T(1)*T(11)-T(2)*T(8)+T(4)*T(6), 520 T(1)*T(12)-T(2)*T(9)+T(5)*T(6), 521 T(1)*T(13)-T(3)*T(8)+T(4)*T(7), 522 T(1)*T(14)-T(3)*T(9)+T(5)*T(7), 523 T(1)*T(15)-T(4)*T(9)+T(5)*T(8), 524 T(2)*T(13)-T(3)*T(11)+T(4)*T(10), 525 T(2)*T(14)-T(3)*T(12)+T(5)*T(10); 2725 526 2726 527 bigintmat Q[5][15] = … … 2731 532 0,0,0,0,1,0,0,0,1,0,0,1,0,1,1; 2732 533 534 int t = rtimer; 535 fan F = gitFan(a, Q); 536 t = rtimer - t; 537 538 int tt = rtimer; 539 fan F = gitFan(a, Q, 4); 540 tt = rtimer - tt; 541 542 t; 543 tt; 544 545 } 546 547 ///////////////////////////////////// 548 // Computes all simplicial orbit cones 549 // w.r.t. the 0-ideal: 550 551 static proc simplicialToricOrbitCones(bigintmat Q) 552 { 553 intvec gam0; 554 list OC; 555 cone c; 556 int r = ncols(Q); 557 int j; 558 559 for(int i = 1; i < 2^r; i++ ) 560 { 561 gam0 = int2face(i,r); 562 563 // each simplicial cone is generated by 564 // exactly nrows(Q) many columns of Q: 565 if(size(gam0) == nrows(Q)) 566 { 567 bigintmat M[size(gam0)][nrows(Q)]; 568 569 for(j = 1; j <= size(gam0); j++) 570 { 571 M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]]; 572 } 573 574 c = coneViaPoints(M); 575 576 if((dimension(c) == nrows(Q)) and (!(listContainsCone(OC, c)))) 577 { 578 OC[size(OC)+1] = c; 579 } 580 581 kill M; 582 } 583 } 584 585 return(OC); 586 } 587 588 ///////////////////////////////////// 589 590 proc gkzFan(bigintmat Q) 591 "USAGE: gkzFan(Q); a: ideal, Q:bigintmat 592 PURPOSE: Returns the GKZ-fan of the matrix Q. 593 RETURN: a fan. 594 EXAMPLE: example gkzFan; shows an example 595 " 596 { 597 // only difference to gitFan: 598 // it suffices to consider all faces 599 // that are simplicial: 2733 600 list OC = simplicialToricOrbitCones(Q); 2734 print(size(OC)); 2735 } 2736 2737 2738 proc G25Action() 2739 "USAGE: G25Action(Q); 2740 PURPOSE: Returns a representation of S5 as a subgroup of S10 with the action on the Grassmannian G25. 2741 RETURN: list of with elements of type permutation. 2742 EXAMPLE: example G25Action; shows an example 2743 " 2744 { 2745 list simplexSymmetryGroup = permutationFromIntvec(intvec( 1 .. 10 )), 2746 permutationFromIntvec(intvec( 1, 2, 4, 3, 5, 7, 6, 9, 8, 10 )), 2747 permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )), 2748 permutationFromIntvec(intvec( 1, 3, 4, 2, 6, 7, 5, 10, 8, 9 )), 2749 permutationFromIntvec(intvec( 1, 4, 2, 3, 7, 5, 6, 9, 10, 8 )), 2750 permutationFromIntvec(intvec( 1, 4, 3, 2, 7, 6, 5, 10, 9, 8 )), 2751 permutationFromIntvec(intvec( 1, 5, 6, 7, 2, 3, 4, 8, 9, 10 )), 2752 permutationFromIntvec(intvec( 1, 5, 7, 6, 2, 4, 3, 9, 8, 10 )), 2753 permutationFromIntvec(intvec( 1, 6, 5, 7, 3, 2, 4, 8, 10, 9 )), 2754 permutationFromIntvec(intvec( 1, 6, 7, 5, 3, 4, 2, 10, 8, 9 )), 2755 permutationFromIntvec(intvec( 1, 7, 5, 6, 4, 2, 3, 9, 10, 8 )), 2756 permutationFromIntvec(intvec( 1, 7, 6, 5, 4, 3, 2, 10, 9, 8 )), 2757 permutationFromIntvec(intvec( 2, 1, 3, 4, 5, 8, 9, 6, 7, 10 )), 2758 permutationFromIntvec(intvec( 2, 1, 4, 3, 5, 9, 8, 7, 6, 10 )), 2759 permutationFromIntvec(intvec( 2, 3, 1, 4, 8, 5, 9, 6, 10, 7 )), 2760 permutationFromIntvec(intvec( 2, 3, 4, 1, 8, 9, 5, 10, 6, 7 )), 2761 permutationFromIntvec(intvec( 2, 4, 1, 3, 9, 5, 8, 7, 10, 6 )), 2762 permutationFromIntvec(intvec( 2, 4, 3, 1, 9, 8, 5, 10, 7, 6 )), 2763 permutationFromIntvec(intvec( 2, 5, 8, 9, 1, 3, 4, 6, 7, 10 )), 2764 permutationFromIntvec(intvec( 2, 5, 9, 8, 1, 4, 3, 7, 6, 10 )), 2765 permutationFromIntvec(intvec( 2, 8, 5, 9, 3, 1, 4, 6, 10, 7 )), 2766 permutationFromIntvec(intvec( 2, 8, 9, 5, 3, 4, 1, 10, 6, 7 )), 2767 permutationFromIntvec(intvec( 2, 9, 5, 8, 4, 1, 3, 7, 10, 6 )), 2768 permutationFromIntvec(intvec( 2, 9, 8, 5, 4, 3, 1, 10, 7, 6 )), 2769 permutationFromIntvec(intvec( 3, 1, 2, 4, 6, 8, 10, 5, 7, 9 )), 2770 permutationFromIntvec(intvec( 3, 1, 4, 2, 6, 10, 8, 7, 5, 9 )), 2771 permutationFromIntvec(intvec( 3, 2, 1, 4, 8, 6, 10, 5, 9, 7 )), 2772 permutationFromIntvec(intvec( 3, 2, 4, 1, 8, 10, 6, 9, 5, 7 )), 2773 permutationFromIntvec(intvec( 3, 4, 1, 2, 10, 6, 8, 7, 9, 5 )), 2774 permutationFromIntvec(intvec( 3, 4, 2, 1, 10, 8, 6, 9, 7, 5 )), 2775 permutationFromIntvec(intvec( 3, 6, 8, 10, 1, 2, 4, 5, 7, 9 )), 2776 permutationFromIntvec(intvec( 3, 6, 10, 8, 1, 4, 2, 7, 5, 9 )), 2777 permutationFromIntvec(intvec( 3, 8, 6, 10, 2, 1, 4, 5, 9, 7 )), 2778 permutationFromIntvec(intvec( 3, 8, 10, 6, 2, 4, 1, 9, 5, 7 )), 2779 permutationFromIntvec(intvec( 3, 10, 6, 8, 4, 1, 2, 7, 9, 5 )), 2780 permutationFromIntvec(intvec( 3, 10, 8, 6, 4, 2, 1, 9, 7, 5 )), 2781 permutationFromIntvec(intvec( 4, 1, 2, 3, 7, 9, 10, 5, 6, 8 )), 2782 permutationFromIntvec(intvec( 4, 1, 3, 2, 7, 10, 9, 6, 5, 8 )), 2783 permutationFromIntvec(intvec( 4, 2, 1, 3, 9, 7, 10, 5, 8, 6 )), 2784 permutationFromIntvec(intvec( 4, 2, 3, 1, 9, 10, 7, 8, 5, 6 )), 2785 permutationFromIntvec(intvec( 4, 3, 1, 2, 10, 7, 9, 6, 8, 5 )), 2786 permutationFromIntvec(intvec( 4, 3, 2, 1, 10, 9, 7, 8, 6, 5 )), 2787 permutationFromIntvec(intvec( 4, 7, 9, 10, 1, 2, 3, 5, 6, 8 )), 2788 permutationFromIntvec(intvec( 4, 7, 10, 9, 1, 3, 2, 6, 5, 8 )), 2789 permutationFromIntvec(intvec( 4, 9, 7, 10, 2, 1, 3, 5, 8, 6 )), 2790 permutationFromIntvec(intvec( 4, 9, 10, 7, 2, 3, 1, 8, 5, 6 )), 2791 permutationFromIntvec(intvec( 4, 10, 7, 9, 3, 1, 2, 6, 8, 5 )), 2792 permutationFromIntvec(intvec( 4, 10, 9, 7, 3, 2, 1, 8, 6, 5 )), 2793 permutationFromIntvec(intvec( 5, 1, 6, 7, 2, 8, 9, 3, 4, 10 )), 2794 permutationFromIntvec(intvec( 5, 1, 7, 6, 2, 9, 8, 4, 3, 10 )), 2795 permutationFromIntvec(intvec( 5, 2, 8, 9, 1, 6, 7, 3, 4, 10 )), 2796 permutationFromIntvec(intvec( 5, 2, 9, 8, 1, 7, 6, 4, 3, 10 )), 2797 permutationFromIntvec(intvec( 5, 6, 1, 7, 8, 2, 9, 3, 10, 4 )), 2798 permutationFromIntvec(intvec( 5, 6, 7, 1, 8, 9, 2, 10, 3, 4 )), 2799 permutationFromIntvec(intvec( 5, 7, 1, 6, 9, 2, 8, 4, 10, 3 )), 2800 permutationFromIntvec(intvec( 5, 7, 6, 1, 9, 8, 2, 10, 4, 3 )), 2801 permutationFromIntvec(intvec( 5, 8, 2, 9, 6, 1, 7, 3, 10, 4 )), 2802 permutationFromIntvec(intvec( 5, 8, 9, 2, 6, 7, 1, 10, 3, 4 )), 2803 permutationFromIntvec(intvec( 5, 9, 2, 8, 7, 1, 6, 4, 10, 3 )), 2804 permutationFromIntvec(intvec( 5, 9, 8, 2, 7, 6, 1, 10, 4, 3 )), 2805 permutationFromIntvec(intvec( 6, 1, 5, 7, 3, 8, 10, 2, 4, 9 )), 2806 permutationFromIntvec(intvec( 6, 1, 7, 5, 3, 10, 8, 4, 2, 9 )), 2807 permutationFromIntvec(intvec( 6, 3, 8, 10, 1, 5, 7, 2, 4, 9 )), 2808 permutationFromIntvec(intvec( 6, 3, 10, 8, 1, 7, 5, 4, 2, 9 )), 2809 permutationFromIntvec(intvec( 6, 5, 1, 7, 8, 3, 10, 2, 9, 4 )), 2810 permutationFromIntvec(intvec( 6, 5, 7, 1, 8, 10, 3, 9, 2, 4 )), 2811 permutationFromIntvec(intvec( 6, 7, 1, 5, 10, 3, 8, 4, 9, 2 )), 2812 permutationFromIntvec(intvec( 6, 7, 5, 1, 10, 8, 3, 9, 4, 2 )), 2813 permutationFromIntvec(intvec( 6, 8, 3, 10, 5, 1, 7, 2, 9, 4 )), 2814 permutationFromIntvec(intvec( 6, 8, 10, 3, 5, 7, 1, 9, 2, 4 )), 2815 permutationFromIntvec(intvec( 6, 10, 3, 8, 7, 1, 5, 4, 9, 2 )), 2816 permutationFromIntvec(intvec( 6, 10, 8, 3, 7, 5, 1, 9, 4, 2 )), 2817 permutationFromIntvec(intvec( 7, 1, 5, 6, 4, 9, 10, 2, 3, 8 )), 2818 permutationFromIntvec(intvec( 7, 1, 6, 5, 4, 10, 9, 3, 2, 8 )), 2819 permutationFromIntvec(intvec( 7, 4, 9, 10, 1, 5, 6, 2, 3, 8 )), 2820 permutationFromIntvec(intvec( 7, 4, 10, 9, 1, 6, 5, 3, 2, 8 )), 2821 permutationFromIntvec(intvec( 7, 5, 1, 6, 9, 4, 10, 2, 8, 3 )), 2822 permutationFromIntvec(intvec( 7, 5, 6, 1, 9, 10, 4, 8, 2, 3 )), 2823 permutationFromIntvec(intvec( 7, 6, 1, 5, 10, 4, 9, 3, 8, 2 )), 2824 permutationFromIntvec(intvec( 7, 6, 5, 1, 10, 9, 4, 8, 3, 2 )), 2825 permutationFromIntvec(intvec( 7, 9, 4, 10, 5, 1, 6, 2, 8, 3 )), 2826 permutationFromIntvec(intvec( 7, 9, 10, 4, 5, 6, 1, 8, 2, 3 )), 2827 permutationFromIntvec(intvec( 7, 10, 4, 9, 6, 1, 5, 3, 8, 2 )), 2828 permutationFromIntvec(intvec( 7, 10, 9, 4, 6, 5, 1, 8, 3, 2 )), 2829 permutationFromIntvec(intvec( 8, 2, 5, 9, 3, 6, 10, 1, 4, 7 )), 2830 permutationFromIntvec(intvec( 8, 2, 9, 5, 3, 10, 6, 4, 1, 7 )), 2831 permutationFromIntvec(intvec( 8, 3, 6, 10, 2, 5, 9, 1, 4, 7 )), 2832 permutationFromIntvec(intvec( 8, 3, 10, 6, 2, 9, 5, 4, 1, 7 )), 2833 permutationFromIntvec(intvec( 8, 5, 2, 9, 6, 3, 10, 1, 7, 4 )), 2834 permutationFromIntvec(intvec( 8, 5, 9, 2, 6, 10, 3, 7, 1, 4 )), 2835 permutationFromIntvec(intvec( 8, 6, 3, 10, 5, 2, 9, 1, 7, 4 )), 2836 permutationFromIntvec(intvec( 8, 6, 10, 3, 5, 9, 2, 7, 1, 4 )), 2837 permutationFromIntvec(intvec( 8, 9, 2, 5, 10, 3, 6, 4, 7, 1 )), 2838 permutationFromIntvec(intvec( 8, 9, 5, 2, 10, 6, 3, 7, 4, 1 )), 2839 permutationFromIntvec(intvec( 8, 10, 3, 6, 9, 2, 5, 4, 7, 1 )), 2840 permutationFromIntvec(intvec( 8, 10, 6, 3, 9, 5, 2, 7, 4, 1 )), 2841 permutationFromIntvec(intvec( 9, 2, 5, 8, 4, 7, 10, 1, 3, 6 )), 2842 permutationFromIntvec(intvec( 9, 2, 8, 5, 4, 10, 7, 3, 1, 6 )), 2843 permutationFromIntvec(intvec( 9, 4, 7, 10, 2, 5, 8, 1, 3, 6 )), 2844 permutationFromIntvec(intvec( 9, 4, 10, 7, 2, 8, 5, 3, 1, 6 )), 2845 permutationFromIntvec(intvec( 9, 5, 2, 8, 7, 4, 10, 1, 6, 3 )), 2846 permutationFromIntvec(intvec( 9, 5, 8, 2, 7, 10, 4, 6, 1, 3 )), 2847 permutationFromIntvec(intvec( 9, 7, 4, 10, 5, 2, 8, 1, 6, 3 )), 2848 permutationFromIntvec(intvec( 9, 7, 10, 4, 5, 8, 2, 6, 1, 3 )), 2849 permutationFromIntvec(intvec( 9, 8, 2, 5, 10, 4, 7, 3, 6, 1 )), 2850 permutationFromIntvec(intvec( 9, 8, 5, 2, 10, 7, 4, 6, 3, 1 )), 2851 permutationFromIntvec(intvec( 9, 10, 4, 7, 8, 2, 5, 3, 6, 1 )), 2852 permutationFromIntvec(intvec( 9, 10, 7, 4, 8, 5, 2, 6, 3, 1 )), 2853 permutationFromIntvec(intvec( 10, 3, 6, 8, 4, 7, 9, 1, 2, 5 )), 2854 permutationFromIntvec(intvec( 10, 3, 8, 6, 4, 9, 7, 2, 1, 5 )), 2855 permutationFromIntvec(intvec( 10, 4, 7, 9, 3, 6, 8, 1, 2, 5 )), 2856 permutationFromIntvec(intvec( 10, 4, 9, 7, 3, 8, 6, 2, 1, 5 )), 2857 permutationFromIntvec(intvec( 10, 6, 3, 8, 7, 4, 9, 1, 5, 2 )), 2858 permutationFromIntvec(intvec( 10, 6, 8, 3, 7, 9, 4, 5, 1, 2 )), 2859 permutationFromIntvec(intvec( 10, 7, 4, 9, 6, 3, 8, 1, 5, 2 )), 2860 permutationFromIntvec(intvec( 10, 7, 9, 4, 6, 8, 3, 5, 1, 2 )), 2861 permutationFromIntvec(intvec( 10, 8, 3, 6, 9, 4, 7, 2, 5, 1 )), 2862 permutationFromIntvec(intvec( 10, 8, 6, 3, 9, 7, 4, 5, 2, 1 )), 2863 permutationFromIntvec(intvec( 10, 9, 4, 7, 8, 3, 6, 2, 5, 1 )), 2864 permutationFromIntvec(intvec( 10, 9, 7, 4, 8, 6, 3, 5, 2, 1 )); 2865 return(simplexSymmetryGroup); 601 602 fan f = refineCones(OC, Q); 603 return(f); 2866 604 } 2867 605 example 2868 606 { 2869 echo = 2; 2870 G25Action(); 2871 } 2872 2873 2874 proc findOrbits(list G, int #) 2875 "USAGE: findOrbits(G [,d]); G list of permutations in a subgroup of the symmetric group; d int minimum cardinality of simplices to be considered; if d is not specified all orbits are computed. 2876 PURPOSE: Computes the orbit decomposition of the action of G. 2877 RETURN: list of intvec. 2878 EXAMPLE: example findOrbits; shows an example 2879 " 2880 { 2881 int d; 2882 if (size(#)>0){d=#;} 2883 int n = size(permutationToIntvec(G[1])); 2884 list listOrbits; 2885 list finished; 2886 int tst; 2887 bigint startel; 2888 intvec startelbin; 2889 list neworbit,neworbitint; 2890 int i,posToInsert; 2891 bigint nn; 2892 while (size(finished)<2^n){ 2893 startel=0; 2894 if (size(finished)>0){ 2895 tst=0; 2896 while ((tst==0) and (startel+1<=size(finished))){ 2897 if (finished[int(startel+1)]<>startel){ 2898 tst=1; 2899 } else { 2900 startel=startel+1; 2901 } 2902 } 2903 } 2904 if (startel==0){ 2905 neworbit[1]= list(); 2906 neworbitint[1]=0; 2907 } else { 2908 startelbin=bigintToBinary(startel,n); 2909 neworbitint=list(startel); 2910 for (i=2;i<=size(G);i++){ 2911 nn=binaryToBigint(applyPermutationToIntvec(startelbin,G[i])); 2912 //nn;neworbitint; 2913 //"place to insert"; 2914 posToInsert = findPlaceToInsert(neworbitint,nn); 2915 //"pos";posToInsert; 2916 if(posToInsert > 0) 2917 { 2918 //"vorher";neworbitint; 2919 neworbitint = insertToList(neworbitint,nn,posToInsert); 2920 //"nachher";neworbitint; 2921 } 2922 } 2923 for (i=1;i<=size(neworbitint);i++){ 2924 neworbit[i]=bigintToBinary(neworbitint[i],n); 2925 } 2926 } 2927 if (size(neworbit[1])>=d){ 2928 listOrbits[size(listOrbits)+1]=neworbit; 2929 } 2930 finished=sort(finished+neworbitint)[1]; 2931 } 2932 return(listOrbits);} 2933 example 2934 { 2935 list G = G25Action(); 2936 list orb = findOrbits(G); 2937 for (int i=1;i<=size(orb);i++){orb[i][1];} 2938 } 2939 2940 2941 2942 static proc GITfanWrapperWithSymmetry(ideal J, intmat Q, list simplexSymmetryGroup){ 2943 list orb = findOrbits(simplexSymmetryGroup,nrows(Q)); 2944 list simplexOrbitRepresentatives; 2945 for (int i=1;i<=size(orb);i++){simplexOrbitRepresentatives[i]=orb[i][1];} 2946 list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives); 2947 list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q); 2948 list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup); 2949 list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits); 2950 list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q); 2951 list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits); 2952 list Asigma = groupActionOnQImage(simplexSymmetryGroup,Q); 2953 list actionOnOrbitconeIndices = groupActionOnHashes(Asigma,listOfOrbitConeOrbits); 2954 list OClist = listOfOrbitConeOrbits[1]; 2955 for (i =2;i<=size(listOfOrbitConeOrbits);i++){ 2956 OClist = OClist + listOfOrbitConeOrbits[i]; 2957 } 2958 cone mov = coneViaPoints(transpose(Q)); 2959 mov = canonicalizeCone(mov); 2960 list Sigma = GITfanParallelSymmetric(OClist, Q, mov, actionOnOrbitconeIndices); 2961 return(Sigma, OClist); 2962 } 607 echo=2; 608 intmat Q[3][4] = 609 1,0,1,0, 610 0,1,0,1, 611 0,0,1,1; 612 613 gkzFan(Q); 614 }
Note: See TracChangeset
for help on using the changeset viewer.