decomposer.cc

   1 #include <vector>
   2 #include <assert.h>
   3 #include <gmp.h>
   4 #include <NTL/mat_ZZ.h>
   5 #include <NTL/LLL.h>
   6 #include <barvinok/barvinok.h>
   7 #include <barvinok/util.h>
   8 #include "conversion.h"
   9 #include "decomposer.h"
  10
  11 #ifdef NTL_STD_CXX
  12 using namespace NTL;
  13 #endif
  14 using std::vector;
  15
  16 /*
  17  * Returns the largest absolute value in the vector
  18  */
  19 static ZZ max(vec_ZZ& v)
  20 {
  21     ZZ max = abs(v[0]);
  22     for (int i = 1; i < v.length(); ++i)
  23         if (abs(v[i]) > max)
  24             max = abs(v[i]);
  25     return max;
  26 }
  27
  28 class cone {
  29 public:
  30     cone(Matrix *M) {
  31         Cone = 0;
  32         Rays = Matrix_Copy(M);
  33         set_det();
  34     }
  35     cone(Polyhedron *C) {
  36         Cone = Polyhedron_Copy(C);
  37         Rays = rays(C);
  38         set_det();
  39     }
  40     void set_det() {
  41         mat_ZZ A;
  42         matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
  43         det = determinant(A);
  44     }
  45
  46     Vector* short_vector(vec_ZZ& lambda, barvinok_options *options) {
  47         Matrix *M = Matrix_Copy(Rays);
  48         Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
  49         int ok = Matrix_Inverse(M, inv);
  50         assert(ok);
  51         Matrix_Free(M);
  52
  53         ZZ det2;
  54         mat_ZZ B;
  55         mat_ZZ U;
  56         matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
  57         long r = LLL(det2, B, U, options->LLL_a, options->LLL_b);
  58
  59         ZZ min = max(B[0]);
  60         int index = 0;
  61         for (int i = 1; i < B.NumRows(); ++i) {
  62             ZZ tmp = max(B[i]);
  63             if (tmp < min) {
  64                 min = tmp;
  65                 index = i;
  66             }
  67         }
  68
  69         Matrix_Free(inv);
  70
  71         lambda = B[index];
  72
  73         Vector *z = Vector_Alloc(U[index].length()+1);
  74         assert(z);
  75         zz2values(U[index], z->p);
  76         value_set_si(z->p[U[index].length()], 0);
  77
  78         Polyhedron *C = poly();
  79         int i;
  80         for (i = 0; i < lambda.length(); ++i)
  81             if (lambda[i] > 0)
  82                 break;
  83         if (i == lambda.length()) {
  84             Value tmp;
  85             value_init(tmp);
  86             value_set_si(tmp, -1);
  87             Vector_Scale(z->p, z->p, tmp, z->Size-1);
  88             value_clear(tmp);
  89         }
  90         return z;
  91     }
  92
  93     ~cone() {
  94         Polyhedron_Free(Cone);
  95         Matrix_Free(Rays);
  96     }
  97
  98     Polyhedron *poly() {
  99         if (!Cone) {
 100             Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
 101             for (int i = 0; i < Rays->NbRows; ++i) {
 102                 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
 103                 value_set_si(M->p[i][0], 1);
 104             }
 105             Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
 106             value_set_si(M->p[Rays->NbRows][0], 1);
 107             value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
 108             Cone = Rays2Polyhedron(M, M->NbRows+1);
 109             assert(Cone->NbConstraints == Cone->NbRays);
 110             Matrix_Free(M);
 111         }
 112         return Cone;
 113     }
 114
 115     ZZ det;
 116     Polyhedron *Cone;
 117     Matrix *Rays;
 118 };
 119
 120 static void primal_decompose(Polyhedron *C, signed_cone_consumer& scc,
 121                              barvinok_options *options)
 122 {
 123     vector<cone *> nonuni;
 124     cone * c = new cone(C);
 125     ZZ det = c->det;
 126     int s = sign(det);
 127     assert(det != 0);
 128     if (abs(det) > 1) {
 129         nonuni.push_back(c);
 130     } else {
 131         try {
 132             options->stats.unimodular_cones++;
 133             scc.handle(signed_cone(C, 1));
 134             delete c;
 135         } catch (...) {
 136             delete c;
 137             throw;
 138         }
 139     }
 140     vec_ZZ lambda;
 141     while (!nonuni.empty()) {
 142         c = nonuni.back();
 143         nonuni.pop_back();
 144         Vector* v = c->short_vector(lambda, options);
 145         for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
 146             if (lambda[i] == 0)
 147                 continue;
 148             Matrix* M = Matrix_Copy(c->Rays);
 149             Vector_Copy(v->p, M->p[i], v->Size);
 150             cone * pc = new cone(M);
 151             assert (pc->det != 0);
 152             if (abs(pc->det) > 1) {
 153                 assert(abs(pc->det) < abs(c->det));
 154                 nonuni.push_back(pc);
 155             } else {
 156                 try {
 157                     options->stats.unimodular_cones++;
 158                     scc.handle(signed_cone(pc->poly(), sign(pc->det) * s));
 159                     delete pc;
 160                 } catch (...) {
 161                     delete c;
 162                     delete pc;
 163                     while (!nonuni.empty()) {
 164                         c = nonuni.back();
 165                         nonuni.pop_back();
 166                         delete c;
 167                     }
 168                     Matrix_Free(M);
 169                     Vector_Free(v);
 170                     throw;
 171                 }
 172             }
 173             Matrix_Free(M);
 174         }
 175         Vector_Free(v);
 176         delete c;
 177     }
 178 }
 179
 180 struct polar_signed_cone_consumer : public signed_cone_consumer {
 181     signed_cone_consumer& scc;
 182     polar_signed_cone_consumer(signed_cone_consumer& scc) : scc(scc) {}
 183     virtual void handle(const signed_cone& sc) {
 184         Polyhedron_Polarize(sc.C);
 185         scc.handle(sc);
 186     }
 187 };
 188
 189 static void polar_decompose(Polyhedron *cone, signed_cone_consumer& scc,
 190                             barvinok_options *options)
 191 {
 192     POL_ENSURE_VERTICES(cone);
 193     Polyhedron_Polarize(cone);
 194     if (cone->NbRays - 1 != cone->Dimension) {
 195         Polyhedron *tmp = cone;
 196         cone = triangulate_cone(cone, options->MaxRays);
 197         Polyhedron_Free(tmp);
 198     }
 199     polar_signed_cone_consumer pssc(scc);
 200     try {
 201         for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
 202             primal_decompose(Polar, pssc, options);
 203         Domain_Free(cone);
 204     } catch (...) {
 205         Domain_Free(cone);
 206         throw;
 207     }
 208 }
 209
 210 void barvinok_decompose(Polyhedron *C, signed_cone_consumer& scc,
 211                         barvinok_options *options)
 212 {
 213     polar_decompose(C, scc, options);
 214 }
 215
 216 void vertex_decomposer::decompose_at_vertex(Param_Vertices *V, int _i,
 217                                             barvinok_options *options)
 218 {
 219     Polyhedron *C = supporting_cone_p(P, V);
 220     vert = _i;
 221     this->V = V;
 222
 223     barvinok_decompose(C, scc, options);
 224 }
 225
 226 /*
 227  * Barvinok's Decomposition of a simplicial cone
 228  *
 229  * Returns two lists of polyhedra
 230  */
 231 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
 232 {
 233     barvinok_options *options = barvinok_options_new_with_defaults();
 234     Polyhedron *pos = *ppos, *neg = *pneg;
 235     vector<cone *> nonuni;
 236     cone * c = new cone(C);
 237     ZZ det = c->det;
 238     int s = sign(det);
 239     assert(det != 0);
 240     if (abs(det) > 1) {
 241         nonuni.push_back(c);
 242     } else {
 243         Polyhedron *p = Polyhedron_Copy(c->Cone);
 244         p->next = pos;
 245         pos = p;
 246         delete c;
 247     }
 248     vec_ZZ lambda;
 249     while (!nonuni.empty()) {
 250         c = nonuni.back();
 251         nonuni.pop_back();
 252         Vector* v = c->short_vector(lambda, options);
 253         for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
 254             if (lambda[i] == 0)
 255                 continue;
 256             Matrix* M = Matrix_Copy(c->Rays);
 257             Vector_Copy(v->p, M->p[i], v->Size);
 258             cone * pc = new cone(M);
 259             assert (pc->det != 0);
 260             if (abs(pc->det) > 1) {
 261                 assert(abs(pc->det) < abs(c->det));
 262                 nonuni.push_back(pc);
 263             } else {
 264                 Polyhedron *p = pc->poly();
 265                 pc->Cone = 0;
 266                 if (sign(pc->det) == s) {
 267                     p->next = pos;
 268                     pos = p;
 269                 } else {
 270                     p->next = neg;
 271                     neg = p;
 272                 }
 273                 delete pc;
 274             }
 275             Matrix_Free(M);
 276         }
 277         Vector_Free(v);
 278         delete c;
 279     }
 280     *ppos = pos;
 281     *pneg = neg;
 282     free(options);
 283 }
 284