range.cc

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