topcom.c

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