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