1 /*
2  * gfanlib_polyhedralfan.cpp
3  *
4  *  Created on: Nov 16, 2010
5  *      Author: anders
6  */
7 
8 #include <stddef.h>
9 #include <sstream>
10 #include "gfanlib_polyhedralfan.h"
11 #include "gfanlib_polymakefile.h"
12 
13 using namespace std;
14 namespace gfan{
15 
PolyhedralFan(int ambientDimension)16 PolyhedralFan::PolyhedralFan(int ambientDimension):
17   n(ambientDimension),
18   symmetries(n)
19 {
20 }
21 
PolyhedralFan(SymmetryGroup const & sym)22 PolyhedralFan::PolyhedralFan(SymmetryGroup const &sym):
23   n(sym.sizeOfBaseSet()),
24   symmetries(sym)
25 {
26 
27 }
28 
fullSpace(int n)29 PolyhedralFan PolyhedralFan::fullSpace(int n)
30 {
31   PolyhedralFan ret(n);
32 
33   ZCone temp(n);
34   temp.canonicalize();
35   ret.cones.insert(temp);
36 
37   return ret;
38 }
39 
40 
facetsOfCone(ZCone const & c)41 PolyhedralFan PolyhedralFan::facetsOfCone(ZCone const &c)
42 {
43   ZCone C(c);
44   C.canonicalize();
45   PolyhedralFan ret(C.ambientDimension());
46 
47   ZMatrix halfSpaces=C.getFacets();
48 
49   for(int i=0;i<halfSpaces.getHeight();i++)
50     {
51       ZMatrix a(0,C.ambientDimension());
52       ZMatrix b(0,C.ambientDimension());
53       b.appendRow(halfSpaces[i]);
54       ZCone N=intersection(ZCone(a,b),c);
55       N.canonicalize();
56       ret.cones.insert(N);
57     }
58   return ret;
59 }
60 
61 
getAmbientDimension() const62 int PolyhedralFan::getAmbientDimension()const
63 {
64   return n;
65 }
66 
isEmpty() const67 bool PolyhedralFan::isEmpty()const
68 {
69   return cones.empty();
70 }
71 
getMaxDimension() const72 int PolyhedralFan::getMaxDimension()const
73 {
74   assert(!cones.empty());
75 
76   return cones.begin()->dimension();
77 }
78 
getMinDimension() const79 int PolyhedralFan::getMinDimension()const
80 {
81   assert(!cones.empty());
82 
83   return cones.rbegin()->dimension();
84 }
85 
86 /*
87 PolyhedralFan refinement(const PolyhedralFan &a, const PolyhedralFan &b, int cutOffDimension, bool allowASingleConeOfCutOffDimension)
88 {
89   //  fprintf(Stderr,"PolyhedralFan refinement: #A=%i #B=%i\n",a.cones.size(),b.cones.size());
90   int conesSkipped=0;
91   int numberOfComputedCones=0;
92   bool linealityConeFound=!allowASingleConeOfCutOffDimension;
93   assert(a.getAmbientDimension()==b.getAmbientDimension());
94 
95   PolyhedralFan ret(a.getAmbientDimension());
96 
97   for(PolyhedralConeList::const_iterator A=a.cones.begin();A!=a.cones.end();A++)
98     for(PolyhedralConeList::const_iterator B=b.cones.begin();B!=b.cones.end();B++)
99       {
100         PolyhedralCone c=intersection(*A,*B);
101         int cdim=c.dimension();
102         //      assert(cdim>=linealitySpaceDimension);
103         bool thisIsLinealityCone=(cutOffDimension>=cdim);
104         if((!thisIsLinealityCone)||(!linealityConeFound))
105           {
106             numberOfComputedCones++;
107             c.canonicalize();
108             ret.cones.insert(c);
109             linealityConeFound=linealityConeFound || thisIsLinealityCone;
110           }
111         else
112           {
113             conesSkipped++;
114           }
115       }
116   //  fprintf(Stderr,"Number of skipped cones: %i, lineality cone found: %i\n",conesSkipped,linealityConeFound);
117   //  fprintf(Stderr,"Number of computed cones: %i, number of unique cones: %i\n",numberOfComputedCones,ret.cones.size());
118 
119   return ret;
120 }
121 */
122 
123 /*
124 PolyhedralFan product(const PolyhedralFan &a, const PolyhedralFan &b)
125 {
126   PolyhedralFan ret(a.getAmbientDimension()+b.getAmbientDimension());
127 
128   for(PolyhedralConeList::const_iterator A=a.cones.begin();A!=a.cones.end();A++)
129     for(PolyhedralConeList::const_iterator B=b.cones.begin();B!=b.cones.end();B++)
130       ret.insert(product(*A,*B));
131 
132   return ret;
133 }
134 */
135 
136 /*IntegerVectorList PolyhedralFan::getRays(int dim)
137 {
138   IntegerVectorList ret;
139   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
140     {
141       if(i->dimension()==dim)
142         ret.push_back(i->getRelativeInteriorPoint());
143     }
144   return ret;
145 }
146 */
147 
148 /*IntegerVectorList PolyhedralFan::getRelativeInteriorPoints()
149 {
150   IntegerVectorList ret;
151   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
152     {
153       ret.push_back(i->getRelativeInteriorPoint());
154     }
155   return ret;
156 }
157 */
158 
159 /*PolyhedralCone const& PolyhedralFan::highestDimensionalCone()const
160 {
161   return *cones.rbegin();
162 }
163 */
insert(ZCone const & c)164 void PolyhedralFan::insert(ZCone const &c)
165 {
166   ZCone temp=c;
167   temp.canonicalize();
168   cones.insert(temp);
169 }
170 
remove(ZCone const & c)171 void PolyhedralFan::remove(ZCone const &c)
172 {
173   cones.erase(c);
174 }
175 
176 /*
177 void PolyhedralFan::removeAllExcept(int a)
178 {
179   PolyhedralConeList::iterator i=cones.begin();
180   while(a>0)
181     {
182       if(i==cones.end())return;
183       i++;
184       a--;
185     }
186   cones.erase(i,cones.end());
187 }
188 */
189 
removeAllLowerDimensional()190 void PolyhedralFan::removeAllLowerDimensional()
191 {
192   if(!cones.empty())
193     {
194       int d=getMaxDimension();
195       PolyhedralConeList::iterator i=cones.begin();
196       while(i!=cones.end() && i->dimension()==d)i++;
197       cones.erase(i,cones.end());
198     }
199 }
200 
201 
facetComplex() const202 PolyhedralFan PolyhedralFan::facetComplex()const
203 {
204   //  fprintf(Stderr,"Computing facet complex...\n");
205   PolyhedralFan ret(n);
206 
207   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
208     {
209       PolyhedralFan a=facetsOfCone(*i);
210       for(PolyhedralConeList::const_iterator j=a.cones.begin();j!=a.cones.end();j++)
211         ret.insert(*j);
212     }
213   //  fprintf(Stderr,"Done computing facet complex.\n");
214   return ret;
215 }
216 
217 
218 /*
219 PolyhedralFan PolyhedralFan::fullComplex()const
220 {
221   PolyhedralFan ret=*this;
222 
223   while(1)
224     {
225       log2 debug<<"looping";
226       bool doLoop=false;
227       PolyhedralFan facets=ret.facetComplex();
228       log2 debug<<"number of facets"<<facets.size()<<"\n";
229       for(PolyhedralConeList::const_iterator i=facets.cones.begin();i!=facets.cones.end();i++)
230         if(!ret.contains(*i))
231           {
232             ret.insert(*i);
233             doLoop=true;
234           }
235       if(!doLoop)break;
236     }
237   return ret;
238 }
239 */
240 
241 
242 #if 0
243 PolyhedralFan PolyhedralFan::facetComplexSymmetry(SymmetryGroup const &sym, bool keepRays, bool dropLinealitySpace)const
244 {
245   log1 fprintf(Stderr,"Computing facet complex...\n");
246   PolyhedralFan ret(n);
247 
248   vector<IntegerVector> relIntTable;
249   vector<int> dimensionTable;
250 
251   if(keepRays)
252     for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
253       if(i->dimension()==i->dimensionOfLinealitySpace()+1)
254         {
255           relIntTable.push_back(i->getRelativeInteriorPoint());
256           dimensionTable.push_back(i->dimension());
257           ret.insert(*i);
258         }
259 
260   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
261     {
262       int iDim=i->dimension();
263       if(dropLinealitySpace && (i->dimension()==i->dimensionOfLinealitySpace()+1))continue;
264 
265       //      i->findFacets();
266       IntegerVectorList normals=i->getHalfSpaces();
267       for(IntegerVectorList::const_iterator j=normals.begin();j!=normals.end();j++)
268         {
269           bool alreadyInRet=false;
270           for(int l=0;l<relIntTable.size();l++)
271             {
272               if(dimensionTable[l]==iDim-1)
273                 for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
274                   {
275                     IntegerVector u=SymmetryGroup::compose(*k,relIntTable[l]);
276                     if((dotLong(*j,u)==0) && (i->contains(u)))
277                       {
278                         alreadyInRet=true;
279                         goto exitLoop;
280                       }
281                   }
282             }
283         exitLoop:
284           if(!alreadyInRet)
285             {
286               IntegerVectorList equations=i->getEquations();
287               IntegerVectorList inequalities=i->getHalfSpaces();
288               equations.push_back(*j);
289               PolyhedralCone c(inequalities,equations,n);
290               c.canonicalize();
291               ret.insert(c);
292               relIntTable.push_back(c.getRelativeInteriorPoint());
293               dimensionTable.push_back(c.dimension());
294             }
295         }
296     }
297   log1 fprintf(Stderr,"Done computing facet complex.\n");
298   return ret;
299 }
300 #endif
301 
302 
303 
getRaysInPrintingOrder(bool upToSymmetry) const304 ZMatrix PolyhedralFan::getRaysInPrintingOrder(bool upToSymmetry)const
305 {
306   /*
307    * The ordering changed in this version. Previously the orbit representatives stored in "rays" were
308    * just the first extreme ray from the orbit that appeared. Now it is gotten using "orbitRepresentative"
309    * which causes the ordering in which the different orbits appear to change.
310    */
311   if(cones.empty())return ZMatrix(0,n);
312 ZMatrix generatorsOfLinealitySpace=cones.begin()->generatorsOfLinealitySpace();//all cones have the same lineality space
313 
314 
315   std::set<ZVector> rays;//(this->getAmbientDimension());
316 //  log1 fprintf(Stderr,"Computing rays of %i cones\n",cones.size());
317   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
318     {
319       ZMatrix temp=i->extremeRays(&generatorsOfLinealitySpace);
320  //     std::cerr<<temp;
321       for(int j=0;j<temp.getHeight();j++)
322         rays.insert(symmetries.orbitRepresentative(temp[j]));
323     }
324   ZMatrix ret(0,getAmbientDimension());
325   if(upToSymmetry)
326     for(set<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++)ret.appendRow(*i);
327   else
328     for(set<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++)
329       {
330         set<ZVector> thisOrbitsRays;
331         for(SymmetryGroup::ElementContainer::const_iterator k=symmetries.elements.begin();k!=symmetries.elements.end();k++)
332           {
333             ZVector temp=k->apply(*i);
334             thisOrbitsRays.insert(temp);
335           }
336 
337         for(set<ZVector>::const_iterator i=thisOrbitsRays.begin();i!=thisOrbitsRays.end();i++)ret.appendRow(*i);
338       }
339   return ret;
340 }
341 
342 
343 
344 /*MARKS CONES AS NONMAXIMAL IN THE SYMMETRIC COMPLEX IN WHICH THEY WILL BE INSERTED -not to be confused with the facet testing in the code
345    */
computeFacets(SymmetricComplex::Cone const & theCone,ZMatrix const & rays,ZMatrix const & facetCandidates,SymmetricComplex const & theComplex)346  static list<SymmetricComplex::Cone> computeFacets(SymmetricComplex::Cone const &theCone, ZMatrix const &rays, ZMatrix const &facetCandidates, SymmetricComplex const &theComplex/*, int linealityDim*/)
347 {
348   set<SymmetricComplex::Cone> ret;
349 
350   for(int i=0;i<facetCandidates.getHeight();i++)
351     {
352       set<int> indices;
353 
354       bool notAll=false;
355       for(unsigned j=0;j<theCone.indices.size();j++)
356         if(dot(rays[theCone.indices[j]].toVector(),facetCandidates[i].toVector()).sign()==0)
357           indices.insert(theCone.indices[j]);
358         else
359           notAll=true;
360 
361       SymmetricComplex::Cone temp(indices,theCone.dimension-1,Integer(),false,theComplex);
362       /*      temp.multiplicity=0;
363       temp.dimension=theCone.dimension-1;
364       temp.setIgnoreSymmetry(true);
365       */
366       if(notAll)ret.insert(temp);
367 
368     }
369   //  fprintf(Stderr,"HEJ!!!!\n");
370 
371   list<SymmetricComplex::Cone> ret2;
372   for(set<SymmetricComplex::Cone>::const_iterator i=ret.begin();i!=ret.end();i++)
373     {
374       bool isMaximal=true;
375 
376       /*      if(i->indices.size()+linealityDim<i->dimension)//#3
377         isMaximal=false;
378         else*/
379         for(set<SymmetricComplex::Cone>::const_iterator j=ret.begin();j!=ret.end();j++)
380           if(i!=j && i->isSubsetOf(*j))
381             {
382               isMaximal=false;
383               break;
384             }
385       if(isMaximal)
386         {
387           SymmetricComplex::Cone temp(i->indexSet(),i->dimension,i->multiplicity,true,theComplex);
388           temp.setKnownToBeNonMaximal(); // THIS IS WHERE WE SET THE CONES NON-MAXIMAL FLAG
389           //      temp.setIgnoreSymmetry(false);
390           ret2.push_back(temp);
391         }
392     }
393   return ret2;
394 }
395 
addFacesToSymmetricComplex(SymmetricComplex & c,ZCone const & cone,ZMatrix const & facetCandidates,ZMatrix const & generatorsOfLinealitySpace)396 void addFacesToSymmetricComplex(SymmetricComplex &c, ZCone const &cone, ZMatrix const &facetCandidates, ZMatrix const &generatorsOfLinealitySpace)
397 {
398   // ZMatrix const &rays=c.getVertices();
399   std::set<int> indices;
400 
401   // for(int j=0;j<rays.getHeight();j++)if(cone.contains(rays[j]))indices.insert(j);
402 
403   ZMatrix l=cone.extremeRays(&generatorsOfLinealitySpace);
404   for(int i=0;i<l.getHeight();i++)indices.insert(c.indexOfVertex(l[i]));
405 
406   addFacesToSymmetricComplex(c,indices,facetCandidates,cone.dimension(),cone.getMultiplicity());
407 }
408 
addFacesToSymmetricComplex(SymmetricComplex & c,std::set<int> const & indices,ZMatrix const & facetCandidates,int dimension,Integer multiplicity)409 void addFacesToSymmetricComplex(SymmetricComplex &c, std::set<int> const &indices, ZMatrix const &facetCandidates, int dimension, Integer multiplicity)
410 {
411   ZMatrix const &rays=c.getVertices();
412   list<SymmetricComplex::Cone> clist;
413   {
414     SymmetricComplex::Cone temp(indices,dimension,multiplicity,true,c);
415     //    temp.dimension=cone.dimension();
416     //   temp.multiplicity=cone.getMultiplicity();
417     clist.push_back(temp);
418   }
419 
420   //  int linealityDim=cone.dimensionOfLinealitySpace();
421 
422   while(!clist.empty())
423     {
424       SymmetricComplex::Cone &theCone=clist.front();
425 
426       if(!c.contains(theCone))
427         {
428           c.insert(theCone);
429           list<SymmetricComplex::Cone> facets=computeFacets(theCone,rays,facetCandidates,c/*,linealityDim*/);
430           clist.splice(clist.end(),facets);
431         }
432       clist.pop_front();
433     }
434 
435 }
436 
437 #if 0
438 /**
439    Produce strings that express the vectors in terms of rays of the fan modulo the lineality space. Symmetry is ignored??
440  */
441 vector<string> PolyhedralFan::renamingStrings(IntegerVectorList const &theVectors, IntegerVectorList const &originalRays, IntegerVectorList const &linealitySpace, SymmetryGroup *sym)const
442 {
443   vector<string> ret;
444   for(IntegerVectorList::const_iterator i=theVectors.begin();i!=theVectors.end();i++)
445     {
446       for(PolyhedralConeList::const_iterator j=cones.begin();j!=cones.end();j++)
447         {
448           if(j->contains(*i))
449             {
450               vector<int> relevantIndices;
451               IntegerVectorList relevantRays=linealitySpace;
452               int K=0;
453               for(IntegerVectorList::const_iterator k=originalRays.begin();k!=originalRays.end();k++,K++)
454                 if(j->contains(*k))
455                   {
456                     relevantIndices.push_back(K);
457                     relevantRays.push_back(*k);
458                   }
459 
460               FieldMatrix LFA(Q,relevantRays.size(),n);
461               int J=0;
462               for(IntegerVectorList::const_iterator j=relevantRays.begin();j!=relevantRays.end();j++,J++)
463                 LFA[J]=integerVectorToFieldVector(*j,Q);
464               FieldVector LFB=concatenation(integerVectorToFieldVector(*i,Q),FieldVector(Q,relevantRays.size()));
465               LFA=LFA.transposed();
466               FieldVector LFX=LFA.solver().canonicalize(LFB);
467               stringstream s;
468               if(LFX.subvector(0,n).isZero())
469                 {
470                   s<<"Was:";
471                   FieldVector S=LFX.subvector(n+linealitySpace.size(),LFX.size());
472                   for(int k=0;k<S.size();k++)
473                     if(!S[k].isZero())
474                       s<<"+"<<S[k].toString()<<"*["<<relevantIndices[k]<<"] ";
475                 }
476               ret.push_back(s.str());
477               break;
478             }
479         }
480     }
481   return ret;
482 }
483 #endif
484 
toSymmetricComplex() const485 SymmetricComplex PolyhedralFan::toSymmetricComplex()const
486 {
487 
488           ZMatrix rays=getRaysInPrintingOrder();
489 
490           ZMatrix generatorsOfLinealitySpace=cones.empty()?ZMatrix::identity(getAmbientDimension()):cones.begin()->generatorsOfLinealitySpace();
491           SymmetricComplex symCom(rays,generatorsOfLinealitySpace,symmetries);
492 
493           for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
494             {
495               addFacesToSymmetricComplex(symCom,*i,i->getFacets(),generatorsOfLinealitySpace);
496             }
497 
498 //          log1 cerr<<"Remapping";
499           symCom.remap();
500 //          log1 cerr<<"Done remapping";
501 
502           return symCom;
503 }
504 
toString(int flags) const505 std::string PolyhedralFan::toString(int flags)const
506 //void PolyhedralFan::printWithIndices(class Printer *p, bool printMultiplicities, SymmetryGroup *sym, bool group, bool ignoreCones, bool xml, bool tPlaneSort, vector<string> const *comments)const
507 {
508   stringstream ret;
509 
510   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
511     {
512       ret<<"Cone\n"<<endl;
513     ret<<*i;
514     }  return ret.str();
515 #if 0
516   PolymakeFile polymakeFile;
517   polymakeFile.create("NONAME","PolyhedralFan","PolyhedralFan",flags&FPF_xml);
518 
519 //  if(!sym)sym=&symm;
520 
521   if(cones.empty())
522     {
523 //      p->printString("Polyhedral fan is empty. Printing not supported.\n");
524       ret<<"Polyhedral fan is empty. Printing not supported.\n";
525       return ret.str();
526     }
527 
528   int h=cones.begin()->dimensionOfLinealitySpace();
529 
530 //  log1 fprintf(Stderr,"Computing rays.\n");
531   ZMatrix rays=getRaysInPrintingOrder();
532 
533   SymmetricComplex symCom(rays,cones.begin()->generatorsOfLinealitySpace(),symmetries);
534 
535   polymakeFile.writeCardinalProperty("AMBIENT_DIM",n);
536   polymakeFile.writeCardinalProperty("DIM",getMaxDimension());
537   polymakeFile.writeCardinalProperty("LINEALITY_DIM",h);
538   polymakeFile.writeMatrixProperty("RAYS",rays,true,comments);
539   polymakeFile.writeCardinalProperty("N_RAYS",rays.size());
540   IntegerVectorList linealitySpaceGenerators=highestDimensionalCone().linealitySpace().dualCone().getEquations();
541   polymakeFile.writeMatrixProperty("LINEALITY_SPACE",rowsToIntegerMatrix(linealitySpaceGenerators,n));
542   polymakeFile.writeMatrixProperty("ORTH_LINEALITY_SPACE",rowsToIntegerMatrix(highestDimensionalCone().linealitySpace().getEquations(),n));
543 
544   if(flags & FPF_primitiveRays)
545   {
546          ZMatrix primitiveRays;
547          for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
548                  for(PolyhedralConeList::const_iterator j=cones.begin();j!=cones.end();j++)
549                          if(j->contains(*i)&&(j->dimensionOfLinealitySpace()+1==j->dimension()))
550                                          primitiveRays.push_back(j->semiGroupGeneratorOfRay());
551 
552           polymakeFile.writeMatrixProperty("PRIMITIVE_RAYS",rowsToIntegerMatrix(primitiveRays,n));
553   }
554 
555 
556   ZMatrix generatorsOfLinealitySpace=cones.begin()->generatorsOfLinealitySpace();
557 
558   log1 fprintf(Stderr,"Building symmetric complex.\n");
559   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
560     {
561       {
562         static int t;
563 //        log1 fprintf(Stderr,"Adding faces of cone %i\n",t++);
564       }
565 //      log2 fprintf(Stderr,"Dim: %i\n",i->dimension());
566 
567       addFacesToSymmetricComplex(symCom,*i,i->getHalfSpaces(),generatorsOfLinealitySpace);
568     }
569 
570 //  log1 cerr<<"Remapping";
571   symCom.remap();
572 //  log1 cerr<<"Done remapping";
573 
574 
575   PolyhedralFan f=*this;
576 
577 //  log1 fprintf(Stderr,"Computing f-vector.\n");
578   ZVector fvector=symCom.fvector();
579   polymakeFile.writeCardinalVectorProperty("F_VECTOR",fvector);
580 //  log1 fprintf(Stderr,"Done computing f-vector.\n");
581 
582   if(flags&FPF_boundedInfo)
583     {
584 //      log1 fprintf(Stderr,"Computing bounded f-vector.\n");
585       ZVector fvectorBounded=symCom.fvector(true);
586       polymakeFile.writeCardinalVectorProperty("F_VECTOR_BOUNDED",fvectorBounded);
587 //      log1 fprintf(Stderr,"Done computing bounded f-vector.\n");
588     }
589   {
590     Integer euler;
591     int mul=-1;
592     for(int i=0;i<fvector.size();i++,mul*=-1)euler+=Integer(mul)*fvector[i];
593     polymakeFile.writeCardinalProperty("MY_EULER",euler);
594   }
595 
596 //  log1 fprintf(Stderr,"Checking if complex is simplicial and pure.\n");
597   polymakeFile.writeCardinalProperty("SIMPLICIAL",symCom.isSimplicial());
598   polymakeFile.writeCardinalProperty("PURE",symCom.isPure());
599 //  log1 fprintf(Stderr,"Done checking.\n");
600 
601 
602   if(flags&FPF_conesCompressed)
603   {
604 //    log1 fprintf(Stderr,"Producing list of cones up to symmetry.\n");
605     polymakeFile.writeStringProperty("CONES_ORBITS",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),false,flags&FPF_group,0,true,flags&FPF_tPlaneSort));
606 //    log1 fprintf(Stderr,"Done producing list of cones up to symmetry.\n");
607 //    log1 fprintf(Stderr,"Producing list of maximal cones up to symmetry.\n");
608     stringstream multiplicities;
609     polymakeFile.writeStringProperty("MAXIMAL_CONES_ORBITS",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,flags&FPF_group, &multiplicities,true,flags&FPF_tPlaneSort));
610     if(flags&FPF_multiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES_ORBITS",multiplicities.str());
611 //    log1 fprintf(Stderr,"Done producing list of maximal cones up to symmetry.\n");
612   }
613 
614   if(flags&FPF_conesExpanded)
615     {
616       if(flags&FPF_cones)
617         {
618 //          log1 fprintf(Stderr,"Producing list of cones.\n");
619           polymakeFile.writeStringProperty("CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),false,flags&FPF_group,0,false,flags&FPF_tPlaneSort));
620 //          log1 fprintf(Stderr,"Done producing list of cones.\n");
621         }
622       if(flags&FPF_maximalCones)
623         {
624 //          log1 fprintf(Stderr,"Producing list of maximal cones.\n");
625           stringstream multiplicities;
626           polymakeFile.writeStringProperty("MAXIMAL_CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,flags&FPF_group, &multiplicities,false,flags&FPF_tPlaneSort));
627           if(flags&FPF_multiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES",multiplicities.str());
628 //          log1 fprintf(Stderr,"Done producing list of maximal cones.\n");
629         }
630     }
631 #endif
632   #if 0
633   if(flags&FPF_values)
634     {
635       {
636         ZMatrix values;
637         for(int i=0;i<linealitySpaceGenerators.getHeight();i++)
638           {
639             ZVector v(1);
640             v[0]=evaluatePiecewiseLinearFunction(linealitySpaceGenerators[i]);
641             values.appendRow(v);
642           }
643         polymakeFile.writeMatrixProperty("LINEALITY_VALUES",rowsToIntegerMatrix(values,1));
644       }
645       {
646         ZMatrix values;
647         for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
648           {
649             ZVector v(1);
650             v[0]=evaluatePiecewiseLinearFunction(*i);
651             values.push_back(v);
652           }
653         polymakeFile.writeMatrixProperty("RAY_VALUES",rowsToIntegerMatrix(values,1));
654       }
655     }
656 #endif
657 
658 
659 //  log1 fprintf(Stderr,"Producing final string for output.\n");
660 /*  stringstream s;
661   polymakeFile.writeStream(s);
662   string S=s.str();
663 //  log1 fprintf(Stderr,"Printing string.\n");
664   p->printString(S.c_str());
665 *///  log1 fprintf(Stderr,"Done printing string.\n");
666 }
667 
668 #if 0
669 PolyhedralFan PolyhedralFan::readFan(string const &filename, bool onlyMaximal, IntegerVector *w, set<int> const *coneIndices, SymmetryGroup const *sym, bool readCompressedIfNotSym)
670 {
671     PolymakeFile inFile;
672     inFile.open(filename.c_str());
673 
674     int n=inFile.readCardinalProperty("AMBIENT_DIM");
675     int nRays=inFile.readCardinalProperty("N_RAYS");
676     IntegerMatrix rays=inFile.readMatrixProperty("RAYS",nRays,n);
677     int linealityDim=inFile.readCardinalProperty("LINEALITY_DIM");
678     IntegerMatrix linealitySpace=inFile.readMatrixProperty("LINEALITY_SPACE",linealityDim,n);
679 
680 
681     const char *sectionName=0;
682     const char *sectionNameMultiplicities=0;
683     if(sym || readCompressedIfNotSym)
684     {
685       sectionName=(onlyMaximal)?"MAXIMAL_CONES_ORBITS":"CONES_ORBITS";
686       sectionNameMultiplicities=(onlyMaximal)?"MULTIPLICITIES_ORBITS":"DUMMY123";
687     }
688       else
689       {  sectionName=(onlyMaximal)?"MAXIMAL_CONES":"CONES";
690       sectionNameMultiplicities=(onlyMaximal)?"MULTIPLICITIES":"DUMMY123";
691       }
692 
693 
694     IntegerVector w2(n);
695     if(w==0)w=&w2;
696 
697     SymmetryGroup sym2(n);
698     if(sym==0)sym=&sym2;
699 
700     vector<list<int> > cones=inFile.readMatrixIncidenceProperty(sectionName);
701     IntegerVectorList r;
702 
703     bool hasMultiplicities=inFile.hasProperty(sectionNameMultiplicities);
704     IntegerMatrix multiplicities(0,0);
705     if(hasMultiplicities)multiplicities=inFile.readMatrixProperty(sectionNameMultiplicities,cones.size(),1);
706 
707 
708     PolyhedralFan ret(n);
709 
710     log2 cerr<< "Number of orbits to expand "<<cones.size()<<endl;
711     for(int i=0;i<cones.size();i++)
712       if(coneIndices==0 || coneIndices->count(i))
713         {
714           log2 cerr<<"Expanding symmetries of cone"<<endl;
715           {
716             IntegerVectorList coneRays;
717             for(list<int>::const_iterator j=cones[i].begin();j!=cones[i].end();j++)
718               coneRays.push_back((rays[*j]));
719             PolyhedralCone C=PolyhedralCone::givenByRays(coneRays,linealitySpace.getRows(),n);
720             if(hasMultiplicities)C.setMultiplicity(multiplicities[i][0]);
721             for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
722               {
723                 if(C.contains(SymmetryGroup::composeInverse(*perm,*w)))
724                   {
725                     PolyhedralCone C2=C.permuted(*perm);
726                     C2.canonicalize();
727                     ret.insert(C2);
728                   }
729               }
730           }
731         }
732     return ret;
733 }
734 #endif
735 
736 #if 0
737 IncidenceList PolyhedralFan::getIncidenceList(SymmetryGroup *sym)const //fan must be pure
738 {
739   IncidenceList ret;
740   SymmetryGroup symm(n);
741   if(!sym)sym=&symm;
742   assert(!cones.empty());
743   int h=cones.begin()->dimensionOfLinealitySpace();
744   IntegerVectorList rays=getRaysInPrintingOrder(sym);
745   PolyhedralFan f=*this;
746 
747   while(f.getMaxDimension()!=h)
748     {
749       IntegerVectorList l;
750       PolyhedralFan done(n);
751       IntegerVector rayIncidenceCounter(rays.size());
752       int faceIndex =0;
753       for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
754         {
755           if(!done.contains(*i))
756             {
757               for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
758                 {
759                   PolyhedralCone cone=i->permuted(*k);
760                   if(!done.contains(cone))
761                     {
762                       int rayIndex=0;
763                       IntegerVector indices(0);
764                       for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
765                         {
766                           if(cone.contains(*j))
767                             {
768                               indices.grow(indices.size()+1);
769                               indices[indices.size()-1]=rayIndex;
770                               rayIncidenceCounter[rayIndex]++;
771                             }
772                           rayIndex++;
773                         }
774                       l.push_back(indices);
775                       faceIndex++;
776                       done.insert(cone);
777                     }
778                 }
779             }
780         }
781       ret[f.getMaxDimension()]=l;
782       f=f.facetComplex();
783     }
784   return ret;
785 }
786 #endif
787 
makePure()788 void PolyhedralFan::makePure()
789 {
790   if(getMaxDimension()!=getMinDimension())removeAllLowerDimensional();
791 }
792 
contains(ZCone const & c) const793 bool PolyhedralFan::contains(ZCone const &c)const
794 {
795   return cones.count(c);
796 }
797 
798 
799 #if 0
800 PolyhedralCone PolyhedralFan::coneContaining(ZVector const &v)const
801 {
802   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
803     if(i->contains(v))return i->faceContaining(v);
804   debug<<"Vector "<<v<<" not contained in support of fan\n";
805   assert(0);
806 }
807 #endif
808 
conesBegin() const809 PolyhedralFan::coneIterator PolyhedralFan::conesBegin()const
810 {
811   return cones.begin();
812 }
813 
814 
conesEnd() const815 PolyhedralFan::coneIterator PolyhedralFan::conesEnd()const
816 {
817   return cones.end();
818 }
819 
820 
821 
link(ZVector const & w,SymmetryGroup * sym) const822 PolyhedralFan PolyhedralFan::link(ZVector const &w, SymmetryGroup *sym)const
823 {
824   SymmetryGroup symL(n);
825   if(!sym)sym=&symL;
826 
827   PolyhedralFan ret(n);
828 
829   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
830     {
831       for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
832         {
833           ZVector w2=perm->applyInverse(w);
834           if(i->contains(w2))
835             {
836               ret.insert(i->link(w2));
837             }
838         }
839     }
840   return ret;
841 }
842 
link(ZVector const & w) const843 PolyhedralFan PolyhedralFan::link(ZVector const &w)const
844 {
845   PolyhedralFan ret(n);
846 
847   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
848     {
849       if(i->contains(w))
850         {
851           ret.insert(i->link(w));
852         }
853     }
854   return ret;
855 }
856 
857 
size() const858 int PolyhedralFan::size()const
859 {
860   return cones.size();
861 }
862 
dimensionOfLinealitySpace() const863 int PolyhedralFan::dimensionOfLinealitySpace()const
864 {
865   assert(cones.size());//slow!
866   return cones.begin()->dimensionOfLinealitySpace();
867 }
868 
869 
870 
871 
removeNonMaximal()872 void PolyhedralFan::removeNonMaximal()
873 {
874   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();)
875     {
876       ZVector w=i->getRelativeInteriorPoint();
877       bool containedInOther=false;
878       for(PolyhedralConeList::iterator j=cones.begin();j!=cones.end();j++)
879         if(j!=i)
880           {
881             if(j->contains(w)){containedInOther=true;break;}
882           }
883       if(containedInOther)
884         {
885           PolyhedralConeList::iterator k=i;
886           i++;
887           cones.erase(k);
888         }
889       else i++;
890     }
891 }
892 
893 
894 }
895