verif_ehrhart_e.c

   1 /*************************************************/
   2 /*     verif_ehrhart.c                           */
   3 /* program to compare effective number of points */
   4 /* in a polytope with the corresponding          */
   5 /* evaluation of the Ehrhart polynomial.         */
   6 /* Parameters vary in range -RANGE to RANGE      */
   7 /* (define below) by default.                    */
   8 /* Can be overridden by specifying               */
   9 /* -r<RANGE>, or -m<min> and -M<max>             */
  10 /*                                               */
  11 /* written by Vincent Loechner (c) 2000.         */
  12 /*  loechner@icps.u-strasbg.fr                   */
  13 /*************************************************/
  14
  15 #include <stdio.h>
  16 #include <string.h>
  17 #include <stdlib.h>
  18
  19 #include <polylib/polylibgmp.h>
  20 #include "ev_operations.h"
  21 #include "util.h"
  22 #include <barvinok.h>
  23 #define MAXRAYS 1024
  24
  25 /* define this to print all the results */
  26 /* else, only a progress bar is printed */
  27 /* #define PRINT_ALL_RESULTS
  28 /* define this to continue the test after first error found */
  29 /* #define DONT_BREAK_ON_ERROR */
  30
  31 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
  32 #define RANGE 50
  33
  34 /* SRANGE : small range for evalutations */
  35 #define SRANGE 15
  36
  37 /* if dimension >= BIDDIM, use SRANGE */
  38 #define BIGDIM 5
  39
  40 /* VSRANGE : very small range for evalutations */
  41 #define VSRANGE 5
  42
  43 /* if dimension >= VBIDDIM, use VSRANGE */
  44 #define VBIGDIM 8
  45
  46 Value min, max;
  47
  48 char **params;
  49
  50 #ifdef DONT_BREAK_ON_ERROR
  51 #define PRINT_ALL_RESULTS
  52 #endif
  53
  54 #ifndef PRINT_ALL_RESULTS
  55 int st;
  56 #endif
  57
  58 /****************************************************/
  59 /* function check_poly :                            */
  60 /* scans the parameter space from min to max (all   */
  61 /* directions). Computes the number of points in    */
  62 /* the polytope using both methods, and compare them*/
  63 /* returns 1 on success                             */
  64 /****************************************************/
  65
  66 int check_poly(Polyhedron *S, Polyhedron *C, evalue *EP,
  67                int exist, int nparam, int pos, Value *z) {
  68
  69   int k;
  70   Value c,tmp;
  71
  72   value_init(c); value_init(tmp);
  73
  74   if(pos == nparam) {
  75
  76     /* Computes the ehrhart polynomial */
  77     value_set_double(c, compute_evalue(EP,&z[S->Dimension-nparam+1])+.25);
  78     /* if c=0 we may be out of context. */
  79     /* scanning is useless in this case*/
  80     if(!in_domain(C,&z[S->Dimension-nparam+1])) {
  81
  82       /* ok */ ;
  83     }
  84     else {
  85
  86 #ifdef PRINT_ALL_RESULTS
  87       printf("EP( ");
  88       value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
  89       for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  90         printf(", ");
  91         value_print(stdout,VALUE_FMT,z[k]);
  92       }
  93       printf(" ) = ");
  94       value_print(stdout,VALUE_FMT,c);
  95       printf(" ");
  96 #endif
  97
  98       /* Manually count the number of points */
  99       count_points_e(1, S, exist, nparam, z, &tmp);
 100 #ifdef PRINT_ALL_RESULTS
 101         printf(", count = ");
 102         value_print(stdout, P_VALUE_FMT, tmp);
 103         printf(". ");
 104 #endif
 105
 106       if(value_ne(tmp,c)) {
 107         printf("\n");
 108         fflush(stdout);
 109         fprintf(stderr,"Error !\n");
 110         fprintf(stderr,"EP( ");
 111         value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
 112         for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 113           fprintf(stderr,", ");
 114           value_print(stderr,VALUE_FMT,z[k]);
 115         }
 116         fprintf(stderr," ) should be ");
 117         value_print(stderr,VALUE_FMT,tmp);
 118         fprintf(stderr,", while EP eval gives ");
 119         value_print(stderr,VALUE_FMT,c);
 120         fprintf(stderr,".\n");
 121         print_evalue(stderr, EP, params);
 122 #ifndef DONT_BREAK_ON_ERROR
 123         value_clear(c); value_clear(tmp);
 124         return(0);
 125 #endif
 126       }
 127
 128 #ifdef PRINT_ALL_RESULTS
 129       else
 130         printf("OK.\n");
 131 #endif
 132     }
 133   }
 134   else
 135     for(value_assign(tmp,min); value_le(tmp,max); value_increment(tmp,tmp)) {
 136
 137 #ifndef PRINT_ALL_RESULTS
 138       k = VALUE_TO_INT(tmp);
 139       if(!pos && !(k%st)) {
 140         printf("o");
 141         fflush(stdout);
 142       }
 143 #endif
 144
 145       value_assign(z[pos+S->Dimension-nparam+1],tmp);
 146       if(!check_poly(S, C, EP, exist, nparam, pos+1, z)) {
 147         value_clear(c); value_clear(tmp);
 148         return(0);
 149       }
 150     }
 151   value_clear(c); value_clear(tmp);
 152   return(1);
 153 } /* check_poly */
 154
 155 int main(int argc,char *argv[])
 156 {
 157   Matrix *C1, *P1;
 158   Polyhedron *C, *P, *S;
 159   Value *p, tmp;
 160   int i,j;
 161   int m,M;
 162   int exist, nparam;
 163   char s[128];
 164   evalue *EP;
 165   int res;
 166
 167 /******* Read the input *********/
 168   P1 = Matrix_Read();
 169
 170     fgets(s, 128, stdin);
 171     while ((*s=='#') || (sscanf(s, "E %d", &exist)<1))
 172         fgets(s, 128, stdin);
 173
 174     fgets(s, 128, stdin);
 175     while ((*s=='#') || (sscanf(s, "P %d", &nparam)<1))
 176         fgets(s, 128, stdin);
 177
 178   P = Constraints2Polyhedron(P1,MAXRAYS);
 179   params = Read_ParamNames(stdin, nparam);
 180   Matrix_Free(P1);
 181
 182   /******* Read the options: initialize min and max ********/
 183   if(P->Dimension >= VBIGDIM)
 184     M = VSRANGE;
 185   else if(P->Dimension >= BIGDIM)
 186     M = SRANGE;
 187   else
 188     M = RANGE;
 189   m = -M;
 190   if(argc != 1 ) {
 191     for(i=1;i<argc;i++) {
 192       if(!strncmp(argv[i],"-m",2)) {
 193
 194         /* min specified */
 195         m = atoi(&argv[i][2]);
 196       }
 197       else if(!strncmp(argv[i],"-M",2)) {
 198
 199         /* max specified */
 200         M = atoi(&argv[i][2]);
 201       }
 202       else if(!strncmp(argv[i], "-r", 2)) {
 203
 204         /* range specified */
 205         M = atoi(&argv[i][2]);
 206         m = -M;
 207       }
 208       else {
 209         fprintf(stderr,"Unknown option: %s\n",argv[i]);
 210         fprintf(stderr,"Usage: %s [-m<>][-M<>][-r<>]\n",argv[0]);
 211         return(-1);
 212       }
 213     }
 214   }
 215   if(m > M) {
 216     fprintf(stderr,"Nothing to do: min > max !\n");
 217     return(0);
 218   }
 219   value_init(min);
 220   value_init(max);
 221   value_set_si(min,m);
 222   value_set_si(max,M);
 223   value_init(tmp);
 224
 225   /******* Compute true context *******/
 226   C1 = Matrix_Alloc(nparam+1,P->Dimension+1);
 227
 228   for(i=0;i<C1->NbRows;i++)
 229     for(j=0;j<C1->NbColumns;j++)
 230       if(i==j-P->Dimension+nparam)
 231         value_set_si(C1->p[i][j],1);
 232       else
 233         value_set_si(C1->p[i][j],0);
 234   C = Polyhedron_Image(P,C1,MAXRAYS);
 235   Matrix_Free(C1);
 236
 237   /******* Compute EP *********/
 238   EP = barvinok_enumerate_e(P, exist, nparam, MAXRAYS);
 239
 240   /******* Initializations for check *********/
 241   p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
 242   for(i=0;i<=P->Dimension;i++) {
 243     value_init(p[i]);
 244     value_set_si(p[i],0);
 245   }
 246   value_init(p[i]);
 247   value_set_si(p[i],1);
 248
 249   /* S = scanning list of polyhedra */
 250   S = Polyhedron_Scan(P,C,MAXRAYS);
 251
 252 #ifndef PRINT_ALL_RESULTS
 253   if(C->Dimension > 0) {
 254     value_substract(tmp,max,min);
 255     if (VALUE_TO_INT(tmp) > 80)
 256       st = 1+(VALUE_TO_INT(tmp))/80;
 257     else
 258       st=1;
 259     for(i=VALUE_TO_INT(min);i<=VALUE_TO_INT(max);i+=st)
 260       printf(".");
 261     printf( "\r" );
 262     fflush(stdout);
 263   }
 264 #endif
 265
 266   /******* CHECK NOW *********/
 267   res = 0;
 268   if(S && !check_poly(S, C, EP, exist, nparam, 0, p)) {
 269     fprintf(stderr,"Check failed !\n");
 270     res = -1;
 271   }
 272
 273 #ifndef PRINT_ALL_RESULTS
 274   printf( "\n" );
 275 #endif
 276
 277   for(i=0;i<=(P->Dimension+1);i++)
 278     value_clear(p[i]);
 279   free(p);
 280   value_clear(tmp);
 281   Free_ParamNames(params, C->Dimension);
 282   Polyhedron_Free(S);
 283   Polyhedron_Free(C);
 284   Polyhedron_Free(P);
 285   free_evalue_refs(EP);
 286   free(EP);
 287   return res;
 288 } /* main */
 289
 290
 291
 292