?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
move Param_Polyhedron_Scale_Integer to scale.c and rename
[barvinok.git] / scale.c
blobca9741d7ada004e48e916c7abd2f16917485a30a
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 Matrix *expansion;
234
235 value_set_si(*det, 1);
236 if (!PP->nbV)
237 return;
238
239 nb_param = PP->D->Domain->Dimension;
240 nb_vars = PP->V->Vertex->NbRows;
241
242 /* Scan the vertices and make an orthogonal expansion of the variable
243 space */
244 /* a- prepare the array of common denominators */
245 denoms = Vector_Alloc(nb_vars);
246 value_init(global_var_lcm);
247
248 /* b- scan the vertices and compute the variables' global lcms */
249 for (V = PP->V; V; V = V->next)
250 for (i = 0; i < nb_vars; i++)
251 Lcm3(denoms->p[i], V->Vertex->p[i][nb_param+1], &denoms->p[i]);
252
253 value_set_si(global_var_lcm, 1);
254 for (i = 0; i < nb_vars; i++) {
255 value_multiply(*det, *det, denoms->p[i]);
256 Lcm3(global_var_lcm, denoms->p[i], &global_var_lcm);
257 }
258
259 /* scale vertices */
260 for (V = PP->V; V; V = V->next)
261 for (i = 0; i < nb_vars; i++) {
262 Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], denoms->p[i], nb_param+1);
263 Vector_Normalize(V->Vertex->p[i], nb_param+2);
264 }
265
266 /* the expansion can be actually writen as global_var_lcm.L^{-1} */
267 /* this is equivalent to multiply the rows of P by denoms_det */
268 for (i = 0; i < nb_vars; i++)
269 value_division(denoms->p[i], global_var_lcm, denoms->p[i]);
270
271 /* OPT : we could use a vector instead of a diagonal matrix here (c- and d-).*/
272 /* c- make the quick expansion matrix */
273 expansion = Matrix_Alloc(nb_vars+nb_param+1, nb_vars+nb_param+1);
274 for (i = 0; i < nb_vars; i++)
275 value_assign(expansion->p[i][i], denoms->p[i]);
276 for (i = nb_vars; i < nb_vars+nb_param+1; i++)
277 value_assign(expansion->p[i][i], global_var_lcm);
278
279 /* d- apply the variable expansion to the polyhedron */
280 if (P)
281 *P = Polyhedron_Preimage(*P, expansion, MaxRays);
282
283 Matrix_Free(expansion);
284 value_clear(global_var_lcm);
285 Vector_Free(denoms);
286 }
287
288 /* adapted from mpolyhedron_inflate in PolyLib */
289 static Polyhedron *Polyhedron_Inflate(Polyhedron *P, unsigned nparam,
290 unsigned MaxRays)
291 {
292 Value sum;
293 int nvar = P->Dimension - nparam;
294 Matrix *C = Polyhedron2Constraints(P);
295 Polyhedron *P2;
296 int i, j;
297
298 value_init(sum);
299 /* subtract the sum of the negative coefficients of each inequality */
300 for (i = 0; i < C->NbRows; ++i) {
301 value_set_si(sum, 0);
302 for (j = 0; j < nvar; ++j)
303 if (value_neg_p(C->p[i][1+j]))
304 value_addto(sum, sum, C->p[i][1+j]);
305 value_subtract(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], sum);
306 }
307 value_clear(sum);
308 P2 = Constraints2Polyhedron(C, MaxRays);
309 Matrix_Free(C);
310 return P2;
311 }
312
313 /* adapted from mpolyhedron_deflate in PolyLib */
314 static Polyhedron *Polyhedron_Deflate(Polyhedron *P, unsigned nparam,
315 unsigned MaxRays)
316 {
317 Value sum;
318 int nvar = P->Dimension - nparam;
319 Matrix *C = Polyhedron2Constraints(P);
320 Polyhedron *P2;
321 int i, j;
322
323 value_init(sum);
324 /* subtract the sum of the positive coefficients of each inequality */
325 for (i = 0; i < C->NbRows; ++i) {
326 value_set_si(sum, 0);
327 for (j = 0; j < nvar; ++j)
328 if (value_pos_p(C->p[i][1+j]))
329 value_addto(sum, sum, C->p[i][1+j]);
330 value_subtract(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], sum);
331 }
332 value_clear(sum);
333 P2 = Constraints2Polyhedron(C, MaxRays);
334 Matrix_Free(C);
335 return P2;
336 }
337
338 Polyhedron *scale_init(Polyhedron *P, Polyhedron *C, struct scale_data *scaling,
339 struct barvinok_options *options)
340 {
341 unsigned nparam = C->Dimension;
342
343 value_init(scaling->det);
344 value_set_si(scaling->det, 1);
345 scaling->save_approximation = options->polynomial_approximation;
346
347 if (options->polynomial_approximation == BV_APPROX_SIGN_NONE ||
348 options->polynomial_approximation == BV_APPROX_SIGN_APPROX)
349 return P;
350
351 if (options->polynomial_approximation == BV_APPROX_SIGN_UPPER)
352 P = Polyhedron_Inflate(P, nparam, options->MaxRays);
353 if (options->polynomial_approximation == BV_APPROX_SIGN_LOWER)
354 P = Polyhedron_Deflate(P, nparam, options->MaxRays);
355
356 /* Don't deflate/inflate again (on this polytope) */
357 options->polynomial_approximation = BV_APPROX_SIGN_NONE;
358
359 return P;
360 }
361
362 Polyhedron *scale(Param_Polyhedron *PP, Polyhedron *P,
363 struct scale_data *scaling, int free_P,
364 struct barvinok_options *options)
365 {
366 Polyhedron *T = P;
367
368 if (options->scale_flags & BV_APPROX_SCALE_FAST)
369 Param_Polyhedron_Scale_Integer_Fast(PP, &T, &scaling->det, options->MaxRays);
370 else
371 Param_Polyhedron_Scale_Integer_Slow(PP, &T, &scaling->det, options->MaxRays);
372 if (free_P)
373 Polyhedron_Free(P);
374 return T;
375 }
376
377 void scale_finish(evalue *e, struct scale_data *scaling,
378 struct barvinok_options *options)
379 {
380 if (value_notone_p(scaling->det))
381 evalue_div(e, scaling->det);
382 value_clear(scaling->det);
383 /* reset options that may have been changed */
384 options->polynomial_approximation = scaling->save_approximation;
385 }
lattice point enumeration library