1 #include "traverser_sphere.h"
2 #include "determinant.h"
3
4 #include "printer.h"
5
6 //void updatePolyhedralCone();
SphereTraverser(vector<PolyhedralCone> const & cones_,map<IntegerVector,list<int>> const & adjacency_,IntegerVector const & startCone,IntegerVector const & normal)7 SphereTraverser::SphereTraverser(vector<PolyhedralCone> const &cones_, map<IntegerVector,list<int> > const &adjacency_, IntegerVector const &startCone, IntegerVector const &normal):
8 ConeTraverser(normal.size()),
9 adjacency(adjacency_),
10 cones(cones_)
11 {
12 currentConeIndex=0;
13 for(;currentConeIndex<cones.size();currentConeIndex++)
14 {
15 if(cones[currentConeIndex].contains(startCone))
16 break;
17 }
18 assert(currentConeIndex!=cones.size());
19
20 currentNormal=normal;
21 }
22
changeCone(IntegerVector const & ridgeVector,IntegerVector const & rayVector)23 void SphereTraverser::changeCone(IntegerVector const &ridgeVector, IntegerVector const &rayVector)
24 {
25 list<int> const &candidates=adjacency[ridgeVector];
26 for(list<int>::const_iterator i=candidates.begin();i!=candidates.end();i++)
27 {
28 PolyhedralCone l=cones[*i].link(ridgeVector);
29 if(l.contains(rayVector))
30 {
31 IntegerVectorList M;
32 M.push_back(cones[currentConeIndex].getRelativeInteriorPoint());
33 M.push_back(currentNormal);
34 currentConeIndex=*i;
35 IntegerVectorList eq=cones[*i].getEquations();
36 assert(eq.size()==1);
37
38 MatrixTermOrder T(M);
39 if(!T(*eq.begin(),currentNormal-currentNormal))//HERE
40 // if(dotLong(currentNormal,*eq.begin())>0)
41 {currentNormal=*eq.begin();}
42 else
43 {currentNormal=-*eq.begin();}
44
45 break;
46 }
47 // assert(0);
48 }
49 }
50
link(IntegerVector const & ridgeVector)51 IntegerVectorList SphereTraverser::link(IntegerVector const &ridgeVector)
52 {
53 list<int> const &candidates=adjacency[ridgeVector];
54
55 PolyhedralCone ridge=cones[currentConeIndex].faceContaining(ridgeVector);
56 ridge.canonicalize();
57 // IntegerVectorList partialBasis=ridge.getEquations();
58 IntegerVectorList partialBasis=ridge.link(ridge.getRelativeInteriorPoint()).dualCone().getEquations();
59
60 IntegerVector myRay=cones[currentConeIndex].link(ridgeVector).getUniquePoint();
61
62 partialBasis.push_back(currentNormal);
63 partialBasis.push_back(myRay);
64 int sign=determinantSign(partialBasis);
65 partialBasis.pop_back();
66 partialBasis.pop_back();
67
68 IntegerVectorList c2;
69
70 for(list<int>::const_iterator i=candidates.begin();i!=candidates.end();i++)
71 {
72 // debug << "TESTET\n";
73 PolyhedralCone temp=cones[*i].link(ridgeVector);
74 IntegerVector v=temp.getUniquePoint();
75 if(dotLong(v,currentNormal)>=0)c2.push_back(v);
76 }
77 assert(c2.size()>0);
78 IntegerVectorList ret2;
79 ret2.push_back(myRay);
80 IntegerVector ret;
81 for(IntegerVectorList::const_iterator i=c2.begin();i!=c2.end();i++)
82 if(*i!=myRay){
83 ret=*i;
84 break;
85 }
86
87 if(ret.size()==0)
88 {
89 PolyhedralFan f(theCone.ambientDimension());
90 f.insert(theCone);
91 f.printWithIndices(&pout);
92 pout<<currentNormal;
93 }
94 for(IntegerVectorList::const_iterator i=c2.begin();i!=c2.end();i++)
95 {
96 if(*i==myRay) continue;
97 partialBasis.push_back(ret);
98 partialBasis.push_back(*i);
99 int s=determinantSign(partialBasis);
100 partialBasis.pop_back();
101 partialBasis.pop_back();
102 if(sign==s)
103 {
104 ret=*i;
105 }
106 }
107 ret2.push_back(ret);
108
109 // debug<<ret2;
110
111 return ret2;
112 }
113
114
refToPolyhedralCone()115 PolyhedralCone & SphereTraverser::refToPolyhedralCone()
116 {
117 theCone=cones[currentConeIndex];
118 return theCone;//cones[currentCone];
119 }
120