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