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