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