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