test_approx.c

   1 #include <limits.h>
   2 #include <sys/times.h>
   3 #include <barvinok/barvinok.h>
   4 #include "verify.h"
   5
   6 struct {
   7     int     sign;
   8     int     method;
   9     int     flags;
  10 } methods[] = {
  11 { BV_APPROX_SIGN_NONE,      BV_APPROX_NONE,     0 },
  12 { BV_APPROX_SIGN_APPROX,    BV_APPROX_DROP,     0 },
  13 { BV_APPROX_SIGN_APPROX,    BV_APPROX_VOLUME,   0 },
  14 { BV_APPROX_SIGN_APPROX,    BV_APPROX_SCALE,    0 },
  15 { BV_APPROX_SIGN_APPROX,    BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST },
  16 { BV_APPROX_SIGN_LOWER,     BV_APPROX_DROP,     0 },
  17 { BV_APPROX_SIGN_LOWER,     BV_APPROX_VOLUME,   0 },
  18 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    0 },
  19 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST },
  20 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW },
  21 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW },
  22 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW2 },
  23 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 },
  24 { BV_APPROX_SIGN_UPPER,     BV_APPROX_DROP,     0 },
  25 { BV_APPROX_SIGN_UPPER,     BV_APPROX_VOLUME,   0 },
  26 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    0 },
  27 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST },
  28 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW },
  29 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW },
  30 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW2 },
  31 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 },
  32 };
  33
  34 #define nr_methods (sizeof(methods)/sizeof(*methods))
  35
  36 struct argp_option argp_options[] = {
  37     { "quiet",      'q' },
  38     { 0 }
  39 };
  40
  41 struct options {
  42     int quiet;
  43     struct verify_options    verify;
  44 };
  45
  46 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  47 {
  48     struct options *options = (struct options*) state->input;
  49
  50     switch (key) {
  51     case ARGP_KEY_INIT:
  52         state->child_inputs[0] = options->verify.barvinok;
  53         state->child_inputs[1] = &options->verify;
  54         options->quiet = 0;
  55         break;
  56     case 'q':
  57         options->quiet = 1;
  58         break;
  59     default:
  60         return ARGP_ERR_UNKNOWN;
  61     }
  62     return 0;
  63 }
  64
  65 struct result_data {
  66     Value                   n;
  67     double                  RE_sum[nr_methods];
  68
  69     clock_t                 ticks[nr_methods];
  70     size_t                  size[nr_methods];
  71 };
  72
  73 void result_data_init(struct result_data *result)
  74 {
  75     int i;
  76     for (i = 0; i < nr_methods; ++i) {
  77         result->RE_sum[i] = 0;
  78         result->ticks[i] = 0;
  79         result->size[i] = 0;
  80     }
  81     value_init(result->n);
  82 }
  83
  84 void result_data_clear(struct result_data *result)
  85 {
  86     int i;
  87     value_clear(result->n);
  88 }
  89
  90 void result_data_print(struct result_data *result, int n)
  91 {
  92     int i;
  93
  94     fprintf(stderr, "%d", result->ticks[0]/n);
  95     for (i = 1; i < nr_methods; ++i)
  96         fprintf(stderr, ", %d", result->ticks[i]/n);
  97     fprintf(stderr, "\n");
  98
  99     fprintf(stderr, "%d", result->size[0]/n);
 100     for (i = 1; i < nr_methods; ++i)
 101         fprintf(stderr, ", %d", result->size[i]/n);
 102     fprintf(stderr, "\n");
 103
 104     fprintf(stderr, "%g\n", VALUE_TO_DOUBLE(result->n));
 105     fprintf(stderr, "%g", result->RE_sum[0]/VALUE_TO_DOUBLE(result->n));
 106     for (i = 1; i < nr_methods; ++i)
 107         fprintf(stderr, ", %g", result->RE_sum[i]/VALUE_TO_DOUBLE(result->n));
 108     fprintf(stderr, "\n");
 109 }
 110
 111 struct test_approx_data {
 112     struct check_poly_data   cp;
 113     evalue                  **EP;
 114     struct result_data      *result;
 115 };
 116
 117 static void eval(const evalue *EP, Value *z, int sign, Value *v)
 118 {
 119     evalue *res;
 120
 121     res = evalue_eval(EP, z);
 122     if (sign == BV_APPROX_SIGN_LOWER)
 123         mpz_cdiv_q(*v, res->x.n, res->d);
 124     else if (sign == BV_APPROX_SIGN_UPPER)
 125         mpz_fdiv_q(*v, res->x.n, res->d);
 126     else if (sign == BV_APPROX_SIGN_APPROX)
 127         mpz_tdiv_q(*v, res->x.n, res->d);
 128     else {
 129         assert(value_one_p(res->d));
 130         value_assign(*v, res->x.n);
 131     }
 132     free_evalue_refs(res);
 133     free(res);
 134 }
 135
 136 static int test_approx(const struct check_poly_data *data, int nparam, Value *z,
 137                        const struct verify_options *options)
 138 {
 139     struct test_approx_data* ta_data = (struct test_approx_data*) data;
 140     Value exact, approx;
 141     int i;
 142
 143     value_init(exact);
 144     value_init(approx);
 145
 146     eval(ta_data->EP[0], z, BV_APPROX_SIGN_NONE, &exact);
 147
 148     /*
 149     value_print(stderr, VALUE_FMT, exact);
 150     */
 151
 152     value_increment(ta_data->result->n, ta_data->result->n);
 153     for (i = 1; i < nr_methods; ++i) {
 154         double error;
 155         eval(ta_data->EP[i], z, methods[i].sign, &approx);
 156         /*
 157         fprintf(stderr, ", ");
 158         value_print(stderr, VALUE_FMT, approx);
 159         */
 160         if (methods[i].sign == BV_APPROX_SIGN_LOWER)
 161             assert(value_le(approx, exact));
 162         if (methods[i].sign == BV_APPROX_SIGN_UPPER)
 163             assert(value_ge(approx, exact));
 164         value_subtract(approx, approx, exact);
 165         if (value_zero_p(exact))
 166             error = abs(VALUE_TO_DOUBLE(approx));
 167         else
 168             error = abs(VALUE_TO_DOUBLE(approx)) / VALUE_TO_DOUBLE(exact);
 169         ta_data->result->RE_sum[i] += error;
 170     }
 171
 172     /*
 173     fprintf(stderr, "\n");
 174     */
 175
 176     value_clear(exact);
 177     value_clear(approx);
 178     return 1;
 179 }
 180
 181 static void test(Polyhedron *P, Polyhedron *C, evalue **EP,
 182                  struct result_data *result,
 183                  struct verify_options *options)
 184 {
 185     Polyhedron *CS;
 186     Vector *p;
 187     unsigned nparam = C->Dimension;
 188     struct test_approx_data data;
 189
 190     CS = check_poly_context_scan(P, &C, C->Dimension, options);
 191
 192     p = Vector_Alloc(P->Dimension+2);
 193     value_set_si(p->p[P->Dimension+1], 1);
 194
 195     check_poly_init(C, options);
 196
 197     data.cp.z = p->p;
 198     data.cp.check = test_approx;
 199     data.EP = EP;
 200     data.result = result;
 201     check_poly(CS, &data.cp, nparam, 0, p->p+P->Dimension-nparam+1,
 202                options);
 203     if (!options->print_all)
 204         printf("\n");
 205
 206     Vector_Free(p);
 207     if (CS) {
 208         Domain_Free(CS);
 209         Domain_Free(C);
 210     }
 211 }
 212
 213 void Matrix_File_Read_Input(FILE *in, Matrix *Mat)
 214 {
 215   Value *p;
 216   int i,j,n;
 217   char *c, s[1024],str[1024];
 218
 219   p = Mat->p_Init;
 220   for (i=0;i<Mat->NbRows;i++) {
 221     do {
 222       c = fgets(s, 1024, in);
 223       while(isspace(*c) && *c!='\n')
 224         ++c;
 225     } while(c && (*c =='#' || *c== '\n'));
 226
 227     if (!c) {
 228       errormsg1( "Matrix_Read", "baddim", "not enough rows" );
 229       break;
 230     }
 231     for (j=0;j<Mat->NbColumns;j++) {
 232       if(!c || *c=='\n' || *c=='#') {
 233         errormsg1("Matrix_Read", "baddim", "not enough columns");
 234         break;
 235       }
 236       if (sscanf(c,"%s%n",str,&n) == 0) {
 237         errormsg1( "Matrix_Read", "baddim", "not enough columns" );
 238         break;
 239       }
 240       value_read(*(p++),str);
 241       c += n;
 242     }
 243   }
 244 } /* Matrix_Read_Input */
 245
 246 /*
 247  * Read the contents of the matrix 'Mat' from standard input.
 248  * A '#' in the first column is a comment line
 249  */
 250 Matrix *Matrix_File_Read(FILE *in)
 251 {
 252   Matrix *Mat;
 253   unsigned NbRows, NbColumns;
 254   char s[1024];
 255
 256   fgets(s, 1024, in);
 257   while ((*s=='#' || *s=='\n') ||
 258          (sscanf(s, "%d %d", &NbRows, &NbColumns)<2))
 259     fgets(s, 1024, in);
 260   Mat = Matrix_Alloc(NbRows,NbColumns);
 261   if(!Mat) {
 262     errormsg1("Matrix_Read", "outofmem", "out of memory space");
 263     return(NULL);
 264   }
 265   Matrix_File_Read_Input(in, Mat);
 266   return Mat;
 267 } /* Matrix_Read */
 268
 269 void handle(FILE *in, struct result_data *result, struct verify_options *options)
 270 {
 271     int i;
 272     Polyhedron *A, *C;
 273     Matrix *M;
 274     char **param_name;
 275     evalue *EP[nr_methods];
 276
 277     M = Matrix_File_Read(in);
 278     A = Constraints2Polyhedron(M, options->barvinok->MaxRays);
 279     Matrix_Free(M);
 280     M = Matrix_File_Read(in);
 281     C = Constraints2Polyhedron(M, options->barvinok->MaxRays);
 282     Matrix_Free(M);
 283     param_name = Read_ParamNames(in, C->Dimension);
 284
 285     for (i = 0; i < nr_methods; ++i) {
 286         struct tms st_cpu;
 287         struct tms en_cpu;
 288         options->barvinok->polynomial_approximation = methods[i].sign;
 289         options->barvinok->approximation_method = methods[i].method;
 290         options->barvinok->scale_flags = methods[i].flags;
 291
 292         times(&st_cpu);
 293         EP[i] = barvinok_enumerate_with_options(A, C, options->barvinok);
 294         times(&en_cpu);
 295         result->ticks[i] = en_cpu.tms_utime - st_cpu.tms_utime;
 296         /*
 297         print_evalue(stdout, EP[i], param_name);
 298         */
 299     }
 300     for (i = 0; i < nr_methods; ++i)
 301         result->size[i] = evalue_size(EP[i])/4;
 302     test(A, C, EP, result, options);
 303     for (i = 0; i < nr_methods; ++i) {
 304         free_evalue_refs(EP[i]);
 305         free(EP[i]);
 306     }
 307
 308     Free_ParamNames(param_name, C->Dimension);
 309     Polyhedron_Free(A);
 310     Polyhedron_Free(C);
 311 }
 312
 313 int main(int argc, char **argv)
 314 {
 315     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 316     char path[PATH_MAX+1];
 317     struct result_data all_result;
 318     int n = 0;
 319     static struct argp_child argp_children[] = {
 320         { &barvinok_argp,       0,      0,              0 },
 321         { &verify_argp,         0,      "verification",         BV_GRP_LAST+1 },
 322         { 0 }
 323     };
 324     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 325     struct options options;
 326
 327     options.verify.barvinok = bv_options;
 328     argp_parse(&argp, argc, argv, 0, 0, &options);
 329
 330     if (options.verify.M == INT_MIN)
 331         options.verify.M = 100;
 332     if (options.verify.m == INT_MAX)
 333         options.verify.m = -100;
 334
 335     result_data_init(&all_result);
 336
 337     while (fgets(path, sizeof(path), stdin)) {
 338         struct result_data result;
 339         FILE *in;
 340         int i;
 341
 342         ++n;
 343         result_data_init(&result);
 344         fprintf(stderr, "%s", path);
 345         *strchr(path, '\n') = '\0';
 346         in = fopen(path, "r");
 347         assert(in);
 348         handle(in, &result, &options.verify);
 349         fclose(in);
 350
 351         if (!options.quiet)
 352             result_data_print(&result, 1);
 353
 354         value_addto(all_result.n, all_result.n, result.n);
 355         for (i = 0; i < nr_methods; ++i) {
 356             all_result.RE_sum[i] += result.RE_sum[i];
 357             all_result.ticks[i] += result.ticks[i];
 358             all_result.size[i] += result.size[i];
 359         }
 360
 361         result_data_clear(&result);
 362
 363         if (!options.quiet) {
 364             fprintf(stderr, "average\n");
 365             result_data_print(&all_result, n);
 366         }
 367     }
 368
 369     result_data_clear(&all_result);
 370
 371     barvinok_options_free(bv_options);
 372
 373     return 0;
 374 }