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