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