decomposer.cc

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