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