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