test_approx.c

   1 #include <assert.h>
   2 #include <ctype.h>
   3 #include <limits.h>
   4 #include <math.h>
   5 #include <stdlib.h>
   6 #include <string.h>
   7 #include <isl/ctx.h>
   8 #include <isl/val.h>
   9 #include <isl/space.h>
  10 #include <isl/set.h>
  11 #include <isl/polynomial.h>
  12 #include <barvinok/polylib.h>
  13 #include <barvinok/isl.h>
  14 #include <barvinok/options.h>
  15 #include "verify.h"
  16 #include "config.h"
  17
  18 #ifdef HAVE_SYS_TIMES_H
  19
  20 #include <sys/times.h>
  21
  22 typedef clock_t         my_clock_t;
  23
  24 static my_clock_t time_diff(struct tms *before, struct tms *after)
  25 {
  26         return after->tms_utime - before->tms_utime;
  27 }
  28
  29 #else
  30
  31 typedef int             my_clock_t;
  32
  33 struct tms { int dummy; };
  34 static void times(struct tms* time)
  35 {
  36 }
  37 static my_clock_t time_diff(struct tms *before, struct tms *after)
  38 {
  39         return 0;
  40 }
  41
  42 #endif
  43
  44 struct {
  45     int     sign;
  46     int     method;
  47     int     flags;
  48 } methods[] = {
  49 { BV_APPROX_SIGN_NONE,      BV_APPROX_NONE,     0 },
  50 { BV_APPROX_SIGN_APPROX,    BV_APPROX_BERNOULLI,        0 },
  51 { BV_APPROX_SIGN_APPROX,    BV_APPROX_DROP,     0 },
  52 { BV_APPROX_SIGN_APPROX,    BV_APPROX_VOLUME,   BV_VOL_LIFT },
  53 { BV_APPROX_SIGN_APPROX,    BV_APPROX_VOLUME,   BV_VOL_VERTEX },
  54 { BV_APPROX_SIGN_APPROX,    BV_APPROX_VOLUME,   BV_VOL_BARYCENTER },
  55 { BV_APPROX_SIGN_APPROX,    BV_APPROX_SCALE,    0 },
  56 { BV_APPROX_SIGN_APPROX,    BV_APPROX_SCALE,    BV_APPROX_SCALE_CHAMBER },
  57 { BV_APPROX_SIGN_APPROX,    BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST },
  58 { BV_APPROX_SIGN_APPROX,    BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
  59 { BV_APPROX_SIGN_LOWER,     BV_APPROX_DROP,     0 },
  60 { BV_APPROX_SIGN_LOWER,     BV_APPROX_VOLUME,   BV_VOL_LIFT },
  61 { BV_APPROX_SIGN_LOWER,     BV_APPROX_VOLUME,   BV_VOL_VERTEX },
  62 { BV_APPROX_SIGN_LOWER,     BV_APPROX_VOLUME,   BV_VOL_BARYCENTER },
  63 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    0 },
  64 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_CHAMBER },
  65 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST },
  66 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
  67 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW },
  68 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW  | BV_APPROX_SCALE_CHAMBER},
  69 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW },
  70 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
  71 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW2 },
  72 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
  73 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 },
  74 { BV_APPROX_SIGN_LOWER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
  75 { BV_APPROX_SIGN_UPPER,     BV_APPROX_DROP,     0 },
  76 { BV_APPROX_SIGN_UPPER,     BV_APPROX_VOLUME,   BV_VOL_LIFT },
  77 { BV_APPROX_SIGN_UPPER,     BV_APPROX_VOLUME,   BV_VOL_VERTEX },
  78 { BV_APPROX_SIGN_UPPER,     BV_APPROX_VOLUME,   BV_VOL_BARYCENTER },
  79 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    0 },
  80 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_CHAMBER },
  81 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST },
  82 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
  83 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW },
  84 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
  85 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW },
  86 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
  87 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW2 },
  88 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
  89 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 },
  90 { BV_APPROX_SIGN_UPPER,     BV_APPROX_SCALE,    BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
  91 };
  92
  93 #define nr_methods (sizeof(methods)/sizeof(*methods))
  94
  95 struct options {
  96     int quiet;
  97     struct verify_options    *verify;
  98 };
  99
 100 ISL_ARGS_START(struct options, options_args)
 101 ISL_ARG_CHILD(struct options, verify, NULL, &verify_options_args, NULL)
 102 ISL_ARG_BOOL(struct options, quiet, 'q', "quiet", 0, NULL)
 103 ISL_ARGS_END
 104
 105 ISL_ARG_DEF(options, struct options, options_args)
 106
 107 struct result_data {
 108     Value                   n;
 109     double                  RE_sum[nr_methods];
 110
 111     my_clock_t              ticks[nr_methods];
 112     size_t                  size[nr_methods];
 113 };
 114
 115 void result_data_init(struct result_data *result)
 116 {
 117     int i;
 118     for (i = 0; i < nr_methods; ++i) {
 119         result->RE_sum[i] = 0;
 120         result->ticks[i] = 0;
 121         result->size[i] = 0;
 122     }
 123     value_init(result->n);
 124 }
 125
 126 void result_data_clear(struct result_data *result)
 127 {
 128     value_clear(result->n);
 129 }
 130
 131 void result_data_print(struct result_data *result, int n)
 132 {
 133     int i;
 134
 135     fprintf(stderr, "%g", (double)result->ticks[0]/n);
 136     for (i = 1; i < nr_methods; ++i)
 137         fprintf(stderr, ", %g", (double)result->ticks[i]/n);
 138     fprintf(stderr, "\n");
 139
 140     fprintf(stderr, "%zd", result->size[0]/n);
 141     for (i = 1; i < nr_methods; ++i)
 142         fprintf(stderr, ", %zd", result->size[i]/n);
 143     fprintf(stderr, "\n");
 144
 145     fprintf(stderr, "%g\n", VALUE_TO_DOUBLE(result->n));
 146     fprintf(stderr, "%g", result->RE_sum[0]/VALUE_TO_DOUBLE(result->n));
 147     for (i = 1; i < nr_methods; ++i)
 148         fprintf(stderr, ", %g", result->RE_sum[i]/VALUE_TO_DOUBLE(result->n));
 149     fprintf(stderr, "\n");
 150 }
 151
 152 struct test_approx_data {
 153     struct verify_point_data vpd;
 154     isl_pw_qpolynomial      **pwqp;
 155     struct result_data      *result;
 156 };
 157
 158 static __isl_give isl_val *eval(__isl_keep isl_pw_qpolynomial *pwqp,
 159         __isl_keep isl_point *pnt, int sign)
 160 {
 161     isl_val *res;
 162
 163     pwqp = isl_pw_qpolynomial_copy(pwqp);
 164     res = isl_pw_qpolynomial_eval(pwqp, isl_point_copy(pnt));
 165     if (sign == BV_APPROX_SIGN_LOWER)
 166         res = isl_val_ceil(res);
 167     else if (sign == BV_APPROX_SIGN_UPPER)
 168         res = isl_val_floor(res);
 169     else if (sign == BV_APPROX_SIGN_APPROX)
 170         res = isl_val_trunc(res);
 171     else
 172         assert(isl_val_is_int(res));
 173
 174     return res;
 175 }
 176
 177 static isl_stat test_approx(__isl_take isl_point *pnt, void *user)
 178 {
 179     struct test_approx_data *ta_data = (struct test_approx_data *) user;
 180     isl_val *exact, *approx;
 181     int i;
 182
 183     ta_data->vpd.n--;
 184
 185     exact = eval(ta_data->pwqp[0], pnt, BV_APPROX_SIGN_NONE);
 186
 187     value_increment(ta_data->result->n, ta_data->result->n);
 188     for (i = 1; i < nr_methods; ++i) {
 189         double error;
 190         approx = eval(ta_data->pwqp[i], pnt, methods[i].sign);
 191         assert(approx);
 192         if (methods[i].sign == BV_APPROX_SIGN_LOWER)
 193             assert(isl_val_le(approx, exact));
 194         if (methods[i].sign == BV_APPROX_SIGN_UPPER)
 195             assert(isl_val_ge(approx, exact));
 196         approx = isl_val_sub(approx, isl_val_copy(exact));
 197         if (isl_val_is_zero(exact))
 198             error = fabs(isl_val_get_d(approx));
 199         else
 200             error = fabs(isl_val_get_d(approx)) / isl_val_get_d(exact);
 201         isl_val_free(approx);
 202         ta_data->result->RE_sum[i] += error;
 203     }
 204
 205     if (!ta_data->vpd.options->print_all &&
 206         (ta_data->vpd.n % ta_data->vpd.s) == 0) {
 207             printf("o");
 208             fflush(stdout);
 209     }
 210
 211     isl_val_free(exact);
 212     isl_point_free(pnt);
 213
 214     return (ta_data->vpd.n >= 1) ? isl_stat_ok : isl_stat_error;
 215 }
 216
 217 static int test(__isl_keep isl_set *context,
 218         __isl_keep isl_pw_qpolynomial **pwqp, struct result_data *result,
 219         struct verify_options *options)
 220 {
 221         struct test_approx_data data = { { options } };
 222         int r;
 223
 224         r = verify_point_data_init(&data.vpd, context);
 225
 226         data.pwqp = pwqp;
 227         data.result = result;
 228         if (r == 0)
 229                 isl_set_foreach_point(context, &test_approx, &data);
 230         if (data.vpd.error)
 231                 r = -1;
 232
 233         verify_point_data_fini(&data.vpd);
 234
 235         return r;
 236 }
 237
 238 /* Turn the set dimensions of "context" into parameters and return
 239  * the corresponding parameter domain.
 240  */
 241 static __isl_give isl_set *to_parameter_domain(__isl_take isl_set *context)
 242 {
 243         context = isl_set_move_dims(context, isl_dim_param, 0, isl_dim_set, 0,
 244                                        isl_set_dim(context, isl_dim_set));
 245         return isl_set_params(context);
 246 }
 247
 248 static int handle(isl_ctx *ctx, FILE *in, struct result_data *result,
 249         struct verify_options *options)
 250 {
 251     int i;
 252     int r;
 253     int nparam;
 254     isl_space *space;
 255     isl_set *set;
 256     isl_set *context;
 257     isl_pw_qpolynomial *pwqp[nr_methods];
 258
 259     set = isl_set_read_from_file(ctx, in);
 260     context = isl_set_read_from_file(ctx, in);
 261
 262     context = to_parameter_domain(context);
 263     nparam = isl_set_dim(context, isl_dim_param);
 264     if (nparam != isl_set_dim(set, isl_dim_param)) {
 265             int dim = isl_set_dim(set, isl_dim_set);
 266             set = isl_set_move_dims(set, isl_dim_param, 0,
 267                                         isl_dim_set, dim - nparam, nparam);
 268     }
 269
 270     set = isl_set_intersect_params(set, context);
 271     context = isl_set_params(isl_set_copy(set));
 272     space = isl_set_get_space(context);
 273
 274     context = verify_context_set_bounds(context, options);
 275
 276     for (i = 0; i < nr_methods; ++i) {
 277         struct tms st_cpu;
 278         struct tms en_cpu;
 279         options->barvinok->approx->approximation = methods[i].sign;
 280         options->barvinok->approx->method = methods[i].method;
 281         if (methods[i].method == BV_APPROX_SCALE)
 282             options->barvinok->approx->scale_flags = methods[i].flags;
 283         else if (methods[i].method == BV_APPROX_VOLUME)
 284             options->barvinok->approx->volume_triangulate = methods[i].flags;
 285
 286         times(&st_cpu);
 287         pwqp[i] = isl_set_card(isl_set_copy(set));
 288         times(&en_cpu);
 289         result->ticks[i] = time_diff(&en_cpu, &st_cpu);
 290     }
 291     for (i = 0; i < nr_methods; ++i)
 292         result->size[i] = 0;
 293     r = test(context, pwqp, result, options);
 294     for (i = 0; i < nr_methods; ++i)
 295         isl_pw_qpolynomial_free(pwqp[i]);
 296
 297     isl_space_free(space);
 298     isl_set_free(context);
 299     isl_set_free(set);
 300
 301     return r;
 302 }
 303
 304 int main(int argc, char **argv)
 305 {
 306     isl_ctx *ctx;
 307     char path[PATH_MAX+1];
 308     struct result_data all_result;
 309     int n = 0;
 310     int r = EXIT_SUCCESS;
 311     struct options *options = options_new_with_defaults();
 312
 313     argc = options_parse(options, argc, argv, ISL_ARG_ALL);
 314     ctx = isl_ctx_alloc_with_options(&options_args, options);
 315
 316     if (options->verify->M == INT_MIN)
 317         options->verify->M = 100;
 318     if (options->verify->m == INT_MAX)
 319         options->verify->m = -100;
 320
 321     result_data_init(&all_result);
 322
 323     while (fgets(path, sizeof(path), stdin)) {
 324         struct result_data result;
 325         FILE *in;
 326         int i;
 327
 328         ++n;
 329         result_data_init(&result);
 330         fprintf(stderr, "%s", path);
 331         *strchr(path, '\n') = '\0';
 332         in = fopen(path, "r");
 333         assert(in);
 334         if (handle(ctx, in, &result, options->verify) < 0)
 335             r = EXIT_FAILURE;
 336         fclose(in);
 337
 338         if (!options->quiet)
 339             result_data_print(&result, 1);
 340
 341         value_addto(all_result.n, all_result.n, result.n);
 342         for (i = 0; i < nr_methods; ++i) {
 343             all_result.RE_sum[i] += result.RE_sum[i];
 344             all_result.ticks[i] += result.ticks[i];
 345             all_result.size[i] += result.size[i];
 346         }
 347
 348         result_data_clear(&result);
 349
 350         if (!options->quiet) {
 351             fprintf(stderr, "average\n");
 352             result_data_print(&all_result, n);
 353         }
 354     }
 355
 356     result_data_clear(&all_result);
 357
 358     isl_ctx_free(ctx);
 359
 360     return r;
 361 }