range.cc

   1 #include <assert.h>
   2 #include <bernstein/bernstein.h>
   3 #include <bernstein/piecewise_lst.h>
   4 #include <barvinok/bernstein.h>
   5 #include <barvinok/options.h>
   6 #include <barvinok/polylib.h>
   7 #include <barvinok/util.h>
   8 #include <barvinok/evalue.h>
   9 #include "range.h"
  10
  11 using namespace GiNaC;
  12 using namespace bernstein;
  13 using namespace barvinok;
  14
  15 using std::cerr;
  16 using std::endl;
  17
  18 #define ALLOC(type) (type*)malloc(sizeof(type))
  19 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  20
  21 struct range_data {
  22     struct barvinok_options *options;
  23     unsigned             nparam;
  24     const evalue        *poly;
  25     int                 *signs;
  26     int                  sign;
  27     int                  test_monotonicity;
  28     int                  monotonicity;
  29     piecewise_lst       *pl;
  30     piecewise_lst       *pl_all;
  31     exvector             params;
  32 };
  33
  34 static int type_offset(enode *p)
  35 {
  36    return p->type == fractional ? 1 :
  37           p->type == flooring ? 1 :
  38           p->type == relation ? 1 : 0;
  39 }
  40
  41 /* Remove all terms from e that are not positive (sign > 0)
  42  * or not negative (sign < 0)
  43  */
  44 static void evalue_fixed_sign_terms(evalue *e, int* signs, int sign)
  45 {
  46     int i, offset;
  47     int sign_odd = sign;
  48
  49     assert(!EVALUE_IS_NAN(*e));
  50     if (value_notzero_p(e->d)) {
  51         if (sign * value_sign(e->x.n) < 0) {
  52             value_set_si(e->x.n, 0);
  53             value_set_si(e->d, 1);
  54         }
  55         return;
  56     }
  57
  58     if (e->x.p->type == relation) {
  59         for (i = e->x.p->size-1; i >= 1; --i)
  60             evalue_fixed_sign_terms(&e->x.p->arr[i], signs, sign);
  61         return;
  62     }
  63
  64     if (e->x.p->type == polynomial)
  65         sign_odd *= signs[e->x.p->pos-1];
  66
  67     offset = type_offset(e->x.p);
  68     for (i = offset; i < e->x.p->size; ++i)
  69         evalue_fixed_sign_terms(&e->x.p->arr[i], signs,
  70                                      (i - offset) % 2 ? sign_odd : sign);
  71 }
  72
  73 static void propagate_on_domain(Polyhedron *D, const evalue *poly,
  74                                 struct range_data *data);
  75
  76 static evalue *bound2evalue(Value *bound, unsigned dim)
  77 {
  78     Value denom;
  79     evalue *b;
  80
  81     value_init(denom);
  82
  83     if (!bound) {
  84         b = evalue_nan();
  85     } else if (value_pos_p(bound[1])) {
  86         value_assign(denom, bound[1]);
  87         b = affine2evalue(bound+2, denom, dim);
  88         evalue_negate(b);
  89     } else {
  90         value_oppose(denom, bound[1]);
  91         b = affine2evalue(bound+2, denom, dim);
  92     }
  93
  94     value_clear(denom);
  95
  96     return b;
  97 }
  98
  99 static int evalue_is_constant(const evalue *e)
 100 {
 101     return EVALUE_IS_NAN(*e) || value_notzero_p(e->d);
 102 }
 103
 104 static void range_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
 105 {
 106     struct range_data *data = (struct range_data *)cb_data;
 107     evalue *l, *u;
 108     unsigned dim = M->NbColumns-2;
 109     evalue *pos, *neg;
 110     evalue **subs;
 111     Matrix *M2;
 112     Polyhedron *T;
 113
 114     M2 = Matrix_Copy(M);
 115     T = Constraints2Polyhedron(M2, data->options->MaxRays);
 116     Matrix_Free(M2);
 117
 118     POL_ENSURE_VERTICES(T);
 119     if (emptyQ(T)) {
 120         Polyhedron_Free(T);
 121         return;
 122     }
 123
 124     l = bound2evalue(lower, dim);
 125     u = bound2evalue(upper, dim);
 126
 127     subs = ALLOCN(evalue *, 1+dim);
 128     for (int i = 0; i < dim; ++i)
 129         subs[1+i] = evalue_var(i);
 130
 131     if (evalue_is_constant(data->poly)) {
 132         pos = evalue_dup(data->poly);
 133     } else if (data->monotonicity) {
 134         pos = evalue_dup(data->poly);
 135         if (data->monotonicity * data->sign > 0)
 136             subs[0] = u;
 137         else
 138             subs[0] = l;
 139         evalue_substitute(pos, subs);
 140     } else {
 141         pos = evalue_dup(data->poly);
 142         neg = evalue_dup(data->poly);
 143
 144         evalue_fixed_sign_terms(pos, data->signs, data->sign);
 145         evalue_fixed_sign_terms(neg, data->signs, -data->sign);
 146
 147         subs[0] = u;
 148         evalue_substitute(pos, subs);
 149
 150         subs[0] = l;
 151         evalue_substitute(neg, subs);
 152
 153         eadd(neg, pos);
 154         evalue_free(neg);
 155     }
 156
 157     for (int i = 0; i < dim; ++i)
 158         evalue_free(subs[1+i]);
 159     free(subs);
 160     evalue_free(l);
 161     evalue_free(u);
 162
 163     if (dim == data->nparam) {
 164         data->pl->add_guarded_lst(T, lst(evalue2ex(pos, data->params)));
 165     } else {
 166         propagate_on_domain(T, pos, data);
 167         Polyhedron_Free(T);
 168     }
 169
 170     evalue_free(pos);
 171 }
 172
 173 static int has_sign(Polyhedron *D, evalue *poly, int sign,
 174                     int *signs, struct barvinok_options *options)
 175 {
 176     struct range_data data_m;
 177     int i;
 178     lst l;
 179
 180     data_m.options = options;
 181     data_m.nparam = 0;
 182     data_m.test_monotonicity = 0;
 183     data_m.signs = signs;
 184
 185     data_m.pl_all = NULL;
 186     data_m.sign = -sign;
 187     propagate_on_domain(D, poly, &data_m);
 188
 189     assert(data_m.pl_all->list.size() == 1);
 190     assert(universeQ(data_m.pl_all->list[0].first));
 191     l = data_m.pl_all->list[0].second;
 192
 193     for (i = 0; i < l.nops(); ++i) {
 194         if (is_exactly_a<fail>(l.op(i)))
 195             break;
 196         if (sign * l.op(i) < 0)
 197             break;
 198     }
 199     delete data_m.pl_all;
 200
 201     return i == l.nops();
 202 }
 203
 204 /* Returns 1 if poly is monotonically increasing,
 205  *        -1 if poly is monotonically decreasing,
 206  *         0 if no conclusion.
 207  */
 208 static int monotonicity(Polyhedron *D, const evalue *poly,
 209                         struct range_data *data)
 210 {
 211     evalue **subs;
 212     evalue *diff;
 213     Value one;
 214     int result = 0;
 215
 216     value_init(one);
 217     value_set_si(one, 1);
 218
 219     subs = ALLOCN(evalue *, D->Dimension);
 220     for (int i = 0; i < D->Dimension; ++i)
 221         subs[i] = evalue_var(i);
 222
 223     evalue_add_constant(subs[0], one);
 224
 225     diff = evalue_dup(poly);
 226     evalue_substitute(diff, subs);
 227
 228     for (int i = 0; i < D->Dimension; ++i)
 229         evalue_free(subs[i]);
 230     free(subs);
 231
 232     evalue_negate(diff);
 233     eadd(poly, diff);
 234     reduce_evalue(diff);
 235     evalue_negate(diff);
 236
 237     value_clear(one);
 238
 239     if (has_sign(D, diff, 1, data->signs, data->options))
 240         result = 1;
 241     else if (has_sign(D, diff, -1, data->signs, data->options))
 242         result = -1;
 243
 244     evalue_free(diff);
 245     return result;
 246 }
 247
 248 static void propagate_on_domain(Polyhedron *D, const evalue *poly,
 249                                 struct range_data *data)
 250 {
 251     const evalue *save_poly = data->poly;
 252     int save_monotonicity = data->monotonicity;
 253
 254     assert(D->Dimension > data->nparam);
 255
 256     if (data->test_monotonicity)
 257         data->monotonicity = monotonicity(D, poly, data);
 258     else
 259         data->monotonicity = 0;
 260
 261     if (D->Dimension == data->nparam+1)
 262         data->pl = new piecewise_lst(data->params);
 263
 264     data->poly = poly;
 265     for_each_lower_upper_bound(D, NULL, range_cb, data);
 266
 267     if (D->Dimension == data->nparam+1) {
 268         if (!data->pl_all)
 269             data->pl_all = data->pl;
 270         else {
 271             data->pl_all->combine(*data->pl);
 272             delete data->pl;
 273         }
 274         data->pl = NULL;
 275     }
 276
 277     data->poly = save_poly;
 278     data->monotonicity = save_monotonicity;
 279 }
 280
 281 /*
 282  * Determines the sign of each variable in D.
 283  * Copied from evalue.c
 284  */
 285 static void domain_signs(Polyhedron *D, int *signs)
 286 {
 287     int j, k;
 288
 289     POL_ENSURE_VERTICES(D);
 290     for (j = 0; j < D->Dimension; ++j) {
 291         signs[j] = 0;
 292         for (k = 0; k < D->NbRays; ++k) {
 293             signs[j] = value_sign(D->Ray[k][1+j]);
 294             if (signs[j])
 295                 break;
 296         }
 297     }
 298 }
 299
 300 static evalue *ex2evalue(const ex& poly, const exvector& params, int pos)
 301 {
 302     if (pos >= params.size()) {
 303         evalue *c = ALLOC(evalue);
 304         value_init(c->d);
 305         value_init(c->x.n);
 306         assert(is_a<numeric>(poly));
 307         numeric2value(ex_to<numeric>(poly).numer(), c->x.n);
 308         numeric2value(ex_to<numeric>(poly).denom(), c->d);
 309         return c;
 310     }
 311
 312     evalue *v = evalue_var(pos);
 313     evalue *sum = ex2evalue(poly.coeff(params[pos], poly.degree(params[pos])),
 314                             params, pos+1);
 315     for (int i = poly.degree(params[pos])-1; i >= 0; --i) {
 316         evalue *t = ex2evalue(poly.coeff(params[pos], i), params, pos+1);
 317         emul(v, sum);
 318         eadd(t, sum);
 319         evalue_free(t);
 320     }
 321     evalue_free(v);
 322     return sum;
 323 }
 324
 325 /* Returns true is poly is no better than any of those from begin to end */
 326 static int is_no_better(const ex& poly, const lst::const_iterator& begin,
 327                         const lst::const_iterator& end, const exvector& params,
 328                         int sign, Polyhedron *D,
 329                         int *signs, barvinok_options *options)
 330 {
 331     lst::const_iterator k;
 332     int no_better = 0;
 333
 334     for (k = begin; k != end; ++k) {
 335         ex diff = *k - poly;
 336         diff = diff.expand();
 337         evalue *e = ex2evalue(diff, params, 0);
 338         no_better = has_sign(D, e, sign, signs, options);
 339         evalue_free(e);
 340         if (no_better)
 341             break;
 342     }
 343     return no_better;
 344 }
 345
 346 static void remove_redundants(piecewise_lst *pl, const exvector& params,
 347                                 barvinok_options *options)
 348 {
 349     if (pl->sign == 0)
 350         return;
 351
 352     int *signs = (int *)alloca(sizeof(int) * params.size());
 353
 354     for (int i = 0; i < pl->list.size(); ++i) {
 355         Polyhedron *D = pl->list[i].first;
 356         lst todo = pl->list[i].second;
 357         lst newlist;
 358         lst::const_iterator j, k;
 359
 360         domain_signs(D, signs);
 361
 362         for (j = todo.begin(); j != todo.end(); ++j) {
 363             k = j; ++k;
 364             if (is_no_better(*j, k, todo.end(), params,
 365                                 pl->sign, D, signs, options))
 366                 continue;
 367
 368             if (is_no_better(*j, newlist.begin(), newlist.end(),
 369                                 params, pl->sign, D, signs, options))
 370                 continue;
 371
 372             newlist.append(*j);
 373         }
 374         pl->list[i].second = newlist;
 375     }
 376 }
 377
 378 piecewise_lst *evalue_range_propagation(piecewise_lst *pl_all, const evalue *e,
 379                                       const exvector& params,
 380                                       barvinok_options *options)
 381 {
 382     evalue *e2;
 383     struct range_data data;
 384
 385     if (EVALUE_IS_ZERO(*e))
 386         return pl_all;
 387
 388     assert(value_zero_p(e->d));
 389     assert(e->x.p->type == partition);
 390     assert(e->x.p->size >= 2);
 391     unsigned nparam = params.size();
 392     unsigned dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
 393     unsigned nvars = dim - nparam;
 394
 395     if (nvars == 0) {
 396         assert(0);
 397         return pl_all;
 398     }
 399
 400     data.nparam = nparam;
 401     data.params = params;
 402     data.options = options;
 403     data.pl_all = pl_all;
 404     if (options->bernstein_optimize == BV_BERNSTEIN_MIN)
 405         data.sign = -1;
 406     else
 407         data.sign = 1;
 408     data.test_monotonicity = 1;
 409
 410     e2 = evalue_dup(e);
 411     evalue_split_domains_into_orthants(e2, options->MaxRays);
 412
 413     data.signs = (int *)alloca(sizeof(int) * dim);
 414
 415     for (int i = 0; i < e2->x.p->size/2; ++i) {
 416         Polyhedron *D = EVALUE_DOMAIN(e2->x.p->arr[2*i]);
 417         for (Polyhedron *P = D; P; P = P->next) {
 418             domain_signs(P, data.signs);
 419             propagate_on_domain(P, &e2->x.p->arr[2*i+1], &data);
 420         }
 421     }
 422     evalue_free(e2);
 423     if (data.pl_all) {
 424         data.pl_all->sign = options->bernstein_optimize;
 425         remove_redundants(data.pl_all, params, options);
 426     }
 427     return data.pl_all;
 428 }