1 #include "symmetriccomplex.h"
2
3 #include <sstream>
4 #include "polyhedralcone.h"
5 #include "printer.h"
6 #include "lp.h"
7 #include "linalg.h"
8 #include "determinant.h"
9 #include "log.h"
10 #include <iostream>
11
12 // {ab+aa+cc+dd,ba+ca+db}
13
Cone(set<int> const & indices_,int dimension_,int multiplicity_,bool sortWithSymmetry,SymmetricComplex const & complex)14 SymmetricComplex::Cone::Cone(set<int> const &indices_, int dimension_, int multiplicity_, bool sortWithSymmetry, SymmetricComplex const &complex):
15 dimension(dimension_),
16 multiplicity(multiplicity_),
17 isKnownToBeNonMaximalFlag(false)
18 {
19 indices=vector<int>(indices_.size());
20 int j=0;
21 for(set<int>::const_iterator i=indices_.begin();i!=indices_.end();i++,j++)
22 indices[j]=*i;
23
24 IntegerMatrix const &vertices=complex.getVertices();
25 IntegerVector sum(vertices.getWidth());
26 for(int i=0;i<indices.size();i++)
27 sum+=vertices[indices[i]];
28
29 if(sortWithSymmetry)
30 {
31 // sortKey=complex.sym.orbitRepresentative(sum,&sortKeyPermutation);
32 sortKey=complex.sym.orbitRepresentative(sum,sortKeyPermutationIndex);
33 // sortKey=sum;
34
35 /*
36 LexicographicTermOrder myOrder;
37
38 for(SymmetryGroup::ElementContainer::const_iterator k=complex.sym.elements.begin();k!=complex.sym.elements.end();k++)
39 if(myOrder(SymmetryGroup::compose(*k,sum),sortKey))
40 sortKey=SymmetryGroup::compose(*k,sum);
41 */
42
43
44 /*
45 int n=sum.size();
46 for(SymmetryGroup::ElementContainer::const_iterator k=complex.sym.elements.begin();k!=complex.sym.elements.end();k++)
47 {
48 bool isBetter=true;
49 for(int i=0;i<n;i++)
50 {
51 if(sum[(*k)[i]]>sortKey[i]){isBetter=false;break;}
52 if(sum[(*k)[i]]<sortKey[i])break;
53 }
54 if(isBetter)
55 {
56 sortKey=SymmetryGroup::compose(*k,sum);
57 }
58 }
59 */
60
61 }
62 else
63 {
64 sortKey=sum;
65 }
66 }
67
68
indexOfVertex(IntegerVector const & v) const69 int SymmetricComplex::indexOfVertex(IntegerVector const &v)const
70 {
71 map<IntegerVector,int>::const_iterator it=indexMap.find(v);
72 assert(it!=indexMap.end());
73 return it->second;
74 }
75
76
remap(SymmetricComplex & complex)77 void SymmetricComplex::Cone::remap(SymmetricComplex &complex)
78 {
79 IntegerMatrix const &vertices=complex.getVertices();
80 IntegerVector sum(vertices.getWidth());
81 for(int i=0;i<indices.size();i++)
82 sum+=vertices[indices[i]];
83
84 int n=sum.size();
85 /* IntegerVector bestPermutation;
86 for(SymmetryGroup::ElementContainer::const_iterator k=complex.sym.elements.begin();k!=complex.sym.elements.end();k++)
87 {
88 if(SymmetryGroup::compose(*k,sum)==sortKey)
89 bestPermutation=*k;
90 }
91 */
92 IntegerVector const &bestPermutation=complex.sym.referenceToIthVector(sortKeyPermutationIndex);
93
94 assert(bestPermutation.size()==n);
95
96 vector<int> indicesNew(indices.size());
97 int I=0;
98 for(vector<int>::const_iterator i=indices.begin();i!=indices.end();i++,I++)
99 {
100 IntegerVector ny=SymmetryGroup::compose(bestPermutation,complex.vertices[*i]);
101 map<IntegerVector,int>::const_iterator it=complex.indexMap.find(ny);
102 assert(it!=complex.indexMap.end());
103 indicesNew[I]=it->second;
104 }
105 indices=indicesNew;
106 }
107
108
indexSet() const109 set<int> SymmetricComplex::Cone::indexSet()const
110 {
111 set<int> ret;
112 for(vector<int>::const_iterator i=indices.begin();i!=indices.end();i++)
113 ret.insert(*i);
114
115 return ret;
116 }
117
isSubsetOf(Cone const & c) const118 bool SymmetricComplex::Cone::isSubsetOf(Cone const &c)const
119 {
120 int next=0;
121 for(int i=0;i<indices.size();i++)
122 {
123 while(1)
124 {
125 if(next>=c.indices.size())return false;
126 if(indices[i]==c.indices[next])break;
127 next++;
128 }
129 }
130 return true;
131 /*
132 set<int> b=c.indexSet();
133
134 for(vector<int>::const_iterator i=indices.begin();i!=indices.end();i++)
135 if(b.count(*i)==0)return false;
136 return true;
137 */
138 }
139
140
permuted(IntegerVector const & permutation,SymmetricComplex const & complex,bool withSymmetry) const141 SymmetricComplex::Cone SymmetricComplex::Cone::permuted(IntegerVector const &permutation, SymmetricComplex const &complex, bool withSymmetry)const
142 {
143 /* Cone ret;
144 ret.dimension=dimension;
145 ret.multiplicity=multiplicity;*/
146
147 set<int> r;
148 for(vector<int>::const_iterator i=indices.begin();i!=indices.end();i++)
149 {
150 IntegerVector ny=SymmetryGroup::compose(permutation,complex.vertices[*i]);
151 map<IntegerVector,int>::const_iterator it=complex.indexMap.find(ny);
152 if(it==complex.indexMap.end())
153 {
154 AsciiPrinter(Stderr).printVector(complex.vertices[*i]);
155 AsciiPrinter(Stderr).printVector(ny);
156
157 assert(0);
158 }
159 r.insert(it->second);
160 }
161
162
163 return Cone(r,dimension,multiplicity,withSymmetry,complex);
164 }
165
166
167 /*
168 void SymmetricComplex::Cone::computeRelativeInteriorPoint(SymmetricComplex const &complex)
169 {
170 IntegerMatrix const &vertices=complex.getVertices();
171 IntegerVector sum(vertices.getWidth());
172 for(const_iterator i=begin();i!=end();i++)
173 sum+=vertices[*i];
174 relativeInteriorPoint=sum;
175 sum.sort();
176 summary=sum;
177 }
178 */
179
180 /*void SymmetricComplex::Cone::computeSmallestRepresentative(SymmetricComplex const &complex)
181 {
182 if(relativeInteriorPoint.size()==0)computeRelativeInteriorPoint(complex);
183
184 LexicographicTermOrder myOrder;
185
186 smallestRepresentative=relativeInteriorPoint;
187
188 for(SymmetryGroup::ElementContainer::const_iterator k=complex.sym.elements.begin();k!=complex.sym.elements.end();k++)
189 if(myOrder(SymmetryGroup::compose(*k,relativeInteriorPoint),smallestRepresentative))
190 smallestRepresentative=SymmetryGroup::compose(*k,relativeInteriorPoint);
191 }
192 */
193
operator <(Cone const & b) const194 bool SymmetricComplex::Cone::operator<(Cone const & b)const
195 {
196 /* if(ignoreSymmetry)return ((set<int>)*this)<(b);*/
197 return sortKey<b.sortKey;
198 }
199
200
isSimplicial(int linealityDim) const201 bool SymmetricComplex::Cone::isSimplicial(int linealityDim)const
202 {
203 return (indices.size()+linealityDim)==dimension;
204 }
205
206
orthogonalComplement(SymmetricComplex & complex) const207 IntegerVectorList SymmetricComplex::Cone::orthogonalComplement(SymmetricComplex &complex)const
208 {
209 IntegerVectorList l;
210 for(int i=0;i<indices.size();i++)
211 l.push_back(complex.vertices[indices[i]]);
212
213 FieldMatrix m=integerMatrixToFieldMatrix(rowsToIntegerMatrix(l,complex.n),Q);
214 return fieldMatrixToIntegerMatrixPrimitive(m.reduceAndComputeKernel()).getRows();
215 }
216
217
SymmetricComplex(int n_,IntegerVectorList const & v,SymmetryGroup const & sym_)218 SymmetricComplex::SymmetricComplex(int n_, IntegerVectorList const &v, SymmetryGroup const &sym_):
219 n(n_),
220 sym(sym_),
221 dimension(-1)
222 {
223 vertices=rowsToIntegerMatrix(v,n);
224 for(int i=0;i<vertices.getHeight();i++)indexMap[vertices[i]]=i;
225 }
226
227
contains(Cone const & c) const228 bool SymmetricComplex::contains(Cone const &c)const
229 {
230 Cone temp=c;//#1
231 /* temp.computeRelativeInteriorPoint(*this);
232
233 temp.computeSmallestRepresentative(*this);
234 */
235 return cones.find(temp)!=cones.end();///////////////////!!!!!!!!!!!!!!!!!!!!!!!
236
237
238 /*
239 set<IntegerVector> possibleMatches;
240 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
241 {
242 */
243 /* AsciiPrinter(Stderr).printVector(i->summary);
244 AsciiPrinter(Stderr).printVector(temp.summary);
245 fprintf(stderr,"\n");*/
246 /* if(i->dimension==temp.dimension)
247 if(i->summary==temp.summary)
248 {
249 possibleMatches.insert(i->relativeInteriorPoint);
250 }
251 }
252 for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
253 {
254 if(possibleMatches.find(SymmetryGroup::compose(*k,temp.relativeInteriorPoint))!=possibleMatches.end())
255 return true;
256 }
257 */
258 /* for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
259 {
260 Cone c2=c.permuted(*k,*this);
261 for(list<Cone>::const_iterator i=cones.begin();i!=cones.end();i++)
262 {
263 if(c2==*i)return true;
264 }
265 }*/
266
267 // return false;
268 }
269
270
insert(Cone const & c)271 void SymmetricComplex::insert(Cone const &c)
272 {
273 if(c.dimension>dimension)dimension=c.dimension;
274 // Cone temp=c;
275 /* temp.computeRelativeInteriorPoint(*this);
276 temp.computeSmallestRepresentative(*this);
277 */
278 if(!contains(c))//#2
279 {
280 cones.insert(c);
281 //cones.push_back(temp);
282 // fprintf(Stderr,"INSERTING\n");
283 // cones.back().computeRelativeInteriorPoint(*this);
284 }
285 else
286 {
287 // if(c.isKnownToBeNonMaximal())cones.find(c)->setKnownToBeNonMaximal();
288 if(c.isKnownToBeNonMaximal()){cones.erase(c);cones.insert(c);}// mark as non-maximal
289 }
290 }
291
292
getMaxDim() const293 int SymmetricComplex::getMaxDim()const
294 {/*
295 int ret=-1;
296 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
297 {
298 if(i->dimension>ret)ret=i->dimension;
299 }
300 return ret;*/
301 return dimension;
302 }
303
304
getMinDim() const305 int SymmetricComplex::getMinDim()const
306 {
307 int ret=100000;
308 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
309 {
310 if(i->dimension<ret)ret=i->dimension;
311 }
312 return ret;
313 }
314
315
isMaximal(Cone const & c) const316 bool SymmetricComplex::isMaximal(Cone const &c)const
317 {
318 if(c.isKnownToBeNonMaximal())return false;
319 if(c.dimension==dimension)return true;
320 for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
321 {
322 Cone c2=c.permuted(*k,*this,false);
323 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
324 {
325 if(i->dimension>c.dimension)
326 if(c2.isSubsetOf(*i) && !i->isSubsetOf(c2))return false;
327 }
328 }
329 return true;
330 }
331
332
dimensionsAtInfinity() const333 IntegerVector SymmetricComplex::dimensionsAtInfinity()const
334 {
335 /* Using a double description like method this routine computes the
336 dimension of the intersection of each cone in the complex with
337 the plane x_0=0 */
338
339 IntegerVector ret(cones.size());
340
341 int I=0;
342 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++,I++)
343 {
344 IntegerVectorList raysAtInfinity;
345 for(vector<int>::const_iterator j=i->indices.begin();j!=i->indices.end();j++)
346 {
347 if(vertices[*j][0]==0)raysAtInfinity.push_back(vertices[*j]);
348 for(vector<int>::const_iterator k=j;k!=i->indices.end();k++)
349 if(vertices[*j][0]*vertices[*k][0]<0)
350 raysAtInfinity.push_back(((vertices[*j][0]>0)?1:-1)*(vertices[*j][0])*vertices[*k].toVector()+
351 ((vertices[*k][0]>0)?1:-1)*(vertices[*k][0])*vertices[*j].toVector());
352 }
353 ret[I]=rankOfMatrix(raysAtInfinity);
354 }
355 return ret;
356 }
357
358
toString(int dimLow,int dimHigh,bool onlyMaximal,bool group,ostream * multiplicities,bool compressed,bool tPlaneSort,bool xml) const359 string SymmetricComplex::toString(int dimLow, int dimHigh, bool onlyMaximal, bool group, ostream *multiplicities, bool compressed, bool tPlaneSort, bool xml)const
360 {
361 stringstream ret;
362
363 if(!onlyMaximal)
364 {
365 if(xml)ret<<"<m>\n";
366 if(xml)if(multiplicities)*multiplicities<<"<m>\n";
367 }
368 IntegerVector additionalSortKeys(cones.size());
369 if(tPlaneSort)additionalSortKeys=dimensionsAtInfinity();
370 int lowKey=additionalSortKeys.min();
371 int highKey=additionalSortKeys.max();
372
373 if(xml)if(onlyMaximal)
374 {
375 if(compressed)
376 ret<<"<m>\n";
377 else
378 ret<<"<m cols=\""<<this->vertices.getHeight()<<"\">\n";
379 if(multiplicities)*multiplicities<<"<m>\n";
380 }
381 for(int d=dimLow;d<=dimHigh;d++)
382 {
383 log1 cerr << "Processing dimension "<<d<<".\n";
384 int numberOfOrbitsOutput=0;
385 int numberOfOrbitsOfThisDimension=0;
386 bool newDimension=true;
387 for(int key=lowKey;key<=highKey;key++)
388 {
389 if(xml)
390 {
391 if(!onlyMaximal)
392 {
393 if(compressed)
394 ret<<"<m>\n";
395 else
396 ret<<"<m cols=\""<<this->vertices.getHeight()<<"\">\n";
397 if(multiplicities)*multiplicities<<"<m>\n";
398 }
399 }
400 int I=0;
401 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++,I++)
402 if(additionalSortKeys[I]==key)
403 if(i->dimension==d)
404 {
405 numberOfOrbitsOfThisDimension++;
406 if(!onlyMaximal || isMaximal(*i))
407 {
408 numberOfOrbitsOutput++;
409 bool isMax=isMaximal(*i);
410 bool newOrbit=true;
411 set<set<int> > temp;
412 for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
413 {
414 Cone temp1=i->permuted(*k,*this,false);
415 temp.insert(temp1.indexSet());
416 if(compressed)break;
417 }
418 for(set<set<int> >::const_iterator j=temp.begin();j!=temp.end();j++)
419 {
420 if(!xml)ret << "{";else ret<<"<v>";
421 for(set<int>::const_iterator a=j->begin();a!=j->end();a++)
422 {
423 if(a!=j->begin())ret<<" ";
424 ret << *a;
425 }
426 if(!xml)ret << "}";else ret<<"</v>";
427 if(!xml)
428 {
429 if(group)if(newOrbit)ret << "\t# New orbit";
430 if(newDimension)ret << "\t# Dimension "<<d;
431 }
432 ret <<endl;
433 if(isMax)if(multiplicities)
434 {
435 if(xml)*multiplicities<<"<v>";
436 *multiplicities << i->multiplicity;
437 if(xml)*multiplicities<<"</v>";
438 if(!xml)
439 {
440 if(group)if(newOrbit)*multiplicities << "\t# New orbit";
441 if(newDimension)*multiplicities << "\t# Dimension "<<d;
442 }
443 *multiplicities << endl;
444 }
445 newOrbit=false;
446 newDimension=false;
447 }
448 }
449 }
450 if(xml)
451 {
452 if(!onlyMaximal)
453 {
454 ret<<"</m>\n";
455 if(multiplicities)*multiplicities<<"</m>\n";
456 }
457 }
458 }
459 log1 cerr<<"Number of orbits of this dimension: " << numberOfOrbitsOfThisDimension << endl;
460 log1 cerr<<"Number of orbits output: " << numberOfOrbitsOutput << endl;
461 }
462 if(xml)if(onlyMaximal)
463 {
464 ret<<"</m>\n";
465 if(multiplicities)*multiplicities<<"</m>\n";
466 }
467
468 if(!onlyMaximal)
469 {
470 if(xml)ret<<"</m>\n";
471 if(xml)if(multiplicities)*multiplicities<<"</m>\n";
472 }
473
474 return ret.str();
475 }
476
477
fvector(bool boundedPart) const478 IntegerVector SymmetricComplex::fvector(bool boundedPart)const
479 {
480 int min=getMinDim();
481 IntegerVector ret(getMaxDim()-min+1);
482
483 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
484 {
485 /* set<Cone> temp;
486 for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
487 temp.insert(i->permuted(*k,*this));*/
488 /* set<IntegerVector> temp;
489 for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
490 temp.insert(SymmetryGroup::compose(*k,i->sortKey));
491 */
492 bool doAdd=!boundedPart;
493 if(boundedPart)
494 {
495 bool isBounded=true;
496 for(vector<int>::const_iterator j=i->indices.begin();j!=i->indices.end();j++)
497 if(vertices[*j][0]==0)isBounded=false;
498 doAdd=isBounded;
499 }
500 if(doAdd)
501 ret[i->dimension-min]+=sym.orbitSize(i->sortKey);
502 }
503 return ret;
504 }
505
506
isPure() const507 bool SymmetricComplex::isPure()const
508 {
509 int dim=-1;
510 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
511 {
512 gfan_log2{static int a;if(!((a++)&63))fprintf(Stderr,"%i\n",a);}//log0
513 if(isMaximal(*i))
514 {
515 int dim2=i->dimension;
516 if(dim==-1)dim=dim2;
517 if(dim!=dim2)return false;
518 }
519 }
520 return true;
521 }
522
523
isSimplicial() const524 bool SymmetricComplex::isSimplicial()const
525 {
526 int linealityDim=getMinDim();
527 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
528 if(!i->isSimplicial(linealityDim))
529 return false;
530 return true;
531 }
532
533
remap()534 void SymmetricComplex::remap()
535 {
536 for(ConeContainer::iterator i=cones.begin();i!=cones.end();i++)
537 {
538 Cone const&j=*i;
539 Cone &j2=const_cast<Cone&>(j);//DANGER: cast away const. This does not change the sort key in the container, so should be OK.
540 j2.remap(*this);
541 }
542 }
543
544
numberOfConesOfDimension(int d) const545 int SymmetricComplex::numberOfConesOfDimension(int d)const
546 {
547 assert(sym.isTrivial());
548
549 int ret=0;
550 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
551 if(d==i->dimension)
552 {
553 ret++;
554 }
555 return ret;
556 }
557
558
dimensionIndex(Cone const & c)559 int SymmetricComplex::dimensionIndex(Cone const &c)
560 {
561 assert(sym.isTrivial());
562
563 int ret=0;
564 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
565 if(c.dimension==i->dimension)
566 {
567 if(!(c<*i)&&!(*i<c))
568 return ret;
569 else
570 ret++;
571 }
572 return ret;
573 }
574
575
boundary(Cone const & c,vector<int> & indices_,vector<int> & signs)576 void SymmetricComplex::boundary(Cone const &c, vector<int> &indices_, vector<int> &signs)
577 {
578 indices_=vector<int>();
579 signs=vector<int>();
580 int d=c.dimension;
581
582
583 IntegerVectorList l;
584 for(int i=0;i<c.indices.size();i++)
585 l.push_back(vertices[c.indices[i]]);
586 IntegerVectorList facetNormals=PolyhedralCone(l,IntegerVectorList(),n).extremeRays();
587 IntegerVectorList complementBasis=c.orthogonalComplement(*this);
588 for(IntegerVectorList::const_iterator i=facetNormals.begin();i!=facetNormals.end();i++)
589 {
590 IntegerVectorList complementBasis1=complementBasis;
591 complementBasis1.push_back(*i);
592 FieldMatrix m=integerMatrixToFieldMatrix(rowsToIntegerMatrix(complementBasis1,n),Q);
593 IntegerVectorList completion=fieldMatrixToIntegerMatrixPrimitive(m.reduceAndComputeKernel()).getRows();
594 for(IntegerVectorList::const_iterator j=completion.begin();j!=completion.end();j++)complementBasis1.push_back(*j);
595 int sign=determinantSign(complementBasis1);
596
597
598
599 set<int> indices;
600 for(vector<int>::const_iterator j=c.indices.begin();j!=c.indices.end();j++)if(dotLong(vertices[*j],*i)==0)indices.insert(*j);
601 Cone facet(indices,d-1,1,true,*this);
602 IntegerVectorList complementBasis2=facet.orthogonalComplement(*this);
603 for(IntegerVectorList::const_iterator j=completion.begin();j!=completion.end();j++)complementBasis2.push_back(*j);
604 indices_.push_back(dimensionIndex(facet));
605 signs.push_back(sign*determinantSign(complementBasis2));
606 }
607 }
608
609
boundaryMap(int d)610 IntegerMatrix SymmetricComplex::boundaryMap(int d)
611 {
612 assert(sym.isTrivial());
613
614 IntegerMatrix ret(numberOfConesOfDimension(d-1),numberOfConesOfDimension(d));
615
616 for(ConeContainer::const_iterator i=cones.begin();i!=cones.end();i++)
617 if(d==i->dimension)
618 {
619 int I=dimensionIndex(*i);
620 vector<int> indices;
621 vector<int> signs;
622 boundary(*i,indices,signs);
623 for(int j=0;j<indices.size();j++)
624 {
625 ret[indices[j]][I]+=signs[j];
626 }
627 }
628 return ret;
629 }
630
631