source: git/dyn_modules/callgfanlib/gitfan.cc @ 62de8de

spielwiese
Last change on this file since 62de8de was 69b2c1, checked in by Yue Ren <ren@…>, 10 years ago
chg: updated callgfanlib to newest version
  • Property mode set to 100644
File size: 10.2 KB
Line 
1/***************************************************************
2 *
3 * File:       gitfan.cc
4 * Purpose:    Computationally intensive procedures for gitfan.lib,
5 *             outsourced to improve the performance.
6 * Authors:    Janko Boehm    boehm@mathematik.uni-kl.de
7 *             Simon Keicher  keicher@mail.mathematik.uni-tuebingen.de
8 *             Yue Ren        ren@mathematik.uni-kl.de
9 *
10 ***************************************************************/
11
12#include <kernel/mod2.h>
13
14#include <Singular/ipid.h>
15#include <Singular/lists.h>
16#include <Singular/ipshell.h>
17#include <libpolys/coeffs/bigintmat.h>
18
19#include <bbcone.h>
20#include <bbfan.h>
21#include <gitfan.h>
22
23namespace gitfan
24{
25
26  facet::facet():
27    eta(gfan::ZCone()),
28    interiorPoint(gfan::ZVector()),
29    facetNormal(gfan::ZVector())
30  {
31  }
32
33  facet::facet(const gitfan::facet &f):
34    eta(f.eta),
35    interiorPoint(f.interiorPoint),
36    facetNormal(f.facetNormal)
37  {
38#ifndef NDEBUG
39    gfan::ZCone c = f.eta;
40    gfan::ZVector v = f.interiorPoint;
41    gfan::ZVector w = f.facetNormal;
42    assume(c.ambientDimension() == (int)v.size());
43    assume(c.ambientDimension() == (int)w.size());
44    assume(c.contains(v));
45    assume(!c.contains(w));
46#endif
47  }
48
49  facet::facet(const gfan::ZCone &c, const gfan::ZVector &v, const gfan::ZVector &w):
50    eta(c),
51    interiorPoint(v),
52    facetNormal(w)
53  {
54#ifndef NDEBUG
55    assume(c.ambientDimension() == (int)v.size());
56    assume(c.ambientDimension() == (int)w.size());
57    assume(c.contains(v));
58    assume(!c.contains(w));
59#endif
60  }
61
62  facet::~facet()
63  {
64#ifndef NDEBUG
65    gfan::ZCone c = this->eta;
66    gfan::ZVector v = this->interiorPoint;
67    gfan::ZVector w = this->facetNormal;
68    assume(c.ambientDimension() == (int)v.size());
69    assume(c.ambientDimension() == (int)w.size());
70    assume(c.contains(v));
71    assume(!c.contains(w));
72#endif
73  }
74
75  void mergeFacets(facets &F, const facets &newFacets)
76  {
77    std::pair<facets::iterator,bool> check(newFacets.begin(),false);
78    for(facets::iterator p=newFacets.begin(); p!=newFacets.end(); p++)
79    {
80      check = F.insert(*p);
81      if(!check.second)
82        F.erase(check.first);
83    }
84  }
85
86}
87
88
89static gfan::ZCone subcone(const lists &cones, const gfan::ZVector &point)
90{
91  gfan::ZCone sigma = gfan::ZCone(gfan::ZMatrix(1,point.size()), gfan::ZMatrix(1,point.size()));
92  gfan::ZCone* zc;
93  for (int i=0; i<=cones->nr; i++)
94  {
95    zc = (gfan::ZCone*) cones->m[i].Data();
96    if (zc->contains(point))
97      sigma = gfan::intersection(sigma,*zc);
98  }
99  return(sigma);
100}
101
102static gitfan::facets interiorFacets(const gfan::ZCone &zc, const gfan::ZCone &bound)
103{
104  gfan::ZMatrix inequalities = zc.getFacets();
105  gfan::ZMatrix equations = zc.getImpliedEquations();
106  int r = inequalities.getHeight();
107  int c = inequalities.getWidth();
108  gitfan::facets F;
109  if (r*c == 0)
110    /***
111     * this is the trivial case where either we are in a zerodimensional ambient space,
112     * or the cone has no facets.
113     **/
114    return F;
115
116  // int index = 0;
117  /* next we iterate over each of the r facets, build the respective cone and add it to the list */
118  /* this is the i=0 case */
119  gfan::ZMatrix newInequalities = inequalities.submatrix(1,0,r,c);
120  gfan::ZMatrix newEquations = equations;
121  newEquations.appendRow(inequalities[0]);
122  gfan::ZCone eta = gfan::ZCone(newInequalities,newEquations);
123  eta.canonicalize();
124  gfan::ZVector v = eta.getRelativeInteriorPoint();
125  gfan::ZVector w = inequalities[0];
126
127  if (bound.containsRelatively(v))
128    F.insert(gitfan::facet(eta,v,w));
129
130  /* these are the cases i=1,...,r-2 */
131  for (int i=1; i<r-1; i++)
132  {
133    newInequalities = inequalities.submatrix(0,0,i,c);
134    newInequalities.append(inequalities.submatrix(i+1,0,r,c));
135    newEquations = equations;
136    newEquations.appendRow(inequalities[i]);
137    eta = gfan::ZCone(newInequalities,newEquations);
138    eta.canonicalize();
139    v = eta.getRelativeInteriorPoint();
140    w = inequalities[i];
141    if (bound.containsRelatively(v))
142      F.insert(gitfan::facet(eta,v,w));
143  }
144
145  /* this is the i=r-1 case */
146  newInequalities = inequalities.submatrix(0,0,r-1,c);
147  newEquations = equations;
148  newEquations.appendRow(inequalities[r-1]);
149  eta = gfan::ZCone(newInequalities,newEquations);
150  eta.canonicalize();
151
152  v = eta.getRelativeInteriorPoint();
153  w = inequalities[r-1];
154  if (bound.containsRelatively(v))
155    F.insert(gitfan::facet(eta,v,w));
156
157  return F;
158}
159
160BOOLEAN refineCones(leftv res, leftv args)
161{
162  leftv u=args;
163  if ((u != NULL) && (u->Typ() == LIST_CMD))
164  {
165    leftv v=u->next;
166    if ((v != NULL) && (v->Typ() == BIGINTMAT_CMD))
167    {
168      lists cones = (lists) u->Data();
169      bigintmat* bim = (bigintmat*) v->Data();
170      gfan::ZMatrix* zm = bigintmatToZMatrix(bim->transpose());
171      gfan::ZCone support = gfan::ZCone::givenByRays(*zm, gfan::ZMatrix(0, zm->getWidth()));
172      delete zm;
173
174      /***
175       * Randomly compute a first full-dimensional cone and insert it into the fan.
176       * Compute a list of facets and relative interior points.
177       * The relative interior points are unique, assuming the cone is stored in canonical form,
178       * which is the case in our algorithm, as we supply no redundant inequalities.
179       * Hence we can decide whether a facet need to be traversed by crosschecking
180       * its relative interior point with this list.
181       **/
182      gfan::ZCone lambda; gfan::ZVector point;
183      do
184      {
185        point = randomPoint(&support);
186        lambda = subcone(cones, point);
187      }
188      while (lambda.dimension() < lambda.ambientDimension());
189      int iterationNumber = 1;
190      std::cout << "cones found: " << iterationNumber++ << std::endl;
191
192      lambda.canonicalize();
193      gfan::ZFan* Sigma = new gfan::ZFan(lambda.ambientDimension());
194      Sigma->insert(lambda);
195      gitfan::facets F = interiorFacets(lambda, support);
196      if (F.empty())
197      {
198        res->rtyp = fanID;
199        res->data = (void*) Sigma;
200        return FALSE;
201      }
202      int mu = 1024;
203
204      gitfan::facet f;
205      gfan::ZCone eta;
206      gfan::ZVector interiorPoint;
207      gfan::ZVector facetNormal;
208      gitfan::facets newFacets;
209      while (!F.empty())
210      {
211        /***
212         * Extract a facet to traverse and its relative interior point.
213         **/
214        f = *(F.begin());
215        eta = f.getEta();
216        interiorPoint = f.getInteriorPoint();
217        facetNormal = f.getFacetNormal();
218
219        /***
220         * construct a point, which lies on the other side of the facet.
221         * make sure it lies in the known support of our fan
222         * and that the cone around the point is maximal, containing eta.
223         **/
224        point = mu * interiorPoint - facetNormal;
225        while (!support.containsRelatively(point))
226        {
227          mu = mu * 16;
228          point = mu * interiorPoint - facetNormal;
229        }
230
231        lambda = subcone(cones,point);
232        while ((lambda.dimension() < lambda.ambientDimension()) && !(lambda.contains(interiorPoint)))
233        {
234          mu = mu * 16;
235          point = mu * interiorPoint - facetNormal;
236          lambda = subcone(cones,point);
237        }
238        std::cout << "cones found: " << iterationNumber++ << std::endl;
239
240        /***
241         * insert lambda into Sigma, and create a list of facets of lambda.
242         * merge the two lists of facets
243         **/
244        lambda.canonicalize();
245        Sigma->insert(lambda);
246        newFacets = interiorFacets(lambda, support);
247        mergeFacets(F,newFacets);
248        newFacets.clear();
249      }
250      res->rtyp = fanID;
251      res->data = (void*) Sigma;
252      return FALSE;
253    }
254  }
255  WerrorS("refineCones: unexpected parameters");
256  return TRUE;
257}
258
259
260static int binomial(int n, int k)
261{
262  if (n<k)
263    return(0);
264  gfan::Integer num = 1;
265  gfan::Integer den = 1;
266  for (int i=1; i<=k; i++)
267    den = den*i;
268  for (int j=n-k+1; j<=n; j++)
269    num = num*j;
270  gfan::Integer bin = num/den;
271  return(bin.toInt());
272}
273
274
275intvec* intToAface(unsigned int v0, int n, int k)
276{
277  intvec* v = new intvec(k);
278  int j = 0;
279  for (int i=0; i<n; i++)
280  {
281    if (v0 & (1<<i))
282      (*v)[j++] = i+1;
283  }
284  return v;
285}
286
287
288BOOLEAN listOfAfacesToCheck(leftv res, leftv args)
289{
290  leftv u = args;
291  if ((u != NULL) && (u->Typ() == INT_CMD))
292  {
293    leftv v = u->next;
294    if ((v != NULL) && (v->Typ() == INT_CMD))
295    {
296      int n = (int)(long) u->Data();
297      int k = (int)(long) v->Data();
298      unsigned int v = 0;
299      for (int i=0; i<k; i++)
300        v |= 1<<i;  // sets the first k bits of v as 1
301
302      lists L = (lists)omAllocBin(slists_bin);
303      int count = (int) binomial(n,k); L->Init(count);
304      unsigned int t;
305      while (!(v & (1<<n)))
306      {
307        L->m[--count].rtyp = INTVEC_CMD;
308        L->m[count].data = (void*) intToAface(v,n,k);
309
310        // t gets v's least significant 0 bits set to 1
311        t = v | (v - 1);
312        // Next set to 1 the most significant bit to change,
313        // set to 0 the least significant ones, and add the necessary 1 bits.
314        v = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1));
315      }
316      res->rtyp = LIST_CMD;
317      res->data = (void*) L;
318      return FALSE;
319    }
320  }
321  WerrorS("listOfAfacesToCheck: unexpected parameter");
322  return TRUE;
323}
324
325
326BOOLEAN nextAfaceToCheck(leftv res, leftv args)
327{
328  leftv u = args;
329  if ((u != NULL) && (u->Typ() == INTVEC_CMD))
330  {
331    leftv v = u->next;
332    if ((v != NULL) && (v->Typ() == INT_CMD))
333    {
334      leftv w = v->next;
335      if ((w != NULL) && (w->Typ() == INT_CMD))
336      {
337        intvec* aface = (intvec*) u->Data();
338        int ambientDimension = (int)(long) v->Data();
339        int dimension = (int)(long) w->Data();
340
341        unsigned int af = 0;
342        for (int i=0; i<aface->length(); i++)
343          af |= 1<<((*aface)[i]-1);
344
345        unsigned int t = af | (af - 1);
346        af = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(af) + 1));
347
348        if (af & (1<<ambientDimension))
349        {
350          res->rtyp = INTVEC_CMD;
351          res->data = (void*) new intvec(1);
352          return FALSE;
353        }
354
355        res->rtyp = INTVEC_CMD;
356        res->data = (void*) intToAface(af,ambientDimension,dimension);
357        return FALSE;
358      }
359    }
360  }
361  WerrorS("nextAfaceToCheck: unexpected parameter");
362  return TRUE;
363}
364
365
366void gitfan_setup(SModulFunctions* p)
367{
368  p->iiAddCproc("","refineCones",FALSE,refineCones);
369  p->iiAddCproc("","listOfAfacesToCheck",FALSE,listOfAfacesToCheck);
370  p->iiAddCproc("","nextAfaceToCheck",FALSE,nextAfaceToCheck);
371}
Note: See TracBrowser for help on using the repository browser.