← Back to branch summary

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

« back to all changes in this revision

Viewing changes to 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:
polyhedralfan.cpp
8
8
#include "symmetry.h"
9
9
#include "polymakefile.h"
10
10
#include "symmetriccomplex.h"
 
11
#include "linalg.h"
 
12
#include "lp.h"
 
13
#include "codimoneconnectedness.h"
 
14
#include "symmetrictraversal.h"
 
15
#include "traverser_groebnerfan.h"
11
16
#include "log.h"
12
17
 
13
18
static Timer polyhedralFanRefinementTimer("Polyhedral fan refinement",1);
22
27
{
23
28
  PolyhedralFan ret(n);
24
29
 
25
 
  ret.cones.insert(PolyhedralCone(n));
 
30
  PolyhedralCone temp(n);
 
31
  temp.canonicalize();
 
32
  ret.cones.insert(temp);
26
33
 
27
34
  return ret;
28
35
}
65
72
  return ret;
66
73
}
67
74
 
 
75
PolyhedralFan PolyhedralFan::complementOfCone(PolyhedralCone const &c, bool includec)
 
76
{
 
77
  PolyhedralCone C=c;
 
78
  C.canonicalize();
 
79
  IntegerVectorList inequalities=C.getHalfSpaces();
 
80
  IntegerVectorList equations=C.getEquations();
 
81
 
 
82
  for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
 
83
    inequalities.push_back(*i);
 
84
 
 
85
  int n=C.ambientDimension();
 
86
 
 
87
  PolyhedralFan ret=PolyhedralFan::fullSpace(n);
 
88
  for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
 
89
    {
 
90
      PolyhedralFan temp(C.ambientDimension());
 
91
      for(int j=-1;j<2;j+=2)
 
92
        {
 
93
          IntegerVectorList inequality;
 
94
          inequality.push_back(j* *i);
 
95
          IntegerVectorList empty;
 
96
          PolyhedralCone tempC(inequality,empty,n);
 
97
          tempC.canonicalize();
 
98
          temp.insert(tempC);
 
99
        }
 
100
      ret=refinement(ret,temp);
 
101
    }
 
102
  if(!includec)ret.remove(C);
 
103
  return ret;
 
104
}
68
105
 
69
106
PolyhedralFan PolyhedralFan::bergmanOfPrincipalIdeal(Polynomial const &p1)
70
107
{
99
136
              equalities.push_back(*j);
100
137
              PolyhedralCone c=PolyhedralCone(inequalities,equalities).withLastCoordinateRemoved();
101
138
              c.canonicalize();
 
139
              c.setLinearForm(i->begin()->getMarked().m.exponent.subvector(0,p1.getNumberOfVariables()));
 
140
              c.setMultiplicity(gcdOfVector(j->subvector(0,j->size()-1)));
102
141
              ret.cones.insert(c);
103
142
            }
104
143
        }
114
153
  PolynomialRing theSecondRing=theRing.withVariablesAppended("H");
115
154
  Polynomial p=p1.homogenization(theSecondRing);
116
155
  PolyhedralFan ret(p1.getNumberOfVariables());
117
 
  PolynomialSet g(theRing);
 
156
  PolynomialSet g(theSecondRing);
 
157
  //  PolynomialSet g(theRing);//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
118
158
  g.push_back(p);
119
159
  buchberger(&g,LexicographicTermOrder());
120
160
 
121
161
  EnumerationTargetCollector gfan;
 
162
/*  {//old enumeration strategy
122
163
  LexicographicTermOrder myTermOrder;
123
164
  ReverseSearch rs(myTermOrder);
124
165
 
125
166
  rs.setEnumerationTarget(&gfan);
126
167
 
127
168
  rs.enumerate(g);
 
169
  }
 
170
*/
 
171
  {//new enumeration strategy
 
172
    GroebnerFanTraverser traverser(g);
 
173
    traverser.setIsKnownToBeComplete(true);
 
174
    TargetGlue target(gfan);
 
175
    symmetricTraverse(traverser,target);
 
176
  }
128
177
 
129
178
  PolynomialSetList theList=gfan.getList();
130
179
  for(PolynomialSetList::const_iterator i=theList.begin();i!=theList.end();i++)
134
183
 
135
184
      PolyhedralCone c=PolyhedralCone(inequalities,equalities).withLastCoordinateRemoved();
136
185
      c.canonicalize();
 
186
      c.setLinearForm(i->begin()->getMarked().m.exponent.subvector(0,c.ambientDimension()));
137
187
      ret.cones.insert(c);
138
188
    }
139
189
 
157
207
      p->printPolyhedralCone(*i);
158
208
    }
159
209
  p->printString("Done printing PolyhedralFan.");
160
 
  p->printNewLine(); 
 
210
  p->printNewLine();
161
211
}
162
212
 
163
213
int PolyhedralFan::getAmbientDimension()const
183
233
 
184
234
  return cones.rbegin()->dimension();
185
235
}
186
 
 
 
236
 
187
237
PolyhedralFan refinement(const PolyhedralFan &a, const PolyhedralFan &b, int cutOffDimension, bool allowASingleConeOfCutOffDimension)
188
238
{
189
239
  TimerScope ts(&polyhedralFanRefinementTimer);
220
270
  return ret;
221
271
}
222
272
 
 
273
 
 
274
PolyhedralFan product(const PolyhedralFan &a, const PolyhedralFan &b)
 
275
{
 
276
  PolyhedralFan ret(a.getAmbientDimension()+b.getAmbientDimension());
 
277
 
 
278
  for(PolyhedralConeList::const_iterator A=a.cones.begin();A!=a.cones.end();A++)
 
279
    for(PolyhedralConeList::const_iterator B=b.cones.begin();B!=b.cones.end();B++)
 
280
      ret.insert(product(*A,*B));
 
281
 
 
282
  return ret;
 
283
}
 
284
 
 
285
 
223
286
IntegerVectorList PolyhedralFan::getRays(int dim)
224
287
{
225
288
  IntegerVectorList ret;
253
316
  cones.insert(c);
254
317
}
255
318
 
 
319
 
 
320
void PolyhedralFan::insertFacetsOfCone(PolyhedralCone const &c)
 
321
{
 
322
  PolyhedralFan facets=facetsOfCone(c);
 
323
  for(PolyhedralConeList::const_iterator i=facets.cones.begin();i!=facets.cones.end();i++)insert(*i);
 
324
}
 
325
 
 
326
 
256
327
void PolyhedralFan::remove(PolyhedralCone const &c)
257
328
{
258
329
  cones.erase(c);
298
369
}
299
370
 
300
371
 
 
372
PolyhedralFan PolyhedralFan::fullComplex()const
 
373
{
 
374
  PolyhedralFan ret=*this;
 
375
 
 
376
  while(1)
 
377
    {
 
378
      log2 debug<<"looping";
 
379
      bool doLoop=false;
 
380
      PolyhedralFan facets=ret.facetComplex();
 
381
      log2 debug<<"number of facets"<<facets.size()<<"\n";
 
382
      for(PolyhedralConeList::const_iterator i=facets.cones.begin();i!=facets.cones.end();i++)
 
383
        if(!ret.contains(*i))
 
384
          {
 
385
            ret.insert(*i);
 
386
            doLoop=true;
 
387
          }
 
388
      if(!doLoop)break;
 
389
    }
 
390
  return ret;
 
391
}
 
392
 
 
393
 
 
394
/*
301
395
PolyhedralFan PolyhedralFan::facetComplexSymmetry(SymmetryGroup const &sym, bool keepRays, bool dropLinealitySpace)const
302
396
{
303
397
  log1 fprintf(Stderr,"Computing facet complex...\n");
332
426
  log1 fprintf(Stderr,"Done computing facet complex.\n");
333
427
  return ret;
334
428
}
 
429
*/
 
430
 
 
431
PolyhedralFan PolyhedralFan::facetComplexSymmetry(SymmetryGroup const &sym, bool keepRays, bool dropLinealitySpace)const
 
432
{
 
433
  log1 fprintf(Stderr,"Computing facet complex...\n");
 
434
  PolyhedralFan ret(n);
 
435
 
 
436
  vector<IntegerVector> relIntTable;
 
437
  vector<int> dimensionTable;
 
438
 
 
439
  if(keepRays)
 
440
    for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
 
441
      if(i->dimension()==i->dimensionOfLinealitySpace()+1)
 
442
        {
 
443
          relIntTable.push_back(i->getRelativeInteriorPoint());
 
444
          dimensionTable.push_back(i->dimension());
 
445
          ret.insert(*i);
 
446
        }
 
447
 
 
448
  for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
 
449
    {
 
450
      int iDim=i->dimension();
 
451
      if(dropLinealitySpace && (i->dimension()==i->dimensionOfLinealitySpace()+1))continue;
 
452
 
 
453
      //      i->findFacets();
 
454
      IntegerVectorList normals=i->getHalfSpaces();
 
455
      for(IntegerVectorList::const_iterator j=normals.begin();j!=normals.end();j++)
 
456
        {
 
457
          bool alreadyInRet=false;
 
458
          for(int l=0;l<relIntTable.size();l++)
 
459
            {
 
460
              if(dimensionTable[l]==iDim-1)
 
461
                for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
 
462
                  {
 
463
                    IntegerVector u=SymmetryGroup::compose(*k,relIntTable[l]);
 
464
                    if((dotLong(*j,u)==0) && (i->contains(u)))
 
465
                      {
 
466
                        alreadyInRet=true;
 
467
                        goto exitLoop;
 
468
                      }
 
469
                  }
 
470
            }
 
471
        exitLoop:
 
472
          if(!alreadyInRet)
 
473
            {
 
474
              IntegerVectorList equations=i->getEquations();
 
475
              IntegerVectorList inequalities=i->getHalfSpaces();
 
476
              equations.push_back(*j);
 
477
              PolyhedralCone c(inequalities,equations,n);
 
478
              c.canonicalize();
 
479
              ret.insert(c);
 
480
              relIntTable.push_back(c.getRelativeInteriorPoint());
 
481
              dimensionTable.push_back(c.dimension());
 
482
            }
 
483
        }
 
484
    }
 
485
  log1 fprintf(Stderr,"Done computing facet complex.\n");
 
486
  return ret;
 
487
}
335
488
 
336
489
/*
337
490
IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup *sym)const
363
516
*/
364
517
 
365
518
 
366
 
static IntegerVector stableRay(PolyhedralCone const &c, SymmetryGroup *sym)
 
519
IntegerVector PolyhedralFan::stableRay(PolyhedralCone const &c, SymmetryGroup const *sym)
367
520
{
368
521
  PolyhedralCone C=c;//cast away const instead?
 
522
 
369
523
  IntegerVector v=C.getRelativeInteriorPoint();
 
524
 
370
525
  IntegerVector ret(v.size());
371
526
  for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
372
527
    {
376
531
  return normalized(ret);
377
532
}
378
533
 
 
534
 
 
535
IntegerVector PolyhedralFan::stableRay(IntegerVector const &v, IntegerVectorList const &equations, IntegerVectorList const &inequalities, SymmetryGroup const *sym)
 
536
{
 
537
  IntegerVector ret(v.size());
 
538
  for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
 
539
    {
 
540
      IntegerVector v2=SymmetryGroup::compose(*k,v);
 
541
      bool containsV2=true;
 
542
 
 
543
      for(IntegerVectorList::const_iterator l=equations.begin();l!=equations.end();l++)
 
544
        if(dotLong(*l,v2)!=0)
 
545
          {
 
546
            containsV2=false;
 
547
            goto leave;
 
548
          }
 
549
      for(IntegerVectorList::const_iterator l=inequalities.begin();l!=inequalities.end();l++)
 
550
        if(dotLong(*l,v2)<0)
 
551
          {
 
552
            containsV2=false;
 
553
            goto leave;
 
554
          }
 
555
    leave:
 
556
      if(containsV2)ret+=v2;
 
557
    }
 
558
  return normalized(ret);
 
559
}
 
560
 
 
561
/* Slow version using facetComplexSymmetry()
379
562
IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup *sym)const
380
563
{
381
564
  assert(!cones.empty());
395
578
          bool found=false;
396
579
          for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
397
580
            if(i->contains(*j))found=true;
398
 
          
 
581
 
399
582
          if(!found)
400
583
            {
401
584
              //IntegerVector interiorPointForOrbit=i->getRelativeInteriorPoint();
415
598
        }
416
599
    }
417
600
  return rays;
418
 
}
 
601
  }*/
 
602
 
 
603
PolyhedralFan PolyhedralFan::rayComplexSymmetry(SymmetryGroup const &sym)const
 
604
{
 
605
  //  log0 fprintf(Stderr,"rayComplexSymmetry - begin\n");
 
606
  PolyhedralFan ret(n);
 
607
  log1 fprintf(Stderr,"Computing rays of %i cones\n",cones.size());
 
608
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
609
    {
 
610
      {
 
611
        static int t;
 
612
        if(!((t++)%10))log1 fprintf(Stderr,"%i\n",t);
 
613
      }
 
614
      //  log0 fprintf(Stderr,"calling\n");
 
615
      IntegerVectorList rays=i->extremeRays();
 
616
      //log0 fprintf(Stderr,"returning\n");
 
617
      for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
 
618
        {
 
619
          bool alreadyInRet=false;
 
620
          for(PolyhedralConeList::const_iterator I=ret.cones.begin();I!=ret.cones.end();I++)
 
621
            for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
 
622
              {
 
623
                IntegerVector u=SymmetryGroup::compose(*k,*j);
 
624
                if((I->contains(u)))
 
625
                  {
 
626
                    alreadyInRet=true;
 
627
                    goto exitLoop;
 
628
                  }
 
629
              }
 
630
        exitLoop:
 
631
          IntegerVectorList equations=i->getEquations();
 
632
          IntegerVectorList inequalities1=i->getHalfSpaces();
 
633
          IntegerVectorList inequalities2;
 
634
          for(IntegerVectorList::const_iterator k=inequalities1.begin();k!=inequalities1.end();k++)
 
635
            {
 
636
              if(dotLong(*j,*k))
 
637
                inequalities2.push_back(*k);
 
638
              else
 
639
                equations.push_back(*k);
 
640
            }
 
641
          PolyhedralCone C(inequalities2,equations,n);
 
642
          C.canonicalize();
 
643
          ret.insert(C);
 
644
        }
 
645
    }
 
646
  //  log0 fprintf(Stderr,"rayComplexSymmetry - end\n");
 
647
  return ret;
 
648
}
 
649
 
 
650
 
 
651
#if 0
 
652
IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup const *sym)const
 
653
{
 
654
  assert(!cones.empty());
 
655
  int h=cones.begin()->dimensionOfLinealitySpace();
 
656
 
 
657
  /*
 
658
  PolyhedralFan f=*this;
 
659
  if(f.getMaxDimension()==h)return IntegerVectorList();
 
660
  while(f.getMaxDimension()>h+1)
 
661
    {
 
662
      f=f.facetComplexSymmetry(*sym,true,true);
 
663
    }
 
664
  */
 
665
  PolyhedralFan f=rayComplexSymmetry(*sym);
 
666
  IntegerVectorList rays;
 
667
 
 
668
  log1 fprintf(Stderr,"Number of cones in RayComplex: %i\n",f.cones.size());
 
669
 
 
670
  for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
 
671
    {
 
672
      static int t;
 
673
      log1 fprintf(Stderr,"%i\n",t++);
 
674
      if(i->dimension()!=i->dimensionOfLinealitySpace())//This check is needed since the above while loop may not be run and therefore the lineality space may not have been removed.
 
675
        {
 
676
          bool found=false;
 
677
          for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
 
678
            if(i->contains(*j))
 
679
              {
 
680
                found=true;
 
681
                break;
 
682
              }
 
683
          if(!found)
 
684
            {
 
685
              //IntegerVector interiorPointForOrbit=i->getRelativeInteriorPoint();
 
686
              IntegerVector interiorPointForOrbit=stableRay(*i,sym);
 
687
              //    PolyhedralFan done(n);
 
688
 
 
689
              //Check that this works:
 
690
              set<IntegerVector> thisOrbitsRays;
 
691
 
 
692
              for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
 
693
                {
 
694
                  IntegerVector temp=SymmetryGroup::compose(*k,interiorPointForOrbit);
 
695
                  thisOrbitsRays.insert(temp);
 
696
                }
 
697
              for(set<IntegerVector>::const_iterator i=thisOrbitsRays.begin();i!=thisOrbitsRays.end();i++)rays.push_back(*i);
 
698
              //Instead of this:
 
699
              /*              for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
 
700
                {
 
701
                  bool found=false;
 
702
                  IntegerVector temp=SymmetryGroup::compose(*k,interiorPointForOrbit);
 
703
                  for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)//REWRITE WITH LOGARITHMIC SEARCH
 
704
                    if(*j==temp)
 
705
                      {
 
706
                        found=true;
 
707
                        break;
 
708
                      }
 
709
                  if(!found)
 
710
                    {
 
711
                      PolyhedralCone cone=i->permuted(*k);
 
712
                      rays.push_back(SymmetryGroup::compose(*k,interiorPointForOrbit));
 
713
                      //      done.insert(cone);
 
714
                    }
 
715
                }
 
716
              */
 
717
            }
 
718
        }
 
719
    }
 
720
  return rays;
 
721
}
 
722
 
 
723
#elseif 0
 
724
//version used until Sep 2010
 
725
IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup const *sym, bool upToSymmetry)const
 
726
{
 
727
        SymmetryGroup localsym(n);
 
728
        if(!sym)sym=&localsym;
 
729
  IntegerVectorList rays;
 
730
  log1 fprintf(Stderr,"Computing rays of %i cones\n",cones.size());
 
731
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
732
    {
 
733
      {
 
734
        static int t;
 
735
        if(!((t++)%10))log1 fprintf(Stderr,"%i\n",t);
 
736
      }
 
737
      IntegerVectorList temp=i->extremeRays();
 
738
      //      AsciiPrinter(Stderr).printVectorList(temp);
 
739
 
 
740
      for(IntegerVectorList::const_iterator j=temp.begin();j!=temp.end();j++)
 
741
        {
 
742
          IntegerVectorList equations=i->getEquations();
 
743
          IntegerVectorList inequalities1=i->getHalfSpaces();
 
744
          IntegerVectorList inequalities2;
 
745
          for(IntegerVectorList::const_iterator k=inequalities1.begin();k!=inequalities1.end();k++)
 
746
            {
 
747
              if(dotLong(*j,*k))
 
748
                inequalities2.push_back(*k);
 
749
              else
 
750
                equations.push_back(*k);
 
751
            }
 
752
          bool isFound=false;
 
753
          for(IntegerVectorList::const_iterator j2=rays.begin();j2!=rays.end();j2++)
 
754
            for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
 
755
              {
 
756
                bool isInCone=true;
 
757
                IntegerVector v=SymmetryGroup::compose(*k,*j2);
 
758
                for(IntegerVectorList::const_iterator l=equations.begin();l!=equations.end();l++)
 
759
                  if(dotLong(*l,v)!=0)
 
760
                    {
 
761
                      isInCone=false;
 
762
                      goto leave;
 
763
                    }
 
764
                for(IntegerVectorList::const_iterator l=inequalities2.begin();l!=inequalities2.end();l++)
 
765
                  if(dotLong(*l,v)<0)
 
766
                    {
 
767
                      isInCone=false;
 
768
                      goto leave;
 
769
                    }
 
770
              leave:
 
771
                if(isInCone)
 
772
                  {
 
773
                    isFound=true;
 
774
                    goto leave2;
 
775
                  }
 
776
              }
 
777
        leave2:
 
778
          if(!isFound)
 
779
            {
 
780
              IntegerVector ray=stableRay(*j,equations,inequalities2,sym);
 
781
              rays.push_back(ray);
 
782
            }
 
783
        }
 
784
    }
 
785
  rays.sort();
 
786
  if(upToSymmetry)return rays;
 
787
  IntegerVectorList ret;
 
788
  for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
 
789
    {
 
790
      set<IntegerVector> thisOrbitsRays;
 
791
      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
 
792
        {
 
793
          IntegerVector temp=SymmetryGroup::compose(*k,*i);
 
794
          thisOrbitsRays.insert(temp);
 
795
        }
 
796
      for(set<IntegerVector>::const_iterator i=thisOrbitsRays.begin();i!=thisOrbitsRays.end();i++)ret.push_back(*i);
 
797
    }
 
798
  return ret;
 
799
}
 
800
#else
 
801
IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup const *sym, bool upToSymmetry)const
 
802
{
 
803
  /*
 
804
   * The ordering changed in this version. Previously the orbit representatives stored in "rays" were
 
805
   * just the first extreme ray from the orbit that appeared. Now it is gotten using "orbitRepresentative"
 
806
   * which causes the ordering in which the different orbits appear to change.
 
807
   */
 
808
 
 
809
  if(cones.empty())return IntegerVectorList();
 
810
  IntegerVectorList generatorsOfLinealitySpace=cones.begin()->generatorsOfLinealitySpace();//all cones have the same lineality space
 
811
 
 
812
        SymmetryGroup localsym(n);
 
813
        if(!sym)sym=&localsym;
 
814
  set<IntegerVector> rays;
 
815
  log1 fprintf(Stderr,"Computing rays of %i cones\n",cones.size());
 
816
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
817
    {
 
818
      {
 
819
        static int t;
 
820
        if(!((t++)%10))log1 fprintf(Stderr,"%i\n",t);
 
821
      }
 
822
      IntegerVectorList temp=i->extremeRays(&generatorsOfLinealitySpace);
 
823
      for(IntegerVectorList::const_iterator j=temp.begin();j!=temp.end();j++)
 
824
        rays.insert(sym->orbitRepresentative(*j));
 
825
    }
 
826
  IntegerVectorList ret;
 
827
  if(upToSymmetry)
 
828
    for(set<IntegerVector>::const_iterator i=rays.begin();i!=rays.end();i++)ret.push_back(*i);
 
829
  else
 
830
    for(set<IntegerVector>::const_iterator i=rays.begin();i!=rays.end();i++)
 
831
      {
 
832
        set<IntegerVector> thisOrbitsRays;
 
833
        for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
 
834
          {
 
835
            IntegerVector temp=SymmetryGroup::compose(*k,*i);
 
836
            thisOrbitsRays.insert(temp);
 
837
          }
 
838
        for(set<IntegerVector>::const_iterator i=thisOrbitsRays.begin();i!=thisOrbitsRays.end();i++)ret.push_back(*i);
 
839
      }
 
840
  return ret;
 
841
}
 
842
 
 
843
 
 
844
#endif
419
845
 
420
846
 
421
847
/*void PolyhedralFan::printWithIndices(class Printer *p, SymmetryGroup *sym)const //fan must be pure
514
940
  PolymakeFile polymakeFile;
515
941
 
516
942
 
517
 
  
 
943
 
518
944
 
519
945
 
520
946
 
549
975
      polymakeFile.writeMatrixProperty("ORTH_LINEALITY_SPACE",rowsToIntegerMatrix(highestDimensionalCone().linealitySpace().getEquations()));
550
976
      polymakeFile.writeCardinalProperty("PURE",1);
551
977
    }
552
 
    
 
978
 
553
979
 
554
980
  PolyhedralFan f=*this;
555
981
  stringstream s;
683
1109
 
684
1110
      IntegerVectorList m;
685
1111
      m.push_back(multiplicities);
686
 
      polymakeFile.writeMatrixProperty("MULTIPLICITIES",rowsToIntegerMatrix(m).transposed());      
 
1112
      polymakeFile.writeMatrixProperty("MULTIPLICITIES",rowsToIntegerMatrix(m).transposed());
687
1113
    }
688
1114
  if(polymakeFileName)
689
1115
    polymakeFile.close();
694
1120
*/
695
1121
 
696
1122
 
697
 
void PolyhedralFan::printWithIndices(class Printer *p, bool printMultiplicities, SymmetryGroup *sym, bool group)const 
 
1123
 
 
1124
  /*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
 
1125
   */
 
1126
 static list<SymmetricComplex::Cone> computeFacets(SymmetricComplex::Cone const &theCone, IntegerMatrix const &rays, IntegerVectorList const &facetCandidates, SymmetricComplex const &theComplex/*, int linealityDim*/)
 
1127
{
 
1128
  set<SymmetricComplex::Cone> ret;
 
1129
 
 
1130
  for(IntegerVectorList::const_iterator i=facetCandidates.begin();i!=facetCandidates.end();i++)
 
1131
    {
 
1132
      set<int> indices;
 
1133
 
 
1134
      bool notAll=false;
 
1135
      for(vector<int>::const_iterator j=theCone.indices.begin();j!=theCone.indices.end();j++)
 
1136
        if(dotLong(rays[*j],*i)==0)
 
1137
          indices.insert(*j);
 
1138
        else
 
1139
          notAll=true;
 
1140
 
 
1141
      SymmetricComplex::Cone temp(indices,theCone.dimension-1,0,false,theComplex);
 
1142
      /*      temp.multiplicity=0;
 
1143
      temp.dimension=theCone.dimension-1;
 
1144
      temp.setIgnoreSymmetry(true);
 
1145
      */
 
1146
      if(notAll)ret.insert(temp);
 
1147
 
 
1148
    }
 
1149
  //  fprintf(Stderr,"HEJ!!!!\n");
 
1150
 
 
1151
  list<SymmetricComplex::Cone> ret2;
 
1152
  for(set<SymmetricComplex::Cone>::const_iterator i=ret.begin();i!=ret.end();i++)
 
1153
    {
 
1154
      bool isMaximal=true;
 
1155
 
 
1156
      /*      if(i->indices.size()+linealityDim<i->dimension)//#3
 
1157
        isMaximal=false;
 
1158
        else*/
 
1159
        for(set<SymmetricComplex::Cone>::const_iterator j=ret.begin();j!=ret.end();j++)
 
1160
          if(i!=j && i->isSubsetOf(*j))
 
1161
            {
 
1162
              isMaximal=false;
 
1163
              break;
 
1164
            }
 
1165
      if(isMaximal)
 
1166
        {
 
1167
          SymmetricComplex::Cone temp(i->indexSet(),i->dimension,i->multiplicity,true,theComplex);
 
1168
          temp.setKnownToBeNonMaximal(); // THIS IS WHERE WE SET THE CONES NON-MAXIMAL FLAG
 
1169
          //      temp.setIgnoreSymmetry(false);
 
1170
          ret2.push_back(temp);
 
1171
        }
 
1172
    }
 
1173
  return ret2;
 
1174
}
 
1175
 
 
1176
 
 
1177
void addFacesToSymmetricComplex(SymmetricComplex &c, PolyhedralCone const &cone, IntegerVectorList const &facetCandidates, IntegerVectorList const &generatorsOfLinealitySpace)
 
1178
{
 
1179
  IntegerMatrix const &rays=c.getVertices();
 
1180
  set<int> indices;
 
1181
 
 
1182
//  for(int j=0;j<rays.getHeight();j++)if(cone.contains(rays[j]))indices.insert(j);
 
1183
 
 
1184
  IntegerVectorList l=cone.extremeRays(&generatorsOfLinealitySpace);
 
1185
  for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)indices.insert(c.indexOfVertex(*i));
 
1186
 
 
1187
  addFacesToSymmetricComplex(c,indices,facetCandidates,cone.dimension(),cone.getMultiplicity());
 
1188
}
 
1189
 
 
1190
void addFacesToSymmetricComplex(SymmetricComplex &c, set<int> const &indices, IntegerVectorList const &facetCandidates, int dimension, int multiplicity)
 
1191
{
 
1192
  IntegerMatrix const &rays=c.getVertices();
 
1193
  list<SymmetricComplex::Cone> clist;
 
1194
  {
 
1195
 
 
1196
    SymmetricComplex::Cone temp(indices,dimension,multiplicity,true,c);
 
1197
    //    temp.dimension=cone.dimension();
 
1198
    //   temp.multiplicity=cone.getMultiplicity();
 
1199
    clist.push_back(temp);
 
1200
  }
 
1201
 
 
1202
  //  int linealityDim=cone.dimensionOfLinealitySpace();
 
1203
 
 
1204
  while(!clist.empty())
 
1205
    {
 
1206
      SymmetricComplex::Cone &theCone=clist.front();
 
1207
 
 
1208
      if(!c.contains(theCone))
 
1209
        {
 
1210
          log2
 
1211
          {
 
1212
            static int t;
 
1213
            if((t&1023)==0)
 
1214
              {
 
1215
                fprintf(Stderr,"clist size:%i\n",clist.size());
 
1216
              }
 
1217
            t++;
 
1218
          }
 
1219
 
 
1220
          c.insert(theCone);
 
1221
          //      log0 fprintf(Stderr,"ADDING\n");
 
1222
          list<SymmetricComplex::Cone> facets=computeFacets(theCone,rays,facetCandidates,c/*,linealityDim*/);
 
1223
          clist.splice(clist.end(),facets);
 
1224
        }
 
1225
      clist.pop_front();
 
1226
    }
 
1227
 
 
1228
}
 
1229
 
 
1230
 
 
1231
/**
 
1232
   Produce strings that express the vectors in terms of rays of the fan modulo the lineality space. Symmetry is ignored??
 
1233
 */
 
1234
vector<string> PolyhedralFan::renamingStrings(IntegerVectorList const &theVectors, IntegerVectorList const &originalRays, IntegerVectorList const &linealitySpace, SymmetryGroup *sym)const
 
1235
{
 
1236
  vector<string> ret;
 
1237
  for(IntegerVectorList::const_iterator i=theVectors.begin();i!=theVectors.end();i++)
 
1238
    {
 
1239
      for(PolyhedralConeList::const_iterator j=cones.begin();j!=cones.end();j++)
 
1240
        {
 
1241
          if(j->contains(*i))
 
1242
            {
 
1243
              vector<int> relevantIndices;
 
1244
              IntegerVectorList relevantRays=linealitySpace;
 
1245
              int K=0;
 
1246
              for(IntegerVectorList::const_iterator k=originalRays.begin();k!=originalRays.end();k++,K++)
 
1247
                if(j->contains(*k))
 
1248
                  {
 
1249
                    relevantIndices.push_back(K);
 
1250
                    relevantRays.push_back(*k);
 
1251
                  }
 
1252
 
 
1253
              FieldMatrix LFA(Q,relevantRays.size(),n);
 
1254
              int J=0;
 
1255
              for(IntegerVectorList::const_iterator j=relevantRays.begin();j!=relevantRays.end();j++,J++)
 
1256
                LFA[J]=integerVectorToFieldVector(*j,Q);
 
1257
              FieldVector LFB=concatenation(integerVectorToFieldVector(*i,Q),FieldVector(Q,relevantRays.size()));
 
1258
              LFA=LFA.transposed();
 
1259
              FieldVector LFX=LFA.solver().canonicalize(LFB);
 
1260
              stringstream s;
 
1261
              if(LFX.subvector(0,n).isZero())
 
1262
                {
 
1263
                  s<<"Was:";
 
1264
                  FieldVector S=LFX.subvector(n+linealitySpace.size(),LFX.size());
 
1265
                  for(int k=0;k<S.size();k++)
 
1266
                    if(!S[k].isZero())
 
1267
                      s<<"+"<<S[k].toString()<<"*["<<relevantIndices[k]<<"] ";
 
1268
                }
 
1269
              ret.push_back(s.str());
 
1270
              break;
 
1271
            }
 
1272
        }
 
1273
    }
 
1274
  return ret;
 
1275
}
 
1276
 
 
1277
SymmetricComplex PolyhedralFan::toSymmetricComplex(SymmetryGroup *sym)
 
1278
{
 
1279
  SymmetryGroup symm(n);
 
1280
          if(!sym)sym=&symm;
 
1281
 
 
1282
          IntegerVectorList rays=getRaysInPrintingOrder(sym);
 
1283
 
 
1284
          SymmetricComplex symCom(n,rays,*sym);
 
1285
 
 
1286
          if(cones.empty())return symCom;
 
1287
          IntegerVectorList generatorsOfLinealitySpace=cones.begin()->generatorsOfLinealitySpace();
 
1288
 
 
1289
          for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1290
            {
 
1291
              {
 
1292
                static int t;
 
1293
                log1 fprintf(Stderr,"Adding faces of cone %i\n",t++);
 
1294
              }
 
1295
              log2 fprintf(Stderr,"Dim: %i\n",i->dimension());
 
1296
 
 
1297
              addFacesToSymmetricComplex(symCom,*i,i->getHalfSpaces(),generatorsOfLinealitySpace);
 
1298
            }
 
1299
 
 
1300
          log1 cerr<<"Remapping";
 
1301
          symCom.remap();
 
1302
          log1 cerr<<"Done remapping";
 
1303
          return symCom;
 
1304
}
 
1305
 
 
1306
void PolyhedralFan::printWithIndices(class Printer *p, int flags, SymmetryGroup *sym, vector<string> const *comments)const
 
1307
//void PolyhedralFan::printWithIndices(class Printer *p, bool printMultiplicities, SymmetryGroup *sym, bool group, bool ignoreCones, bool xml, bool tPlaneSort, vector<string> const *comments)const
698
1308
{
699
1309
  assert(p);
700
1310
  //  IntegerVector multiplicities;
701
1311
  SymmetryGroup symm(n);
702
1312
 
703
1313
  PolymakeFile polymakeFile;
704
 
  polymakeFile.create("NONAME","PolyhedralFan","PolyhedralFan");
 
1314
//  polymakeFile.create("NONAME","fan","PolyhedralFan",flags&FPF_xml);
 
1315
  polymakeFile.create("NONAME","fan","SymmetricFan",flags&FPF_xml);
 
1316
 
 
1317
  bool produceXml=polymakeFile.isXmlFormat();
 
1318
 
705
1319
 
706
1320
  if(!sym)sym=&symm;
707
1321
 
712
1326
    }
713
1327
 
714
1328
  int h=cones.begin()->dimensionOfLinealitySpace();
 
1329
 
 
1330
  log1 fprintf(Stderr,"Computing rays.\n");
715
1331
  IntegerVectorList rays=getRaysInPrintingOrder(sym);
716
1332
 
717
1333
  SymmetricComplex symCom(n,rays,*sym);
719
1335
  polymakeFile.writeCardinalProperty("AMBIENT_DIM",n);
720
1336
  polymakeFile.writeCardinalProperty("DIM",getMaxDimension());
721
1337
  polymakeFile.writeCardinalProperty("LINEALITY_DIM",h);
722
 
  polymakeFile.writeMatrixProperty("RAYS",rowsToIntegerMatrix(rays,n),true);
 
1338
  polymakeFile.writeMatrixProperty("RAYS",rowsToIntegerMatrix(rays,n),true,comments);
723
1339
  polymakeFile.writeCardinalProperty("N_RAYS",rays.size());
724
 
  polymakeFile.writeMatrixProperty("LINEALITY_SPACE",rowsToIntegerMatrix(highestDimensionalCone().linealitySpace().dualCone().getEquations(),n));
 
1340
  IntegerVectorList linealitySpaceGenerators=highestDimensionalCone().linealitySpace().dualCone().getEquations();
 
1341
  polymakeFile.writeMatrixProperty("LINEALITY_SPACE",rowsToIntegerMatrix(linealitySpaceGenerators,n));
725
1342
  polymakeFile.writeMatrixProperty("ORTH_LINEALITY_SPACE",rowsToIntegerMatrix(highestDimensionalCone().linealitySpace().getEquations(),n));
726
 
    
 
1343
 
 
1344
  if(flags & FPF_primitiveRays)
 
1345
  {
 
1346
         IntegerVectorList primitiveRays;
 
1347
         for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
 
1348
                 for(PolyhedralConeList::const_iterator j=cones.begin();j!=cones.end();j++)
 
1349
                         if(j->contains(*i)&&(j->dimensionOfLinealitySpace()+1==j->dimension()))
 
1350
                                         primitiveRays.push_back(j->semiGroupGeneratorOfRay());
 
1351
 
 
1352
          polymakeFile.writeMatrixProperty("PRIMITIVE_RAYS",rowsToIntegerMatrix(primitiveRays,n));
 
1353
  }
 
1354
 
 
1355
 
 
1356
  IntegerVectorList generatorsOfLinealitySpace=cones.begin()->generatorsOfLinealitySpace();
 
1357
 
 
1358
  log1 fprintf(Stderr,"Building symmetric complex.\n");
 
1359
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1360
    {
 
1361
      {
 
1362
        static int t;
 
1363
        log1 fprintf(Stderr,"Adding faces of cone %i\n",t++);
 
1364
      }
 
1365
      log2 fprintf(Stderr,"Dim: %i\n",i->dimension());
 
1366
 
 
1367
      addFacesToSymmetricComplex(symCom,*i,i->getHalfSpaces(),generatorsOfLinealitySpace);
 
1368
    }
 
1369
 
 
1370
  log1 cerr<<"Remapping";
 
1371
  symCom.remap();
 
1372
  log1 cerr<<"Done remapping";
 
1373
 
727
1374
 
728
1375
  PolyhedralFan f=*this;
729
1376
 
730
1377
  //  IntegerVector fvector(f.getMaxDimension()-h);
731
1378
 
 
1379
 
 
1380
 
732
1381
  //fprintf(Stderr,"maxdim %i h %i\n",f.getMaxDimension(),h);
733
 
  while(!f.cones.empty())
 
1382
  /*  while(!f.cones.empty())
734
1383
    {
735
1384
      int currentDimension=f.getMaxDimension()-h;
736
1385
      IntegerVector rayIncidenceCounter(rays.size());
756
1405
      //      fvector[f.getMaxDimension()-h-1]=faceIndex;
757
1406
      f=f.facetComplexSymmetry(*sym);
758
1407
    }
759
 
 
760
 
  polymakeFile.writeCardinalVectorProperty("F_VECTOR",symCom.fvector());
761
 
  
762
 
  polymakeFile.writeStringProperty("CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),false,group));
763
 
  stringstream multiplicities;
764
 
  polymakeFile.writeStringProperty("MAXIMAL_CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,group, &multiplicities));
765
 
 
766
 
  polymakeFile.writeCardinalProperty("PURE",symCom.isPure());
767
 
 
768
 
  //  IntegerVectorList m;
769
 
  //  m.push_back(multiplicities);
770
 
  //  polymakeFile.writeMatrixProperty("MULTIPLICITIES",rowsToIntegerMatrix(m).transposed());      
771
 
  if(printMultiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES",multiplicities.str());
772
 
 
 
1408
  */
 
1409
  log1 fprintf(Stderr,"Computing f-vector.\n");
 
1410
  IntegerVector fvector=symCom.fvector();
 
1411
  polymakeFile.writeCardinalVectorProperty("F_VECTOR",fvector);
 
1412
  log1 fprintf(Stderr,"Done computing f-vector.\n");
 
1413
 
 
1414
  if(flags&FPF_boundedInfo)
 
1415
    {
 
1416
      log1 fprintf(Stderr,"Computing bounded f-vector.\n");
 
1417
      IntegerVector fvectorBounded=symCom.fvector(true);
 
1418
      polymakeFile.writeCardinalVectorProperty("F_VECTOR_BOUNDED",fvectorBounded);
 
1419
      log1 fprintf(Stderr,"Done computing bounded f-vector.\n");
 
1420
    }
 
1421
/*
 
1422
 * Removed to make the Polymake people happy.
 
1423
 *    {
 
1424
    int euler=0;
 
1425
    int mul=-1;
 
1426
    for(int i=0;i<fvector.size();i++,mul*=-1)euler+=mul*fvector[i];
 
1427
    polymakeFile.writeCardinalProperty("MY_EULER",euler);
 
1428
  }
 
1429
*/
 
1430
  log1 fprintf(Stderr,"Checking if complex is simplicial and pure.\n");
 
1431
  polymakeFile.writeBooleanProperty("SIMPLICIAL",symCom.isSimplicial());
 
1432
  polymakeFile.writeBooleanProperty("PURE",symCom.isPure());
 
1433
  log1 fprintf(Stderr,"Done checking.\n");
 
1434
 
 
1435
 
 
1436
  if(flags&FPF_conesCompressed)
 
1437
  {
 
1438
    polymakeFile.writeArrayArrayIntProperty("SYMMETRY_GENERATORS",rowsToIntegerMatrix(sym->getUniqueGenerators(),n));
 
1439
    log1 fprintf(Stderr,"Producing list of cones up to symmetry.\n");
 
1440
    polymakeFile.writeStringProperty("CONES_ORBITS",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),false,flags&FPF_group,0,true,flags&FPF_tPlaneSort,produceXml));
 
1441
    log1 fprintf(Stderr,"Done producing list of cones up to symmetry.\n");
 
1442
    log1 fprintf(Stderr,"Producing list of maximal cones up to symmetry.\n");
 
1443
    stringstream multiplicities;
 
1444
    polymakeFile.writeStringProperty("MAXIMAL_CONES_ORBITS",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,flags&FPF_group, &multiplicities,true,flags&FPF_tPlaneSort,produceXml));
 
1445
    if(flags&FPF_multiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES_ORBITS",multiplicities.str());
 
1446
    log1 fprintf(Stderr,"Done producing list of maximal cones up to symmetry.\n");
 
1447
  }
 
1448
 
 
1449
  if(flags&FPF_conesExpanded)
 
1450
    {
 
1451
      if(flags&FPF_cones)
 
1452
        {
 
1453
          log1 fprintf(Stderr,"Producing list of cones.\n");
 
1454
          polymakeFile.writeStringProperty("CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),false,flags&FPF_group,0,false,flags&FPF_tPlaneSort,produceXml));
 
1455
          log1 fprintf(Stderr,"Done producing list of cones.\n");
 
1456
        }
 
1457
      if(flags&FPF_maximalCones)
 
1458
        {
 
1459
          log1 fprintf(Stderr,"Producing list of maximal cones.\n");
 
1460
          stringstream multiplicities;
 
1461
          polymakeFile.writeStringProperty("MAXIMAL_CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,flags&FPF_group, &multiplicities,false,flags&FPF_tPlaneSort,produceXml));
 
1462
          if(flags&FPF_multiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES",multiplicities.str());
 
1463
          log1 fprintf(Stderr,"Done producing list of maximal cones.\n");
 
1464
        }
 
1465
    }
 
1466
 
 
1467
  if(flags&FPF_values)
 
1468
    {
 
1469
      {
 
1470
        IntegerVectorList values;
 
1471
        for(IntegerVectorList::const_iterator i=linealitySpaceGenerators.begin();i!=linealitySpaceGenerators.end();i++)
 
1472
          {
 
1473
            IntegerVector v(1);
 
1474
            v[0]=evaluatePiecewiseLinearFunction(*i);
 
1475
            values.push_back(v);
 
1476
          }
 
1477
        polymakeFile.writeMatrixProperty("LINEALITY_VALUES",rowsToIntegerMatrix(values,1));
 
1478
      }
 
1479
      {
 
1480
        IntegerVectorList values;
 
1481
        for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
 
1482
          {
 
1483
            IntegerVector v(1);
 
1484
            v[0]=evaluatePiecewiseLinearFunction(*i);
 
1485
            values.push_back(v);
 
1486
          }
 
1487
        polymakeFile.writeMatrixProperty("RAY_VALUES",rowsToIntegerMatrix(values,1));
 
1488
      }
 
1489
    }
 
1490
 
 
1491
 
 
1492
  log1 fprintf(Stderr,"Producing final string for output.\n");
773
1493
  stringstream s;
774
1494
  polymakeFile.writeStream(s);
775
1495
  string S=s.str();
 
1496
  log1 fprintf(Stderr,"Printing string.\n");
776
1497
  p->printString(S.c_str());
777
 
}
 
1498
  log1 fprintf(Stderr,"Done printing string.\n");
 
1499
}
 
1500
 
 
1501
 
 
1502
PolyhedralFan PolyhedralFan::readFan(string const &filename, bool onlyMaximal, IntegerVector *w, set<int> const *coneIndices, SymmetryGroup const *sym, bool readCompressedIfNotSym)
 
1503
{
 
1504
    PolymakeFile inFile;
 
1505
    inFile.open(filename.c_str());
 
1506
 
 
1507
    int n=inFile.readCardinalProperty("AMBIENT_DIM");
 
1508
    int nRays=inFile.readCardinalProperty("N_RAYS");
 
1509
    IntegerMatrix rays=inFile.readMatrixProperty("RAYS",nRays,n);
 
1510
    int linealityDim=inFile.readCardinalProperty("LINEALITY_DIM");
 
1511
    IntegerMatrix linealitySpace=inFile.readMatrixProperty("LINEALITY_SPACE",linealityDim,n);
 
1512
 
 
1513
 
 
1514
    const char *sectionName=0;
 
1515
    const char *sectionNameMultiplicities=0;
 
1516
    if(sym || readCompressedIfNotSym)
 
1517
    {
 
1518
      sectionName=(onlyMaximal)?"MAXIMAL_CONES_ORBITS":"CONES_ORBITS";
 
1519
      sectionNameMultiplicities=(onlyMaximal)?"MULTIPLICITIES_ORBITS":"DUMMY123";
 
1520
    }
 
1521
      else
 
1522
      {  sectionName=(onlyMaximal)?"MAXIMAL_CONES":"CONES";
 
1523
      sectionNameMultiplicities=(onlyMaximal)?"MULTIPLICITIES":"DUMMY123";
 
1524
      }
 
1525
 
 
1526
 
 
1527
    IntegerVector w2(n);
 
1528
    if(w==0)w=&w2;
 
1529
 
 
1530
    SymmetryGroup sym2(n);
 
1531
    if(sym==0)sym=&sym2;
 
1532
 
 
1533
    vector<list<int> > cones=inFile.readMatrixIncidenceProperty(sectionName);
 
1534
    IntegerVectorList r;
 
1535
 
 
1536
    bool hasMultiplicities=inFile.hasProperty(sectionNameMultiplicities);
 
1537
    IntegerMatrix multiplicities(0,0);
 
1538
    if(hasMultiplicities)multiplicities=inFile.readMatrixProperty(sectionNameMultiplicities,cones.size(),1);
 
1539
 
 
1540
 
 
1541
    PolyhedralFan ret(n);
 
1542
 
 
1543
    log2 cerr<< "Number of orbits to expand "<<cones.size()<<endl;
 
1544
    for(int i=0;i<cones.size();i++)
 
1545
      if(coneIndices==0 || coneIndices->count(i))
 
1546
        {
 
1547
          log2 cerr<<"Expanding symmetries of cone"<<endl;
 
1548
          /*      for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
 
1549
            {
 
1550
              IntegerVectorList coneRays;
 
1551
              for(list<int>::const_iterator j=cones[i].begin();j!=cones[i].end();j++)
 
1552
                coneRays.push_back(SymmetryGroup::compose(*perm,rays[*j]));
 
1553
              if(isInNonNegativeSpan(*w,coneRays,linealitySpace.getRows()))
 
1554
                {
 
1555
                  PolyhedralCone C=PolyhedralCone::givenByRays(coneRays,linealitySpace.getRows(),n);
 
1556
                  C.canonicalize();
 
1557
                  ret.insert(C);
 
1558
                }
 
1559
            }
 
1560
          */
 
1561
          {
 
1562
            IntegerVectorList coneRays;
 
1563
            for(list<int>::const_iterator j=cones[i].begin();j!=cones[i].end();j++)
 
1564
              coneRays.push_back((rays[*j]));
 
1565
            PolyhedralCone C=PolyhedralCone::givenByRays(coneRays,linealitySpace.getRows(),n);
 
1566
            if(hasMultiplicities)C.setMultiplicity(multiplicities[i][0]);
 
1567
            for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
 
1568
              {
 
1569
                if(C.contains(SymmetryGroup::composeInverse(*perm,*w)))
 
1570
                  {
 
1571
                    PolyhedralCone C2=C.permuted(*perm);
 
1572
                    C2.canonicalize();
 
1573
                    ret.insert(C2);
 
1574
                  }
 
1575
              }
 
1576
          }
 
1577
        }
 
1578
    return ret;
 
1579
}
 
1580
 
778
1581
 
779
1582
IncidenceList PolyhedralFan::getIncidenceList(SymmetryGroup *sym)const //fan must be pure
780
1583
{
838
1641
}
839
1642
 
840
1643
 
 
1644
PolyhedralCone PolyhedralFan::coneContaining(IntegerVector const &v)const
 
1645
{
 
1646
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1647
    if(i->contains(v))return i->faceContaining(v);
 
1648
  debug<<"Vector "<<v<<" not contained in support of fan\n";
 
1649
  assert(0);
 
1650
}
 
1651
 
 
1652
 
841
1653
PolyhedralFan::coneIterator PolyhedralFan::conesBegin()const
842
1654
{
843
1655
  return cones.begin();
861
1673
          PolyhedralCone c=intersection(*i,*j);
862
1674
          if(c.dimension()==n)
863
1675
            {
864
 
             
 
1676
 
865
1677
              if(!j->contains(*i))
866
1678
                return false;
867
1679
            }
886
1698
      t+=t1;
887
1699
    }
888
1700
  fprintf(Stderr,"%i\n",t);
889
 
  
 
1701
 
890
1702
      /*      bool found=false;
891
1703
      for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
892
1704
          if(j->contains(*i)){
902
1714
                AsciiPrinter(Stderr).printVector(c.getRelativeInteriorPoint());
903
1715
              }
904
1716
            }
905
 
          
 
1717
 
906
1718
          return false;
907
1719
          }*/
908
1720
 
909
1721
  return true;
910
1722
}
 
1723
 
 
1724
/*          for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
 
1725
              {
 
1726
                if(C.contains(SymmetryGroup::composeInverse(*perm,*w)))
 
1727
                  {
 
1728
                    PolyhedralCone C2=C.permuted(*perm);
 
1729
                    C2.canonicalize();
 
1730
                    ret.insert(C2);
 
1731
                  }
 
1732
              }
 
1733
*/
 
1734
PolyhedralFan PolyhedralFan::link(IntegerVector const &w, SymmetryGroup *sym)const
 
1735
{
 
1736
  SymmetryGroup symL(n);
 
1737
  if(!sym)sym=&symL;
 
1738
 
 
1739
  PolyhedralFan ret(n);
 
1740
 
 
1741
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1742
    {
 
1743
      for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
 
1744
        {
 
1745
          IntegerVector w2=SymmetryGroup::composeInverse(*perm,w);
 
1746
          if(i->contains(w2))
 
1747
            {
 
1748
              IntegerVectorList equations=i->getEquations();
 
1749
              IntegerVectorList inequalities1=i->getHalfSpaces();
 
1750
              IntegerVectorList inequalities2;
 
1751
              for(IntegerVectorList::const_iterator j=inequalities1.begin();j!=inequalities1.end();j++)
 
1752
                if(dotLong(w2,*j)==0)inequalities2.push_back(*j);
 
1753
              PolyhedralCone C(inequalities2,equations,n);
 
1754
              C.canonicalize();
 
1755
              C.setLinearForm(i->getLinearForm());
 
1756
              PolyhedralCone C2=C.permuted(*perm);
 
1757
              C2.canonicalize();
 
1758
              C2.setMultiplicity(i->getMultiplicity());
 
1759
              ret.insert(C2);
 
1760
            }
 
1761
        }
 
1762
    }
 
1763
  return ret;
 
1764
}
 
1765
 
 
1766
PolyhedralFan PolyhedralFan::link(IntegerVector const &w)const
 
1767
{
 
1768
  PolyhedralFan ret(n);
 
1769
 
 
1770
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1771
    {
 
1772
      if(i->contains(w))
 
1773
        {
 
1774
          IntegerVectorList equations=i->getEquations();
 
1775
          IntegerVectorList inequalities1=i->getHalfSpaces();
 
1776
          IntegerVectorList inequalities2;
 
1777
          for(IntegerVectorList::const_iterator j=inequalities1.begin();j!=inequalities1.end();j++)
 
1778
            if(dotLong(w,*j)==0)inequalities2.push_back(*j);
 
1779
          PolyhedralCone C(inequalities2,equations,n);
 
1780
          C.canonicalize();
 
1781
          C.setLinearForm(i->getLinearForm());
 
1782
      C.setMultiplicity(i->getMultiplicity());
 
1783
          ret.insert(C);
 
1784
        }
 
1785
    }
 
1786
  return ret;
 
1787
}
 
1788
 
 
1789
 
 
1790
int64 PolyhedralFan::evaluatePiecewiseLinearFunction(IntegerVector const &x)const
 
1791
{
 
1792
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1793
    {
 
1794
      if(i->contains(x))return dotLong(i->getLinearForm(),x);
 
1795
    }
 
1796
  assert(0);
 
1797
  return 0;
 
1798
}
 
1799
 
 
1800
 
 
1801
FieldElement PolyhedralFan::volume(int d, SymmetryGroup *sym)const
 
1802
{
 
1803
  SymmetryGroup symL(n);
 
1804
  if(!sym)sym=&symL;
 
1805
 
 
1806
  FieldElement ret(Q);
 
1807
 
 
1808
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1809
    {
 
1810
      if(i->dimension()==d)
 
1811
        {
 
1812
          IntegerVector w=stableRay(*i,sym);
 
1813
          ret=ret+Q.zHomomorphism(sym->orbitSize(w))*i->volume();
 
1814
        }
 
1815
    }
 
1816
  return ret;
 
1817
}
 
1818
 
 
1819
 
 
1820
bool PolyhedralFan::isConnected(SymmetryGroup *sym)const
 
1821
{
 
1822
  SymmetryGroup symL(n);
 
1823
  if(!sym)sym=&symL;
 
1824
 
 
1825
  CodimOneConnectednessTester ct;
 
1826
 
 
1827
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1828
    {
 
1829
      log2 cerr<<"Computing ridges of facet." << endl;
 
1830
      PolyhedralFan ridges=facetsOfCone(*i);
 
1831
      IntegerVectorList interiorPoints;
 
1832
      for(PolyhedralConeList::const_iterator j=ridges.cones.begin();j!=ridges.cones.end();j++)
 
1833
        interiorPoints.push_back(sym->orbitRepresentative(j->getUniquePoint()));
 
1834
      ct.insertFacetOrbit(interiorPoints);
 
1835
    }
 
1836
  return ct.isConnected();
 
1837
}
 
1838
 
 
1839
 
 
1840
int PolyhedralFan::size()const
 
1841
{
 
1842
  return cones.size();
 
1843
}
 
1844
 
 
1845
int PolyhedralFan::dimensionOfLinealitySpace()const
 
1846
{
 
1847
  assert(cones.size());//slow!
 
1848
  return cones.begin()->dimensionOfLinealitySpace();
 
1849
}
 
1850
 
 
1851
PolyhedralFan PolyhedralFan::negated()const
 
1852
{
 
1853
  PolyhedralFan ret(n);
 
1854
 
 
1855
  for(coneIterator i=conesBegin();i!=conesEnd();i++)
 
1856
    ret.insert(i->negated());
 
1857
  return ret;
 
1858
}
 
1859
 
 
1860
 
 
1861
bool PolyhedralFan::isCompatible(PolyhedralCone const &c)const
 
1862
{
 
1863
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1864
    {
 
1865
      PolyhedralCone C=intersection(c,*i);
 
1866
      C.canonicalize();
 
1867
      if(!c.hasFace(C))return false;
 
1868
      if(!i->hasFace(C))return false;
 
1869
    }
 
1870
  return true;
 
1871
}
 
1872
 
 
1873
void PolyhedralFan::merge(PolyhedralCone const &c)
 
1874
{
 
1875
  AsciiPrinter P(Stderr);
 
1876
 
 
1877
  if(isCompatible(c))insert(c);
 
1878
  else
 
1879
    {
 
1880
      assert(0);//Does not work in general.
 
1881
    }
 
1882
  /*
 
1883
  //  P<<"BEFORE MERGE-------------------------" <<*this;
 
1884
 
 
1885
  PolyhedralFan ret=complementOfCone(c,false);
 
1886
  //  P<<"COMPLEMENT OF CONE-------------------------" <<ret;
 
1887
  ret=refinement(*this,ret);
 
1888
 
 
1889
  PolyhedralFan C(c.ambientDimension());
 
1890
  C.insert(c);
 
1891
 
 
1892
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
 
1893
    C=refinement(complementOfCone(*i,true),C);
 
1894
 
 
1895
  for(PolyhedralConeList::const_iterator i=C.cones.begin();i!=C.cones.end();i++)
 
1896
    ret.insert(*i);
 
1897
  *this=ret;
 
1898
 
 
1899
  // P<<"AFTER MERGE" <<*this;
 
1900
  */
 
1901
}
 
1902
 
 
1903
 
 
1904
 
 
1905
void PolyhedralFan::removeNonMaximal()
 
1906
{
 
1907
  for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();)
 
1908
    {
 
1909
      IntegerVector w=i->getRelativeInteriorPoint();
 
1910
      bool containedInOther=false;
 
1911
      for(PolyhedralConeList::iterator j=cones.begin();j!=cones.end();j++)
 
1912
        if(j!=i)
 
1913
          {
 
1914
            if(j->contains(w)){containedInOther=true;break;}
 
1915
          }
 
1916
      if(containedInOther)
 
1917
        {
 
1918
          PolyhedralConeList::iterator k=i;
 
1919
          i++;
 
1920
          cones.erase(k);
 
1921
        }
 
1922
      else i++;
 
1923
    }
 
1924
}
 
1925
 
 
1926