File: tropicalbasis.cpp

package info (click to toggle)
gfan 0.3dfsg-1.1
  • links: PTS
  • area: main
  • in suites: wheezy
  • size: 2,016 kB
  • sloc: cpp: 17,728; makefile: 251
file content (190 lines) | stat: -rw-r--r-- 6,377 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
#include "tropicalbasis.h"
#include "buchberger.h"

#include "tropical.h"
#include "tropical2.h"
#include "division.h"
#include "wallideal.h"
#include "log.h"

#include "timer.h"

static Timer iterativeTropicalBasisTimer("Iterative tropical basis",1);

typedef set<IntegerVector> IntegerVectorSet;

PolynomialSet tropicalBasisOfCurve(int n, PolynomialSet g, PolyhedralFan *intersectionFan, int linealitySpaceDimension) //Assuming g is homogeneous
{
  int homog=linealitySpaceDimension;
  assert(homog>0 || n==0);
  TimerScope ts(&iterativeTropicalBasisTimer);
  PolyhedralFan f(n);
  if(!intersectionFan)intersectionFan=&f;
     
  *intersectionFan=tropicalPrincipalIntersection(n,g,linealitySpaceDimension);

  IntegerVectorSet containsNoMonomialCache;

  while(1)
    {
      PolyhedralFan::coneIterator i;
      for(i=intersectionFan->conesBegin();i!=intersectionFan->conesEnd();i++)
	{
	  IntegerVector relativeInteriorPoint=i->getRelativeInteriorPoint();
	  if(containsNoMonomialCache.count(relativeInteriorPoint)>0)
	    {
	      log2 fprintf(Stderr,"Weight vector found in cache.... contains no monomial.\n");
	    }
	  else
	    {
	      WeightReverseLexicographicTermOrder t(relativeInteriorPoint);
	      log2 fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
	      log2 AsciiPrinter(Stderr).printVector(relativeInteriorPoint);
	      log2 fprintf(Stderr,"\n");
	      PolynomialSet h2=g;
	      buchberger(&h2,t);
	      log2 fprintf(Stderr,"Done computing Gr\"obner basis.\n");
	      
	      log3 AsciiPrinter(Stderr).printPolynomialSet(h2);

	      PolynomialSet wall=initialFormsAssumeMarked(h2,relativeInteriorPoint);
	      
	      if(containsMonomial(wall))
		{
		  log2 fprintf(Stderr,"Initial ideal contains a monomial.\n");
		  Polynomial m(computeTermInIdeal(wall));
		  log2 fprintf(Stderr,"Done computing term in ideal\n");
		  
		  Polynomial temp=m-division(m,h2,LexicographicTermOrder());
		  g.push_back(temp);
		  
		  log2 fprintf(Stderr,"Adding element to basis:\n");
		  log2 AsciiPrinter(Stderr).printPolynomial(temp);
		  log2 fprintf(Stderr,"\n");
		  
		  *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
		  break;
		}
	      else
		{
		  if(i->dimension()<=1+homog)
		    //if(!containsMonomial(wall) && i->dimension()<=1+homog)//line for testing perturbation code
		    {
		      log2 fprintf(Stderr,"Initial ideal contains no monomial... caching weight vector.\n");
		      containsNoMonomialCache.insert(relativeInteriorPoint);
		    }
		  else
		    {
		      /* We need to compute the initial ideal with
			 respect to "relativeInteriorPoint" perturbed
			 with a basis of the span of the cone. Instead
			 of perturbing we may as well compute initial
			 ideal successively. We have already computed
			 the initial ideal with respect to
			 "relativeInteriorPoint". To get the perturbed
			 initial ideal we take initial ideal with
			 respect to each vector in the basis of the
			 span.*/
		      IntegerVectorList empty;
		      PolyhedralCone dual=PolyhedralCone(empty,i->getEquations(),i->ambientDimension()).dualCone();
		      dual.canonicalize();
		      IntegerVectorList basis=dual.getEquations();
		      PolynomialSet witnessLiftBasis=h2;//basis with respect to relativeInteriorPoint
		      for(IntegerVectorList::const_iterator j=basis.begin();j!=basis.end();j++)
			{
			  WeightReverseLexicographicTermOrder t(*j);
			  PolynomialSet h3=wall;
			  buchberger(&h3,t);
			  wall=initialFormsAssumeMarked(h3,*j);
			  witnessLiftBasis=liftBasis(h3,witnessLiftBasis);
			}
		      if(containsMonomial(wall))
			{
			  Polynomial m(computeTermInIdeal(wall));
			  Polynomial temp=m-division(m,witnessLiftBasis,LexicographicTermOrder());
			  g.push_back(temp);
			  *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
			  break;
			}
		      else
			{
			  assert(0);
			}
		    }
		}
	    }
	}
      if(i==intersectionFan->conesEnd())break;
    }
  return g;
}

/*
PolynomialSet iterativeTropicalBasisNoPerturbation(int n, PolynomialSet g, PolyhedralFan *intersectionFan, int linealitySpaceDimension, bool doPrint) //Assuming g is homogeneous
{
  TimerScope ts(&iterativeTropicalBasisTimer);
  PolyhedralFan f(n);
  if(!intersectionFan)intersectionFan=&f;
     
  *intersectionFan=tropicalPrincipalIntersection(n,g,linealitySpaceDimension);

  IntegerVectorSet containsNoMonomialCache;

  while(1)
    {
      //      AsciiPrinter(Stderr).printPolyhedralFan(*intersectionFan);
      //      assert(f.getMaxDimension()==1);

      IntegerVectorList l=intersectionFan->getRelativeInteriorPoints();

      IntegerVectorList::const_iterator i;
      for(i=l.begin();i!=l.end();i++)
	{
	  if(containsNoMonomialCache.count(*i)>0)
	    {
	      if(doPrint)fprintf(Stderr,"Weight vector found in cache.... contains no monomial.\n");
	    }
	  else
	    {
	      WeightReverseLexicographicTermOrder t(*i);
	      if(doPrint)fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
	      if(doPrint)AsciiPrinter(Stderr).printVector(*i);
	      if(doPrint)fprintf(Stderr,"\n");
	      PolynomialSet h2=g;
	      buchberger(&h2,t);
	      if(doPrint)fprintf(Stderr,"Done computing Gr\"obner basis.\n");
	      
	      //	      AsciiPrinter(Stderr).printPolynomialSet(h2);
	      PolynomialSet wall=initialFormsAssumeMarked(h2,*i);
	      //fprintf(Stderr,"Wall ideal:\n");
	      //AsciiPrinter(Stderr).printPolynomialSet(wall);
	      
	      if(containsMonomial(wall))
		{
		  if(doPrint)fprintf(Stderr,"Initial ideal contains a monomial.\n");
		  Polynomial m(computeTermInIdeal(wall));
		  if(doPrint)fprintf(Stderr,"Done computing term in ideal\n");
		  
		  Polynomial temp=m-division(m,h2,LexicographicTermOrder());
		  g.push_back(temp);
		  
		  if(doPrint)fprintf(Stderr,"Adding element to basis:\n");
		  if(doPrint)AsciiPrinter(Stderr).printPolynomial(temp);
		  if(doPrint)fprintf(Stderr,"\n");
		  
		  *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
		  break;
		}
	      else
		{
		  if(doPrint)fprintf(Stderr,"Initial ideal contains no monomial... caching weight vector.\n");
		  containsNoMonomialCache.insert(*i);
		}
	    }
	}
      if(i==l.end())break;
    }
  return g;
}

*/