1 #include "parser.h"
2 #include "printer.h"
3 #include "polynomial.h"
4 #include "division.h"
5 #include "buchberger.h"
6 #include "wallideal.h"
7 #include "lp.h"
8 #include "reversesearch.h"
9 #include "breadthfirstsearch.h"
10 #include "termorder.h"
11 #include "ep_standard.h"
12 #include "ep_xfig.h"
13 #include "gfanapplication.h"
14 #include "timer.h"
15 #include "log.h"
16 #include "matrix.h"
17 #include "lll.h"
18 #include "polyhedralfan.h"
19 #include "linalg.h"
20 #include "determinant.h"
21 #include "triangulation.h"
22 #include "intsinpolytope.h"
23 #include "graph.h"
24 
25 #include "triangulation2.h"
26 
27 #include "traverser_secondaryfan.h"
28 #include "symmetrictraversal.h"
29 
30 #include <iostream>
31 #include <algorithm>
32 
33 #define NOOUTPUT 0
34 
35 class SecondaryFanApplication : public GFanApplication
36 {
37   SimpleOption hirschOption;
38   SimpleOption searchOption;
39   IntegerOption scaleOption;
40   StringOption optionRestrictingFan;
41   SimpleOption symmetryOption;
42   SimpleOption optionIgnoreCones;
43   InterruptOption optionInterrupt;
44 public:
helpText()45   const char *helpText()
46   {
47     return "This program computes the secondary fan of a vector configuration. The configuration is given as an ordered list of vectors. In order to compute the secondary fan of a point configuration an additional coordinate of ones must be added. For example {(1,0),(1,1),(1,2),(1,3)}.\n";
48   }
SecondaryFanApplication()49   SecondaryFanApplication():
50     searchOption("--unimodular","Use heuristics to search for unimodular triangulation rather than computing the complete secondary fan"),
51     scaleOption("--scale","Assuming that the first coordinate of each vector is 1, this option will take the polytope in the 1 plane and scale it. The point configuration will be all lattice points in that scaled polytope. The polytope must have maximal dimension. When this option is used the vector configuration must have full rank. This option may be removed in the future."),
52     symmetryOption("--symmetry","Tells the program to read in generators for a group of symmetries (subgroup of $S_n$) after having read in the vector configuration. The program checks that the configuration stays fixed when permuting the variables with respect to elements in the group. The output is grouped according to the symmetry.\n"),
53     optionRestrictingFan("--restrictingfan","Specify the name of a file containing a polyhedral fan in Polymake format. The computation of the Secondary fan will be restricted to this fan. If the --symmetry option is used then this restricting fan must be invariant under the symmetry and the orbits in the file must be with respect to the specified group of symmetries. The orbits of maximal cones of the file are then read in rather than the maximal cones.\n",0),
54     optionIgnoreCones("--nocones","Tells the program not to output the CONES and MAXIMAL_CONES sections, but still output CONES_COMPRESSED and MAXIMAL_CONES_COMPRESSED if --symmetry is used."),
55     hirschOption("--hirsch","")
56   {
57     hirschOption.hide();
58     registerOptions();
59   }
60 
name()61   const char *name()
62   {
63     return "_secondaryfan";
64   }
65 
enumerate(Triangulation2 const & t)66   PolyhedralFan enumerate(Triangulation2 const &t)
67   {
68     PolyhedralFan ret(t.getN());
69     list<Triangulation2> active;
70     active.push_back(t);
71     IntegerVectorList interiorPoints;
72     interiorPoints.push_back(t.interiorPoint());
73     while(!active.empty())
74       {
75 	Triangulation2 a=active.front();
76 	PolyhedralCone C=a.secondaryCone();
77 
78 
79 	//	if(active.size()>100)break;//SLETMIGGGGG
80 	//log0 fprintf(stderr,"a\n");
81 	/*	{
82 	  PolyhedralCone C2=C;
83 	  C2.canonicalize();
84 	  }*/
85 	C.canonicalize();
86 	//log0 fprintf(stderr,"b\n");
87 	ret.insert(C);
88 	AsciiPrinter P(Stderr);
89 	//	C.print(&P);
90 	active.pop_front();
91 	//	fprintf(stderr,"pop\n");
92 	IntegerVectorList flips=a.facets();
93 	for(IntegerVectorList::const_iterator i=flips.begin();i!=flips.end();i++)
94 	  {
95 	    {
96 	      IntegerVectorList t=C.getEquations();
97 	      t.push_back(*i);
98 	      PolyhedralCone CF(C.getHalfSpaces(),t);
99 	      CF.findFacets();
100 	      //	      CF.canonicalize();
101 	    }
102 
103 	    if(!i->isNonNegative()) //is this the right condition or should i be negated?
104 	    //   if(!(-*i).isNonNegative()) //is this the right condition or should i be negated?
105 	      {
106 		Triangulation2 b=a;
107 		log1 AsciiPrinter(Stderr)<<*i;
108 		/*fprintf(stderr,"Before:");
109 		  b.print(P);*/
110 		//		b.flip(*i);
111 		b.flipNew(-*i);
112 		/*fprintf(stderr,"After:");
113 		  b.print(P);*/
114 		if(!b.isEmpty())
115 		  {
116 		    IntegerVectorList inequalities=b.inequalities();
117 		    bool isKnown=false;
118 		    for(IntegerVectorList::const_iterator j=interiorPoints.begin();j!=interiorPoints.end();j++)
119 		      {
120 			bool match=true;
121 			for(IntegerVectorList::const_iterator k=inequalities.begin();k!=inequalities.end();k++)
122 			  {
123 			    if(dotLong(-*k,*j)<0)
124 			      {
125 				match=false;
126 				break;
127 			      }
128 			  }
129 			if(match)isKnown=true;
130 		      }
131 		    if(!isKnown)
132 		      {
133 			active.push_back(b);
134 			interiorPoints.push_back(b.interiorPoint());
135 		      }
136 		  }
137 	      }
138 	  }
139       }
140     return ret;
141   }
interactive(Triangulation2 const & t)142   PolyhedralFan interactive(Triangulation2 const &t)
143   {
144     Triangulation2 a=t;
145 
146     while(1)
147       {
148 	fprintf(stdout,"Triangles in current triangulation:%i\n",(int)a.bases.size());
149 	PolyhedralCone C=a.secondaryCone();
150 
151 	C.canonicalize();
152 
153 
154 	AsciiPrinter Pstd(Stderr);
155 	IntegerVectorList flips=a.facets();
156 	int I=0;
157 	for(IntegerVectorList::const_iterator i=flips.begin();i!=flips.end();i++,I++)
158 	  {
159 	    fprintf(stdout,"%i:\n",I);
160 	    Pstd.printVector(*i);
161 
162 
163 	    if(!i->isNonNegative())
164 	      {
165 		Triangulation2 b=a;
166 		//fprintf(stderr,"Before:");
167 		//b.print(P);
168 		//		b.flip(*i);
169 		/*		b.flipNew(*i);
170 		//fprintf(stderr,"After:");
171 		//b.print(P);
172 		fprintf(stdout,"Triangles in new triangulation:%i\n",b.bases.size());
173 		*/
174 	      }
175 
176 	    fprintf(stdout,"\n");
177 	  }
178 	int s;
179 	//cin >> s;
180 	int err=fscanf(stdin,"%i",&s);
181 	assert(err==1);
182 	if((s>=0)&&(s<I))
183 	  {
184 	    IntegerVectorList::const_iterator i=flips.begin();
185 	    for(int I=0;I<s;I++)i++;
186 	    a.flipNew(*i);
187 	  }
188 	else
189 	  {
190 	    a.print(Pstd);
191 	  }
192       }
193     PolyhedralFan ret(0);
194 
195     return ret;
196   }
automatic(Triangulation2 const & t,int abortAtSize)197   PolyhedralFan automatic(Triangulation2 const &t, int abortAtSize)
198   {
199     Triangulation2 a=t;
200 
201     cout << "Looking for triangulation with "<<abortAtSize<<" simplices."<<endl;
202 
203     while(1)
204       {
205 	fprintf(stdout,"Triangles in current triangulation:%i\n",(int)a.bases.size());
206 	//	PolyhedralCone C=a.secondaryCone();
207 
208 	//	C.canonicalize();
209 
210 
211 	AsciiPrinter Pstd(Stderr);
212 	IntegerVectorList flips=a.facets();
213 	int I=0;
214 	cerr << "Number of facets of secondary cone: " << flips.size() <<endl;
215 	for(IntegerVectorList::const_iterator i=flips.begin();i!=flips.end();i++,I++)
216 	  {
217 	    if(!i->isNonNegative())
218 	      {
219 		Triangulation2 b=a;
220 		b.flipNew(-*i);
221 		fprintf(stdout,"Triangles in new triangulation:%i\n",(int)b.bases.size());
222 
223 		if(b.bases.size()==abortAtSize)
224 		  {
225 		    b.print(Pstd);
226 
227 		    exit(0);
228 		  }
229 
230 		if((b.bases.size()>a.bases.size())||((rand()&127)==0))
231 		  {
232 		    a=b;
233 		    break;
234 		  }
235 	      }
236 	  }
237       }
238     PolyhedralFan ret(0);
239 
240     return ret;
241   }
automaticHirsch(Triangulation2 const & t)242   PolyhedralFan automaticHirsch(Triangulation2 const &t)
243   {
244     Triangulation2 a=t;
245     while(1)
246       {
247 	int nVertices=a.bases.size();
248 	int nEdges=a.coDimensionOneTriangles().size();
249 	int diameter=a.edgeGraph().diameter();
250 	int dimension=a.getD();
251 	int nFacets=a.usedRays().size();
252 	fprintf(stdout,"NVER: %i NEDGES: %i DIAMETER:%i DIMENSION:%i NFACETS:%i\n",nVertices,nEdges,diameter,dimension,nFacets);
253 
254 	AsciiPrinter Pstd(Stderr);
255 	IntegerVectorList flips=a.facets();
256 	int I=0;
257 	float currentScore=a.hirschScore();
258 	cerr << "Current score: " << currentScore <<endl;
259 	for(IntegerVectorList::const_iterator i=flips.begin();i!=flips.end();i++,I++)
260 	  {
261 	    if(!i->isNonNegative())
262 	      {
263 		Triangulation2 b=a;
264 		b.flipNew(-*i);
265 		float bScore=b.hirschScore();
266 		fprintf(stdout,"New score:%f\n",bScore);
267 
268 		if((bScore>currentScore)||((rand()&31)==0))
269 		  {
270 		    a=b;
271 		    break;
272 		  }
273 	      }
274 	  }
275       }
276     PolyhedralFan ret(0);
277 
278     return ret;
279   }
main()280   int main()
281   {
282 
283 
284 //    debug<<(int)sizeof(SymmetricComplex::Cone);
285 //    debug<<(int)sizeof(IntegerVector);
286 //    debug<<(int)sizeof(vector<int>);
287     IntegerMatrix A=rowsToIntegerMatrix(FileParser(Stdin).parseIntegerVectorList()).transposed();
288 
289     int n=A.getWidth();
290 
291     SymmetryGroup s(n);
292     if(symmetryOption.getValue())
293       {
294 	IntegerVectorList generators=FileParser(Stdin).parseIntegerVectorList();
295 	for(IntegerVectorList::const_iterator i=generators.begin();i!=generators.end();i++)
296 	  {
297 	    assert(i->size()==n);
298 	    FieldMatrix M1=integerMatrixToFieldMatrix(A,Q);
299 	    FieldMatrix M2=integerMatrixToFieldMatrix(rowsToIntegerMatrix(SymmetryGroup::permuteIntegerVectorList(A.getRows(),*i)),Q);
300 	    M1.reduce();
301 	    M1.REformToRREform(true);
302 	    M2.reduce();
303 	    M2.REformToRREform(true);
304 
305 	    if(!(M1==M2))
306 	      {
307 		AsciiPrinter(Stderr) << "Permutation "<< *i <<
308 		  " does not keep the configuration fixed.\n";
309 		assert(0);
310 	      }
311 	  }
312 	s.computeClosure(generators);
313         s.createTrie();
314       }
315 
316 
317     if(scaleOption.getValue())
318       {
319 	if(::rank(A)!=A.getHeight())
320 	  {
321 	    cerr << "The vector configuration must have full rank in order to use the scale option.\n";
322 	    assert(0);
323 	  }
324 
325 	int s=scaleOption.getValue();
326 
327 	cout << "Input configuration:" << endl;
328 	AsciiPrinter(Stdout)<<A.transposed().getRows();
329 
330 	IntegerVectorList empty;
331 	PolyhedralCone dual(A.transposed().getRows(),empty,n);
332 	dual.canonicalize();
333 	assert(dual.dimensionOfLinealitySpace()==0);
334 	IntegerVectorList inequalities=dual.extremeRays();
335 
336 	IntegerMatrix M1=rowsToIntegerMatrix(inequalities);
337 	IntegerMatrix M2=-1*M1.submatrix(0,1,M1.getHeight(),M1.getWidth());
338 	IntegerVector rightHandSide=s*M1.transposed().submatrix(0,0,1,M1.getWidth())[0].toVector();
339 	IntegerVector v=s*A.submatrix(1,0,A.getHeight(),1).transposed()[0].toVector();
340 
341 	IntegerVectorList l=intsInPolytopeGivenIneqAndPt(M2,rightHandSide,v);
342 
343 	IntegerVectorList lT=rowsToIntegerMatrix(l).transposed().getRows();
344 	lT.push_front(IntegerVector::allOnes(l.size()));
345 	l=rowsToIntegerMatrix(lT).transposed().getRows();
346 
347 	cout << "New configuration:" << endl;
348 	AsciiPrinter(Stdout)<<l;
349 	A=rowsToIntegerMatrix(l).transposed();
350       }
351 
352     /* If the vector configuration A does not have full rank then
353        change coordinates. */
354     Triangulation2::makeConfigurationFullDimensional(A);
355 
356     Triangulation2 t(A);
357     t.triangulate();
358 
359     if(searchOption.getValue())
360       {
361 	PolyhedralFan f=automatic(t,t.totalVolume());
362       }
363     else if(hirschOption.getValue())
364       {
365 	PolyhedralFan f=automaticHirsch(t);
366       }
367     else
368 {
369 #if !NOOUTPUT
370     	SymmetricTargetFanBuilder target(n,s);
371 #else
372     	SymmetricTargetNull target(n,s);
373 #endif
374 
375 		if(!optionRestrictingFan.getValue())
376 		{
377 			SecondaryFanTraverser traverser(t);
378 			SymmetricTargetCounterInterrupted target2(target,optionInterrupt.getValue());
379 			symmetricTraverse(traverser,target2,&s);
380 		}
381 		else
382 		{
383 		    PolyhedralFan f1=PolyhedralFan::readFan(optionRestrictingFan.getValue(),true,0,0,/*optionSymmetry.getValue()?&s:0*/0,false/*true*/);
384 
385 		    for(PolyhedralFan::coneIterator i=f1.conesBegin();i!=f1.conesEnd();i++)
386 			{
387 		    	static int a;
388 		    	log2 cerr<<"Processing Cone "<<a++<<" which has dimension "<<i->dimension()<<endl;
389 		    	SecondaryFanTraverser traverser(triangulationWithFullDimensionalIntersection(t,*i),*i);
390 				SymmetricTargetCounterInterrupted target2(target,optionInterrupt.getValue());
391 				symmetricTraverse(traverser,target2,&s);
392 			}
393 		}
394 
395 #if !NOOUTPUT
396 		target.getFanRef().printWithIndices(&pout,
397 	                                    (symmetryOption.getValue()?FPF_group|FPF_conesCompressed:0)|
398 	                                    (optionIgnoreCones.getValue()?0:FPF_conesExpanded)|
399 	                                    FPF_maximalCones|FPF_cones,
400 	                                    &s);
401 #endif
402 /*		target.getFanRef().printWithIndices(&p,
403 											FPF_default|
404 											(symmetryOption.getValue()?FPF_group|FPF_conesCompressed:0),
405 											&s);
406 */
407 /*    	PolyhedralFan f=enumerate(t).negated();//Changing sign
408 
409         f.printWithIndices(&p,
410 			   FPF_default|
411 			   (symmetryOption.getValue()?FPF_group|FPF_conesCompressed:0),
412 			   &s);
413   */
414 		}
415     return 0;
416   }
417 };
418 
419 static SecondaryFanApplication theApplication;
420 
421