1 #include "minkowskidual.h"
2 #include "field_rationals.h"
3 #include "linalg.h"
4
5
6
7 class MyCone{
8 public:
9 struct Pair{
10 int rayIndex;
11 vector<IntegerVector> const &rays;
12 IntegerVector const &permutation;
PairMyCone::Pair13 Pair(vector<IntegerVector> const &rays_, int rayIndex_, IntegerVector const & permutation_):
14 rays(rays_),
15 rayIndex(rayIndex_),
16 permutation(permutation_)
17 {
18 }
operator <MyCone::Pair19 bool operator<(Pair const &b)const
20 {
21 return SymmetryGroup::compose(permutation,rays[rayIndex])<SymmetryGroup::compose(b.permutation,b.rays[b.rayIndex]);
22 }
23 };
24 set<Pair> rays;
print()25 void print()
26 {
27 for(set<Pair>::const_iterator i=rays.begin();i!=rays.end();i++)
28 {
29 fprintf(Stderr,"--:%i ",i->rayIndex);
30 AsciiPrinter(Stderr).printVector((i->permutation));
31 fprintf(Stderr,"\n");
32 }
33 }
34 };
35
computeFacets(SymmetricComplex::Cone const & theCone,IntegerMatrix const & rays,IntegerVectorList const & facetCandidates,SymmetricComplex const & theComplex)36 static list<SymmetricComplex::Cone> computeFacets(SymmetricComplex::Cone const &theCone, IntegerMatrix const &rays, IntegerVectorList const &facetCandidates, SymmetricComplex const &theComplex/*, int linealityDim*/)
37 {
38 set<SymmetricComplex::Cone> ret;
39
40 for(IntegerVectorList::const_iterator i=facetCandidates.begin();i!=facetCandidates.end();i++)
41 {
42 set<int> indices;
43
44 bool notAll=false;
45 for(vector<int>::const_iterator j=theCone.indices.begin();j!=theCone.indices.end();j++)
46 if(dotLong(rays[*j],*i)==0)
47 indices.insert(*j);
48 else
49 notAll=true;
50
51 SymmetricComplex::Cone temp(indices,theCone.dimension-1,0,false,theComplex);
52 /* temp.multiplicity=0;
53 temp.dimension=theCone.dimension-1;
54 temp.setIgnoreSymmetry(true);
55 */
56 if(notAll)ret.insert(temp);
57
58 }
59 // fprintf(Stderr,"HEJ!!!!\n");
60
61 list<SymmetricComplex::Cone> ret2;
62 for(set<SymmetricComplex::Cone>::const_iterator i=ret.begin();i!=ret.end();i++)
63 {
64 bool isMaximal=true;
65
66 /* if(i->indices.size()+linealityDim<i->dimension)//#3
67 isMaximal=false;
68 else*/
69 for(set<SymmetricComplex::Cone>::const_iterator j=ret.begin();j!=ret.end();j++)
70 if(i!=j && i->isSubsetOf(*j))
71 {
72 isMaximal=false;
73 break;
74 }
75 if(isMaximal)
76 {
77 SymmetricComplex::Cone temp(i->indexSet(),i->dimension,i->multiplicity,true,theComplex);
78 // temp.setIgnoreSymmetry(false);
79 ret2.push_back(temp);
80 }
81 }
82 return ret2;
83 }
84
85
dualMinkowskiMixed(PolynomialSet const & g,SymmetryGroup const & sym,PolyhedralFan const & cones)86 SymmetricComplex dualMinkowskiMixed(PolynomialSet const &g, SymmetryGroup const &sym, PolyhedralFan const &cones)
87 {
88 /* The set of rays of the dual fan is the vertices of the Minkowski
89 sum after centering around zero. We only generate them up to symmetry. */
90 int n=cones.getAmbientDimension();
91 int nMaxOrbits=0;
92 for(PolyhedralFan::coneIterator i=cones.conesBegin();i!=cones.conesEnd();i++)nMaxOrbits++;
93 vector<IntegerVector> rays(nMaxOrbits);
94 FieldVector translationVector(Q,n);
95 for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
96 translationVector+=Q.zHomomorphism(i->numberOfTerms()).inverse()*integerVectorToFieldVector(i->exponentsSum(),Q);
97 {
98 int j=0;
99 for(PolyhedralFan::coneIterator i=cones.conesBegin();i!=cones.conesEnd();i++,j++)
100 {
101 FieldVector temp=integerVectorToFieldVector(i->getRelativeInteriorPoint(),Q);
102 temp+=(Q.zHomomorphism(-1)*translationVector);
103 rays[j]=temp.primitive();
104 }
105 }
106
107 fprintf(Stderr,"The new rays:\n");
108 for(int i=0;i<nMaxOrbits;i++)
109 {
110 fprintf(Stderr,"%i:",i);
111 AsciiPrinter(Stderr).printVector(rays[i]);
112 fprintf(Stderr,"\n");
113 }
114
115 /* The new cones of the dual fan comes from the rays of the origanal
116 Minkowski sum's normal fan. For every primary ray we add the set of dual rays whose primary cone contains the primary ray.
117
118 That is, add the vertices to the facets of the minkowskisum.
119 */
120 IntegerVectorList primaryRays=cones.getRaysInPrintingOrder(&sym,true);
121
122 vector<MyCone> newCones(primaryRays.size());
123 IntegerVectorList::const_iterator K=primaryRays.begin();
124 for(int k=0;k<newCones.size();k++,K++)
125 {
126 for(SymmetryGroup::ElementContainer::const_iterator l=sym.elements.begin();l!=sym.elements.end();l++)
127 {
128 IntegerVector v=SymmetryGroup::compose(*l,*K);
129 int j=0;
130 for(PolyhedralFan::coneIterator i=cones.conesBegin();i!=cones.conesEnd();i++,j++)
131 {
132 if(i->contains(v))
133 {
134 newCones[k].rays.insert(MyCone::Pair(rays,j,*l));
135 }
136 }
137 }
138 }
139
140 fprintf(Stderr,"Facets with their vertices:\n");
141 for(int k=0;k<newCones.size();k++)
142 {
143 fprintf(Stderr,"Printing facet:%i\n\n",k);
144 // newCones[k].print();
145 }
146
147
148
149 /* To make things easier we store for each primary cone the indices of rays it contains (among all rays)*/
150 IntegerVectorList allPrimaryRays=cones.getRaysInPrintingOrder(&sym,false);
151 vector<set<int> > allRaysIndices(nMaxOrbits);
152
153 {
154 int j=0;
155 for(PolyhedralFan::coneIterator i=cones.conesBegin();i!=cones.conesEnd();i++,j++)
156 {
157 int k=0;
158 for(IntegerVectorList::const_iterator l=allPrimaryRays.begin();l!=allPrimaryRays.end();l++,k++)
159 if(i->contains(*l))allRaysIndices[j].insert(k);
160 }
161 }
162
163 fprintf(Stderr,"Facets with their vertices (as indices):\n");
164 for(int k=0;k<nMaxOrbits;k++)
165 {
166
167 fprintf(Stderr,"Cone %i\n",k);
168 for(set<int>::const_iterator j=allRaysIndices[k].begin();j!=allRaysIndices[k].end();j++)
169 fprintf(Stderr,"%i ",*j);
170 fprintf(Stderr,"\n\n");
171 // newCones[k].print();
172 }
173
174
175 /* We prepare the complex to return. */
176
177 SymmetricComplex ret(n,allPrimaryRays,sym);
178
179 /* For every dual cone we need to compute facet normal. This is done
180 by cddlib. Afterward we use these facet normals to add the faces
181 of the dual cone to our collection of cones (pre-complex) - each
182 time checking that the corresponding face is mixed. */
183
184 fprintf(Stderr,"Adding mixed faces for cones:\n");
185 for(int k=0;k<newCones.size();k++)
186 {
187 IntegerVectorList a;
188 for(set<MyCone::Pair>::const_iterator i=newCones[k].rays.begin();i!=newCones[k].rays.end();i++)
189 {
190 a.push_back(SymmetryGroup::compose(i->permutation,rays[i->rayIndex]));
191 }
192
193 AsciiPrinter(Stderr).printVectorList(a);
194
195 IntegerVectorList empty;
196 PolyhedralCone C(a,empty,n);
197 IntegerVectorList facetCandidates=C.extremeRays();
198
199
200 fprintf(Stderr,"Cone %i\nFacetNormals:\n",k);
201 AsciiPrinter(Stderr).printVectorList(facetCandidates);
202 /* list<SymmetricComplex::Cone> clist;
203 {
204 set<int> indices;
205 for(int j=0;j<rays.getHeight();j++)if(cone.contains(rays[j]))indices.insert(j);
206 SymmetricComplex::Cone temp(indices,cone.dimension(),cone.getMultiplicity(),true,c);
207 clist.push_back(temp);
208 }
209
210 while(!clist.empty())
211 {
212 SymmetricComplex::Cone &theCone=clist.front();
213
214 if()//Check mixed face property
215 if(!ret.contains(theCone))
216 {
217 c.insert(theCone);
218
219 list<SymmetricComplex::Cone> facets=computeFacets(theCone,rays,facetCandidates,c);
220 clist.splice(clist.end(),facets);
221 }
222 clist.pop_front();
223 }
224 */
225 }
226
227
228
229 /* For every cone in the final "pre-complex" we translate it into a cone in the primary */
230
231 return ret;
232 }
233