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 "reduce_domain.h"
  11
  12 #ifdef NTL_STD_CXX
  13 using namespace NTL;
  14 #endif
  15 using std::vector;
  16 using std::cerr;
  17 using std::endl;
  18
  19 /*
  20  * Returns the largest absolute value in the vector
  21  */
  22 static ZZ max(vec_ZZ& v)
  23 {
  24     ZZ max = abs(v[0]);
  25     for (int i = 1; i < v.length(); ++i)
  26         if (abs(v[i]) > max)
  27             max = abs(v[i]);
  28     return max;
  29 }
  30
  31 /* Remove common divisor of elements of cols of B */
  32 static void normalize_cols(mat_ZZ& B)
  33 {
  34     ZZ gcd;
  35     for (int i = 0; i < B.NumCols(); ++i) {
  36         gcd = B[0][i];
  37         for (int j = 1 ; gcd != 1 && j < B.NumRows(); ++j)
  38             GCD(gcd, gcd, B[j][i]);
  39         if (gcd != 1)
  40             for (int j = 0; j < B.NumRows(); ++j)
  41                 B[j][i] /= gcd;
  42     }
  43 }
  44
  45 class cone {
  46 public:
  47     cone(const mat_ZZ& r, int row, const vec_ZZ& w) {
  48         rays = r;
  49         rays[row] = w;
  50         set_det();
  51         set_closed(NULL);
  52     }
  53     cone(const signed_cone& sc) {
  54         rays = sc.rays;
  55         set_det();
  56         set_closed(sc.closed);
  57     }
  58     void set_det() {
  59         det = determinant(rays);
  60     }
  61     void set_closed(int *cl) {
  62         closed = NULL;
  63         if (cl) {
  64             closed = new int[rays.NumRows()];
  65             for (int i = 0; i < rays.NumRows(); ++i)
  66                 closed[i] = cl[i];
  67         }
  68     }
  69     bool needs_split(barvinok_options *options) {
  70         index = abs(det);
  71         if (IsOne(index))
  72             return false;
  73         if (options->primal && index <= options->max_index)
  74             return false;
  75
  76         Matrix *M = rays2matrix(rays);
  77         Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
  78         int ok = Matrix_Inverse(M, inv);
  79         assert(ok);
  80         Matrix_Free(M);
  81
  82         matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
  83         Matrix_Free(inv);
  84
  85         if (!options->primal && options->max_index > 1) {
  86             mat_ZZ B2 = B;
  87             normalize_cols(B2);
  88             index = abs(determinant(B2));
  89             if (index <= options->max_index)
  90                 return false;
  91         }
  92
  93         return true;
  94     }
  95
  96     void short_vector(vec_ZZ& v, vec_ZZ& lambda, barvinok_options *options) {
  97         ZZ det2;
  98         mat_ZZ U;
  99         long r = LLL(det2, B, U, options->LLL_a, options->LLL_b);
 100
 101         ZZ min = max(B[0]);
 102         int index = 0;
 103         for (int i = 1; i < B.NumRows(); ++i) {
 104             ZZ tmp = max(B[i]);
 105             if (tmp < min) {
 106                 min = tmp;
 107                 index = i;
 108             }
 109         }
 110
 111         lambda = B[index];
 112
 113         v = U[index];
 114
 115         int i;
 116         for (i = 0; i < lambda.length(); ++i)
 117             if (lambda[i] > 0)
 118                 break;
 119         if (i == lambda.length()) {
 120             v = -v;
 121             lambda = -lambda;
 122         }
 123     }
 124
 125     ~cone() {
 126         if (closed)
 127             delete [] closed;
 128     }
 129
 130     ZZ det;
 131     ZZ index;
 132     mat_ZZ rays;
 133     mat_ZZ B;
 134     int *closed;
 135 };
 136
 137 static void decompose(const signed_cone& sc, signed_cone_consumer& scc,
 138                              barvinok_options *options)
 139 {
 140     vector<cone *> nonuni;
 141     cone * c = new cone(sc);
 142     ZZ det = c->det;
 143     int s = sign(det);
 144     assert(det != 0);
 145     if (c->needs_split(options)) {
 146         nonuni.push_back(c);
 147     } else {
 148         try {
 149             options->stats->base_cones++;
 150             scc.handle(signed_cone(sc.C, sc.rays, sc.sign, to_ulong(c->index),
 151                                    sc.closed), options);
 152             delete c;
 153         } catch (...) {
 154             delete c;
 155             throw;
 156         }
 157         return;
 158     }
 159     vec_ZZ lambda;
 160     vec_ZZ v;;
 161     int closed[c->rays.NumRows()];
 162     while (!nonuni.empty()) {
 163         c = nonuni.back();
 164         nonuni.pop_back();
 165         c->short_vector(v, lambda, options);
 166         for (int i = 0; i < c->rays.NumRows(); ++i) {
 167             if (lambda[i] == 0)
 168                 continue;
 169             cone *pc = new cone(c->rays, i, v);
 170             if (c->closed) {
 171                 for (int j = 0; j < c->rays.NumRows(); ++j) {
 172                     if (lambda[j] == 0)
 173                         closed[j] = c->closed[j];
 174                     else if (j == i) {
 175                         if (lambda[i] > 0)
 176                             closed[j] = c->closed[j];
 177                         else
 178                             closed[j] = !c->closed[j];
 179                     } else if (sign(lambda[i]) == sign(lambda[j])) {
 180                         if (c->closed[i] == c->closed[j])
 181                             closed[j] = i < j;
 182                         else
 183                             closed[j] = c->closed[j];
 184                     } else
 185                         closed[j] = c->closed[i] && c->closed[j];
 186                 }
 187                 pc->set_closed(closed);
 188             }
 189             assert (pc->det != 0);
 190             if (pc->needs_split(options)) {
 191                 assert(abs(pc->det) < abs(c->det));
 192                 nonuni.push_back(pc);
 193             } else {
 194                 try {
 195                     options->stats->base_cones++;
 196                     scc.handle(signed_cone(pc->rays, sign(pc->det) * s,
 197                                            to_ulong(pc->index),
 198                                            pc->closed), options);
 199                     delete pc;
 200                 } catch (...) {
 201                     delete c;
 202                     delete pc;
 203                     while (!nonuni.empty()) {
 204                         c = nonuni.back();
 205                         nonuni.pop_back();
 206                         delete c;
 207                     }
 208                     throw;
 209                 }
 210             }
 211         }
 212         delete c;
 213     }
 214 }
 215
 216 struct polar_signed_cone_consumer : public signed_cone_consumer {
 217     signed_cone_consumer& scc;
 218     mat_ZZ r;
 219     polar_signed_cone_consumer(signed_cone_consumer& scc) : scc(scc) {}
 220     virtual void handle(const signed_cone& sc, barvinok_options *options) {
 221         Polyhedron *C = sc.C;
 222         if (!sc.C) {
 223             Matrix *M = Matrix_Alloc(sc.rays.NumRows()+1, sc.rays.NumCols()+2);
 224             for (int i = 0; i < sc.rays.NumRows(); ++i) {
 225                 zz2values(sc.rays[i], M->p[i]+1);
 226                 value_set_si(M->p[i][0], 1);
 227             }
 228             value_set_si(M->p[sc.rays.NumRows()][0], 1);
 229             value_set_si(M->p[sc.rays.NumRows()][1+sc.rays.NumCols()], 1);
 230             C = Rays2Polyhedron(M, M->NbRows+1);
 231             assert(C->NbConstraints == C->NbRays);
 232             Matrix_Free(M);
 233         }
 234         Polyhedron_Polarize(C);
 235         rays(C, r);
 236         scc.handle(signed_cone(C, r, sc.sign, sc.det), options);
 237         if (!sc.C)
 238             Polyhedron_Free(C);
 239     }
 240 };
 241
 242 /* Remove common divisor of elements of rows of B */
 243 static void normalize_rows(mat_ZZ& B)
 244 {
 245     ZZ gcd;
 246     for (int i = 0; i < B.NumRows(); ++i) {
 247         gcd = B[i][0];
 248         for (int j = 1 ; gcd != 1 && j < B.NumCols(); ++j)
 249             GCD(gcd, gcd, B[i][j]);
 250         if (gcd != 1)
 251             for (int j = 0; j < B.NumCols(); ++j)
 252                 B[i][j] /= gcd;
 253     }
 254 }
 255
 256 static void polar_decompose(Polyhedron *cone, signed_cone_consumer& scc,
 257                             barvinok_options *options)
 258 {
 259     POL_ENSURE_VERTICES(cone);
 260     Polyhedron_Polarize(cone);
 261     if (cone->NbRays - 1 != cone->Dimension) {
 262         Polyhedron *tmp = cone;
 263         cone = triangulate_cone(cone, options->MaxRays);
 264         Polyhedron_Free(tmp);
 265     }
 266     polar_signed_cone_consumer pssc(scc);
 267     mat_ZZ r;
 268     try {
 269         for (Polyhedron *Polar = cone; Polar; Polar = Polar->next) {
 270             rays(Polar, r);
 271             normalize_rows(r);
 272             decompose(signed_cone(Polar, r, 1), pssc, options);
 273         }
 274         Domain_Free(cone);
 275     } catch (...) {
 276         Domain_Free(cone);
 277         throw;
 278     }
 279 }
 280
 281 static void primal_decompose(Polyhedron *cone, signed_cone_consumer& scc,
 282                              barvinok_options *options)
 283 {
 284     POL_ENSURE_VERTICES(cone);
 285     Polyhedron *parts;
 286     if (cone->NbRays - 1 == cone->Dimension)
 287         parts = cone;
 288     else
 289         parts = triangulate_cone(cone, options->MaxRays);
 290     int closed[cone->Dimension];
 291     Vector *average = NULL;
 292     Value tmp;
 293     if (parts != cone) {
 294         value_init(tmp);
 295         average = inner_point(cone);
 296     }
 297     mat_ZZ ray;
 298     try {
 299         for (Polyhedron *simple = parts; simple; simple = simple->next) {
 300             for (int i = 0, r = 0; r < simple->NbRays; ++r) {
 301                 if (value_notzero_p(simple->Ray[r][1+simple->Dimension]))
 302                     continue;
 303                 if (simple == cone) {
 304                     closed[i] = 1;
 305                 } else {
 306                     int f;
 307                     for (f = 0; f < simple->NbConstraints; ++f) {
 308                         Inner_Product(simple->Ray[r]+1, simple->Constraint[f]+1,
 309                                       simple->Dimension, &tmp);
 310                         if (value_notzero_p(tmp))
 311                             break;
 312                     }
 313                     assert(f < simple->NbConstraints);
 314                     closed[i] = is_internal(average, simple->Constraint[f]);
 315                 }
 316                 ++i;
 317             }
 318             rays(simple, ray);
 319             decompose(signed_cone(simple, ray, 1, 0, closed), scc, 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     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 = supporting_cone_p(P, V);
 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, options);
 386     *ppos = pc.pos;
 387     *pneg = pc.neg;
 388     barvinok_options_free(options);
 389 }
 390