1 /*
2 * gfanlib_tropicalhomotopy.h
3 *
4 * Created on: Apr 10, 2016
5 * Author: anders
6 */
7
8 #ifndef GFANLIB_TROPICALHOMOTOPY_H_
9 #define GFANLIB_TROPICALHOMOTOPY_H_
10
11 #include "gfanlib_paralleltraverser.h"
12 #include "gfanlib_matrix.h"
13
14 namespace gfan{
15
16 class Flags{
17 public:
18 static const bool computeDotProductInMatrix=true;
19 };
20
21 /**
22 * We identify six possibly different types needed with possibly varying precission:
23 * 1) The entries of the circuits (or possibly their packed representation)
24 * 2) The mixed volume contribution of a single cell. This is obtained from an entry of a circuit and therefore can be represented by the above type.
25 * 3) The accumulated mixed volume. This will exceed the bound of the above type in many cases. Overflows are easily checked.
26 * 4) The type that dotVector uses as a result when dotting with the target. (Also used in campareInequalities)
27 * 5) The intermediate type for dotVector.
28 * 6) The type used in compareRevLexInverted
29 *
30 *
31 * Type 1 and 2 are the same.
32 * Type 3 is typically different.
33 *
34 * To simplify our design:
35 * we assume that type 4 is the same as 1 and 2. This is reasonable, as we need some bound to make type 6 efficient.
36 * we use a special (longer) type for 5, as that allows to do overflow checks at the end, assuming some bound on the target.
37 * In 6, we observe that there is no accumulation taking place. Moreover, with the assumption that 4 and 1 are the same, we only need a type with double precission to do the comparisons here, and now overflow check will be required.
38 *
39 *
40 * To conclude, we make two types. A single precision type for 1,2,4 and a double precision type for 3,5,6
41 * We further need to make assumptions on the absolute value of the entries of the target vector and the number of entries in combination to ensure that dot product computations do not overflow.
42 * Overflow checks are then only needed:
43 * when casting the return value of dotVector
44 * when doing the dotDivVector operation. But since bounds are known, in most cases checks are not needed
45 * when accumulating the mixed volume
46 *
47 *
48 * Suggested implementations:
49 * a pair of 32/64 bit integers with only the overflow checking listed above
50 * a pair of gmp integers for exact precision.
51 */
52
53
54
55
56
57 template<class mvtyp>
simplex(int n,mvtyp d)58 static Matrix<mvtyp> simplex(int n, mvtyp d)
59 {
60 Matrix<mvtyp> ret(n,n+1);
61 for(int i=0;i<n;i++)ret[i][i+1]=d;
62 return ret;
63 }
64
65
66
67 template<class mvtyp, class mvtypDouble, class mvtypDivisor>
68 class SingleTropicalHomotopyTraverser{
69 class InequalityComparisonResult{//actual comparison functions were moved to the InequalityTable
70 public:
71 bool empty;
72 int configurationIndex;//used for finding best
73 int columnIndex;//used for finding best
74 };
75 class InequalityTable //circuit table // This table has been moved inside the IntegersectionTraverser simply because it is used nowhere else and is specific to mixed cells in Cayley configurations.
76 {
77 /* All methods are marked to show if they can overflow without throwing/asserting.
78 * Overflowing methods at the moment are:
79 * replaceFirstOrSecond: subroutine calls may overflow (dotDivVector)
80 * compareInequalities: only if target entries are too big
81 * dotVector: only if target entries are too big
82 * setChoicesFromEarlierHomotopy: only if tuple entries are too big
83 */
84 std::vector<Matrix<mvtyp> > tuple;
85 std::vector<int> offsets;
86 std::vector<std::pair<int,int> > choices;
87 Matrix<mvtyp> A;//one row for each potential inequality, to store entries with indices in chosen
88 Vector<mvtyp> tempA;
89 Vector<mvtyp> Abounds;// a negative bound for each row of A, bounding the absolute value of the rows;
90 std::vector<mvtyp> svec;//used locally
91 int subconfigurationIndex;
92 mvtyp denominator;
93 int m;
94 int k;
isLegalIndex(int subconfigurationIndex,int columnIndex)95 bool isLegalIndex(int subconfigurationIndex, int columnIndex)const
96 {// Cannot overflow
97 return choices[subconfigurationIndex].first!=columnIndex && choices[subconfigurationIndex].second!=columnIndex;
98 }
99 mvtyp dotVector(int subconfigurationIndex, int columnIndex, Vector<mvtyp> const &target, int onlyK=-1)const
100 { // May overflow if entries of target are too big.
101 //if onlyK!=-1 then only the onlyKth subconfiguration is considered
102 mvtypDouble ret;
103 if(onlyK!=-1)
104 {
105 if(onlyK==subconfigurationIndex)
106 {
107 int i=subconfigurationIndex;
108 ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
109 ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
110 ret-=extendedMultiplication(denominator,target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));// the multiplication can be merged with multiplication above except that that could cause and overflow.
111 ret+=extendedMultiplication(denominator,target.UNCHECKEDACCESS(columnIndex+offsets[i]));
112 return ret.castToSingle();
113 }
114 else
115 {
116 int i=onlyK;
117 if(target.UNCHECKEDACCESS((choices)[i].first+offsets[i]).isNonZero())
118 {
119 ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
120 ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
121 }
122 return ret.castToSingle();
123 }
124 }
125 for(int i=0;i<tuple.size();i++)
126 {
127 ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
128 if(i==subconfigurationIndex)
129 {
130 ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
131 ret-=extendedMultiplication(denominator,target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));// the multiplication can be merged with multiplication above except that that could cause and overflow.
132 ret+=extendedMultiplication(denominator,target.UNCHECKEDACCESS(columnIndex+offsets[i]));
133 }
134 else
135 {
136 ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
137 }
138 }
139 return ret.castToSingle();
140 }
141 void assignDotProducts(Vector<mvtyp> const &target, int onlyK=-1)
142 {// Cannot overflow
143 int J=0;
144 for(int i=0;i<k;i++)
145 for(int j=0;j<tuple[i].getWidth();j++,J++)
146 A[k][J]=dotVector(i,j,target,onlyK);
147 }
isReverseLexInvertedLessThanZero(int subconfigurationIndex,int columnIndex)148 bool isReverseLexInvertedLessThanZero(int subconfigurationIndex, int columnIndex)const __attribute__ ((always_inline))//As in ReverseLexicographicInvertedTermOrder. Compare against zero
149 {// Cannot overflow
150 int i;
151 int index=columnIndex+offsets[subconfigurationIndex];
152 for(i=0;i<subconfigurationIndex;i++)
153 if(A.UNCHECKEDACCESS(i,index).isNonZero())
154 {
155 if(choices[i].first<choices[i].second)
156 return A.UNCHECKEDACCESS(i,index).isNegative();
157 else
158 return A.UNCHECKEDACCESS(i,index).isPositive();
159 }
160
161 mvtyp a=A.UNCHECKEDACCESS(i,index);
162 {
163 int firstIndex=choices[i].first;
164 int secondIndex=choices[i].second;
165 int thirdIndex=columnIndex;
166 mvtyp firstValue=-a-denominator;
167 mvtyp secondValue=a;
168 mvtyp thirdValue=denominator;
169
170 // Bubble sort
171 if(secondIndex<firstIndex)
172 {
173 std::swap(secondIndex,firstIndex);
174 std::swap(secondValue,firstValue);
175 }
176 if(thirdIndex<secondIndex)
177 {
178 std::swap(secondIndex,thirdIndex);
179 std::swap(secondValue,thirdValue);
180 }
181 if(secondIndex<firstIndex)
182 {
183 std::swap(secondIndex,firstIndex);
184 std::swap(secondValue,firstValue);
185 }
186
187 if(firstValue.isNonZero())
188 return firstValue.isPositive();
189 if(secondValue.isNonZero())
190 return secondValue.isPositive();
191 if(thirdValue.isNonZero())
192 return thirdValue.isPositive();
193 }
194 i++;
195 for(;i<k;i++)
196 if(A.UNCHECKEDACCESS(i,index).isNonZero())
197 {
198 if(choices[i].first<choices[i].second)
199 return A.UNCHECKEDACCESS(i,index).isNegative();
200 else
201 return A.UNCHECKEDACCESS(i,index).isPositive();
202 }
203
204
205 return false;
206 }
207 public:
computeABounds()208 void computeABounds()
209 {//Cannot overflow
210 for(int i=0;i<A.getHeight();i++)
211 Abounds[i]=mvtyp::computeNegativeBound(&(A[i][0]),A.getWidth());
212 }
checkABounds()213 void checkABounds()const//remove?
214 {//Cannot overflow
215 for(int i=0;i<A.getHeight();i++)
216 {
217 mvtyp M=0;
218 mvtyp m=0;
219 for(int j=0;j<A.getWidth();j++)
220 {
221 if(M<=A[i][j])M=A[i][j];
222 if(A[i][j]<=m)m=A[i][j];
223 }
224 assert(Abounds[i]<=m);
225 assert(Abounds[i]<=-M);
226 }
227 }
getCoordinateOfInequality(int subconfigurationIndex,int columnIndex,int i,int j)228 mvtypDouble getCoordinateOfInequality(int subconfigurationIndex, int columnIndex, int i, int j)const
229 {// Cannot overflows
230 //get (i,j)th coordinate of (subconfigurationIndex,columnIndex)th inequality
231 if(i==subconfigurationIndex)
232 {
233 if(choices[i].first==j)return -A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend()-denominator.extend();
234 else if(choices[i].second==j)return A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
235 else if(j==columnIndex)return denominator.extend();
236 else return mvtypDouble(0);
237 }
238 else
239 if(choices[i].first==j)return -A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
240 else if(choices[i].second==j)return A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
241 else return mvtypDouble(0);
242 }
sort2uniquely(int * v,int a,int b)243 int sort2uniquely(int *v, int a, int b)const//a and b different
244 {// Cannot overflow
245 v[a>b]=a;
246 v[b>a]=b;
247 return 2;
248 }
sort3uniquely(int * v,int a,int b,int c)249 int sort3uniquely(int *v, int a, int b, int c)const//a and b and c different
250 {// Cannot overflow
251 v[(a>b)+int(a>c)]=a;
252 v[(b>a)+int(b>c)]=b;
253 v[(c>a)+int(c>b)]=c;
254 return 3;
255 }
sort4uniquely(int * v,int a,int b,int c,int d)256 int sort4uniquely(int *v, int a, int b, int c, int d)const// a and b different and different from c and d, but c may equal d
257 {// Cannot overflow
258 if(c!=d)
259 {
260 v[(a>b)+int(a>c)+int(a>d)]=a;
261 v[(b>a)+int(b>c)+int(b>d)]=b;
262 v[(c>a)+int(c>b)+int(c>d)]=c;
263 v[(d>a)+int(d>b)+int(d>c)]=d;
264 return 4;
265 }
266 else return sort3uniquely(v,a,b,c);
267 }
compareReverseLexicographicInverted(int i1,int j1,int i2,int j2,mvtyp s1,mvtyp s2)268 bool compareReverseLexicographicInverted(int i1, int j1, int i2, int j2, mvtyp s1, mvtyp s2)const//s1 and s2 are always negative
269 {// cannot overflow
270 for(int i=0;i<k;i++)
271 {
272 if(i1!=i && i2!=i)
273 {
274 int temp=determinantSign1(A.UNCHECKEDACCESS(i,offsets[i1]+j1),s1,A.UNCHECKEDACCESS(i,offsets[i2]+j2),s2);
275 if(temp)
276 {
277 if(choices[i].first<choices[i].second)
278 return temp<0;
279 else
280 return temp>0;
281 }
282 }
283 int indices[4];
284 int F=choices[i].first;
285 int S=choices[i].second;
286 int toCheck;
287 if(i1==i)
288 if(i2==i)
289 toCheck=sort4uniquely(indices,F,S,j1,j2);
290 else
291 toCheck=sort3uniquely(indices,F,S,j1);
292 else
293 if(i2==i)
294 toCheck=sort3uniquely(indices,F,S,j2);
295 else
296 toCheck=sort2uniquely(indices,F,S);
297
298 for(int J=0;J<toCheck;J++)
299 {
300 int j=indices[J];
301 int temp=determinantSign2(getCoordinateOfInequality(i1,j1,i,j),s1,getCoordinateOfInequality(i2,j2,i,j),s2);
302 if(temp>0)
303 return true;
304 else if(temp<0)
305 return false;
306 }
307 }
308 return false;
309 }
getVolume()310 mvtyp getVolume()
311 {// Cannot overflow
312 return denominator;
313 }
replaceFirstOrSecond(bool first,int subconfigurationIndex,int newIndex,Vector<mvtyp> const & target)314 void replaceFirstOrSecond(bool first, int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target)__attribute__ ((always_inline))//updates the inequality table according to replacing first or second by newIndex in the subconfigurationIndex'th configuration
315 {// Cannot overflow
316 int newIndex2=newIndex;for(int i=0;i<subconfigurationIndex;i++)newIndex2+=tuple[i].getWidth();
317 int oldIndex=first?choices[subconfigurationIndex].first:choices[subconfigurationIndex].second;
318 (first?choices[subconfigurationIndex].first:choices[subconfigurationIndex].second)=newIndex;
319
320 mvtyp nextDenominator=(first?A[subconfigurationIndex][newIndex2].extend()+denominator.extend():-A[subconfigurationIndex][newIndex2].extend()).castToSingle();
321 mvtyp t=nextDenominator;
322 mvtypDivisor denominatorDivisor(denominator);
323 for(int c=0;c<k+Flags::computeDotProductInMatrix;c++)tempA.UNCHECKEDACCESS(c)=A.UNCHECKEDACCESS(c,newIndex2);
324
325 (-Abounds[subconfigurationIndex].extend()).castToSingle();//This conversion will fail if the following fails
326 for(int u=0;u<m;u++)
327 svec[u]=first?-A.UNCHECKEDACCESS(subconfigurationIndex,u):A.UNCHECKEDACCESS(subconfigurationIndex,u);
328 mvtyp svecBound=Abounds[subconfigurationIndex];
329
330 if(first)
331 for(int j=0;j<tuple[subconfigurationIndex].getWidth();j++)
332 svec[offsets[subconfigurationIndex]+j].subWithOverflowChecking(denominator);
333 svecBound.subWithOverflowChecking(denominator);
334
335
336 for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
337 if(first||a!=subconfigurationIndex)
338 Abounds.UNCHECKEDACCESS(a)=mvtyp::dotDivVector(&A.UNCHECKEDACCESS(a,0),&(svec[0]),t,tempA.UNCHECKEDACCESS(a),denominatorDivisor,m,Abounds[a],svecBound);
339 // for(int u=0;u<m;u++)
340 // *A.UNCHECKEDACCESSRESTRICT(a,u)=dotDiv(svec[u],tempA.UNCHECKEDACCESS(a),t,A.UNCHECKEDACCESS(a,u),denominatorDivisor);
341
342
343 {
344 for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
345 {
346 A[a][offsets[subconfigurationIndex]+oldIndex]=tempA[a].negatedWithOverflowChecking();
347 Abounds[a]=min(Abounds[a],negabs(tempA[a]));
348 }
349
350 if(!first)
351 {
352 A[subconfigurationIndex][offsets[subconfigurationIndex]+oldIndex]=denominator.negatedWithOverflowChecking();
353 Abounds[subconfigurationIndex].subWithOverflowChecking(denominator);
354 }
355 }
356 denominator=nextDenominator;
357
358 // We clear these unused entries of A to ensure that these columns are not chosen when comparing inequalities
359 for(int i=0;i<k+Flags::computeDotProductInMatrix;i++)
360 {
361 A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+choices[subconfigurationIndex].first)=0;
362 A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+choices[subconfigurationIndex].second)=0;
363 }
364 }
replaceFirst(int subconfigurationIndex,int newIndex,Vector<mvtyp> const & target)365 void replaceFirst(int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target){replaceFirstOrSecond(true,subconfigurationIndex,newIndex,target);}
replaceSecond(int subconfigurationIndex,int newIndex,Vector<mvtyp> const & target)366 void replaceSecond(int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target){replaceFirstOrSecond(false,subconfigurationIndex,newIndex,target);}
367
InequalityTable(std::vector<Matrix<mvtyp>> const & tuple_,int subconfigurationIndex_)368 InequalityTable(std::vector<Matrix<mvtyp> > const &tuple_, int subconfigurationIndex_):
369 tempA(tuple_.size()+Flags::computeDotProductInMatrix),
370 tuple(tuple_),
371 choices(tuple_.size()),
372 subconfigurationIndex(subconfigurationIndex_),
373 offsets(tuple_.size())
374 {// Cannot overflow
375 k=tuple.size();
376 m=0;
377 for(int i=0;i<tuple.size();i++)m+=tuple[i].getWidth();
378 svec.resize(m);
379 A=Matrix<mvtyp>(k+Flags::computeDotProductInMatrix,m);
380 {int offset=0;for(int i=0;i<tuple.size();i++){offsets[i]=offset;offset+=tuple[i].getWidth();}}
381 Abounds=Vector<mvtyp>(k+Flags::computeDotProductInMatrix);
382 }
setChoicesInitially()383 void setChoicesInitially()
384 {// Cannot overflow
385 //THIS WILL ONLY WORK FOR THE STARTING CONFIGURATION
386 //sets denominators,A and choices (these must have been initialized to the right sizes)
387 for(int i=0;i<k;i++)
388 choices[i]=std::pair<int, int> (i+0,i+1);
389 for(int i=0;i<m;i++)
390 for(int j=0;j<k;j++)
391 A[j][i]=0;
392 //we follow the lemma in the article. Why does the proof of the lemma split into 3 cases and not just 2?
393 int a=0;
394 for(int i=0;i<k;i++)
395 for(int gamma=0;gamma<tuple[i].getWidth();gamma++,a++)
396 if(gamma>i+1)
397 for(int ii=i;ii<gamma;ii++)
398 A[ii][a]=-1;
399 else
400 for(int ii=gamma;ii<i;ii++)
401 A[ii][a]=1;
402 denominator=1;
403 for(int i=0;i<k;i++)Abounds[i]=-1;
404 // assignDotProducts(target);
405 }
406 void compareInequalities(InequalityComparisonResult &result, Vector<mvtyp> const &target, int onlyK=-1)
407 {// Can only overflow if target entries are too big. Actually it seems that it is not this function that will overflow but dotVector.
408 bool empty=true;
409 int bestConfigurationIndex=-1;
410 int bestColumnIndex=-1;
411 mvtyp targetDotBest=0;
412
413 for(int i=0;i<k;i++)
414 {
415 // TO DO: exchange the line with the template parameter version:
416 // gfan::Matrix<mvtyp>::const_RowRef Ak=const_cast<const gfan::Matrix<mvtyp>&>(A)[k];
417 gfan::Matrix<CircuitTableInt32>::const_RowRef Ak=const_cast<const gfan::Matrix<mvtyp>&>(A)[k];
418 int offsetsi=offsets[i];
419 int tupleiwidth=tuple[i].getWidth();
420
421 if(onlyK!=-1)if(i!=onlyK)continue;
422 for(int j=0;j<tupleiwidth;j++)
423 if(Flags::computeDotProductInMatrix || isLegalIndex(i,j))//unused inequalities will have value 0. Therefore isLegalIndex(i,j) is not required if values are stored.
424 {
425 mvtyp ineqDotTarget=Flags::computeDotProductInMatrix?Ak.UNCHECKEDACCESS(offsetsi+j):dotVector(i,j,target,onlyK);
426 if(ineqDotTarget.isNegative())
427 {
428 if(!isReverseLexInvertedLessThanZero(i,j))
429 {
430 if(empty||compareReverseLexicographicInverted(bestConfigurationIndex,bestColumnIndex,i,j,ineqDotTarget,targetDotBest))
431 {
432 targetDotBest=ineqDotTarget;
433 empty=false;
434 bestConfigurationIndex=i;
435 bestColumnIndex=j;
436 }
437 }
438 // assert(!ineq.isZero());
439 }
440 }
441 }
442 result.empty=empty;
443 result.configurationIndex=bestConfigurationIndex;
444 result.columnIndex=bestColumnIndex;
445 // assert(denominator>0);
446 }
setChoicesFromEarlierHomotopy(InequalityTable const & parent,mvtyp degreeScaling,Vector<mvtyp> const & target)447 void setChoicesFromEarlierHomotopy(InequalityTable const &parent, mvtyp degreeScaling, Vector<mvtyp> const &target)
448 { // Overflows may happen, but only if the tuple entries are too big.
449 // Notice that the code below has overflow checks everywhere, except in the loop marked with a comment.
450
451 //sets denominators,A and choices (these must have been initialized to the right sizes
452 //columns have been added to configuration this->subconfigurationIndex
453 //for that reason we need to introduce new circuits. Old circuits are preserved.
454 //chioices are "relative" so no update is needed.
455
456 choices=parent.choices;
457 int numberToDrop=(subconfigurationIndex!=0) ? numberToDrop=k+1 : 0;
458
459 choices[subconfigurationIndex-1].first-=numberToDrop;
460 choices[subconfigurationIndex-1].second-=numberToDrop;
461
462 denominator=parent.denominator;
463 int offsetOld=0;
464 int offsetNew=0;
465 for(int i=0;i<k;i++)
466 {
467 int localNumberToDrop=0;
468 if(i==subconfigurationIndex-1)
469 localNumberToDrop=numberToDrop;
470 for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
471 if(a==subconfigurationIndex)
472 for(int j=0;j<parent.tuple[i].getWidth()-localNumberToDrop;j++)
473 A.UNCHECKEDACCESS(a,offsetNew+j)=parent.A.UNCHECKEDACCESS(a,offsetOld+j+localNumberToDrop);
474 else
475 {
476 // We can get an estimated 3 percent speed up by using known bounds to check if any of the multiplications below can overflow.
477 // Moreover, these multiplications can be moved outside the outer loop, to give better vectorisation.
478 for(int j=0;j<parent.tuple[i].getWidth()-localNumberToDrop;j++)
479 A.UNCHECKEDACCESS(a,offsetNew+j)=extendedMultiplication(degreeScaling,parent.A.UNCHECKEDACCESS(a,offsetOld+j+localNumberToDrop)).castToSingle();
480 }
481 if(i==subconfigurationIndex)
482 {
483 for(int j=parent.tuple[i].getWidth();j<tuple[i].getWidth();j++)
484 {
485 for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
486 A.UNCHECKEDACCESS(a,offsetNew+j)=0;
487 {
488 int b=choices[subconfigurationIndex].second-1;
489 if(b>=0)
490 A.UNCHECKEDACCESS(subconfigurationIndex,offsetNew+j).msubWithOverflowChecking(tuple[i].UNCHECKEDACCESS(b,j),denominator);
491 }
492 for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
493 { // If tuple entries are not bounded, overflow can happen in the following loop
494 mvtypDouble tempDouble=A.UNCHECKEDACCESS(a,offsetNew+j).extend();
495 for(int b=0;b<k;b++)
496 if(choices[subconfigurationIndex].second!=b+1 &&choices[subconfigurationIndex].first!=b+1)
497 tempDouble+=extendedMultiplication(tuple[i].UNCHECKEDACCESS(b,j),A.UNCHECKEDACCESS(a,offsetNew+b+1));
498 A.UNCHECKEDACCESS(a,offsetNew+j)=tempDouble.castToSingle();
499 }
500 mvtypDouble left=degreeScaling.extend();
501 for(int b=0;b<k;b++)
502 left-=mvtyp(tuple[i].UNCHECKEDACCESS(b,j)).extend();
503
504 mvtyp leftS=left.castToSingle();
505
506 if(choices[subconfigurationIndex].second==0)
507 A.UNCHECKEDACCESS(subconfigurationIndex,offsetNew+j).msubWithOverflowChecking(leftS,denominator);
508 else if(choices[subconfigurationIndex].first!=0)
509 for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
510 A.UNCHECKEDACCESS(a,offsetNew+j).maddWithOverflowChecking(leftS,mvtyp(A.UNCHECKEDACCESS(a,offsetNew)));
511 }
512 for(int j=0;j<parent.tuple[i].getWidth();j++)// <-- this loop has side effects on addExpressionOfCeb() above. Therefore it is placed after the loop above.
513 for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
514 A.UNCHECKEDACCESS(a,offsetNew+j)=extendedMultiplication(A.UNCHECKEDACCESS(a,offsetNew+j),degreeScaling).castToSingle();
515 }
516 offsetOld+=parent.tuple[i].getWidth();
517 offsetNew+=tuple[i].getWidth();
518 }
519 denominator=extendedMultiplication(denominator,degreeScaling).castToSingle();
520 if(Flags::computeDotProductInMatrix)assignDotProducts(target,subconfigurationIndex);
521 computeABounds();
522 for(int a=0;a<k;a++)
523 for(int i=0;i<k+Flags::computeDotProductInMatrix;i++)
524 {
525 A[i][offsets[a]+choices[a].first]=0;
526 A[i][offsets[a]+choices[a].second]=0;
527 }
528 }
529 };
530 struct StackItem{
531 int columnIndex;
532 int configurationIndex;
533 bool b;
534 int choice;
535 bool useFirstChanged,useSecondChanged;
StackItemStackItem536 StackItem(int columnIndex_, int configurationIndex_, bool b_, int choice_, bool useFirstChanged_, bool useSecondChanged_):
537 columnIndex(columnIndex_),
538 configurationIndex(configurationIndex_),
539 b(b_),
540 choice(choice_),
541 useFirstChanged(useFirstChanged_),
542 useSecondChanged(useSecondChanged_)
543 {
544 }
545 };
546 public:
547 std::vector<std::pair<int,int> > choices;
548 Vector<mvtyp> target;
549 bool useFirstChanged;
550 bool useSecondChanged;
551 std::vector<StackItem> stack;
552 int eliminatedK;
553 int eliminatedKOffset;
554 std::vector<Matrix<mvtyp> > tuple;
555 std::vector<int> offsets;
556 int m;
557 InequalityComparisonResult result;
558 InequalityTable inequalityTable;
constructInequalityTableFromParent(InequalityTable const & parentTable,mvtyp degreeScaling)559 void constructInequalityTableFromParent(InequalityTable const &parentTable, mvtyp degreeScaling)
560 {
561 inequalityTable.setChoicesFromEarlierHomotopy(parentTable, degreeScaling, target);
562 }
constructInequalityTableInitially(mvtyp degreeScaling)563 void constructInequalityTableInitially(mvtyp degreeScaling)
564 {
565 std::vector<Matrix<mvtyp> > tempTuple;for(int i=0;i<tuple.size();i++)tempTuple.push_back(simplex<mvtyp>(tuple.size(),1));
566 InequalityTable tempTable(tempTuple,-1);
567 tempTable.setChoicesInitially();
568 constructInequalityTableFromParent(tempTable,degreeScaling);
569 }
SingleTropicalHomotopyTraverser(std::vector<Matrix<mvtyp>> const & tuple_,int m_,std::vector<std::pair<int,int>> const & choices_,Vector<mvtyp> const & target_,int eliminatedK_)570 SingleTropicalHomotopyTraverser(std::vector<Matrix<mvtyp> > const &tuple_, int m_, std::vector<std::pair<int,int> > const &choices_, Vector<mvtyp> const &target_, int eliminatedK_):
571 choices(choices_),
572 target(target_),
573 eliminatedK(eliminatedK_),
574 tuple(tuple_),
575 m(m_),
576 inequalityTable(tuple,eliminatedK),
577 offsets(tuple_.size())
578 {
579 eliminatedKOffset=0;for(int i=0;i<eliminatedK;i++)eliminatedKOffset+=tuple_[i].getWidth();
580 {int offset=0;for(int i=0;i<tuple.size();i++){offsets[i]=offset;offset+=tuple[i].getWidth();}}
581
582 /* We also need to check that the input is within the required bounds
583 * This is important to not have overflows in:
584 * dotVector
585 * setChoicesFromEarlierHomotopy
586 * The number of summands is bounded by 3k+4*/
587
588 bool isOK=mvtyp::isSuitableSummationNumber(3*tuple.size()+4);
589 for(int i=0;i<target.size();i++)isOK&=mvtyp(target[i]).fitsInHalf();
590 for(int i=0;i<tuple.size();i++)
591 for(int j=0;j<tuple[i].getHeight();j++)
592 for(int k=0;k<tuple[i].getWidth();k++)
593 isOK&=mvtyp(tuple[i][j][k]).fitsInHalf();
594 if(!isOK)throw MVMachineIntegerOverflow;
595 }
process()596 virtual void process()
597 {
598 }
findOutgoingAndProcess(bool doProcess)599 bool findOutgoingAndProcess(bool doProcess)//sets up useFirstChanged and useSecondChanged and processes if process is true
600 {//returns true if we are at a leaf
601 // inequalityTable.checkABounds();
602 useFirstChanged=false;
603 useSecondChanged=false;
604 int onlyK=-1;
605 #if 1
606 if(eliminatedK!=-1)
607 if(target[choices[eliminatedK].first+eliminatedKOffset]==target[choices[eliminatedK].second+eliminatedKOffset])
608 onlyK=eliminatedK;
609 #endif
610
611 inequalityTable.compareInequalities(result,target,onlyK);
612
613 if(result.empty)
614 {
615 if(doProcess)process();
616 return true;
617 }
618
619 // reverse search rule:
620
621 mvtypDouble bestAtFirst=inequalityTable.getCoordinateOfInequality(result.configurationIndex,result.columnIndex,result.configurationIndex,choices[result.configurationIndex].first);
622 mvtypDouble bestAtSecond=inequalityTable.getCoordinateOfInequality(result.configurationIndex,result.columnIndex,result.configurationIndex,choices[result.configurationIndex].second);
623 if(bestAtFirst.isNegative())
624 {
625 if(bestAtSecond.isNegative())
626 {
627 useFirstChanged=true;
628 useSecondChanged=true;
629 }
630 else
631 {
632 if(bestAtSecond.isZero())
633 useFirstChanged=true;
634 else
635 if(choices[result.configurationIndex].second<result.columnIndex)
636 useFirstChanged=true;
637 }
638 }
639 else
640 {
641 if(bestAtSecond.isNegative())
642 {
643 if(bestAtFirst.isZero())
644 useSecondChanged=true;
645 else
646 if(choices[result.configurationIndex].first<result.columnIndex)
647 useSecondChanged=true;
648 }
649 else
650 {
651 assert(0);
652 }
653 }
654 return false;
655 }
goToFirstChild()656 void goToFirstChild()
657 {
658 // debug<<"First: configuration index:"<<data.configurationIndex<<" column index:"<<data.columnIndex<<"\n";
659 assert(useFirstChanged);
660 {
661 stack.push_back(StackItem(
662 result.columnIndex,
663 result.configurationIndex,
664 false,
665 choices[result.configurationIndex].first,
666 useFirstChanged,
667 useSecondChanged));
668 choices[result.configurationIndex].first=result.columnIndex;
669 inequalityTable.replaceFirst(result.configurationIndex,result.columnIndex,target);
670 }
671 }
goToSecondChild()672 void goToSecondChild()
673 {
674 // debug<<"Second: configuration index:"<<data.configurationIndex<<" column index:"<<data.columnIndex<<"\n";
675 assert(useSecondChanged);
676 {
677 stack.emplace_back(StackItem(
678 result.columnIndex,
679 result.configurationIndex,
680 true,
681 choices[result.configurationIndex].second,
682 useFirstChanged,
683 useSecondChanged));
684 choices[result.configurationIndex].second=result.columnIndex;
685 inequalityTable.replaceSecond(result.configurationIndex,result.columnIndex,target);
686 }
687 }
numberOfChildren()688 int numberOfChildren()
689 {
690 return int(useFirstChanged)+int(useSecondChanged);
691 }
goToNthChild(int n)692 void goToNthChild(int n)
693 {
694 if(n==0)
695 if(useFirstChanged)
696 goToFirstChild();
697 else
698 goToSecondChild();
699 else
700 goToSecondChild();
701 }
goBack()702 void goBack()
703 {
704 StackItem &B=stack.back();
705 result.columnIndex=B.columnIndex;
706 result.configurationIndex=B.configurationIndex;
707 if(B.b)
708 {
709 choices[result.configurationIndex].second=B.choice;
710 inequalityTable.replaceSecond(result.configurationIndex,B.choice,target);
711 }
712 else
713 {
714 choices[result.configurationIndex].first=B.choice;
715 inequalityTable.replaceFirst(result.configurationIndex,B.choice,target);
716 }
717 useFirstChanged=B.useFirstChanged;
718 useSecondChanged=B.useSecondChanged;
719 stack.pop_back();
720 }
atRoot()721 bool atRoot()
722 {
723 return stack.empty();
724 }
725 };
726
727
728 /*
729 * This class concatenates several homotopies to a single tree.
730 */
731 template<class mvtyp, class mvtypDouble, class mvtypDivisor>
732 class TropicalRegenerationTraverser{
733 // The following class is an attempt to separate homotopy data from the traversal logic.
734 class Data{
degree(Matrix<mvtyp> const & m)735 static mvtyp degree(Matrix<mvtyp> const &m)//assumes entries of m are non-negative
736 {
737 mvtyp ret=0;
738 for(int i=0;i<m.getWidth();i++)
739 {
740 mvtyp s(0);
741 for(int j=0;j<m.getHeight();j++)
742 s.addWithOverflowChecking(m[j][i]);
743 ret=max(s,ret);
744 }
745 return ret;
746 }
747 public:
748 std::vector<Vector<mvtyp> > targets;
749 std::vector<Matrix<mvtyp> > tuple;
750 std::vector<std::vector<Matrix<mvtyp> > > tuples;
751 Vector<mvtyp> degrees;
isFiniteIndex(int level,int index)752 bool isFiniteIndex(int level, int index)
753 {
754 return index>=tuple[0].getHeight()+1;
755 }
756
produceIthTuple(int i)757 std::vector<Matrix<mvtyp> > produceIthTuple(int i)
758 {
759 int n=tuple[0].getHeight();
760 std::vector<Matrix<mvtyp> > ret;
761 for(int j=0;j<tuple.size();j++)
762 {
763 if(j<i)ret.push_back(tuple[j]);
764 if(j==i)ret.push_back(combineLeftRight(simplex<mvtyp>(n,degree(tuple[j])),tuple[j]));
765 if(j>i)ret.push_back(simplex<mvtyp>(n,1));
766 }
767 return ret;
768 }
Data(std::vector<Matrix<mvtyp>> const & tuple_)769 Data(std::vector<Matrix<mvtyp> > const &tuple_):tuple(tuple_)
770 {
771 int n=tuple[0].getHeight();
772
773 {//adjusting to positive orthant
774 for(int i=0;i<tuple.size();i++)
775 for(int j=0;j<tuple[i].getHeight();j++)
776 {
777 mvtyp m;
778 if(tuple[i].getWidth()==0)
779 m=0;
780 else
781 m=tuple[i][j][0];
782 for(int k=0;k<tuple[i].getWidth();k++)m=min(m,tuple[i][j][k]);
783 for(int k=0;k<tuple[i].getWidth();k++)tuple[i][j][k].subWithOverflowChecking(m);
784 }
785 }
786
787 for(int i=0;i<tuple.size();i++)
788 degrees.push_back(degree(tuple[i]));
789
790 for(int i=0;i<tuple.size();i++)
791 tuples.push_back(produceIthTuple(i));
792
793 for(int i=0;i<tuple.size();i++)
794 {
795 Vector<mvtyp> targ;
796 for(int j=0;j<tuple.size();j++)
797 {
798 if(j==i)
799 targ=concatenation(targ,concatenation(Vector<mvtyp>::allOnes(n+1),Vector<mvtyp>(tuple[i].getWidth())));
800 else
801 targ=concatenation(targ,Vector<mvtyp>(tuples[i][j].getWidth()));
802 }
803 targets.push_back(targ);
804 }
805 };
806
firstIntersection()807 std::vector<std::pair<int,int> > firstIntersection()
808 {
809 std::vector<std::pair<int,int> > ret;
810 for(int i=0;i<tuple.size();i++)
811 ret.push_back(std::pair<int,int>(i+0,i+1));
812 return ret;
813 }
814
castToNextLevel(std::vector<std::pair<int,int>> const & choices,int i,int S,std::vector<std::pair<int,int>> & ret)815 void castToNextLevel(std::vector<std::pair<int,int> > const &choices, int i, int S, std::vector<std::pair<int,int> > &ret)
816 {
817 assert(ret.size()==choices.size());
818 for(int j=0;j<choices.size();j++)
819 ret[j]=choices[j];
820
821 assert(ret[i].first>=S);
822 assert(ret[i].second>=S);
823 ret[i].first-=S;
824 ret[i].second-=S;
825 }
826 };
cayleyConfigurationWidth(std::vector<Matrix<mvtyp>> const & tuple)827 static int cayleyConfigurationWidth(std::vector<Matrix<mvtyp> > const &tuple)
828 {
829 int m2=0;
830 for(int i=0;i<tuple.size();i++)
831 m2+=tuple[i].getWidth();
832 return m2;
833 }
834 public:
835 int depth;
836 int counter;
837 typedef SingleTropicalHomotopyTraverser<mvtyp,mvtypDouble,mvtypDivisor> MySingleTropicalHomotopyTraverser;
838 std::vector<MySingleTropicalHomotopyTraverser> traversers;
839 Data fullData;
840 int level;
841 bool deadEnd;
842 bool isLevelLeaf;
843 bool isSolutionVertex;
844 std::vector<bool> isLevelLeafStack;
TropicalRegenerationTraverser(std::vector<Matrix<mvtyp>> const & tuple_)845 TropicalRegenerationTraverser(std::vector<Matrix<mvtyp> > const &tuple_):
846 fullData(tuple_),counter(0),depth(0)
847 {
848 assert(tuple_.size());
849 for(int i=0;i<tuple_.size();i++)
850 traversers.push_back(MySingleTropicalHomotopyTraverser(fullData.tuples[i],cayleyConfigurationWidth(fullData.tuples[i]),fullData.firstIntersection(),fullData.targets[i],i));
851 traversers[0].constructInequalityTableInitially(fullData.degrees[0]);
852 level=0;
853 }
process()854 virtual void process()
855 {
856 }
findOutgoingAndProcess(bool doProcess)857 bool findOutgoingAndProcess(bool doProcess)
858 {
859 isSolutionVertex=false;
860 deadEnd=false;
861 isLevelLeaf=traversers[level].findOutgoingAndProcess(false);
862 if(isLevelLeaf)
863 {//leaf
864 bool isFinite=fullData.isFiniteIndex(level,traversers[level].choices[level].first)&&fullData.isFiniteIndex(level,traversers[level].choices[level].second);
865 deadEnd=!isFinite;
866 if(isFinite && (level==fullData.tuple.size()-1))
867 {
868 isSolutionVertex=true;
869 if(doProcess){process();}
870 return true;
871 }
872 }
873 return false;
874 }
numberOfChildren()875 int numberOfChildren()
876 {
877 if(isLevelLeaf&&(level==fullData.tuple.size()-1))return 0;
878 if(!isLevelLeaf)
879 return traversers[level].numberOfChildren();
880 else
881 return 1-deadEnd;
882 }
goToNthChild(int n)883 void goToNthChild(int n)
884 {
885 depth++;
886 isLevelLeafStack.push_back(isLevelLeaf);
887 if(!isLevelLeaf)
888 traversers[level].goToNthChild(n);
889 else
890 {
891 fullData.castToNextLevel(traversers[level].choices,level,fullData.tuples[level][level].getWidth()-fullData.tuples[level+1][level].getWidth(),traversers[level+1].choices);
892 traversers[level+1].constructInequalityTableFromParent(traversers[level].inequalityTable,fullData.degrees[level+1]);
893 level++;
894 }
895 }
print()896 void print()
897 {
898 }
goBack()899 void goBack()
900 {
901 depth--;
902 counter++;
903 deadEnd=false;
904 if(traversers[level].atRoot())
905 level--;
906 else
907 traversers[level].goBack();
908 isLevelLeaf=isLevelLeafStack.back();
909 isLevelLeafStack.pop_back();
910 }
911 };
912
913 /*
914 * This class glues Bjarne Knudsen's parallel traverser to our MultiLevelIntersectionTraverser.
915 * This class should eventually be written so that it works on any homotopy.
916 *
917 * Since we are not sure about how exceptions would be handled by the parallel traverser,
918 * we write ugly aborting code, which may actually work.
919 */
920 template<class mvtyp, class mvtypDouble, class mvtypDivisor>
921 class SpecializedRTraverser: public Traverser
922 {
923 public:
924 typedef TropicalRegenerationTraverser<mvtyp,mvtypDouble,mvtypDivisor> MyTropicalRegenerationTraverser;
925 MyTropicalRegenerationTraverser T;
926 mvtypDouble mixedVolume;
927 int numberOfExpensiveSteps;
SpecializedRTraverser(std::vector<Matrix<mvtyp>> const & tuple_)928 SpecializedRTraverser(std::vector<Matrix<mvtyp> > const &tuple_):
929 mixedVolume(),
930 numberOfExpensiveSteps(0),
931 T(tuple_) //Constructor my throw if input is too big.
932 {
933 numberOfExpensiveSteps++;
934 T.findOutgoingAndProcess(false);
935 }
getEdgeCountNext(void)936 int getEdgeCountNext( void )
937 {
938 if(!aborting)
939 {
940 try{
941 return T.numberOfChildren();
942 }
943 catch(...){abort();}
944 }
945 return 0;
946 }
947
moveToNext(int index,bool collect_info)948 int moveToNext( int index,
949 bool collect_info )
950 {
951 if(!aborting)
952 {
953 try{
954 T.goToNthChild(index);
955 numberOfExpensiveSteps++;
956 T.findOutgoingAndProcess(false);
957 }
958 catch(...){abort();}
959 }
960 return 0;
961 }
962
moveToPrev(int index)963 void moveToPrev( int index )
964 {
965 if(!aborting)
966 {
967 try{
968 T.goBack(); //index ignored
969 }
970 catch(...){abort();}
971 }
972 }
973
collectInfo(void)974 void collectInfo( void )
975 {
976 if(!aborting)
977 {
978 try{
979 if(T.isSolutionVertex)
980 mixedVolume.addWithOverflowCheck(T.traversers[T.level].inequalityTable.getVolume().extend());
981 }
982 catch(...){abort();}
983 }
984 }
985
printState(void)986 void printState( void )
987 {
988 T.print();
989 }
990 };
991 }
992
993 #endif /* GFANLIB_TROPICALHOMOTOPY_H_ */
994