decomposer.cc

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