← Back to branch summary

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

« back to all changes in this revision

Viewing changes to app_groebnercone.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:
app_groebnercone.cpp
7
7
#include "gfanapplication.h"
8
8
#include "polyhedralfan.h"
9
9
#include "halfopencone.h"
 
10
#include "linalg.h"
 
11
#include "log.h"
10
12
 
11
13
class GroebnerConeApplication : public GFanApplication
12
14
{
13
15
  SimpleOption optionRestrict;
14
16
  SimpleOption optionPair;
15
17
  SimpleOption optionAsFan;
 
18
  SimpleOption optionVectorInput;
16
19
public:
17
20
  const char *helpText()
18
21
  {
21
24
  GroebnerConeApplication():
22
25
    optionRestrict("--restrict","Add an inequality for each coordinate, so that the the cone is restricted to the non-negative orthant."),
23
26
    optionPair("--pair","The Groebner cone is given by a pair of compatible Groebner bases. The first basis is for the initial ideal and the second for the ideal itself. See the tropical section of the manual."),
24
 
    optionAsFan("--asfan","Writes the cone as a polyhedral fan with all its faces instead. In this way the extreme rays of the cone are also computed.")
 
27
    optionAsFan("--asfan","Writes the cone as a polyhedral fan with all its faces instead. In this way the extreme rays of the cone are also computed."),
 
28
    optionVectorInput("--vectorinput","Compute a cone given list of inequalities rather than a Groebner cone. The input is an integer which specifies the dimension of the ambient space, a list of inequalities given as vectors and a list of equations.")
25
29
  {
26
30
    registerOptions();
27
31
  }
28
32
 
29
 
  char *name()
 
33
  const char *name()
30
34
  {
31
35
    return "_groebnercone";
32
36
  }
33
37
 
34
38
  int main()
35
39
  {
36
 
    PolynomialSet m=FileParser(Stdin).parsePolynomialSetWithRing();
37
 
    PolynomialSet g(m.getRing());
38
 
    if(optionPair.getValue())
 
40
    IntegerVectorList equalities;
 
41
    IntegerVectorList normals;
 
42
    int n;
 
43
 
 
44
    if(optionVectorInput.getValue())
39
45
      {
40
 
        g=FileParser(Stdin).parsePolynomialSet(m.getRing());
 
46
        FileParser P(Stdin);
 
47
        n=P.parseInt();
 
48
        normals=P.parseIntegerVectorList();
 
49
        equalities=P.parseIntegerVectorList();
41
50
      }
42
51
    else
43
 
      {
44
 
        g=m;
45
 
        m=g.markedTermIdeal();
46
 
      }
47
 
    int n=g.getRing().getNumberOfVariables();
48
 
        
49
 
    IntegerVectorList equalities=wallInequalities(m);
50
 
    IntegerVectorList normals=algebraicTest(wallInequalities(g),g);
51
 
    if(optionRestrict.getValue())
52
 
      {
53
 
        for(int i=0;i<n;i++)
54
 
          normals.push_back(IntegerVector::standardVector(n,i));
55
 
      }
56
 
 
57
 
    // AsciiPrinter(Stderr).printVectorList(normals);
 
52
    {
 
53
      PolynomialSet m=FileParser(Stdin).parsePolynomialSetWithRing();
 
54
      PolynomialSet g(m.getRing());
 
55
      if(optionPair.getValue())
 
56
      {
 
57
          g=FileParser(Stdin).parsePolynomialSet(m.getRing());
 
58
      }
 
59
      else
 
60
      {
 
61
          g=m;
 
62
          m=g.markedTermIdeal();
 
63
      }
 
64
 
 
65
      {
 
66
          //Check that markings are consistent
 
67
          PolynomialSet M1=g.markedTermIdeal();
 
68
          PolynomialSet M2=m.markedTermIdeal();
 
69
          assert(M1.size()==M2.size());
 
70
          PolynomialSet::const_iterator i1=M1.begin();
 
71
          for(PolynomialSet::const_iterator i2=M2.begin();i2!=M2.end();i1++,i2++)
 
72
          {
 
73
                  assert((*i1-*i2).isZero());
 
74
          }
 
75
      }
 
76
 
 
77
      n=g.getRing().getNumberOfVariables();
 
78
 
 
79
      equalities=wallInequalities(m);
 
80
      normals=(wallInequalities(g));
 
81
      if(optionRestrict.getValue())
 
82
      {
 
83
          for(int i=0;i<n;i++)
 
84
                  normals.push_back(IntegerVector::standardVector(n,i));
 
85
      }
 
86
    }
 
87
 
 
88
 
 
89
 
 
90
    {
 
91
      FieldMatrix A=integerMatrixToFieldMatrix(rowsToIntegerMatrix(equalities,n),Q);
 
92
      A.reduce();
 
93
      A.removeZeroRows();
 
94
      equalities=IntegerVectorList();
 
95
      for(int i=0;i<A.getHeight();i++)
 
96
        equalities.push_back(A[i].primitive());
 
97
 
 
98
      set<IntegerVector> newInequalities;
 
99
      for(IntegerVectorList::const_iterator i=normals.begin();i!=normals.end();i++)
 
100
        newInequalities.insert(A.canonicalize(integerVectorToFieldVector(*i,Q)).primitive());
 
101
 
 
102
      normals=IntegerVectorList();
 
103
      for(set<IntegerVector>::const_iterator i=newInequalities.begin();i!=newInequalities.end();i++)
 
104
        normals.push_back(*i);
 
105
    }
 
106
 
 
107
 
 
108
    log1 fprintf(Stderr,"Inequalities");
 
109
    log1 AsciiPrinter(Stderr).printVectorList(normals);
 
110
    log1 fprintf(Stderr,"Equations");
 
111
    log1 AsciiPrinter(Stderr).printVectorList(equalities);
58
112
 
59
113
    PolyhedralCone c(normals,equalities,n);
60
114
    c.canonicalize();
66
120
        IntegerVectorList empty;
67
121
        HalfOpenCone C(n,c.getEquations(),c.getHalfSpaces(),empty,true);
68
122
        //      PolyhedralFan F=faceComplexOfCone(C);
69
 
        PolyhedralFan F(n);F.insert(C.closure());
70
 
        F.printWithIndices(&P,false,0);
 
123
        PolyhedralCone C2=C.closure();
 
124
        C2.canonicalize();
 
125
        PolyhedralFan F(n);F.insert(C2);
 
126
        F.printWithIndices(&P,FPF_default);
 
127
        //      F.printWithIndices(&P,false,0,false,false,optionXml.getValue());
71
128
      }
72
129
 
73
130
    return 0;