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