← Back to branch summary

~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

« back to all changes in this revision

Viewing changes to gfanlib_polyhedralfan.cpp

  • Committer: Package Import Robot
  • Author(s): Cédric Boutillier
  • Date: 2013-07-09 10:44:01 UTC
  • mfrom: (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130709104401-5q66ozz5j5af0dak
Tags: 0.5+dfsg-3
* Upload to unstable.
* modify remove_failing_tests_on_32bits.patch to replace command of
  0009RenderStairCase test with an empty one instead of deleting it.
* remove lintian override about spelling error

Show diffs side-by-side

added added

removed removed

Lines of Context:
gfanlib_polyhedralfan.cpp
 
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
 
 
12
using namespace std;
 
13
namespace gfan{
 
14
 
 
15
PolyhedralFan::PolyhedralFan(int ambientDimension):
 
16
  n(ambientDimension),
 
17
  symmetries(n)
 
18
{
 
19
}
 
20
 
 
21
PolyhedralFan::PolyhedralFan(SymmetryGroup const &sym):
 
22
  n(sym.sizeOfBaseSet()),
 
23
  symmetries(sym)
 
24
{
 
25
 
 
26
}
 
27
 
 
28
PolyhedralFan 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
 
 
40
PolyhedralFan 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
 
 
61
int PolyhedralFan::getAmbientDimension()const
 
62
{
 
63
  return n;
 
64
}
 
65
 
 
66
bool PolyhedralFan::isEmpty()const
 
67
{
 
68
  return cones.empty();
 
69
}
 
70
 
 
71
int PolyhedralFan::getMaxDimension()const
 
72
{
 
73
  assert(!cones.empty());
 
74
 
 
75
  return cones.begin()->dimension();
 
76
}
 
77
 
 
78
int PolyhedralFan::getMinDimension()const
 
79
{
 
80
  assert(!cones.empty());
 
81
 
 
82
  return cones.rbegin()->dimension();
 
83
}
 
84
 
 
85
/*
 
86
PolyhedralFan 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
/*
 
123
PolyhedralFan 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
*/
 
163
void PolyhedralFan::insert(ZCone const &c)
 
164
{
 
165
  ZCone temp=c;
 
166
  temp.canonicalize();
 
167
  cones.insert(temp);
 
168
}
 
169
 
 
170
void PolyhedralFan::remove(ZCone const &c)
 
171
{
 
172
  cones.erase(c);
 
173
}
 
174
 
 
175
/*
 
176
void 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
 
 
189
void 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
 
 
201
PolyhedralFan 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
/*
 
218
PolyhedralFan 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
 
242
PolyhedralFan 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
 
 
303
ZMatrix 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);
 
312
ZMatrix 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
 
 
395
void 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
 
 
408
void 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
 */
 
442
vector<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
 
 
486
SymmetricComplex 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
 
 
513
std::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
 
677
PolyhedralFan 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
 
745
IncidenceList 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
 
 
796
void PolyhedralFan::makePure()
 
797
{
 
798
  if(getMaxDimension()!=getMinDimension())removeAllLowerDimensional();
 
799
}
 
800
 
 
801
bool PolyhedralFan::contains(ZCone const &c)const
 
802
{
 
803
  return cones.count(c);
 
804
}
 
805
 
 
806
 
 
807
#if 0
 
808
PolyhedralCone 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
 
 
817
PolyhedralFan::coneIterator PolyhedralFan::conesBegin()const
 
818
{
 
819
  return cones.begin();
 
820
}
 
821
 
 
822
 
 
823
PolyhedralFan::coneIterator PolyhedralFan::conesEnd()const
 
824
{
 
825
  return cones.end();
 
826
}
 
827
 
 
828
 
 
829
 
 
830
PolyhedralFan 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
 
 
851
PolyhedralFan 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
 
 
866
int PolyhedralFan::size()const
 
867
{
 
868
  return cones.size();
 
869
}
 
870
 
 
871
int PolyhedralFan::dimensionOfLinealitySpace()const
 
872
{
 
873
  assert(cones.size());//slow!
 
874
  return cones.begin()->dimensionOfLinealitySpace();
 
875
}
 
876
 
 
877
 
 
878
 
 
879
 
 
880
void 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
}