source: git/gfanlib/gfanlib_polyhedralfan.cpp

spielwiese
Last change on this file was 5443c1, checked in by Hans Schoenemann <hannes@…>, 5 years ago
porting: gcc-5 (Ubuntu 16.04.5 LTS)
  • 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 <stddef.h>
9#include <sstream>
10#include "gfanlib_polyhedralfan.h"
11#include "gfanlib_polymakefile.h"
12
13using namespace std;
14namespace gfan{
15
16PolyhedralFan::PolyhedralFan(int ambientDimension):
17  n(ambientDimension),
18  symmetries(n)
19{
20}
21
22PolyhedralFan::PolyhedralFan(SymmetryGroup const &sym):
23  n(sym.sizeOfBaseSet()),
24  symmetries(sym)
25{
26
27}
28
29PolyhedralFan 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
41PolyhedralFan 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
62int PolyhedralFan::getAmbientDimension()const
63{
64  return n;
65}
66
67bool PolyhedralFan::isEmpty()const
68{
69  return cones.empty();
70}
71
72int PolyhedralFan::getMaxDimension()const
73{
74  assert(!cones.empty());
75
76  return cones.begin()->dimension();
77}
78
79int PolyhedralFan::getMinDimension()const
80{
81  assert(!cones.empty());
82
83  return cones.rbegin()->dimension();
84}
85
86/*
87PolyhedralFan 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/*
124PolyhedralFan 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*/
164void PolyhedralFan::insert(ZCone const &c)
165{
166  ZCone temp=c;
167  temp.canonicalize();
168  cones.insert(temp);
169}
170
171void PolyhedralFan::remove(ZCone const &c)
172{
173  cones.erase(c);
174}
175
176/*
177void 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
190void 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
202PolyhedralFan 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/*
219PolyhedralFan 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
243PolyhedralFan 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
304ZMatrix 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);
312ZMatrix 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   */
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
396void 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
409void 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 */
441vector<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
485SymmetricComplex 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
505std::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
669PolyhedralFan 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
737IncidenceList 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
788void PolyhedralFan::makePure()
789{
790  if(getMaxDimension()!=getMinDimension())removeAllLowerDimensional();
791}
792
793bool PolyhedralFan::contains(ZCone const &c)const
794{
795  return cones.count(c);
796}
797
798
799#if 0
800PolyhedralCone 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
809PolyhedralFan::coneIterator PolyhedralFan::conesBegin()const
810{
811  return cones.begin();
812}
813
814
815PolyhedralFan::coneIterator PolyhedralFan::conesEnd()const
816{
817  return cones.end();
818}
819
820
821
822PolyhedralFan 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
843PolyhedralFan 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
858int PolyhedralFan::size()const
859{
860  return cones.size();
861}
862
863int PolyhedralFan::dimensionOfLinealitySpace()const
864{
865  assert(cones.size());//slow!
866  return cones.begin()->dimensionOfLinealitySpace();
867}
868
869
870
871
872void 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}
Note: See TracBrowser for help on using the repository browser.