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, Enumeration *en,
  67                int exist, int nparam, int pos, Value *z) {
  68
  69   int k;
  70   Value c,tmp,*ctmp;
  71
  72   value_init(c); value_init(tmp);
  73
  74   if(pos == nparam) {
  75
  76     /* Computes the ehrhart polynomial */
  77     value_assign(c,*(ctmp=compute_poly(en,&z[S->Dimension-nparam+1])));
  78     free(ctmp);
  79     /* if c=0 we may be out of context. */
  80     /* scanning is useless in this case*/
  81     if(!in_domain(C,&z[S->Dimension-nparam+1])) {
  82
  83       /* ok */ ;
  84     }
  85     else {
  86
  87 #ifdef PRINT_ALL_RESULTS
  88       printf("EP( ");
  89       value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
  90       for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  91         printf(", ");
  92         value_print(stdout,VALUE_FMT,z[k]);
  93       }
  94       printf(" ) = ");
  95       value_print(stdout,VALUE_FMT,c);
  96       printf(" ");
  97 #endif
  98
  99       /* Manually count the number of points */
 100       count_points_e(1, S, exist, nparam, z, &tmp);
 101 #ifdef PRINT_ALL_RESULTS
 102         printf(", count = ");
 103         value_print(stdout, P_VALUE_FMT, tmp);
 104         printf(". ");
 105 #endif
 106
 107       if(value_ne(tmp,c)) {
 108         printf("\n");
 109         fflush(stdout);
 110         fprintf(stderr,"Error !\n");
 111         fprintf(stderr,"EP( ");
 112         value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
 113         for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 114           fprintf(stderr,", ");
 115           value_print(stderr,VALUE_FMT,z[k]);
 116         }
 117         fprintf(stderr," ) should be ");
 118         value_print(stderr,VALUE_FMT,tmp);
 119         fprintf(stderr,", while EP eval gives ");
 120         value_print(stderr,VALUE_FMT,c);
 121         fprintf(stderr,".\n");
 122         {
 123                  Enumeration *ee;
 124                  Enumeration_Print(stderr, en, params);
 125                  ee = en;
 126                  while (ee) {
 127                 if (in_domain(ee->ValidityDomain,&z[S->Dimension-nparam+1])) {
 128                          Print_Domain(stderr, ee->ValidityDomain, params);
 129                          print_evalue(stderr, &ee->EP, params);
 130                 }
 131                 ee = ee->next;
 132           }
 133         }
 134 #ifndef DONT_BREAK_ON_ERROR
 135         value_clear(c); value_clear(tmp);
 136         return(0);
 137 #endif
 138       }
 139
 140 #ifdef PRINT_ALL_RESULTS
 141       else
 142         printf("OK.\n");
 143 #endif
 144     }
 145   }
 146   else
 147     for(value_assign(tmp,min); value_le(tmp,max); value_increment(tmp,tmp)) {
 148
 149 #ifndef PRINT_ALL_RESULTS
 150       k = VALUE_TO_INT(tmp);
 151       if(!pos && !(k%st)) {
 152         printf("o");
 153         fflush(stdout);
 154       }
 155 #endif
 156
 157       value_assign(z[pos+S->Dimension-nparam+1],tmp);
 158       if(!check_poly(S, C, en, exist, nparam, pos+1, z)) {
 159         value_clear(c); value_clear(tmp);
 160         return(0);
 161       }
 162     }
 163   value_clear(c); value_clear(tmp);
 164   return(1);
 165 } /* check_poly */
 166
 167 int main(int argc,char *argv[])
 168 {
 169   Matrix *C1, *P1;
 170   Polyhedron *C, *P, *S;
 171   Value *p, tmp;
 172   int i,j;
 173   int m,M;
 174   int exist, nparam;
 175   char s[128];
 176   evalue *EP;
 177   Enumeration *en;
 178
 179 /******* Read the input *********/
 180   P1 = Matrix_Read();
 181
 182     fgets(s, 128, stdin);
 183     while ((*s=='#') || (sscanf(s, "E %d", &exist)<1))
 184         fgets(s, 128, stdin);
 185
 186     fgets(s, 128, stdin);
 187     while ((*s=='#') || (sscanf(s, "P %d", &nparam)<1))
 188         fgets(s, 128, stdin);
 189
 190   P = Constraints2Polyhedron(P1,MAXRAYS);
 191   params = Read_ParamNames(stdin, nparam);
 192   Matrix_Free(P1);
 193
 194   /******* Read the options: initialize min and max ********/
 195   if(P->Dimension >= VBIGDIM)
 196     M = VSRANGE;
 197   else if(P->Dimension >= BIGDIM)
 198     M = SRANGE;
 199   else
 200     M = RANGE;
 201   m = -M;
 202   if(argc != 1 ) {
 203     for(i=1;i<argc;i++) {
 204       if(!strncmp(argv[i],"-m",2)) {
 205
 206         /* min specified */
 207         m = atoi(&argv[i][2]);
 208       }
 209       else if(!strncmp(argv[i],"-M",2)) {
 210
 211         /* max specified */
 212         M = atoi(&argv[i][2]);
 213       }
 214       else if(!strncmp(argv[i], "-r", 2)) {
 215
 216         /* range specified */
 217         M = atoi(&argv[i][2]);
 218         m = -M;
 219       }
 220       else {
 221         fprintf(stderr,"Unknown option: %s\n",argv[i]);
 222         fprintf(stderr,"Usage: %s [-m<>][-M<>][-r<>]\n",argv[0]);
 223         return(-1);
 224       }
 225     }
 226   }
 227   if(m > M) {
 228     fprintf(stderr,"Nothing to do: min > max !\n");
 229     return(0);
 230   }
 231   value_init(min);
 232   value_init(max);
 233   value_set_si(min,m);
 234   value_set_si(max,M);
 235   value_init(tmp);
 236
 237   /******* Compute true context *******/
 238   C1 = Matrix_Alloc(nparam+1,P->Dimension+1);
 239
 240   for(i=0;i<C1->NbRows;i++)
 241     for(j=0;j<C1->NbColumns;j++)
 242       if(i==j-P->Dimension+nparam)
 243         value_set_si(C1->p[i][j],1);
 244       else
 245         value_set_si(C1->p[i][j],0);
 246   C = Polyhedron_Image(P,C1,MAXRAYS);
 247
 248   /******* Compute EP *********/
 249   EP = barvinok_enumerate_e(P, exist, nparam, MAXRAYS);
 250   en = partition2enumeration(EP);
 251
 252   /******* Initializations for check *********/
 253   p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
 254   for(i=0;i<=P->Dimension;i++) {
 255     value_init(p[i]);
 256     value_set_si(p[i],0);
 257   }
 258   value_init(p[i]);
 259   value_set_si(p[i],1);
 260
 261   /* S = scanning list of polyhedra */
 262   S = Polyhedron_Scan(P,C,MAXRAYS);
 263
 264 #ifndef PRINT_ALL_RESULTS
 265   if(C->Dimension > 0) {
 266     value_substract(tmp,max,min);
 267     if (VALUE_TO_INT(tmp) > 80)
 268       st = 1+(VALUE_TO_INT(tmp))/80;
 269     else
 270       st=1;
 271     for(i=VALUE_TO_INT(min);i<=VALUE_TO_INT(max);i+=st)
 272       printf(".");
 273     printf( "\r" );
 274     fflush(stdout);
 275   }
 276 #endif
 277
 278   /******* CHECK NOW *********/
 279   if(S && !check_poly(S, C, en, exist, nparam, 0, p)) {
 280     fprintf(stderr,"Check failed !\n");
 281     for(i=0;i<=(P->Dimension+1);i++)
 282       value_clear(p[i]);
 283     value_clear(tmp);
 284     return(-1);
 285   }
 286
 287 #ifndef PRINT_ALL_RESULTS
 288   printf( "\n" );
 289 #endif
 290
 291   for(i=0;i<=(P->Dimension+1);i++)
 292     value_clear(p[i]);
 293   value_clear(tmp);
 294   Free_ParamNames(params, C->Dimension);
 295   return(0);
 296 } /* main */
 297
 298
 299
 300