source: git/gfanlib/gfanlib_polyhedralfan.cpp @ 75fa84

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