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         set_closed(NULL);
  35     }
  36     cone(const signed_cone& sc) {
  37         Cone = Polyhedron_Copy(sc.C);
  38         Rays = rays(sc.C);
  39         set_det();
  40         set_closed(sc.closed);
  41     }
  42     void set_det() {
  43         mat_ZZ A;
  44         matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
  45         det = determinant(A);
  46     }
  47     void set_closed(int *cl) {
  48         closed = NULL;
  49         if (cl) {
  50             closed = new int[Rays->NbRows-1];
  51             for (int i = 0; i < Rays->NbRows-1; ++i)
  52                 closed[i] = cl[i];
  53         }
  54     }
  55
  56     Vector* short_vector(vec_ZZ& lambda, barvinok_options *options) {
  57         Matrix *M = Matrix_Copy(Rays);
  58         Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
  59         int ok = Matrix_Inverse(M, inv);
  60         assert(ok);
  61         Matrix_Free(M);
  62
  63         ZZ det2;
  64         mat_ZZ B;
  65         mat_ZZ U;
  66         matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
  67         long r = LLL(det2, B, U, options->LLL_a, options->LLL_b);
  68
  69         ZZ min = max(B[0]);
  70         int index = 0;
  71         for (int i = 1; i < B.NumRows(); ++i) {
  72             ZZ tmp = max(B[i]);
  73             if (tmp < min) {
  74                 min = tmp;
  75                 index = i;
  76             }
  77         }
  78
  79         Matrix_Free(inv);
  80
  81         lambda = B[index];
  82
  83         Vector *z = Vector_Alloc(U[index].length()+1);
  84         assert(z);
  85         zz2values(U[index], z->p);
  86         value_set_si(z->p[U[index].length()], 0);
  87
  88         Polyhedron *C = poly();
  89         int i;
  90         for (i = 0; i < lambda.length(); ++i)
  91             if (lambda[i] > 0)
  92                 break;
  93         if (i == lambda.length()) {
  94             Value tmp;
  95             value_init(tmp);
  96             value_set_si(tmp, -1);
  97             Vector_Scale(z->p, z->p, tmp, z->Size-1);
  98             value_clear(tmp);
  99         }
 100         return z;
 101     }
 102
 103     ~cone() {
 104         Polyhedron_Free(Cone);
 105         Matrix_Free(Rays);
 106         if (closed)
 107             delete [] closed;
 108     }
 109
 110     Polyhedron *poly() {
 111         if (!Cone) {
 112             Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
 113             for (int i = 0; i < Rays->NbRows; ++i) {
 114                 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
 115                 value_set_si(M->p[i][0], 1);
 116             }
 117             Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
 118             value_set_si(M->p[Rays->NbRows][0], 1);
 119             value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
 120             Cone = Rays2Polyhedron(M, M->NbRows+1);
 121             assert(Cone->NbConstraints == Cone->NbRays);
 122             Matrix_Free(M);
 123         }
 124         return Cone;
 125     }
 126
 127     ZZ det;
 128     Polyhedron *Cone;
 129     Matrix *Rays;
 130     int *closed;
 131 };
 132
 133 static void decompose(const signed_cone& sc, signed_cone_consumer& scc,
 134                              barvinok_options *options)
 135 {
 136     vector<cone *> nonuni;
 137     cone * c = new cone(sc);
 138     ZZ det = c->det;
 139     int s = sign(det);
 140     assert(det != 0);
 141     if (abs(det) > 1) {
 142         nonuni.push_back(c);
 143     } else {
 144         try {
 145             options->stats.unimodular_cones++;
 146             scc.handle(sc);
 147             delete c;
 148         } catch (...) {
 149             delete c;
 150             throw;
 151         }
 152         return;
 153     }
 154     vec_ZZ lambda;
 155     int closed[c->Rays->NbRows-1];
 156     while (!nonuni.empty()) {
 157         c = nonuni.back();
 158         nonuni.pop_back();
 159         Vector* v = c->short_vector(lambda, options);
 160         for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
 161             if (lambda[i] == 0)
 162                 continue;
 163             Matrix* M = Matrix_Copy(c->Rays);
 164             Vector_Copy(v->p, M->p[i], v->Size);
 165             cone * pc = new cone(M);
 166             if (c->closed) {
 167                 bool same_sign = sign(c->det) * sign(pc->det) > 0;
 168                 for (int j = 0; j < c->Rays->NbRows - 1; ++j) {
 169                     if (lambda[j] == 0)
 170                         closed[j] = c->closed[j];
 171                     else if (j == i) {
 172                         if (same_sign)
 173                             closed[j] = c->closed[j];
 174                         else
 175                             closed[j] = !c->closed[j];
 176                     } else if (sign(lambda[i]) * sign(lambda[j]) > 0) {
 177                         if (c->closed[i] == c->closed[j])
 178                             closed[j] = i < j;
 179                         else
 180                             closed[j] = c->closed[j];
 181                     } else
 182                         closed[j] = c->closed[i] && c->closed[j];
 183                 }
 184                 pc->set_closed(closed);
 185             }
 186             assert (pc->det != 0);
 187             if (abs(pc->det) > 1) {
 188                 assert(abs(pc->det) < abs(c->det));
 189                 nonuni.push_back(pc);
 190             } else {
 191                 try {
 192                     options->stats.unimodular_cones++;
 193                     if (pc->closed)
 194                         scc.handle(signed_cone(pc->poly(), sign(pc->det) * s,
 195                                                pc->closed));
 196                     else
 197                         scc.handle(signed_cone(pc->poly(), sign(pc->det) * s));
 198                     delete pc;
 199                 } catch (...) {
 200                     delete c;
 201                     delete pc;
 202                     while (!nonuni.empty()) {
 203                         c = nonuni.back();
 204                         nonuni.pop_back();
 205                         delete c;
 206                     }
 207                     Matrix_Free(M);
 208                     Vector_Free(v);
 209                     throw;
 210                 }
 211             }
 212             Matrix_Free(M);
 213         }
 214         Vector_Free(v);
 215         delete c;
 216     }
 217 }
 218
 219 struct polar_signed_cone_consumer : public signed_cone_consumer {
 220     signed_cone_consumer& scc;
 221     polar_signed_cone_consumer(signed_cone_consumer& scc) : scc(scc) {}
 222     virtual void handle(const signed_cone& sc) {
 223         Polyhedron_Polarize(sc.C);
 224         scc.handle(sc);
 225     }
 226 };
 227
 228 static void polar_decompose(Polyhedron *cone, signed_cone_consumer& scc,
 229                             barvinok_options *options)
 230 {
 231     POL_ENSURE_VERTICES(cone);
 232     Polyhedron_Polarize(cone);
 233     if (cone->NbRays - 1 != cone->Dimension) {
 234         Polyhedron *tmp = cone;
 235         cone = triangulate_cone(cone, options->MaxRays);
 236         Polyhedron_Free(tmp);
 237     }
 238     polar_signed_cone_consumer pssc(scc);
 239     try {
 240         for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
 241             decompose(signed_cone(Polar, 1), pssc, options);
 242         Domain_Free(cone);
 243     } catch (...) {
 244         Domain_Free(cone);
 245         throw;
 246     }
 247 }
 248
 249 static void primal_decompose(Polyhedron *cone, signed_cone_consumer& scc,
 250                              barvinok_options *options)
 251 {
 252     POL_ENSURE_VERTICES(cone);
 253     Polyhedron *parts;
 254     if (cone->NbRays - 1 == cone->Dimension)
 255         parts = cone;
 256     else
 257         parts = triangulate_cone(cone, options->MaxRays);
 258     int closed[cone->Dimension];
 259     Vector *average = NULL;
 260     Value tmp;
 261     if (parts != cone) {
 262         value_init(tmp);
 263         average = Vector_Alloc(cone->Dimension);
 264         for (int i = 0; i < cone->NbRays; ++i) {
 265             if (value_notzero_p(cone->Ray[i][1+cone->Dimension]))
 266                 continue;
 267             Vector_Add(average->p, cone->Ray[i]+1, average->p, cone->Dimension);
 268         }
 269     }
 270     try {
 271         for (Polyhedron *simple = parts; simple; simple = simple->next) {
 272             for (int i = 0, r = 0; r < simple->NbRays; ++r) {
 273                 if (value_notzero_p(simple->Ray[r][1+simple->Dimension]))
 274                     continue;
 275                 if (simple == cone) {
 276                     closed[i] = 1;
 277                 } else {
 278                     int f;
 279                     for (f = 0; f < simple->NbConstraints; ++f) {
 280                         Inner_Product(simple->Ray[r]+1, simple->Constraint[f]+1,
 281                                       simple->Dimension, &tmp);
 282                         if (value_notzero_p(tmp))
 283                             break;
 284                     }
 285                     assert(f < simple->NbConstraints);
 286                     Inner_Product(simple->Constraint[f]+1, average->p,
 287                                   simple->Dimension, &tmp);
 288                     if (value_notzero_p(tmp))
 289                         closed[i] = value_pos_p(tmp);
 290                     else {
 291                         int p = First_Non_Zero(simple->Constraint[f]+1,
 292                                                simple->Dimension);
 293                         closed[i] = value_pos_p(simple->Constraint[f][1+p]);
 294                     }
 295                 }
 296                 ++i;
 297             }
 298             decompose(signed_cone(simple, 1, closed), scc, options);
 299         }
 300         Domain_Free(parts);
 301         if (parts != cone) {
 302             Domain_Free(cone);
 303             value_clear(tmp);
 304             Vector_Free(average);
 305         }
 306     } catch (...) {
 307         Domain_Free(parts);
 308         if (parts != cone) {
 309             Domain_Free(cone);
 310             value_clear(tmp);
 311             Vector_Free(average);
 312         }
 313         throw;
 314     }
 315 }
 316
 317 void barvinok_decompose(Polyhedron *C, signed_cone_consumer& scc,
 318                         barvinok_options *options)
 319 {
 320     if (options->primal)
 321         primal_decompose(C, scc, options);
 322     else
 323         polar_decompose(C, scc, options);
 324 }
 325
 326 void vertex_decomposer::decompose_at_vertex(Param_Vertices *V, int _i,
 327                                             barvinok_options *options)
 328 {
 329     Polyhedron *C = supporting_cone_p(P, V);
 330     vert = _i;
 331     this->V = V;
 332
 333     barvinok_decompose(C, scc, options);
 334 }
 335
 336 struct posneg_collector : public signed_cone_consumer {
 337     posneg_collector(Polyhedron *pos, Polyhedron *neg) : pos(pos), neg(neg) {}
 338     virtual void handle(const signed_cone& sc) {
 339         Polyhedron *p = Polyhedron_Copy(sc.C);
 340         if (sc.sign > 0) {
 341             p->next = pos;
 342             pos = p;
 343         } else {
 344             p->next = neg;
 345             neg = p;
 346         }
 347     }
 348     Polyhedron *pos;
 349     Polyhedron *neg;
 350 };
 351
 352 /*
 353  * Barvinok's Decomposition of a simplicial cone
 354  *
 355  * Returns two lists of polyhedra
 356  */
 357 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
 358 {
 359     barvinok_options *options = barvinok_options_new_with_defaults();
 360     posneg_collector pc(*ppos, *pneg);
 361     POL_ENSURE_VERTICES(C);
 362     decompose(signed_cone(C, 1), pc, options);
 363     *ppos = pc.pos;
 364     *pneg = pc.neg;
 365     free(options);
 366 }
 367