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