source: git/gfanlib/gfanlib_polyhedralfan.cpp @ 43317d

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