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