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