scale.c

   1 #include <barvinok/util.h>
   2 #include <barvinok/options.h>
   3 #include "scale.h"
   4
   5 /* If a vertex is described by A x + B p + c = 0, then
   6  * M = [A B] and we want to compute a linear transformation L such
   7  * that H L = A and H \Z contains both A \Z and B \Z.
   8  * We compute
   9  *             [ A B ] = [ H 0 ] [ U_11  U_12 ]
  10  *                               [ U_21  U_22 ]
  11  *
  12  * U_11 is the required linear transformation.
  13  * Note that this also works if M has more rows than there are variables,
  14  * i.e., if some rows in M are linear combinations of other rows.
  15  * These extra rows only affect and H and not U.
  16  */
  17 static Lattice *extract_lattice(Matrix *M, unsigned nvar)
  18 {
  19     int row;
  20     Matrix *H, *Q, *U, *Li;
  21     Lattice *L;
  22     int ok;
  23
  24     left_hermite(M, &H, &Q, &U);
  25     Matrix_Free(U);
  26
  27     Li = Matrix_Alloc(nvar+1, nvar+1);
  28     L = Matrix_Alloc(nvar+1, nvar+1);
  29     value_set_si(Li->p[nvar][nvar], 1);
  30
  31     for (row = 0; row < nvar; ++row)
  32         Vector_Copy(Q->p[row], Li->p[row], nvar);
  33     Matrix_Free(H);
  34     Matrix_Free(Q);
  35
  36     ok = Matrix_Inverse(Li, L);
  37     assert(ok);
  38     Matrix_Free(Li);
  39
  40     return L;
  41 }
  42
  43 /* Returns the smallest (wrt inclusion) lattice that contains both X and Y */
  44 static Lattice *LatticeJoin(Lattice *X, Lattice *Y)
  45 {
  46     int i;
  47     int dim = X->NbRows-1;
  48     Value lcm;
  49     Value tmp;
  50     Lattice *L;
  51     Matrix *M, *H, *U, *Q;
  52
  53     assert(X->NbColumns-1 == dim);
  54     assert(Y->NbRows-1 == dim);
  55     assert(Y->NbColumns-1 == dim);
  56
  57     value_init(lcm);
  58     value_init(tmp);
  59
  60     M = Matrix_Alloc(dim, 2*dim);
  61     value_lcm(X->p[dim][dim], Y->p[dim][dim], &lcm);
  62
  63     value_division(tmp, lcm, X->p[dim][dim]);
  64     for (i = 0; i < dim; ++i)
  65         Vector_Scale(X->p[i], M->p[i], tmp, dim);
  66     value_division(tmp, lcm, Y->p[dim][dim]);
  67     for (i = 0; i < dim; ++i)
  68         Vector_Scale(Y->p[i], M->p[i]+dim, tmp, dim);
  69
  70     left_hermite(M, &H, &Q, &U);
  71     Matrix_Free(M);
  72     Matrix_Free(Q);
  73     Matrix_Free(U);
  74
  75     L = Matrix_Alloc(dim+1, dim+1);
  76     value_assign(L->p[dim][dim], lcm);
  77     for (i = 0; i < dim; ++i)
  78         Vector_Copy(H->p[i], L->p[i], dim);
  79     Matrix_Free(H);
  80
  81     value_clear(tmp);
  82     value_clear(lcm);
  83     return L;
  84 }
  85
  86 static void Param_Vertex_Common_Denominator(Param_Vertices *V)
  87 {
  88     unsigned dim;
  89     Value lcm;
  90     int i;
  91
  92     assert(V->Vertex->NbRows > 0);
  93     dim = V->Vertex->NbColumns-2;
  94
  95     value_init(lcm);
  96
  97     value_assign(lcm, V->Vertex->p[0][dim+1]);
  98     for (i = 1; i < V->Vertex->NbRows; ++i)
  99         value_lcm(V->Vertex->p[i][dim+1], lcm, &lcm);
 100
 101     for (i = 0; i < V->Vertex->NbRows; ++i) {
 102         if (value_eq(V->Vertex->p[i][dim+1], lcm))
 103             continue;
 104         value_division(V->Vertex->p[i][dim+1], lcm, V->Vertex->p[i][dim+1]);
 105         Vector_Scale(V->Vertex->p[i], V->Vertex->p[i],
 106                      V->Vertex->p[i][dim+1], dim+1);
 107         value_assign(V->Vertex->p[i][dim+1], lcm);
 108     }
 109
 110     value_clear(lcm);
 111 }
 112
 113 static void Param_Vertex_Image(Param_Vertices *V, Matrix *T)
 114 {
 115     unsigned nvar  = V->Vertex->NbRows;
 116     unsigned nparam = V->Vertex->NbColumns - 2;
 117     Matrix *Vertex;
 118     int i;
 119
 120     Param_Vertex_Common_Denominator(V);
 121     Vertex = Matrix_Alloc(V->Vertex->NbRows, V->Vertex->NbColumns);
 122     Matrix_Product(T, V->Vertex, Vertex);
 123     for (i = 0; i < nvar; ++i) {
 124         value_assign(Vertex->p[i][nparam+1], V->Vertex->p[i][nparam+1]);
 125         Vector_Normalize(Vertex->p[i], nparam+2);
 126     }
 127     Matrix_Free(V->Vertex);
 128     V->Vertex = Vertex;
 129 }
 130
 131 /* Scales the parametric polyhedron with constraints *P and vertices PP
 132  * such that the number of integer points can be represented by a polynomial.
 133  * Both *P and P->Vertex are adapted according to the scaling.
 134  * The scaling factor is returned in *det.
 135  * The enumerator of the scaled parametric polyhedron should be divided
 136  * by this number to obtain an approximation of the enumerator of the
 137  * original parametric polyhedron.
 138  *
 139  * The algorithm is described in "Approximating Ehrhart Polynomials using
 140  * affine transformations" by B. Meister.
 141  */
 142 void Param_Polyhedron_Scale_Integer_Slow(Param_Polyhedron *PP, Polyhedron **P,
 143                                          Value *det, unsigned MaxRays)
 144 {
 145     Param_Vertices *V;
 146     unsigned dim = (*P)->Dimension;
 147     unsigned nparam;
 148     unsigned nvar;
 149     Lattice *L = NULL, *Li;
 150     Matrix *T;
 151     Matrix *expansion;
 152     int i;
 153     int ok;
 154
 155     if (!PP->nbV)
 156         return;
 157
 158     nparam = PP->V->Vertex->NbColumns - 2;
 159     nvar = dim - nparam;
 160
 161     for (V = PP->V; V; V = V->next) {
 162         Lattice *L2;
 163         Matrix *M;
 164         int i, j, n;
 165         unsigned char *supporting;
 166
 167         supporting = supporting_constraints(*P, V, &n);
 168         M = Matrix_Alloc(n, (*P)->Dimension);
 169         for (i = 0, j = 0; i < (*P)->NbConstraints; ++i)
 170             if (supporting[i])
 171                 Vector_Copy((*P)->Constraint[i]+1, M->p[j++], (*P)->Dimension);
 172         free(supporting);
 173         L2 = extract_lattice(M, nvar);
 174         Matrix_Free(M);
 175
 176         if (!L)
 177             L = L2;
 178         else {
 179             Lattice *L3 = LatticeJoin(L, L2);
 180             Matrix_Free(L);
 181             Matrix_Free(L2);
 182             L = L3;
 183         }
 184     }
 185
 186     /* apply the variable expansion to the polyhedron (constraints) */
 187     expansion = Matrix_Alloc(nvar + nparam + 1,  nvar + nparam + 1);
 188     for (i = 0; i < nvar; ++i)
 189         Vector_Copy(L->p[i], expansion->p[i], nvar);
 190     for (i = nvar; i < nvar+nparam+1; ++i)
 191         value_assign(expansion->p[i][i], L->p[nvar][nvar]);
 192
 193     *P = Polyhedron_Preimage(*P, expansion, MaxRays);
 194     Matrix_Free(expansion);
 195
 196     /* apply the variable expansion to the parametric vertices */
 197     Li = Matrix_Alloc(nvar+1, nvar+1);
 198     ok = Matrix_Inverse(L, Li);
 199     assert(ok);
 200     Matrix_Free(L);
 201     assert(value_one_p(Li->p[nvar][nvar]));
 202     T = Matrix_Alloc(nvar, nvar);
 203     value_set_si(*det, 1);
 204     for (i = 0; i < nvar; ++i) {
 205         value_multiply(*det, *det, Li->p[i][i]);
 206         Vector_Copy(Li->p[i], T->p[i], nvar);
 207     }
 208     Matrix_Free(Li);
 209     for (V = PP->V; V; V = V->next)
 210         Param_Vertex_Image(V, T);
 211     Matrix_Free(T);
 212 }
 213
 214 /* Scales the parametric polyhedron with constraints *P and vertices PP
 215  * such that the number of integer points can be represented by a polynomial.
 216  * Both *P and P->Vertex are adapted according to the scaling.
 217  * The scaling factor is returned in *det.
 218  * The enumerator of the scaled parametric polyhedron should be divided
 219  * by this number to obtain an approximation of the enumerator of the
 220  * original parametric polyhedron.
 221  *
 222  * The algorithm is described in "Approximating Ehrhart Polynomials using
 223  * affine transformations" by B. Meister.
 224  */
 225 void Param_Polyhedron_Scale_Integer_Fast(Param_Polyhedron *PP, Polyhedron **P,
 226                                          Value *det, unsigned MaxRays)
 227 {
 228   int i;
 229   int nb_param, nb_vars;
 230   Vector *denoms;
 231   Param_Vertices *V;
 232   Value global_var_lcm;
 233   Value tmp;
 234   Matrix *expansion;
 235
 236   value_set_si(*det, 1);
 237   if (!PP->nbV)
 238     return;
 239
 240   nb_param = PP->D->Domain->Dimension;
 241   nb_vars = PP->V->Vertex->NbRows;
 242
 243   /* Scan the vertices and make an orthogonal expansion of the variable
 244      space */
 245   /* a- prepare the array of common denominators */
 246   denoms = Vector_Alloc(nb_vars);
 247   value_init(global_var_lcm);
 248
 249   value_init(tmp);
 250   /* b- scan the vertices and compute the variables' global lcms */
 251   for (V = PP->V; V; V = V->next) {
 252     for (i = 0; i < nb_vars; i++) {
 253       Vector_Gcd(V->Vertex->p[i], nb_param, &tmp);
 254       Gcd(tmp, V->Vertex->p[i][nb_param+1], &tmp);
 255       value_division(tmp, V->Vertex->p[i][nb_param+1], tmp);
 256       Lcm3(denoms->p[i], tmp, &denoms->p[i]);
 257     }
 258   }
 259   value_clear(tmp);
 260
 261   value_set_si(global_var_lcm, 1);
 262   for (i = 0; i < nb_vars; i++) {
 263     value_multiply(*det, *det, denoms->p[i]);
 264     Lcm3(global_var_lcm, denoms->p[i], &global_var_lcm);
 265   }
 266
 267   /* scale vertices */
 268   for (V = PP->V; V; V = V->next)
 269     for (i = 0; i < nb_vars; i++) {
 270       Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], denoms->p[i], nb_param+1);
 271       Vector_Normalize(V->Vertex->p[i], nb_param+2);
 272     }
 273
 274   /* the expansion can be actually writen as global_var_lcm.L^{-1} */
 275   /* this is equivalent to multiply the rows of P by denoms_det */
 276   for (i = 0; i < nb_vars; i++)
 277     value_division(denoms->p[i], global_var_lcm, denoms->p[i]);
 278
 279   /* OPT : we could use a vector instead of a diagonal matrix here (c- and d-).*/
 280   /* c- make the quick expansion matrix */
 281   expansion = Matrix_Alloc(nb_vars+nb_param+1, nb_vars+nb_param+1);
 282   for (i = 0; i < nb_vars; i++)
 283     value_assign(expansion->p[i][i], denoms->p[i]);
 284   for (i = nb_vars; i < nb_vars+nb_param+1; i++)
 285     value_assign(expansion->p[i][i], global_var_lcm);
 286
 287   /* d- apply the variable expansion to the polyhedron */
 288   if (P)
 289     *P = Polyhedron_Preimage(*P, expansion, MaxRays);
 290
 291   Matrix_Free(expansion);
 292   value_clear(global_var_lcm);
 293   Vector_Free(denoms);
 294 }
 295
 296 /* adapted from mpolyhedron_inflate in PolyLib */
 297 static Polyhedron *Polyhedron_Inflate(Polyhedron *P, unsigned nparam,
 298                                       unsigned MaxRays)
 299 {
 300     Value sum;
 301     int nvar = P->Dimension - nparam;
 302     Matrix *C = Polyhedron2Constraints(P);
 303     Polyhedron *P2;
 304     int i, j;
 305
 306     value_init(sum);
 307     /* subtract the sum of the negative coefficients of each inequality */
 308     for (i = 0; i < C->NbRows; ++i) {
 309         value_set_si(sum, 0);
 310         for (j = 0; j < nvar; ++j)
 311             if (value_neg_p(C->p[i][1+j]))
 312                 value_addto(sum, sum, C->p[i][1+j]);
 313         value_subtract(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], sum);
 314     }
 315     value_clear(sum);
 316     P2 = Constraints2Polyhedron(C, MaxRays);
 317     Matrix_Free(C);
 318     return P2;
 319 }
 320
 321 /* adapted from mpolyhedron_deflate in PolyLib */
 322 static Polyhedron *Polyhedron_Deflate(Polyhedron *P, unsigned nparam,
 323                                       unsigned MaxRays)
 324 {
 325     Value sum;
 326     int nvar = P->Dimension - nparam;
 327     Matrix *C = Polyhedron2Constraints(P);
 328     Polyhedron *P2;
 329     int i, j;
 330
 331     value_init(sum);
 332     /* subtract the sum of the positive coefficients of each inequality */
 333     for (i = 0; i < C->NbRows; ++i) {
 334         value_set_si(sum, 0);
 335         for (j = 0; j < nvar; ++j)
 336             if (value_pos_p(C->p[i][1+j]))
 337                 value_addto(sum, sum, C->p[i][1+j]);
 338         value_subtract(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], sum);
 339     }
 340     value_clear(sum);
 341     P2 = Constraints2Polyhedron(C, MaxRays);
 342     Matrix_Free(C);
 343     return P2;
 344 }
 345
 346 Polyhedron *scale_init(Polyhedron *P, Polyhedron *C, struct scale_data *scaling,
 347                        struct barvinok_options *options)
 348 {
 349     unsigned nparam = C->Dimension;
 350
 351     value_init(scaling->det);
 352     value_set_si(scaling->det, 1);
 353     scaling->save_approximation = options->polynomial_approximation;
 354
 355     if (options->polynomial_approximation == BV_APPROX_SIGN_NONE ||
 356         options->polynomial_approximation == BV_APPROX_SIGN_APPROX)
 357         return P;
 358
 359     if (options->polynomial_approximation == BV_APPROX_SIGN_UPPER)
 360         P = Polyhedron_Inflate(P, nparam, options->MaxRays);
 361     if (options->polynomial_approximation == BV_APPROX_SIGN_LOWER)
 362         P = Polyhedron_Deflate(P, nparam, options->MaxRays);
 363
 364     /* Don't deflate/inflate again (on this polytope) */
 365     options->polynomial_approximation = BV_APPROX_SIGN_NONE;
 366
 367     return P;
 368 }
 369
 370 Polyhedron *scale(Param_Polyhedron *PP, Polyhedron *P,
 371                   struct scale_data *scaling, int free_P,
 372                   struct barvinok_options *options)
 373 {
 374     Polyhedron *T = P;
 375     unsigned MaxRays;
 376
 377     MaxRays = options->MaxRays;
 378     POL_UNSET(options->MaxRays, POL_INTEGER);
 379
 380     if (options->scale_flags & BV_APPROX_SCALE_FAST)
 381         Param_Polyhedron_Scale_Integer_Fast(PP, &T, &scaling->det, options->MaxRays);
 382     else
 383         Param_Polyhedron_Scale_Integer_Slow(PP, &T, &scaling->det, options->MaxRays);
 384     if (free_P)
 385         Polyhedron_Free(P);
 386
 387     options->MaxRays = MaxRays;
 388
 389     return T;
 390 }
 391
 392 void scale_finish(evalue *e, struct scale_data *scaling,
 393                   struct barvinok_options *options)
 394 {
 395     if (value_notone_p(scaling->det))
 396         evalue_div(e, scaling->det);
 397     value_clear(scaling->det);
 398     /* reset options that may have been changed */
 399     options->polynomial_approximation = scaling->save_approximation;
 400 }