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