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