1 #include "halfopencone.h"
2 
3 #include "buchberger.h"
4 #include "enumeration.h"
5 #include "reversesearch.h"
6 #include "wallideal.h"
7 
8 #include "printer.h"
9 #include "parser.h"
10 
11 #include "lp.h"
12 
13 
printHalfOpenCone(Printer & P,HalfOpenCone c)14 static void printHalfOpenCone(Printer &P, HalfOpenCone c)
15 {
16   P.printPolyhedralCone(c.closure());
17 }
18 
19 
printHalfOpenConeList(Printer & P,HalfOpenConeList const & l)20 static void printHalfOpenConeList(Printer &P, HalfOpenConeList const &l)
21 {
22   P.printString("Begin HalfOpenConeList\n");
23   for(HalfOpenConeList::const_iterator i=l.begin();i!=l.end();i++)
24     printHalfOpenCone(P,*i);
25   P.printString("End HalfOpenConeList\n");
26 }
27 
28 
contains(IntegerVector const & v) const29 bool HalfOpenCone::contains(IntegerVector const &v)const
30 {
31   IntegerVectorList inequalityList=lifted.getHalfSpaces();
32   IntegerVectorList equationList=lifted.getLinealitySpace();
33 
34   IntegerVectorList strict,nonstrict;
35 
36 
37   for(IntegerVectorList::const_iterator i=inequalityList.begin();i!=inequalityList.end();i++)
38     if((*i)[i->size()-1]<0)
39       strict.push_back(*i);
40     else if((*i)[i->size()-1]==0)
41       nonstrict.push_back(*i);
42     else
43       {//CHANGED
44 	assert(i->subvector(0,i->size()-1).isZero());
45 	strict.push_back(*i);
46       }
47 
48   for(IntegerVectorList::const_iterator i=equationList.begin();i!=equationList.end();i++)
49     if(dotLong(i->subvector(0,i->size()-1),v)!=0)return false;
50 
51   for(IntegerVectorList::const_iterator i=nonstrict.begin();i!=nonstrict.end();i++)
52     if(dotLong(i->subvector(0,i->size()-1),v)<0)return false;
53 
54   for(IntegerVectorList::const_iterator i=strict.begin();i!=strict.end();i++)
55     if(dotLong(i->subvector(0,i->size()-1),v)<=0)return false;
56 
57   return true;
58 }
59 
appendList(IntegerVectorList & to,IntegerVectorList const & from,int appendValue)60 void HalfOpenCone::appendList(IntegerVectorList &to, IntegerVectorList const &from, int appendValue)
61 {
62   for(IntegerVectorList::const_iterator i=from.begin();i!=from.end();i++)
63     {
64       IntegerVector v=*i;
65       v.resize(v.size()+1);
66       v[v.size()-1]=appendValue;
67       to.push_back(v);
68     }
69 }
70 
71 
HalfOpenCone(int dimension_,PolyhedralCone const & lifted_)72 HalfOpenCone::HalfOpenCone(int dimension_, PolyhedralCone const &lifted_):
73   dimension(dimension_),
74   liftedDimension(dimension_+1),
75   lifted(lifted_)
76 {
77   //  lifted.findFacets();
78 }
79 
80 
HalfOpenCone(int dimension_,IntegerVectorList const & equations,IntegerVectorList const & nonstrict,IntegerVectorList const & strict,bool findFacets)81 HalfOpenCone::HalfOpenCone(int dimension_, IntegerVectorList const &equations, IntegerVectorList const &nonstrict, IntegerVectorList const &strict, bool findFacets):
82   dimension(dimension_),
83   liftedDimension(dimension_+1),
84   lifted(dimension_+1)
85 {
86   IntegerVectorList equationList,inequalityList;
87 
88   appendList(equationList,equations,0);
89   appendList(inequalityList,nonstrict,0);
90   appendList(inequalityList,strict,-1);
91   inequalityList.push_back(IntegerVector::standardVector(liftedDimension,dimension)); //CHANGED
92   //      AsciiPrinter(Stderr).printVectorList(inequalityList);
93   //      AsciiPrinter(Stderr).printVectorList(equationList);
94   //      AsciiPrinter(Stderr).printInteger(liftedDimension);
95   lifted=PolyhedralCone(inequalityList,equationList,liftedDimension);
96   if(findFacets)lifted.findFacets();
97 }
98 
99 
swapFirstLast(const IntegerVectorList & l)100 static IntegerVectorList swapFirstLast(const IntegerVectorList &l)
101 {
102   IntegerVectorList ret;
103 
104   for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
105     {
106       IntegerVector v=*i;
107       int t=v[0];
108       v[0]=v[v.size()-1];
109       v[v.size()-1]=t;
110       ret.push_back(v);
111     }
112 
113   return ret;
114 }
115 
116 
isEmpty()117 bool HalfOpenCone::isEmpty()
118 {
119     bool ret1=!hasHomogeneousSolution(liftedDimension,
120 				   swapFirstLast(lifted.getHalfSpaces()),
121 				   swapFirstLast(lifted.getLinealitySpace())
122 				   );
123     /*
124   IntegerVectorList inequalityList;
125   inequalityList.push_back(IntegerVector::standardVector(liftedDimension,dimension));
126   PolyhedralCone temp=intersection(lifted,PolyhedralCone(inequalityList,IntegerVectorList(),liftedDimension));
127   IntegerVector v=temp.getRelativeInteriorPoint();
128   //  AsciiPrinter(Stderr).printVector(v);
129   bool ret2=(v[dimension]==0);
130     */
131 
132   /*  fprintf(Stderr,"Inequalities:\n");
133   AsciiPrinter(Stderr).printVectorList(lifted.getHalfSpaces());
134   fprintf(Stderr,"Equations:\n");
135   AsciiPrinter(Stderr).printVectorList(lifted.getLinealitySpace());
136   fprintf(Stderr,"hasSolution=%i\n",ret1);
137   */
138   //  assert(ret1==ret2);
139 
140   return ret1;
141 }
142 
143 
haveEmptyIntersection(const HalfOpenCone & a,const HalfOpenCone & b)144 bool haveEmptyIntersection(const HalfOpenCone &a, const HalfOpenCone &b)
145 {
146   assert(a.dimension==b.dimension);
147   IntegerVectorList inequalityList=a.lifted.getHalfSpaces();
148   IntegerVectorList equationList=a.lifted.getLinealitySpace();
149 
150   IntegerVectorList inequalityList2=b.lifted.getHalfSpaces();
151   IntegerVectorList equationList2=b.lifted.getLinealitySpace();
152 
153   inequalityList.splice(inequalityList.begin(),inequalityList2);
154   equationList.splice(equationList.begin(),equationList2);
155 
156   bool ret1=!hasHomogeneousSolution(a.liftedDimension,swapFirstLast(inequalityList),swapFirstLast(equationList));
157 
158   /*  HalfOpenCone c=intersection(a,b);
159   if(c.isEmpty()!=ret1)
160     {
161       AsciiPrinter(Stderr).printVectorList(inequalityList);
162       AsciiPrinter(Stderr).printVectorList(equationList);
163 
164       AsciiPrinter(Stderr).printVectorList(c.lifted.getHalfSpaces());
165       AsciiPrinter(Stderr).printVectorList(c.lifted.getLinealitySpace());
166 
167       fprintf(Stderr,"hasHomogeneousSolution siger %i\n",!ret1);
168 
169       assert(0);
170     }
171   */
172 
173   return ret1;
174 }
175 
176 
177 /*bool HalfOpenCone::isEmpty()
178 {
179   IntegerVector v(liftedDimension);
180   v[dimension]=-1;
181 
182   IntegerVectorList equationList,inequalityList;
183   inequalityList.push_back(v);
184   PolyhedralCone c(inequalityList,equationList,liftedDimension);
185   PolyhedralCone c2=intersection(c,lifted);
186   lifted.canonicalize();
187   c2.canonicalize();
188 
189   return !(c2!=lifted);
190   }*/
191 
intersection(const HalfOpenCone & a,const HalfOpenCone & b,bool findFacets)192 HalfOpenCone intersection(const HalfOpenCone &a, const HalfOpenCone &b, bool findFacets)
193 {
194   assert(a.dimension==b.dimension);
195 
196   /*  fprintf(Stderr,"-----------------------------------------------------------\n");
197   fprintf(Stderr,"Intersecting:\n");
198   AsciiPrinter P(Stderr);
199   printHalfOpenCone(P,a);
200   printHalfOpenCone(P,b);
201   */
202 
203   HalfOpenCone ret=HalfOpenCone(a.dimension,intersection(a.lifted,b.lifted));
204 
205   {
206     static int t;
207     t++;
208     if(!(t&7))ret.lifted.findFacets();
209 
210     //1 4:53
211     //3 3:38
212     //7
213   }
214 
215 
216   /*  fprintf(Stderr,"Result:\n");
217   printHalfOpenCone(P,ret);
218   fprintf(Stderr,"Is empty:%i\n",ret.isEmpty());
219   fprintf(Stderr,"-----------------------------------------------------------\n");
220   fprintf(Stderr,"States: %i,%i,%i\n",a.lifted.getState(),b.lifted.getState(),ret.lifted.getState());
221   fprintf(Stderr,"-----------------------------------------------------------\n");
222   */
223 
224   return ret;
225 }
226 
227 
shrink(const IntegerVectorList & l)228 IntegerVectorList HalfOpenCone::shrink(const IntegerVectorList &l)
229 {
230   IntegerVectorList ret;
231   for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
232     ret.push_back(i->subvector(0,i->size()-1));
233   return ret;
234 }
235 
236 
closure()237 PolyhedralCone HalfOpenCone::closure()
238 {
239   lifted.findFacets();
240   return PolyhedralCone(shrink(lifted.getHalfSpaces()),shrink(lifted.getLinealitySpace()),dimension);
241 }
242 
243 
244 /*
245  */
orientedBoundary(PolyhedralCone C,TermOrder const & t)246 HalfOpenConeList orientedBoundary(PolyhedralCone C, TermOrder const &t)
247 {
248   int dimension=C.ambientDimension();
249   HalfOpenConeList ret;
250   C.findFacets();
251   assert(C.dimension()==C.ambientDimension());
252 
253   IntegerVectorList facets=C.getHalfSpaces();
254 
255   IntegerVectorList strictList,nonStrictList;
256   for(IntegerVectorList::const_iterator i=facets.begin();i!=facets.end();i++)
257     {
258       if(t(*i,*i-*i))
259 	strictList.push_back(*i);
260       else
261 	nonStrictList.push_back(*i);
262     }
263 
264   // Let's make the non-strict inequalities strict one at a time and add a cone for each iteration
265   while(!nonStrictList.empty())
266     {
267       IntegerVector v=nonStrictList.front();
268       nonStrictList.pop_front();
269 
270 
271       IntegerVectorList equationList;
272       equationList.push_back(v);
273       ret.push_back(HalfOpenCone(dimension,equationList,nonStrictList,strictList,true));
274       strictList.push_back(v);
275     }
276   return ret;
277 }
278 
279 
tropicalHyperSurface(Polynomial const & p1)280 HalfOpenConeList tropicalHyperSurface(Polynomial const &p1)
281 {
282   Polynomial p=p1.homogenization();
283   HalfOpenConeList ret;
284   PolynomialSet g;
285   g.push_back(p);
286   buchberger(&g,LexicographicTermOrder());
287 
288   EnumerationTargetCollector gfan;
289   LexicographicTermOrder myTermOrder;
290   ReverseSearch rs(myTermOrder);
291 
292   rs.setEnumerationTarget(&gfan);
293 
294   fprintf(Stderr,"Starting enumeratioin\n");
295   rs.enumerate(g);
296   fprintf(Stderr,"Done\n");
297 
298   PolynomialSetList theList=gfan.getList();
299   for(PolynomialSetList::const_iterator i=theList.begin();i!=theList.end();i++)
300     {
301       HalfOpenConeList temp=orientedBoundary(groebnerCone(i->deHomogenization(),false),myTermOrder);
302       ret.splice(ret.begin(),temp);
303     }
304 
305   //  AsciiPrinter P(Stderr);
306   //  printHalfOpenConeList(P,ret);
307 
308   return ret;
309 }
310 
311 
refinement(HalfOpenConeList const & a,HalfOpenConeList const & b)312 HalfOpenConeList refinement(HalfOpenConeList const &a, HalfOpenConeList const &b)
313 {
314   HalfOpenConeList ret;
315   for(HalfOpenConeList::const_iterator i=a.begin();i!=a.end();i++)
316     for(HalfOpenConeList::const_iterator j=b.begin();j!=b.end();j++)
317       if(!haveEmptyIntersection(*i,*j))
318       {
319 	HalfOpenCone c=intersection(*i,*j);
320 	//	c.isEmpty();
321 	//	c.isEmpty();
322 	//	if(!c.isEmpty())
323 	  ret.push_back(c);
324       }
325 
326   return ret;
327 }
328 
329 
tropicalHyperSurfaceIntersection2(int dimension,PolynomialSet const & g)330 HalfOpenConeList tropicalHyperSurfaceIntersection2(int dimension, PolynomialSet const &g)
331 {
332   HalfOpenConeList intersection;
333   intersection.push_back(HalfOpenCone(dimension,IntegerVectorList(),IntegerVectorList(),IntegerVectorList()));
334 
335   for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
336     {
337       HalfOpenConeList surface=tropicalHyperSurface(*i);
338 
339       fprintf(Stderr,"Number of cones in current intersection:%i\n",intersection.size());
340       fprintf(Stderr,"Number of cones in next surface:%i\n",surface.size());
341 
342       fprintf(Stderr,"A\n");
343       intersection=refinement(intersection,surface);
344       fprintf(Stderr,"B\n");
345     }
346   fprintf(Stderr,"%i",intersection.size());
347 
348   return intersection;
349 }
350 
351 
tropicalHyperSurfaceIntersectionClosed(int dimension,PolynomialSet const & g)352 PolyhedralFan tropicalHyperSurfaceIntersectionClosed(int dimension, PolynomialSet const &g)
353 {
354   HalfOpenConeList intersection=tropicalHyperSurfaceIntersection(dimension,g);
355 
356 
357   AsciiPrinter P(Stderr);
358   printHalfOpenConeList(intersection, P);
359 
360   PolyhedralFan ret(dimension);
361   for(HalfOpenConeList::iterator i=intersection.begin();i!=intersection.end();i++)
362     {
363       PolyhedralCone c=i->closure();
364       c.canonicalize();
365       ret.insert(c);
366     }
367 
368   return ret;
369 }
370 
371 
splitIntoRelativelyOpenCones(list<HalfOpenCone> & l)372 void HalfOpenCone::splitIntoRelativelyOpenCones(list<HalfOpenCone> &l)
373 {
374   //  fprintf(Stderr,"BEGIN\n");
375   //  AsciiPrinter P(Stderr);
376   //  print(P);
377   lifted.findFacets();
378   //  print(P);
379 
380   /*	{
381 	  IntegerVector v=StringParser("(3,0,3,2,0,3)").parseIntegerVector();
382 	  if(contains(v))fprintf(stderr,"??????????????????????????????????????????\n");
383 	  }*/
384 
385   IntegerVectorList inequalityList=lifted.getHalfSpaces();
386   IntegerVectorList equationList=lifted.getLinealitySpace();
387 
388   IntegerVectorList strict,nonstrict;
389 
390 
391   for(IntegerVectorList::const_iterator i=inequalityList.begin();i!=inequalityList.end();i++)
392     if((*i)[i->size()-1]<0)
393       strict.push_back(*i);
394     else if((*i)[i->size()-1]==0)
395       nonstrict.push_back(*i);
396     else
397       {//CHANGED
398 	assert(i->subvector(0,i->size()-1).isZero());
399 	strict.push_back(*i);
400       }
401 
402   //  AsciiPrinter(Stderr).printVectorList(nonstrict);
403   //  AsciiPrinter(Stderr).printVectorList(strict);
404   //  AsciiPrinter(Stderr).printVectorList(equationList);
405 
406   if(nonstrict.size()==0)
407     {
408       l.push_back(*this);
409     }
410   else
411     {
412       IntegerVector chosen=*nonstrict.begin();
413       nonstrict.pop_front();
414 
415       strict.push_front(chosen);
416       (*strict.begin())[strict.begin()->size()-1]=-1;
417       IntegerVectorList a=nonstrict;
418       IntegerVectorList tempa=strict;
419       a.splice(a.begin(),tempa);
420 
421       //  fprintf(Stderr,"New inequalities:\n");
422       //  AsciiPrinter(Stderr).printVectorList(a);
423 
424       HalfOpenCone A(dimension,PolyhedralCone(a,equationList,liftedDimension));
425       A.splitIntoRelativelyOpenCones(l);
426 
427 
428       equationList.push_front(chosen);
429       strict.pop_front();
430 
431       IntegerVectorList b=nonstrict;
432       IntegerVectorList tempb=strict;
433       b.splice(b.begin(),tempb);
434       // fprintf(Stderr,"New inequalities:\n");
435       // AsciiPrinter(Stderr).printVectorList(b);
436       // fprintf(Stderr,"New equationList:\n");
437       // AsciiPrinter(Stderr).printVectorList(equationList);
438 
439       HalfOpenCone B(dimension,PolyhedralCone(b,equationList,liftedDimension));
440       B.splitIntoRelativelyOpenCones(l);
441     }
442       //  AsciiPrinter(Stderr).print
443   //  fprintf(Stderr,"END\n");
444 }
445 
446 
print(class Printer & p) const447 void HalfOpenCone::print(class Printer &p)const
448 {
449   p.printString("Printing HalfOpenCone\n");
450   lifted.print(&p);
451   p.printString("Done printing HalfOpenCone\n");
452 }
453 
splitIntoRelativelyOpenCones(HalfOpenConeList const & l)454 HalfOpenConeList splitIntoRelativelyOpenCones(HalfOpenConeList const &l)
455 {
456   AsciiPrinter P(Stderr);
457   HalfOpenConeList ret;
458   for(HalfOpenConeList::const_iterator i=l.begin();i!=l.end();i++)
459     {
460       fprintf(Stderr,"A");
461       HalfOpenCone temp=*i;
462       HalfOpenConeList tempSplit;
463       //  fprintf(Stderr,"---------------------------------------------------------------\n");
464       //  temp.print(P);
465       //  fprintf(Stderr,"---------------------------------------------------------------\n");
466       temp.splitIntoRelativelyOpenCones(tempSplit);
467 
468       //  fprintf(Stderr,"Splits into:");
469       //  for(HalfOpenConeList::const_iterator i=tempSplit.begin();i!=tempSplit.end();i++)
470       //    i->print(P);
471       //  fprintf(Stderr,"Splits into End.");
472 
473       ret.splice(ret.begin(),tempSplit);
474 
475       fprintf(Stderr,"B\n");
476     }
477 
478   return ret;
479 }
480 
isSubsetOf(IntegerVector const & v,IntegerVector const & u)481 static bool isSubsetOf(IntegerVector const &v, IntegerVector const &u)
482 {
483   for(int i=0;i<v.size();i++)
484     {
485       bool found=false;
486       for(int j=0;j<u.size();j++)
487 	if(v[i]==u[j])
488 	  {
489 	    found=true;
490 	    break;
491 	  }
492       if(!found)return false;
493     }
494   return true;
495 }
496 
isSubsetOf(IntegerVector const & v,IntegerVectorList const & l)497 static bool isSubsetOf(IntegerVector const &v, IntegerVectorList const &l)
498 {
499   for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
500     if(isSubsetOf(v,*i))return true;
501 
502   return false;
503 }
504 
505 
printHalfOpenConeList(HalfOpenConeList const & l,class Printer & p)506 void printHalfOpenConeList(HalfOpenConeList const &l, class Printer & p)
507 {
508   HalfOpenConeList L=splitIntoRelativelyOpenCones(l);
509 
510   list<PolyhedralCone> cones;
511   for(HalfOpenConeList::iterator i=L.begin();i!=L.end();i++)
512     cones.push_back(i->closure());
513 
514   int homog=1000000;
515   int largest=0;
516   int ambientDimension=-1;
517   for(list<PolyhedralCone>::const_iterator i=cones.begin();i!=cones.end();i++)
518     {
519       if(i->dimension()<homog)homog=i->dimension();
520       if(i->dimension()>largest)largest=i->dimension();
521       ambientDimension=i->ambientDimension();
522     }
523   assert(homog!=1000000);
524 
525   for(list<PolyhedralCone>::const_iterator i=cones.begin();i!=cones.end();i++)
526     {
527       assert(i->dimensionOfLargestContainedSubspace()==homog);
528     }
529   fprintf(stderr,"Ambient dimension: %i, maximal dimension: %i, dimension of lineality space: %i\n",ambientDimension,largest,homog);
530 
531   IntegerVectorList rays;
532   for(list<PolyhedralCone>::iterator i=cones.begin();i!=cones.end();i++)
533     if(i->dimension()==homog+1)rays.push_back(i->getRelativeInteriorPoint());
534 
535   p.printString("Rays:\n");
536 
537   p.printVectorList(rays,true);
538 
539   list<IntegerVectorList> subsets;
540   for(int d=homog;d<=largest;d++)
541     {
542       IntegerVectorList thisDimension;
543       list<PolyhedralCone> cones2;
544       for(list<PolyhedralCone>::iterator i=cones.begin();i!=cones.end();i++)
545 	if(i->dimension()==d)
546 	  cones2.push_back(*i);
547 
548       for(list<PolyhedralCone>::const_iterator i=cones2.begin();i!=cones2.end();i++)
549 	{
550 	  IntegerVector v(0);
551 	  int J=0;
552 	  for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
553 	    {
554 	      if(i->contains(*j))
555 		{
556 		  v.grow(v.size()+1);
557 		  v[v.size()-1]=J;
558 		}
559 	      J++;
560 	    }
561 	  thisDimension.push_back(v);
562 	}
563       subsets.push_back(thisDimension);
564     }
565 
566 
567   list<IntegerVectorList>::const_iterator subsetIterator=subsets.begin();
568   list<IntegerVectorList>::const_iterator subsetIteratorNext=subsets.begin();
569   for(int d=homog;d<=largest;d++)
570     {
571       subsetIteratorNext++;
572       IntegerVectorList maximal,nonmaximal;
573 
574       if(subsetIteratorNext!=subsets.end())
575 	{
576 	  for(IntegerVectorList::const_iterator i=subsetIterator->begin();i!=subsetIterator->end();i++)
577 	    if(isSubsetOf(*i,*subsetIteratorNext))
578 	      nonmaximal.push_back(*i);
579 	    else
580 	      maximal.push_back(*i);
581 	}
582       else
583 	maximal=*subsetIterator;
584 
585       p.printString("Printing ");p.printInteger(subsetIterator->size());p.printString(" ");p.printInteger(d);p.printString("-dimensional cones (");p.printInteger(maximal.size());p.printString(" maximal cones):\n");
586       p.printString("{");
587       {
588 	bool first=true;
589 	for(IntegerVectorList::const_iterator i=maximal.begin();i!=maximal.end();i++)
590 	  {
591 	    if(!first)p.printString(",\n");
592 	    p.printVector(*i);
593 	    first=false;
594 	  }
595 	if(!first && nonmaximal.size()!=0)p.printString(",\n");
596 	p.printString("\n");
597 	first=true;
598 	for(IntegerVectorList::const_iterator i=nonmaximal.begin();i!=nonmaximal.end();i++)
599 	  {
600 	    if(!first)
601 	      p.printString(",\n");
602 	    p.printVector(*i);
603 	    first=false;
604 	  }
605       }
606       p.printString("}\n");
607       subsetIterator++;
608    }
609   /*  for(int d=homog;d<=largest;d++)
610     {
611       list<PolyhedralCone> cones2;
612       for(list<PolyhedralCone>::iterator i=cones.begin();i!=cones.end();i++)
613 	if(i->dimension()==d)
614 	  cones2.push_back(*i);
615 
616       p.printString("Printing ");p.printInteger(cones2.size());p.printString(" ");p.printInteger(d);p.printString("-dimensional cones:\n");
617       p.printString("{");
618       for(list<PolyhedralCone>::const_iterator i=cones2.begin();i!=cones2.end();i++)
619 	{
620 	  IntegerVector v(0);
621 	  int J=0;
622 	  for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
623 	    {
624 	      if(i->contains(*j))
625 		{
626 		  v.grow(v.size()+1);
627 		  v[v.size()-1]=J;
628 		}
629 	      J++;
630 	    }
631 	  if(i!=cones2.begin())p.printString(",\n");
632 	  p.printVector(v);
633 	}
634       p.printString("}\n");
635     }
636 */
637 }
638 
639 
640 class BitSet
641 {
642   vector<int> v;
643 public:
BitSet()644   BitSet()
645   {
646   }
BitSet(int n)647   BitSet(int n):
648     v(n)
649   {
650     for(int i=0;i<n;i++)v[i]=0;
651   }
operator [](int n)652   int& operator[](int n){assert(n>=0 && n<v.size());return (v[n]);}
operator [](int n) const653   const int& operator[](int n)const{assert(n>=0 && n<v.size());return (v[n]);}
add(BitSet const & b)654   void add(BitSet const &b)
655   {
656     assert(b.size()==v.size());
657     for(int i=0;i<b.size();i++)
658       {
659 	//	fprintf(stderr,"%i\n",i);
660 	v[i]=v[i]||b[i];
661       }
662   }
size() const663   int size()const
664   {
665     return v.size();
666   }
print(Printer & P) const667   void print(Printer &P)const
668   {
669     P.printString("(");
670     for(int i=0;i<v.size();i++)
671       {
672 	if(i!=0)P.printString(", ");
673 	P.printInteger(v[i]);
674       }
675     P.printString(")\n");
676   }
negated() const677   BitSet negated()const
678   {
679     BitSet ret(size());
680     for(int i=0;i<size();i++)ret[i]=1-v[i];
681     return ret;
682   }
sizeOfSubset() const683   int sizeOfSubset()const
684   {
685     int ret=0;
686     for(int i=0;i<size();i++)if(v[i])ret++;
687     return ret;
688   }
689 };
690 
691 class Table
692 {
693   vector<vector<vector<BitSet> > > table;
694 public:
Table(vector<vector<HalfOpenCone>> const & l)695   Table(vector<vector<HalfOpenCone> > const &l):
696     table(l.size())
697   {
698     int N=l.size();
699     for(int i=0;i<N;i++)
700       {
701 	vector<vector<BitSet> > v(N);
702 	for(int j=0;j<N;j++)
703 	  {
704 	    vector<BitSet> w(l[i].size());
705 	    for(int k=0;k<l[i].size();k++)
706 	      w[k]=BitSet(l[j].size());
707 	    v[j]=w;
708 	  }
709 	table[i]=v;
710       }
711   }
lookUp(int fan1,int cone1,int fan2,int cone2)712   bool lookUp(int fan1, int cone1, int fan2, int cone2)
713   {
714     assert(fan1<table.size());
715     assert(fan2<table[fan1].size());
716     assert(cone1<table[fan1][fan2].size());
717     assert(cone2<table[fan1][fan2][cone1].size());
718 
719     return table[fan1][fan2][cone1][cone2];
720   }
set(int fan1,int cone1,int fan2,int cone2)721   void set(int fan1, int cone1, int fan2, int cone2)
722   {
723     assert(fan1<table.size());
724     assert(fan2<table[fan1].size());
725     assert(cone1<table[fan1][fan2].size());
726     assert(cone2<table[fan1][fan2][cone1].size());
727 
728     table[fan1][fan2][cone1][cone2]=true;
729     table[fan2][fan1][cone2][cone1]=true;
730   }
nonCandidates(int fan1,int cone1,int fan2) const731   BitSet const& nonCandidates(int fan1, int cone1, int fan2)const
732   {
733     assert(fan1<table.size());
734     assert(fan2<table[fan1].size());
735     assert(cone1<table[fan1][fan2].size());
736 
737     return table[fan1][fan2][cone1];
738   }
print() const739   void print()const
740   {
741     AsciiPrinter P(Stderr);
742     for(int i=0;i<table.size();i++)
743       for(int j=0;j<table[i].size();j++)
744 	{
745 	  fprintf(Stderr,"Entry (%i,%i)\n",i,j);
746 	  for(int k=0;k<table[i][j].size();k++)
747 	    table[i][j][k].print(P);
748 	}
749   }
750 };
751 
752 class RelationTable
753 {
754   vector<vector<HalfOpenCone> > fanList;
755   Table knownEmptyIntersectionInIntersection;
756   Table knownNonEmptyIntersection;
757 public:
758   int numberOfSolvedLPs;
RelationTable(vector<vector<HalfOpenCone>> const & l)759   RelationTable(vector<vector<HalfOpenCone> > const &l):
760     fanList(l),
761     knownEmptyIntersectionInIntersection(l),
762     knownNonEmptyIntersection(l),
763     numberOfSolvedLPs(0)
764   {
765 
766   }
intersectTriviallyInIntersection(int fan1,int cone1,int fan2,int cone2)767   bool intersectTriviallyInIntersection(int fan1, int cone1, int fan2, int cone2)
768   {
769     assert(fan1<fanList.size());
770     assert(fan2<fanList.size());
771     assert(cone1<fanList[fan1].size());
772     assert(cone2<fanList[fan2].size());
773 
774     if(knownEmptyIntersectionInIntersection.lookUp(fan1,cone1,fan2,cone2))
775       return true;
776     if(knownNonEmptyIntersection.lookUp(fan1,cone1,fan2,cone2))
777       return false;
778 
779     //    fprintf(Stderr,"UPDATING:f1:%i,c1:%i,f2:%i,c2:%i\n",fan1,cone1,fan2,cone2);
780     bool ret=haveEmptyIntersection(fanList[fan1][cone1],fanList[fan2][cone2]);
781     numberOfSolvedLPs++;
782     if(ret)
783       knownEmptyIntersectionInIntersection.set(fan1,cone1,fan2,cone2);
784     else
785       knownNonEmptyIntersection.set(fan1,cone1,fan2,cone2);
786     return ret;
787   }
getNonCandidates(int fan1,int cone1,int fan2)788   const BitSet &getNonCandidates(int fan1, int cone1, int fan2)
789   {
790     for(int c2=0;c2<fanList[fan2].size();c2++)
791       intersectTriviallyInIntersection(fan1,cone1,fan2,c2);
792 
793     return knownEmptyIntersectionInIntersection.nonCandidates(fan1,cone1,fan2);
794   }
markNoIntersectionInIntersection(int fan1,int cone1,int fan2,int cone2)795   void markNoIntersectionInIntersection(int fan1, int cone1, int fan2, int cone2)
796   {
797     knownEmptyIntersectionInIntersection.set(fan1,cone1,fan2,cone2);
798   }
print() const799   void print()const
800   {
801     fprintf(Stderr,"knownEmptyIntersectionInIntersection:");
802     knownEmptyIntersectionInIntersection.print();
803     fprintf(Stderr,"knownNonEmptyIntersection:");
804     knownNonEmptyIntersection.print();
805   }
806 };
807 
808 
809 struct RecursionData
810 {
811   vector<vector<HalfOpenCone> > fans;
812   IntegerVector chosen;
813   IntegerVector chosenFans;
814   IntegerVector iterators; //just used for printing
815   IntegerVector nCandidates; //just used for printing
816   BitSet usedFans;
817   int numberOfUsefulCalls;
818   int totalNumberOfCalls;
819 public:
820   RelationTable table;
RecursionDataRecursionData821   RecursionData(vector<vector<HalfOpenCone> > const &fans_):
822     table(fans_),
823     fans(fans_),
824     chosen(fans_.size()),
825     chosenFans(fans_.size()),
826     usedFans(fans_.size()),
827     iterators(fans_.size()),
828     nCandidates(fans_.size()),
829     numberOfUsefulCalls(0),
830     totalNumberOfCalls(0)
831   {
832   }
833 
834   HalfOpenConeList ret;
835 
computeCandidatesRecursionData836   BitSet computeCandidates(int index, int fanNumber)
837   {
838     BitSet nonCandidates(fans[fanNumber].size());
839     for(int i=0;i<index;i++)
840       {
841 	nonCandidates.add(table.getNonCandidates(chosenFans[i],chosen[i],fanNumber));
842       }
843     return nonCandidates.negated();
844   }
rekRecursionData845   bool rek(int index, HalfOpenCone const &current)
846   {
847     totalNumberOfCalls++;
848 
849     bool success=false;
850 
851     if(index == fans.size())
852       {
853 	fprintf(Stderr,"ADDING CONE\n");
854 	AsciiPrinter P(Stderr);
855 	//  current.print(P);
856 	ret.push_back(current);
857 	numberOfUsefulCalls++;
858 	return true;
859       }
860     else
861       {
862 	AsciiPrinter P(Stderr);
863 
864 	int bestIndex=-1;
865 	int bestNumberOfCandidates=1000000;
866 	for(int i=0;i<fans.size();i++)
867 	  {
868 	    if(!usedFans[i])
869 	      {
870 		int n=computeCandidates(index,i).sizeOfSubset();
871 		//		fprintf(Stderr,"Number of candidates for fan %i: %i\n",i,n);
872 		if(n<=bestNumberOfCandidates)  //we could choose a strict inequality
873 		  {
874 		    bestNumberOfCandidates=n;
875 		    bestIndex=i;
876 		  }
877 	      }
878 	  }
879 	assert(bestIndex!=-1);
880 	BitSet candidates=computeCandidates(index,bestIndex);
881 
882 
883 	chosenFans[index]=bestIndex;
884 	usedFans[chosenFans[index]]=true;
885 
886 
887 	nCandidates[index]=bestNumberOfCandidates;//just for printing
888 
889 	static int iterationNumber;
890 	if(!(iterationNumber++ & 31))
891 	{
892 	  fprintf(Stderr,"Iteration level:%i, Chosen fan:%i, Number of candidates:%i, Iteration Number:%i, Useful (%i/%i)=%f\n",index,bestIndex,bestNumberOfCandidates,iterationNumber,numberOfUsefulCalls,totalNumberOfCalls,float(numberOfUsefulCalls)/totalNumberOfCalls);
893 	  fprintf(Stderr,"Chosen fans vector: ");
894 	  P.printVector(chosenFans,false,2);
895 	  fprintf(Stderr,"\nChosen cone vector: ");
896 	  P.printVector(chosen,false,2);
897 	  fprintf(Stderr,"\nNcandidates vector: ");
898 	  P.printVector(nCandidates,false,2);
899 	  fprintf(Stderr,"\nIterator vector:    ");
900 	  P.printVector(iterators,false,2);
901 	  fprintf(Stderr,"\n\n");
902 	}
903 
904 	//
905 	/*	if(index>1)
906 	  {
907 	    smallest=1000000;
908 	    bx=-1;
909 	    by=-1;
910 	    for(int x=0;x<index;x++)
911 	      for(int y=x;y<index;y++)
912 		{
913 		  BitSet c1=table.getNonCandidates(chosenFans[x],chosen[x],bestIndex);
914 		  c1.add(table.getNonCandidates(chosenFans[y],chosen[y],bestIndex));
915 		  int n=c1.sizeOfSubset();
916 		  if(n<smallest)
917 		    {
918 		      bx=x;
919 		      by=y;
920 		    }
921 		}
922 	    bool skipThisOne=true;
923 	    for(int i=0;i<fans[chosenFans[index]].size();i++)
924 	      {
925 		if(!table.intersectTriviallyInIntersection(chosenFans[x],chosen[x],chosenFans[index],i))skipThisOne=false;
926 		if(!table.intersectTriviallyInIntersection(chosenFans[y],chosen[y],chosenFans[index],i))skipThisOne=false;
927 	      }
928 	      }*/
929 
930 
931 	// P.printInteger(fans[index].size());
932 	for(int i=0;i<fans[chosenFans[index]].size();i++)
933 	  if(candidates[i])
934 	    {
935 	      bool ok=true;
936 	      for(int j=0;j<index;j++)
937 		{
938 		  if(table.intersectTriviallyInIntersection(chosenFans[j],chosen[j],chosenFans[index],i))
939 		    {
940 		      ok=false;
941 		      break;
942 		    }
943 		}
944 	      if(ok && !haveEmptyIntersection(current,fans[chosenFans[index]][i]))
945 		{
946 		  chosen[index]=i;
947 
948 		  HalfOpenCone next=intersection(current,fans[chosenFans[index]][i],false);
949 		  if((fans.size()!=10)&&(fans.size()!=9))
950 		    {
951 		      bool s=rek(index+1,next);
952 		      success|=s;
953 		    }
954 		  else
955 		    {
956 		      vector<vector<HalfOpenCone> > L2;
957 		      {
958 			for(int i=0;i<fans.size();i++)
959 			  if(!usedFans[i])
960 			    {
961 			      vector<HalfOpenCone> L;
962 			      for(vector<HalfOpenCone>::const_iterator j=fans[i].begin();j!=fans[i].end();j++)
963 				{
964 				  if(!haveEmptyIntersection(next,*j))
965 				    {
966 				      L.push_back(intersection(next,*j,true));
967 				    }
968 				}
969 			      fprintf(stderr,"New fan size:%i\n",L.size());
970 			      L2.push_back(L);
971 			    }
972 		      }
973 		      RecursionData data(L2);
974 		      data.completeTable();
975 		      success|=data.rek(0, HalfOpenCone(next.dimension,IntegerVectorList(),IntegerVectorList(),IntegerVectorList()));
976 		      ret.splice(ret.begin(),data.ret);
977 		    }
978 		  chosen[index]=-1;//just for printing
979 		}
980 	      iterators[index]++;//just for printing
981 	    }
982 	/*	if(!success)
983 	  {
984 	    for(int x=0;x<index;x++)
985 	      for(int y=x;y<index;y++)
986 		{
987 		  BitSet c1=table.getNonCandidates(chosenFans[x],chosen[x],bestIndex);
988 		  c1.add(table.getNonCandidates(chosenFans[y],chosen[y],bestIndex));
989 		  int n=c1.negated().sizeOfSubset();
990 		  if(n==0)
991 		    {
992 		      table.markNoIntersectionInIntersection(chosenFans[x],chosen[x],chosenFans[y],chosen[y]);
993 		      fprintf(stderr,"ONE FOR FREE\n");
994 		    }
995 		}
996 		}*/
997 
998 	nCandidates[index]=-1;//just for printing
999 	iterators[index]=0;//just for printing
1000 
1001 	usedFans[chosenFans[index]]=false;
1002 	chosenFans[index]=-1;
1003       }
1004     if(success)numberOfUsefulCalls++;
1005     return success;
1006   }
transitiveClosureRecursionData1007   void transitiveClosure()
1008   {
1009     bool rep;
1010     do{
1011       rep=false;
1012       int a=0;
1013     for(int f1=0;f1<fans.size();f1++)
1014       {
1015 	//	fprintf(stderr,"%i\n",f1);
1016 	for(int f2=f1+1;f2<fans.size();f2++)
1017 	  for(int c1=0;c1<fans[f1].size();c1++)
1018 	    for(int c2=0;c2<fans[f2].size();c2++)
1019 	      {
1020 		if(!table.intersectTriviallyInIntersection(f1,c1,f2,c2))
1021 		  {
1022 		    bool dontintersect=false;
1023 		    for(int f3=0;f3<fans.size();f3++)
1024 		      {
1025 			BitSet c=table.getNonCandidates(f1,c1,f3);
1026 			c.add(table.getNonCandidates(f2,c2,f3));
1027 			if(c.negated().sizeOfSubset()==0)
1028 			  {
1029 			    dontintersect=true;
1030 			    a++;
1031 			    break;
1032 			  }
1033 		      }
1034 		    if(dontintersect)
1035 		      {
1036 			table.markNoIntersectionInIntersection(f1,c1,f2,c2);
1037 			rep=true;
1038 		      }
1039 		  }
1040 	      }
1041       }
1042     fprintf(stderr,"%i FOR FREE\n",a);
1043     }
1044     while(rep);
1045   }
1046 
completeTableRecursionData1047   void completeTable()
1048   {
1049     for(int f1=0;f1<fans.size();f1++)
1050       {
1051 	fprintf(stderr,"%i\n",f1);
1052 	for(int f2=f1+1;f2<fans.size();f2++)
1053 	for(int c1=0;c1<fans[f1].size();c1++)
1054 	  for(int c2=0;c2<fans[f2].size();c2++)
1055 	    table.intersectTriviallyInIntersection(f1,c1,f2,c2);
1056       }
1057     transitiveClosure();
1058   }
1059 };
1060 
tropicalHyperSurfaceIntersection(int dimension,PolynomialSet const & g)1061 HalfOpenConeList tropicalHyperSurfaceIntersection(int dimension, PolynomialSet const &g)
1062 {
1063   vector<vector<HalfOpenCone> > L2;
1064   {
1065     for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
1066       {
1067 	HalfOpenConeList l=tropicalHyperSurface(*i);
1068 	vector<HalfOpenCone> L;
1069 	for(HalfOpenConeList::const_iterator i=l.begin();i!=l.end();i++)
1070 	  {
1071 	    L.push_back(*i);
1072 	  }
1073 	L2.push_back(L);
1074       }
1075   }
1076 
1077   RecursionData data(L2);
1078   //  data.completeTable();
1079   data.rek(0, HalfOpenCone(dimension,IntegerVectorList(),IntegerVectorList(),IntegerVectorList()));
1080   fprintf(stderr,"LPs solved:%i for relation table\n",data.table.numberOfSolvedLPs);
1081   return data.ret;
1082 }
1083