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

« back to all changes in this revision

Viewing changes to tropical_weildivisor.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:
 
1
#include "tropical_weildivisor.h"
 
2
#include "tropical2.h"
 
3
#include "log.h"
 
4
#include "printer.h"
 
5
#include <iostream>
 
6
 
 
7
PolyhedralFan weilDivisor(PolyhedralFan const &F, PolyhedralFan const &G)//, Polynomial const &g)
 
8
{
 
9
  //  PolynomialRing R=g.getRing();
 
10
  //  PolyhedralFan G=PolyhedralFan::bergmanOfPrincipalIdeal(g);
 
11
  int n=G.getAmbientDimension();
 
12
  int d=F.getMaxDimension();
 
13
 
 
14
  PolyhedralFan retTemp(n);
 
15
 
 
16
  log1 cerr<<"Computing refinement"<<endl;
 
17
  for(PolyhedralConeList::const_iterator i=F.conesBegin();i!=F.conesEnd();i++)
 
18
    {
 
19
      log1 cerr<<"*";
 
20
 
 
21
#if 1
 
22
      bool found=false;
 
23
      {
 
24
        IntegerVector v=i->getRelativeInteriorPoint();
 
25
        PolyhedralCone c=G.coneContaining(v);
 
26
        if(!(c!=*i))
 
27
          {
 
28
            retTemp.insert(c);
 
29
            found=true;
 
30
          }
 
31
      }
 
32
      if(!found)
 
33
#endif
 
34
      for(PolyhedralConeList::const_iterator j=G.conesBegin();j!=G.conesEnd();j++)
 
35
        {
 
36
          PolyhedralCone c=intersection(*i,*j);
 
37
          c.canonicalize();
 
38
          retTemp.insert(c);
 
39
        }
 
40
    }
 
41
 
 
42
  log1 cerr<<"Computing full complex"<<endl;
 
43
  retTemp=retTemp.fullComplex();
 
44
  PolyhedralFan ret(n);
 
45
 
 
46
  log1 cerr<<"Computing divisor"<<endl;
 
47
  for(PolyhedralConeList::iterator i=retTemp.conesBegin();i!=retTemp.conesEnd();i++)
 
48
    {
 
49
      log1 cerr<<"*";
 
50
      if(i->dimension()==d-1)
 
51
        {
 
52
          IntegerVector v=i->getRelativeInteriorPoint();
 
53
 
 
54
          AsciiPrinter P(Stderr);
 
55
 
 
56
          log2 P<<v<<v<<"\n";
 
57
 
 
58
          int multiplicity=0;
 
59
          IntegerVector evaluationVector(n);
 
60
 
 
61
          //      Polynomial localg=initialForm(g,v);
 
62
          //      PolyhedralFan localG=PolyhedralFan::normalFanOfNewtonPolytope(localg);
 
63
          PolyhedralFan localG=G.link(v);
 
64
 
 
65
          for(PolyhedralConeList::const_iterator j=F.conesBegin();j!=F.conesEnd();j++)
 
66
            {
 
67
              if(j->contains(v))
 
68
                {
 
69
                  IntegerVectorList equations=j->getEquations();
 
70
                  IntegerVectorList inequalities1=j->getHalfSpaces();
 
71
                  IntegerVectorList inequalities2;
 
72
                  for(IntegerVectorList::const_iterator i=inequalities1.begin();i!=inequalities1.end();i++)
 
73
                    if(dotLong(v,*i)==0)inequalities2.push_back(*i);
 
74
                  PolyhedralCone localJ(inequalities2,equations,n);
 
75
                  localJ.canonicalize();
 
76
 
 
77
                  PolyhedralFan refinement(n);
 
78
                  for(PolyhedralConeList::const_iterator k=localG.conesBegin();k!=localG.conesEnd();k++)
 
79
                    {
 
80
                      PolyhedralCone sigma=intersection(localJ,*k);
 
81
                      sigma.canonicalize();
 
82
                      refinement.insert(sigma);
 
83
                    }
 
84
 
 
85
                  for(PolyhedralConeList::const_iterator k=refinement.conesBegin();k!=refinement.conesEnd();k++)
 
86
                    {
 
87
                      PolyhedralCone const &sigma(*k);
 
88
                      if(sigma.dimension()==d)
 
89
                        {
 
90
                          /*IntegerVectorList rays=sigma.extremeRays();
 
91
                          assert(rays.size()==1);//SHOULD ALWAYS BE TRUE????
 
92
                          {
 
93
                            IntegerVector ray=*rays.begin();
 
94
                        */
 
95
                          IntegerVector ray=sigma.semiGroupGeneratorOfRay();
 
96
                          evaluationVector+=j->getMultiplicity()*ray;
 
97
                          //multiplicity+=j->getMultiplicity()*localg.degree(ray);
 
98
                          multiplicity+=j->getMultiplicity()*localG.evaluatePiecewiseLinearFunction(ray);
 
99
                        }
 
100
                    }
 
101
 
 
102
                }
 
103
            }
 
104
          //multiplicity-=localg.degree(evaluationVector);
 
105
          multiplicity-=localG.evaluatePiecewiseLinearFunction(evaluationVector);
 
106
          if(multiplicity!=0)
 
107
            {
 
108
              PolyhedralCone c=*i;
 
109
              c.setMultiplicity(multiplicity);
 
110
              ret.insert(c);
 
111
            }
 
112
        }
 
113
    }
 
114
 
 
115
  return ret;
 
116
}