Context Navigation

source: git/gfanlib/gfanlib_tropicalhomotopy.h @ 7cc3fd

spielwiese

Last change on this file since 7cc3fd was 15813d, checked in by Hans Schoenemann <hannes@…>, 8 years ago
format
Property mode set to `100644`
File size: 56.3 KB

Line
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>
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;
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	}
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:
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	}
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	}
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	}
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	}
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	}
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	}
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	}
310	mvtyp getVolume()
311	{// Cannot overflow
312	return denominator;
313	}
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	}
365	void replaceFirst(int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target){replaceFirstOrSecond(true,subconfigurationIndex,newIndex,target);}
366	void replaceSecond(int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target){replaceFirstOrSecond(false,subconfigurationIndex,newIndex,target);}
367
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	}
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	}
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;
536	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;
559	void constructInequalityTableFromParent(InequalityTable const &parentTable, mvtyp degreeScaling)
560	{
561	inequalityTable.setChoicesFromEarlierHomotopy(parentTable, degreeScaling, target);
562	}
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	}
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	}
596	virtual void process()
597	{
598	}
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	}
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	}
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	}
688	int numberOfChildren()
689	{
690	return int(useFirstChanged)+int(useSecondChanged);
691	}
692	void goToNthChild(int n)
693	{
694	if(n==0)
695	if(useFirstChanged)
696	goToFirstChild();
697	else
698	goToSecondChild();
699	else
700	goToSecondChild();
701	}
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	}
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{
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;
752	bool isFiniteIndex(int level, int index)
753	{
754	return index>=tuple[0].getHeight()+1;
755	}
756
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	}
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
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
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	};
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;
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	}
854	virtual void process()
855	{
856	}
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	}
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	}
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	}
896	void print()
897	{
898	}
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;
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	}
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
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
963	void moveToPrev( int index )
964	{
965	if(!aborting)
966	{
967	try{
968	T.goBack(); //index ignored
969	}
970	catch(...){abort();}
971	}
972	}
973
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
986	void printState( void )
987	{
988	T.print();
989	}
990	};
991	}
992
993	#endif /* GFANLIB_TROPICALHOMOTOPY_H_ */

Note: See TracBrowser for help on using the repository browser.

Download in other formats: