piputil.c

   1 #include <polylib/polylibgmp.h>
   2 #include <stdio.h>  /* needed for piplib ! */
   3 #include <piplib/piplibMP.h>
   4 #include <assert.h>
   5 #include "config.h"
   6
   7 #ifdef HAVE_GROWING_CHERNIKOVA
   8 #define MAXRAYS    POL_NO_DUAL
   9 #else
  10 #define MAXRAYS  600
  11 #endif
  12
  13 static PipMatrix *poly2pip(Polyhedron *P, int pos, int n, int nparam)
  14 {
  15     PipMatrix *     matrix;
  16     int i, j;
  17     unsigned int extra = P->Dimension - pos - n - nparam;
  18
  19     matrix = pip_matrix_alloc(P->NbConstraints, P->Dimension+2);
  20     for (i = 0; i < P->NbConstraints; ++i) {
  21         value_assign(matrix->p[i][0], P->Constraint[i][0]);
  22         for (j = 0; j < n; ++j)
  23             value_assign(matrix->p[i][1+j], P->Constraint[i][1+pos+j]);
  24         for (j = 0; j < pos; ++j)
  25             value_assign(matrix->p[i][1+n+j], P->Constraint[i][1+j]);
  26         for (j = 0; j < extra; ++j)
  27             value_assign(matrix->p[i][1+n+pos+j], P->Constraint[i][1+n+pos+j]);
  28         for (j = 1+pos+n+extra; j < P->Dimension+2; ++j)
  29             value_assign(matrix->p[i][j], P->Constraint[i][j]);
  30     }
  31
  32     return matrix;
  33 }
  34
  35 static int max_new(PipQuast *q, int max, int d, int *maxd)
  36 {
  37     PipNewparm *p;
  38
  39     for (p = q->newparm; p; p = p->next)
  40         if (p->rank > max)
  41             max = p->rank;
  42     if (q->condition) {
  43         if (++d > *maxd)
  44             *maxd = d;
  45         max = max_new(q->next_else, max, d, maxd);
  46         max = max_new(q->next_then, max, d, maxd);
  47     }
  48     return max;
  49 }
  50
  51 static int vectorpos2cpos(int v, int pos, int n, int e, int d, int j)
  52 {
  53     int np = d-v-e;
  54     int l;
  55
  56     if (j < pos)
  57         l = 1+j;
  58     else if (j < v-n)
  59         l = 1+n+j;
  60     else if (j < v-n+np)
  61         l = 1+n+e+j;
  62     else
  63         l = 1+n+j-np;
  64
  65     return l;
  66 }
  67
  68 /*
  69  * i: next row
  70  * v: total number of variables in original problem
  71  * pos: position of minimized variables
  72  * n: dimension of space in which minimum was taken
  73  *    i.e., number of minimized vars in original problem
  74  * e: number of extra existential variables
  75  * d: total dimension
  76  *
  77  * Dimensions in the resulting domain are:
  78  *      v n e (d-v-n-e)
  79  */
  80 static void add_quast(Polyhedron**D, Matrix* M, PipQuast *q,
  81                       int i, int v, int pos, int n, int e, int d)
  82 {
  83     PipNewparm *p;
  84
  85     for (p = q->newparm; p; p = p->next) {
  86         int j, l;
  87         PipVector *vec = p->vector;
  88
  89         value_set_si(M->p[i][0], 1);
  90         Vector_Set(M->p[i]+1, 0, d);
  91
  92         for (j = 0; j < vec->nb_elements-1; ++j) {
  93             l = vectorpos2cpos(v, pos, n, e, d, j);
  94             value_assign(M->p[i][l], vec->the_vector[j]);
  95         }
  96         l = vectorpos2cpos(v, pos, n, e, d, p->rank);
  97         value_oppose(M->p[i][l], p->deno);
  98         value_assign(M->p[i][1+d], vec->the_vector[j]);
  99         ++i;
 100
 101         value_set_si(M->p[i][0], 1);
 102         Vector_Set(M->p[i]+1, 0, d);
 103
 104         for (j = 0; j < vec->nb_elements-1; ++j) {
 105             l = vectorpos2cpos(v, pos, n, e, d, j);
 106             value_oppose(M->p[i][l], vec->the_vector[j]);
 107         }
 108         l = vectorpos2cpos(v, pos, n, e, d, p->rank);
 109         value_assign(M->p[i][l], p->deno);
 110         value_oppose(M->p[i][1+d], vec->the_vector[j]);
 111         value_addto(M->p[i][1+d], M->p[i][1+d], p->deno);
 112         value_decrement(M->p[i][1+d], M->p[i][1+d]);
 113         ++i;
 114     }
 115
 116     if (q->condition) {
 117         int j;
 118         value_set_si(M->p[i][0], 1);
 119         Vector_Set(M->p[i]+1, 0, d);
 120
 121         for (j = 0; j < q->condition->nb_elements-1; ++j) {
 122             int l = vectorpos2cpos(v, pos, n, e, d, j);
 123             value_assign(M->p[i][l], q->condition->the_vector[j]);
 124         }
 125         value_assign(M->p[i][1+d], q->condition->the_vector[j]);
 126         add_quast(D, M, q->next_then, i+1, v, pos, n, e, d);
 127
 128         for (j = 0; j < q->condition->nb_elements-1; ++j) {
 129             int l = vectorpos2cpos(v, pos, n, e, d, j);
 130             value_oppose(M->p[i][l], q->condition->the_vector[j]);
 131         }
 132         value_oppose(M->p[i][1+d], q->condition->the_vector[j]);
 133         value_decrement(M->p[i][1+d], M->p[i][1+d]);
 134         add_quast(D, M, q->next_else, i+1, v, pos, n, e, d);
 135     } else if (q->list) {
 136         PipList *l;
 137         Matrix *C;
 138         Polyhedron *P;
 139         int j, k;
 140         for (j = 0, l = q->list; l; ++j, l = l->next) {
 141             Vector_Set(M->p[i], 0, d+1);
 142             value_set_si(M->p[i][1+pos+j], -1);
 143
 144             for (k = 0; k < l->vector->nb_elements-1; ++k) {
 145                 int ll = vectorpos2cpos(v, pos, n, e, d, k);
 146                 value_assign(M->p[i][ll], l->vector->the_vector[k]);
 147             }
 148             value_assign(M->p[i][1+d], l->vector->the_vector[k]);
 149
 150             ++i;
 151         }
 152         C = Matrix_Alloc(i, 1+d+1);
 153         Vector_Copy(M->p[0], C->p[0], i * (1+d+1));
 154         P = Constraints2Polyhedron(C, MAXRAYS);
 155         Matrix_Free(C);
 156         P->next = *D;
 157         *D = P;
 158     }
 159 }
 160
 161 /*
 162  * nvar: total number of variables, including the minimized variables,
 163  *       but excluding the parameters
 164  * pos: position of the first minimized variable
 165  * n: number of minimized variables
 166  */
 167 static Polyhedron *quast2poly(PipQuast *q, int nvar, int pos, int n)
 168 {
 169     int                 nparam;
 170     int                 nexist;
 171     PipList*            l;
 172     int                 i, d, rows;
 173     Matrix*             M;
 174     Polyhedron*         D;
 175
 176     assert(n != 0);     /* required ? */
 177     nparam = q->newparm ? q->newparm->rank :
 178              q->condition ? q->condition->nb_elements-1 :
 179                             q->list->vector->nb_elements-1;
 180     d = 0;
 181     nexist = max_new(q, nparam-1, 0, &d) - nparam+1;
 182     rows = nparam + 2 * nexist + d + n;
 183
 184     /* nparam now refers to the number of parameters in the original polyhedron */
 185     nparam -= nvar - n;
 186     M = Matrix_Alloc(rows, 1+nvar+nexist+nparam+1);
 187     /* All parameters are/should be positive */
 188     for (i = 0; i < pos; ++i) {
 189         value_set_si(M->p[i][0], 1);
 190         value_set_si(M->p[i][1+i], 1);
 191     }
 192     for (i = pos+n; i < nvar; ++i) {
 193         value_set_si(M->p[i-n][0], 1);
 194         value_set_si(M->p[i-n][1+i], 1);
 195     }
 196     for (i = 0; i < nparam; ++i) {
 197         value_set_si(M->p[nvar-n+i][0], 1);
 198         value_set_si(M->p[nvar-n+i][1+nvar+nexist+i], 1);
 199     }
 200     D = 0;
 201     add_quast(&D, M, q, nparam, nvar, pos, n, nexist, nvar+nexist+nparam);
 202     Matrix_Free(M);
 203     return D;
 204 }
 205
 206 Polyhedron *pip_lexmin(Polyhedron *P, int pos, int n, int nparam)
 207 {
 208     PipOptions  *options;
 209     PipMatrix * domain;
 210     unsigned int nvar = P->Dimension - nparam;
 211     PipMatrix   *context = pip_matrix_alloc(0, P->Dimension - n + 2);
 212     PipQuast    *sol;
 213     Polyhedron  *min;
 214
 215     POL_ENSURE_INEQUALITIES(P);
 216
 217     domain = poly2pip(P, pos, n, nparam);
 218     /*
 219     pip_matrix_print(stderr, domain);
 220     pip_matrix_print(stderr, context);
 221     */
 222
 223     options = pip_options_init();
 224     sol = pip_solve(domain, context, -1, options);
 225
 226     /*
 227     pip_quast_print(stderr, sol, 0);
 228     */
 229
 230     min = quast2poly(sol, nvar, pos, n);
 231
 232     pip_quast_free(sol);
 233     pip_matrix_free(context);
 234     pip_matrix_free(domain);
 235     pip_options_free(options);
 236
 237     return min;
 238 }