topcom.c

   1 #include <assert.h>
   2 #include <barvinok/util.h>
   3 #include <barvinok/options.h>
   4 #include <unistd.h>
   5 #include "normalization.h"
   6 #include "param_util.h"
   7 #include "topcom.h"
   8 #include "config.h"
   9
  10 #define ALLOC(type) (type*)malloc(sizeof(type))
  11
  12 void run_points2triangs(pid_t *child, int *in, int *out)
  13 {
  14     int in_fd[2], out_fd[2];
  15
  16     if (pipe(in_fd))
  17         assert(0);
  18     if (pipe(out_fd))
  19         assert(0);
  20     *child = fork();
  21     assert(*child >= 0);
  22     if (!*child) {
  23         int rc;
  24
  25         dup2(in_fd[0], 0);
  26         dup2(out_fd[1], 1);
  27         close(in_fd[0]);
  28         close(out_fd[1]);
  29         close(in_fd[1]);
  30         close(out_fd[0]);
  31
  32         rc = execl(POINTS2TRIANGS_PATH, "points2triangs", "--regular", NULL);
  33         assert(0);
  34     }
  35     close(in_fd[0]);
  36     close(out_fd[1]);
  37     *in = in_fd[1];
  38     *out = out_fd[0];
  39 }
  40
  41 struct domain {
  42     Param_Domain    domain;
  43     int             F_len;
  44 };
  45
  46 static Param_Vertices *construct_vertex(unsigned *vertex_facets,
  47                                        Matrix *Constraints,
  48                                        int d, unsigned nparam, unsigned MaxRays)
  49 {
  50     unsigned nvar = Constraints->NbColumns-2 - nparam;
  51     Matrix *A = Matrix_Alloc(nvar+1, nvar+1);
  52     Matrix *inv = Matrix_Alloc(nvar+1, nvar+1);
  53     Matrix *B = Matrix_Alloc(nvar, nparam+2);
  54     Matrix *V = Matrix_Alloc(nvar, nparam+2);
  55     Matrix *Domain = Matrix_Alloc(d-nvar, nparam+2);
  56     Polyhedron *AD;
  57     unsigned bx;
  58     int i, j, ix;
  59     int ok;
  60     Param_Vertices *vertex;
  61
  62     for (j = 0, i = 0, ix = 0, bx = MSB; i < d; ++i) {
  63         if ((vertex_facets[ix] & bx) == bx) {
  64             Vector_Copy(Constraints->p[i]+1, A->p[j], nvar);
  65             Vector_Oppose(Constraints->p[i]+1+nvar, B->p[j++], nparam+1);
  66         }
  67         NEXT(ix, bx);
  68     }
  69     assert(j == nvar);
  70     value_set_si(A->p[nvar][nvar], 1);
  71     ok = Matrix_Inverse(A, inv);
  72     assert(ok);
  73     Matrix_Free(A);
  74     inv->NbRows = nvar;
  75     inv->NbColumns = nvar;
  76     Matrix_Product(inv, B, V);
  77     Matrix_Free(B);
  78     for (i = 0; i < nvar; ++i) {
  79         value_assign(V->p[i][nparam+1], inv->p[nvar][nvar]);
  80         Vector_Normalize(V->p[i], V->NbColumns);
  81     }
  82     Matrix_Free(inv);
  83     for (j = 0, i = 0, ix = 0, bx = MSB; i < d; ++i) {
  84         if ((vertex_facets[ix] & bx) == bx) {
  85             NEXT(ix, bx);
  86             continue;
  87         }
  88         Param_Inner_Product(Constraints->p[i], V, Domain->p[j]);
  89         if (First_Non_Zero(Domain->p[j]+1, nparam+1) == -1)
  90             vertex_facets[ix] |= bx;
  91         else
  92             value_set_si(Domain->p[j++][0], 1);
  93         NEXT(ix, bx);
  94     }
  95     Domain->NbRows = j;
  96     A = Matrix_Copy(Domain);
  97     AD = Constraints2Polyhedron(A, MaxRays);
  98     Matrix_Free(A);
  99     POL_ENSURE_VERTICES(AD);
 100     /* A vertex with a lower-dimensional activity domain
 101      * saturates more facets than derived above and is actually
 102      * the superimposition of two or more vertices.
 103      * We therefore discard the domain and (ultimately)
 104      * the chamber containing it.
 105      * We keep the vertex, though, since it may reappear
 106      * in other chambers, which will then likewise be discarded.
 107      * The same holds if the activity domain is empty.
 108      */
 109     if (AD->NbEq > 0) {
 110         Matrix_Free(Domain);
 111         Domain = NULL;
 112     }
 113     Polyhedron_Free(AD);
 114     vertex = calloc(1, sizeof(Param_Vertices));
 115     vertex->Facets = vertex_facets;
 116     vertex->Vertex = V;
 117     vertex->Domain = Domain;
 118     return vertex;
 119 }
 120
 121 static int add_vertex_to_domain(Param_Vertices **vertices, int words,
 122                                 unsigned *vertex_facets,
 123                                 Matrix *Constraints, int d, unsigned nparam,
 124                                 struct domain *domain,
 125                                 unsigned MaxRays)
 126 {
 127     Param_Vertices *vertex;
 128     unsigned vbx;
 129     int vi, vix;
 130
 131     for (vi = 0, vix = 0, vbx = MSB;
 132             *vertices;
 133             vertices = &(*vertices)->next, ++vi) {
 134         int i;
 135         for (i = 0; i < words; ++i)
 136             if (((*vertices)->Facets[i] & vertex_facets[i]) != vertex_facets[i])
 137                 break;
 138         if (i >= words) {
 139             if (!(*vertices)->Domain)
 140                 domain->F_len = 0;
 141             else
 142                 domain->domain.F[vix] |= vbx;
 143             free(vertex_facets);
 144             return 0;
 145         }
 146         NEXT(vix, vbx);
 147     }
 148     if (domain->F_len <= vix) {
 149         domain->F_len++;
 150         domain->domain.F = realloc(domain->domain.F,
 151                                    domain->F_len * sizeof(unsigned));
 152         domain->domain.F[domain->F_len-1] = 0;
 153     }
 154     vertex = construct_vertex(vertex_facets, Constraints, d, nparam, MaxRays);
 155     if (!vertex->Domain)
 156         domain->F_len = 0;
 157     else
 158         domain->domain.F[vix] |= vbx;
 159     vertex->next = *vertices;
 160     *vertices = vertex;
 161     return 1;
 162 }
 163
 164 static void compute_domain(struct domain *domain, Param_Vertices *vertices,
 165                            Polyhedron *C, unsigned MaxRays)
 166 {
 167     unsigned bx;
 168     int ix, j;
 169     int nbV = bit_vector_count(domain->domain.F, domain->F_len);
 170     unsigned cols;
 171     unsigned rows;
 172     Matrix *Constraints;
 173
 174     for (ix = 0, bx = MSB; !vertices->Domain; vertices = vertices->next)
 175         NEXT(ix, bx);
 176
 177     cols = vertices->Domain->NbColumns;
 178     rows = vertices->Domain->NbRows;
 179     Constraints = Matrix_Alloc(nbV * rows + C->NbConstraints, cols);
 180
 181     for (j = 0; vertices; vertices = vertices->next) {
 182         if ((domain->domain.F[ix] & bx) == bx)
 183             Vector_Copy(vertices->Domain->p[0],
 184                         Constraints->p[(j++)*rows], rows * cols);
 185         NEXT(ix, bx);
 186     }
 187     Vector_Copy(C->Constraint[0], Constraints->p[j*rows], C->NbConstraints * cols);
 188     domain->domain.Domain = Constraints2Polyhedron(Constraints, MaxRays);
 189     Matrix_Free(Constraints);
 190 }
 191
 192 static void add_domain(struct domain **domains, struct domain *domain,
 193                        Param_Vertices *vertices, Polyhedron *C,
 194                        struct barvinok_options *options)
 195 {
 196     options->stats->topcom_chambers++;
 197
 198     for (; *domains; domains = (struct domain **)&(*domains)->domain.next) {
 199         int i;
 200         for (i = 0; i < (*domains)->F_len; ++i)
 201             if (((*domains)->domain.F[i] & domain->domain.F[i])
 202                         != domain->domain.F[i])
 203                 break;
 204         if (i < (*domains)->F_len)
 205             continue;
 206         for (; i < domain->F_len; ++i)
 207             if (domain->domain.F[i])
 208                 break;
 209         if (i >= domain->F_len) {
 210             Param_Domain_Free(&domain->domain);
 211             return;
 212         }
 213     }
 214     options->stats->topcom_distinct_chambers++;
 215     compute_domain(domain, vertices, C, options->MaxRays);
 216     *domains = domain;
 217 }
 218
 219 #define INT_BITS (sizeof(unsigned) * 8)
 220
 221 /* Remove any empty or lower-dimensional chamber.  The latter
 222  * lie on the boundary of the context and are facets of other chambers.
 223  *
 224  * While we are examining each chamber, also extend the F vector
 225  * of each chamber to the maximum.
 226  */
 227 static void remove_empty_chambers(Param_Domain **PD, unsigned vertex_words)
 228 {
 229     while (*PD) {
 230         int remove = 0;
 231         int i;
 232
 233         if ((*PD)->Domain->NbEq > 0)
 234             remove = 1;
 235         else {
 236             POL_ENSURE_FACETS((*PD)->Domain);
 237             if ((*PD)->Domain->NbEq > 0)
 238                 remove = 1;
 239         }
 240         if (remove) {
 241             Param_Domain *D = *PD;
 242             *PD = (*PD)->next;
 243             D->next = NULL;
 244             Param_Domain_Free(D);
 245             continue;
 246         }
 247         if ((i = ((struct domain*)(*PD))->F_len) < vertex_words)
 248             (*PD)->F = realloc((*PD)->F, vertex_words * sizeof(unsigned));
 249         for (; i < vertex_words; ++i)
 250             (*PD)->F[i] = 0;
 251         PD = &(*PD)->next;
 252     }
 253 }
 254
 255 static Param_Polyhedron *points2triangs(Matrix *K, Polyhedron *P, Polyhedron *C,
 256                                         struct barvinok_options *options)
 257 {
 258     int in, out;
 259     int i, j;
 260     pid_t child;
 261     FILE *fin, *fout;
 262     int d = K->NbRows;
 263     int words = (d+INT_BITS-1)/INT_BITS;
 264     Param_Vertices *vertices = NULL;
 265     struct domain *domains = NULL;
 266     int vertex_words = 1;
 267     Param_Polyhedron *PP = ALLOC(Param_Polyhedron);
 268     unsigned MaxRays = options->MaxRays;
 269
 270     PP->Rays = NULL;
 271     PP->nbV = 0;
 272     PP->Constraints = Polyhedron2Constraints(P);
 273     /* We need the exact facets, because we may make some of them open later */
 274     POL_UNSET(options->MaxRays, POL_INTEGER);
 275
 276     run_points2triangs(&child, &in, &out);
 277
 278     fin = fdopen(in, "w");
 279     fprintf(fin, "[\n");
 280     for (i = 0; i < K->NbRows; ++i) {
 281         fprintf(fin, "[");
 282         for (j = 0; j < K->NbColumns; ++j)
 283             value_print(fin, P_VALUE_FMT, K->p[i][j]);
 284         fprintf(fin, "]");
 285     }
 286     fprintf(fin, "]\n");
 287     fclose(fin);
 288
 289     fout = fdopen(out, "r");
 290     while (fscanf(fout, "T[%d]:=", &i) == 1) {
 291         int a, b, c;
 292         struct domain *domain = ALLOC(struct domain);
 293         memset(domain, 0, sizeof(*domain));
 294         domain->domain.F = calloc(vertex_words, sizeof(unsigned));
 295         domain->F_len = vertex_words;
 296
 297         c = fgetc(fout);
 298         if (c == '[')
 299                 fscanf(fout, "%d,%d:{", &a, &b);
 300
 301         while (fgetc(fout) == '{') {    /* '{' or closing '}' */
 302             int c;
 303             unsigned *F = calloc(words, sizeof(unsigned));
 304
 305             for (j = 0; j < K->NbColumns; ++j) {
 306                 unsigned v, shift;
 307                 fscanf(fout, "%d", &v);
 308                 shift = INT_BITS - (v % INT_BITS) - 1;
 309                 F[v / INT_BITS] |= 1u << shift;
 310                 fgetc(fout); /* , or } */
 311             }
 312             if (!domain->F_len)
 313                 free(F);
 314             else if (add_vertex_to_domain(&vertices, words, F, PP->Constraints,
 315                                           d, C->Dimension,
 316                                           domain, options->MaxRays))
 317                 ++PP->nbV;
 318             if ((c = fgetc(fout)) != ',')       /* , or } */
 319                 ungetc(c, fout);
 320         }
 321         if (domain->F_len)
 322             vertex_words = domain->F_len;
 323         if (c == '[')
 324                 fgetc(fout); /* ] */
 325         fgetc(fout); /* ; */
 326         fgetc(fout); /* \n */
 327         if (bit_vector_count(domain->domain.F, domain->F_len) > 0)
 328             add_domain(&domains, domain, vertices, C, options);
 329         else {
 330             options->stats->topcom_empty_chambers++;
 331             Param_Domain_Free(&domain->domain);
 332         }
 333     }
 334     fclose(fout);
 335
 336     PP->V = vertices;
 337     PP->D = &domains->domain;
 338
 339     remove_empty_chambers(&PP->D, vertex_words);
 340
 341     options->MaxRays = MaxRays;
 342
 343     return PP;
 344 }
 345
 346 /* Assumes M is of full row rank */
 347 static Matrix *null_space(Matrix *M)
 348 {
 349     Matrix *H, *Q, *U;
 350     Matrix *N;
 351     int i;
 352
 353     left_hermite(M, &H, &Q, &U);
 354     N = Matrix_Alloc(M->NbColumns, M->NbColumns - M->NbRows);
 355     for (i = 0; i < N->NbRows; ++i)
 356         Vector_Copy(U->p[i] + M->NbRows, N->p[i], N->NbColumns);
 357     Matrix_Free(H);
 358     Matrix_Free(Q);
 359     Matrix_Free(U);
 360     return N;
 361 }
 362
 363 static void SwapColumns(Value **V, int n, int i, int j)
 364 {
 365     int r;
 366
 367     for (r = 0; r < n; ++r)
 368         value_swap(V[r][i], V[r][j]);
 369 }
 370
 371 /* C is assumed to be the "true" context, i.e., it has been intersected
 372  * with the projection of P onto the parameter space.
 373  * Furthermore, P and C are assumed to be full-dimensional.
 374  */
 375 Param_Polyhedron *TC_P2PP(Polyhedron *P, Polyhedron *C,
 376                           struct barvinok_options *options)
 377 {
 378     unsigned nparam = C->Dimension;
 379     unsigned nvar = P->Dimension - C->Dimension;
 380     int i, j;
 381     Matrix *H;
 382     Matrix *A;
 383     Matrix *K;
 384     int rows;
 385     Param_Polyhedron *PP;
 386     Matrix M;
 387
 388     assert(C->NbEq == 0);
 389
 390     Polyhedron_Matrix_View(P, &M, P->NbConstraints);
 391     H = standard_constraints(&M, nparam, &rows, NULL);
 392
 393     A = Matrix_Alloc(rows, nvar+rows);
 394     for (i = nvar; i < H->NbRows; ++i) {
 395         Vector_Oppose(H->p[i], A->p[i-nvar], H->NbColumns);
 396         value_set_si(A->p[i-nvar][i], 1);
 397     }
 398     for (i = 0, j = H->NbRows-nvar; i < nvar; ++i) {
 399         if (First_Non_Zero(H->p[i], i) == -1)
 400             continue;
 401         Vector_Oppose(H->p[i], A->p[j], H->NbColumns);
 402         value_set_si(A->p[j][j+nvar], 1);
 403         SwapColumns(A->p, A->NbRows, j+nvar, i);
 404         ++j;
 405     }
 406     K = null_space(A);
 407     Matrix_Free(A);
 408     /* Ignore extra constraints */
 409     K->NbRows = H->NbRows;
 410     Matrix_Free(H);
 411     PP = points2triangs(K, P, C, options);
 412     Matrix_Free(K);
 413     return PP;
 414 }