← Back to branch summary

~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

« back to all changes in this revision

Viewing changes to tropical2.cpp

  • Committer: Package Import Robot
  • Author(s): Cédric Boutillier
  • Date: 2013-07-09 10:44:01 UTC
  • mfrom: (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130709104401-5q66ozz5j5af0dak
Tags: 0.5+dfsg-3
* Upload to unstable.
* modify remove_failing_tests_on_32bits.patch to replace command of
  0009RenderStairCase test with an empty one instead of deleting it.
* remove lintian override about spelling error

Show diffs side-by-side

added added

removed removed

Lines of Context:
tropical2.cpp
1
1
#include "tropical2.h"
2
2
 
3
3
#include <stdlib.h>
 
4
#include <iostream>
 
5
 
4
6
#include "buchberger.h"
5
7
#include "division.h"
6
8
#include "tropical.h"
7
9
#include "wallideal.h"
8
10
#include "dimension.h"
9
11
#include "halfopencone.h"
 
12
#include "breadthfirstsearch.h"
10
13
 
11
14
#include "timer.h"
12
15
#include "log.h"
13
16
 
14
17
static Timer tropicalPrincipalIntersectionTimer("Tropical principal intersection",1);
15
18
 
 
19
static void startingConeError()
 
20
{
 
21
  fprintf(Stderr,"UNABLE TO COMPUTE STARTING CONE.\n");
 
22
  fprintf(Stderr,"THE STARTING CONE ALGORITHM IN GFAN IS BASED ON HEURISTICS WHICH HAVE FAILED ON THIS EXAMPLE.\n");
 
23
  assert(0);
 
24
}
 
25
 
16
26
PolynomialSet initialIdeal(PolynomialSet const &g, IntegerVector const &weight)
17
27
//Assume homogeneous
18
28
{
90
100
{
91
101
  PolynomialRing theRing=groebnerBasis.getRing();
92
102
  PolynomialSet r(theRing);
 
103
  if(theRing.getNumberOfVariables()!=weight.size())
 
104
    {
 
105
      cerr << "Error: Number of varaibles in polynomial ring "<<theRing.getNumberOfVariables()<< " length of weight vector " << weight.size() <<endl;
 
106
      assert(0);
 
107
    }
93
108
 
94
109
  for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
95
110
    {
130
145
      PolynomialSet h2=groebnerBasis;
131
146
      buchberger(&h2,t);
132
147
      log2 fprintf(Stderr,"Done computing Gr\"obner basis.\n");
133
 
      
 
148
 
134
149
      log3 AsciiPrinter(Stderr).printPolynomialSet(h2);
135
150
      PolynomialSet wall=initialFormsAssumeMarked(h2,*i);
136
151
 
137
152
      log3 AsciiPrinter(Stderr).printString("Initial ideal:\n");
138
153
      log3 AsciiPrinter(Stderr).printPolynomialSet(wall);
139
 
      
 
154
 
140
155
      int hdim2=dimensionOfHomogeneitySpace(wall);
141
156
      if(hdim2>h)
142
157
        {
145
160
              log1 fprintf(Stderr,"Iterating recursively.\n");
146
161
              //PolynomialSet initialIdeal=guessInitialIdealWithoutMonomial(wall,0);
147
162
              PolynomialSet initialIdeal=guessInitialIdealWithoutMonomial(wall,fullNeighbourBasis,onlyCheckRays);
148
 
              
 
163
 
149
164
              if(fullNeighbourBasis)
150
165
                {
151
166
                  //*fullNeighbourBasis=liftBasis(initialIdeal,h2);
152
167
                  *fullNeighbourBasis=liftBasis(*fullNeighbourBasis,h2);
153
168
                }
154
 
                  
155
 
              
 
169
 
 
170
 
156
171
              result=true;
157
172
              return initialIdeal;
158
173
            }
165
180
PolynomialSet guessInitialIdealWithoutMonomial(PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, bool onlyCheckRays) //ideal must be homogeneous
166
181
  // fullNeighbourBasis is set to a Groebner basis of the full ideal. The returned basis and fullNeighbourBasis have at least one termorder in common
167
182
{
 
183
  //  log0 fprintf(Stderr,"A\n");
168
184
  assert(groebnerBasis.isValid());
 
185
  //  log0 fprintf(Stderr,"B\n");
169
186
  if(fullNeighbourBasis)
170
187
    {
171
188
      assert(fullNeighbourBasis->isValid());
172
189
    }
 
190
  //  log0 fprintf(Stderr,"C\n");
173
191
 
174
192
  int n=groebnerBasis.numberOfVariablesInRing();
 
193
  //  log0 fprintf(Stderr,"D\n");
175
194
  int h=dimensionOfHomogeneitySpace(groebnerBasis);
 
195
  //  log0 fprintf(Stderr,"E\n");
176
196
  int d=krullDimension(groebnerBasis);
 
197
  //  log0 fprintf(Stderr,"F\n");
177
198
 
178
199
  if(d==h)
179
200
    {
182
203
    }
183
204
 
184
205
  {
185
 
    IntegerVectorList a;
186
 
    PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
 
206
    log2 fprintf(Stderr,"Computing extreme rays.\n");
 
207
    //IntegerVectorList a;
 
208
    PolyhedralCone p=coneFromMarkedBasis(groebnerBasis);
 
209
    //PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
187
210
    IntegerVectorList extreme=p.extremeRays();
188
211
    log2 fprintf(Stderr,"Extreme rays of Groebner cone:\n");
189
212
    log2 AsciiPrinter(Stderr).printVectorList(extreme);
192
215
    PolynomialSet r=checkList(extreme,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
193
216
    if(result)return r;
194
217
  }
195
 
  if(onlyCheckRays)
196
 
    {
197
 
      fprintf(Stderr,"UNABLE TO COMPUTE STARTING CONE.\n");
198
 
      assert(0); // We could not guess an intial ideal
199
 
    }
200
 
 
 
218
  if(onlyCheckRays)startingConeError();
201
219
  PolyhedralFan f=PolyhedralFan::fullSpace(n);
202
220
  /*  for(int i=0;i<d-1;i++)
203
221
    {
219
237
      fprintf(Stderr,"Max dimension: %i\n",f.getMaxDimension());
220
238
      f=refinement(f,PolyhedralFan::bergmanOfPrincipalIdeal(*i));
221
239
      f.removeAllExcept(3);
222
 
      
 
240
 
223
241
      IntegerVectorList l=f.getRelativeInteriorPoints();
224
 
      
 
242
 
225
243
      bool result;
226
244
      PolynomialSet r=checkList(l,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
227
245
      if(result)return r;
228
246
    }
229
 
  fprintf(Stderr,"UNABLE TO COMPUTE STARTING CONE.\n");
230
 
  assert(0); // We could not guess an intial ideal
 
247
  startingConeError();
 
248
  return groebnerBasis;
 
249
}
 
250
 
 
251
static PolynomialSet checkListStably(IntegerVectorList const &l, PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, int h, bool &result, bool onlyCheckRays)
 
252
{
 
253
  debug<< "Checklist called on"<<groebnerBasis;
 
254
  for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
 
255
    {
 
256
      WeightReverseLexicographicTermOrder t(*i);
 
257
      log2 fprintf(Stderr,"Taking initial forms with respect to:");
 
258
      log2 AsciiPrinter(Stderr).printVector(*i);
 
259
      log2 fprintf(Stderr,"\n");
 
260
      PolynomialSet h2=groebnerBasis;
 
261
      log2 fprintf(Stderr,"Done computing Gr\"obner basis.\n");
 
262
 
 
263
      log3 AsciiPrinter(Stderr).printPolynomialSet(h2);
 
264
      PolynomialSet wall=initialForms(h2,*i);
 
265
 
 
266
      log3 AsciiPrinter(Stderr).printString("Initial ideal:\n");
 
267
      log3 AsciiPrinter(Stderr).printPolynomialSet(wall);
 
268
 
 
269
      int hdim2=dimensionOfHomogeneitySpace(wall);
 
270
      if(hdim2>h)
 
271
        {
 
272
          if(nonEmptyStableIntersection(wall))
 
273
            {
 
274
              log1 fprintf(Stderr,"Iterating recursively.\n");
 
275
              //PolynomialSet initialIdeal=guessInitialIdealWithoutMonomial(wall,0);
 
276
              PolynomialSet initialIdeal=guessInitialIdealWithoutMonomialStably(wall,fullNeighbourBasis,onlyCheckRays);
 
277
 
 
278
              if(fullNeighbourBasis)
 
279
                {
 
280
                  //*fullNeighbourBasis=liftBasis(initialIdeal,h2);
 
281
//                  *fullNeighbourBasis=liftBasis(*fullNeighbourBasis,h2);
 
282
                  *fullNeighbourBasis=groebnerBasis;
 
283
                  fullNeighbourBasis->copyMarkings(initialIdeal);
 
284
                }
 
285
 
 
286
 
 
287
              result=true;
 
288
              return initialIdeal;
 
289
            }
 
290
        }
 
291
    }
 
292
  result=false;
 
293
  return groebnerBasis;
 
294
}
 
295
 
 
296
PolynomialSet guessInitialIdealWithoutMonomialStably(PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, bool onlyCheckRays) //ideal must be homogeneous
 
297
  // fullNeighbourBasis is set to a Groebner basis of the full ideal. The returned basis and fullNeighbourBasis have at least one termorder in common
 
298
{
 
299
  int n=groebnerBasis.numberOfVariablesInRing();
 
300
  int h=dimensionOfHomogeneitySpace(groebnerBasis);
 
301
  int d=n-groebnerBasis.size();//krullDimension(groebnerBasis);
 
302
 
 
303
  debug<</*"d"<<d<<*/"h"<<h<<"n"<<n<<"\n";
 
304
 
 
305
 
 
306
  if(d==h)
 
307
    {
 
308
      if(fullNeighbourBasis)*fullNeighbourBasis=groebnerBasis;
 
309
      return groebnerBasis;
 
310
    }
 
311
 
 
312
  {
 
313
    log2 fprintf(Stderr,"Computing extreme rays.\n");
 
314
    //IntegerVectorList a;
 
315
    PolyhedralCone p=coneFromMarkedBasis(groebnerBasis);
 
316
    //PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
 
317
    IntegerVectorList extreme=p.extremeRays();
 
318
    log2 fprintf(Stderr,"Extreme rays of Groebner cone:\n");
 
319
    log2 AsciiPrinter(Stderr).printVectorList(extreme);
 
320
 
 
321
    bool result;
 
322
    PolynomialSet r=checkListStably(extreme,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
 
323
    if(result)return r;
 
324
  }
 
325
  if(onlyCheckRays)startingConeError();
 
326
 
 
327
  PolyhedralFan f=PolyhedralFan::fullSpace(n);
 
328
 
 
329
  int hypersurfacesToGo=groebnerBasis.size();
 
330
  for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
 
331
    {
 
332
      fprintf(Stderr,"Hypersurfaces to go:%i\n",hypersurfacesToGo--);
 
333
      fprintf(Stderr,"Max dimension: %i\n",f.getMaxDimension());
 
334
      f=refinement(f,PolyhedralFan::bergmanOfPrincipalIdeal(*i));
 
335
      f.removeAllExcept(3);
 
336
 
 
337
      IntegerVectorList l=f.getRelativeInteriorPoints();
 
338
 
 
339
      bool result;
 
340
      PolynomialSet r=checkListStably(l,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
 
341
      if(result)return r;
 
342
    }
 
343
  startingConeError();
231
344
  return groebnerBasis;
232
345
}
233
346