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] = u;
 149         evalue_substitute(pos, subs);
 150
 151         subs[0] = l;
 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         propagate_on_domain(T, pos, data);
 168         Polyhedron_Free(T);
 169     }
 170
 171     evalue_free(pos);
 172 }
 173
 174 static int has_sign(Polyhedron *D, evalue *poly, int sign,
 175                     int *signs, struct barvinok_options *options)
 176 {
 177     struct range_data data_m;
 178     int i;
 179     lst l;
 180
 181     data_m.options = options;
 182     data_m.nparam = 0;
 183     data_m.test_monotonicity = 0;
 184     data_m.signs = signs;
 185
 186     data_m.pl_all = NULL;
 187     data_m.sign = -sign;
 188     propagate_on_domain(D, poly, &data_m);
 189
 190     assert(data_m.pl_all->list.size() == 1);
 191     assert(universeQ(data_m.pl_all->list[0].first));
 192     l = data_m.pl_all->list[0].second;
 193
 194     for (i = 0; i < l.nops(); ++i) {
 195         if (is_exactly_a<fail>(l.op(i)))
 196             break;
 197         if (sign * l.op(i) < 0)
 198             break;
 199     }
 200     delete data_m.pl_all;
 201
 202     return i == l.nops();
 203 }
 204
 205 /* Returns 1 if poly is monotonically increasing,
 206  *        -1 if poly is monotonically decreasing,
 207  *         0 if no conclusion.
 208  */
 209 static int monotonicity(Polyhedron *D, const evalue *poly,
 210                         struct range_data *data)
 211 {
 212     evalue **subs;
 213     evalue *diff;
 214     Value one;
 215     int result = 0;
 216
 217     value_init(one);
 218     value_set_si(one, 1);
 219
 220     subs = ALLOCN(evalue *, D->Dimension);
 221     for (int i = 0; i < D->Dimension; ++i)
 222         subs[i] = evalue_var(i);
 223
 224     evalue_add_constant(subs[0], one);
 225
 226     diff = evalue_dup(poly);
 227     evalue_substitute(diff, subs);
 228
 229     for (int i = 0; i < D->Dimension; ++i)
 230         evalue_free(subs[i]);
 231     free(subs);
 232
 233     evalue_negate(diff);
 234     eadd(poly, diff);
 235     reduce_evalue(diff);
 236     evalue_negate(diff);
 237
 238     value_clear(one);
 239
 240     if (has_sign(D, diff, 1, data->signs, data->options))
 241         result = 1;
 242     else if (has_sign(D, diff, -1, data->signs, data->options))
 243         result = -1;
 244
 245     evalue_free(diff);
 246     return result;
 247 }
 248
 249 static void propagate_on_domain(Polyhedron *D, const evalue *poly,
 250                                 struct range_data *data)
 251 {
 252     const evalue *save_poly = data->poly;
 253     int save_monotonicity = data->monotonicity;
 254
 255     assert(D->Dimension > data->nparam);
 256
 257     if (data->test_monotonicity)
 258         data->monotonicity = monotonicity(D, poly, data);
 259     else
 260         data->monotonicity = 0;
 261
 262     if (D->Dimension == data->nparam+1)
 263         data->pl = new piecewise_lst(data->params);
 264
 265     data->poly = poly;
 266     for_each_lower_upper_bound(D, NULL, range_cb, data);
 267
 268     if (D->Dimension == data->nparam+1) {
 269         if (!data->pl_all)
 270             data->pl_all = data->pl;
 271         else {
 272             data->pl_all->combine(*data->pl);
 273             delete data->pl;
 274         }
 275         data->pl = NULL;
 276     }
 277
 278     data->poly = save_poly;
 279     data->monotonicity = save_monotonicity;
 280 }
 281
 282 /*
 283  * Determines the sign of each variable in D.
 284  * Copied from evalue.c
 285  */
 286 static void domain_signs(Polyhedron *D, int *signs)
 287 {
 288     int j, k;
 289
 290     POL_ENSURE_VERTICES(D);
 291     for (j = 0; j < D->Dimension; ++j) {
 292         signs[j] = 0;
 293         for (k = 0; k < D->NbRays; ++k) {
 294             signs[j] = value_sign(D->Ray[k][1+j]);
 295             if (signs[j])
 296                 break;
 297         }
 298     }
 299 }
 300
 301 /* Returns true is poly is no better than any of those from begin to end */
 302 static int is_no_better(const ex& poly, const lst::const_iterator& begin,
 303                         const lst::const_iterator& end, const exvector& params,
 304                         int sign, Polyhedron *D,
 305                         int *signs, barvinok_options *options)
 306 {
 307     lst::const_iterator k;
 308     int no_better = 0;
 309
 310     for (k = begin; k != end; ++k) {
 311         ex diff = *k - poly;
 312         diff = diff.expand();
 313         evalue *e = ex2evalue(diff, params);
 314         no_better = has_sign(D, e, sign, signs, options);
 315         evalue_free(e);
 316         if (no_better)
 317             break;
 318     }
 319     return no_better;
 320 }
 321
 322 static void remove_redundants(piecewise_lst *pl, const exvector& params,
 323                                 barvinok_options *options)
 324 {
 325     if (pl->sign == 0)
 326         return;
 327
 328     int *signs = (int *)alloca(sizeof(int) * params.size());
 329
 330     for (int i = 0; i < pl->list.size(); ++i) {
 331         Polyhedron *D = pl->list[i].first;
 332         lst todo = pl->list[i].second;
 333         lst newlist;
 334         lst::const_iterator j, k;
 335
 336         domain_signs(D, signs);
 337
 338         for (j = todo.begin(); j != todo.end(); ++j) {
 339             k = j; ++k;
 340             if (is_no_better(*j, k, todo.end(), params,
 341                                 pl->sign, D, signs, options))
 342                 continue;
 343
 344             if (is_no_better(*j, newlist.begin(), newlist.end(),
 345                                 params, pl->sign, D, signs, options))
 346                 continue;
 347
 348             newlist.append(*j);
 349         }
 350         pl->list[i].second = newlist;
 351     }
 352 }
 353
 354 piecewise_lst *evalue_range_propagation(piecewise_lst *pl_all, const evalue *e,
 355                                       const exvector& params,
 356                                       barvinok_options *options)
 357 {
 358     evalue *e2;
 359     struct range_data data;
 360
 361     if (EVALUE_IS_ZERO(*e))
 362         return pl_all;
 363
 364     assert(value_zero_p(e->d));
 365     assert(e->x.p->type == partition);
 366     assert(e->x.p->size >= 2);
 367     unsigned nparam = params.size();
 368     unsigned dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
 369     unsigned nvars = dim - nparam;
 370
 371     if (nvars == 0) {
 372         assert(0);
 373         return pl_all;
 374     }
 375
 376     data.nparam = nparam;
 377     data.params = params;
 378     data.options = options;
 379     data.pl_all = pl_all;
 380     if (options->bernstein_optimize == BV_BERNSTEIN_MIN)
 381         data.sign = -1;
 382     else
 383         data.sign = 1;
 384     data.test_monotonicity = 1;
 385
 386     e2 = evalue_dup(e);
 387     evalue_split_domains_into_orthants(e2, options->MaxRays);
 388
 389     data.signs = (int *)alloca(sizeof(int) * dim);
 390
 391     for (int i = 0; i < e2->x.p->size/2; ++i) {
 392         Polyhedron *D = EVALUE_DOMAIN(e2->x.p->arr[2*i]);
 393         for (Polyhedron *P = D; P; P = P->next) {
 394             domain_signs(P, data.signs);
 395             propagate_on_domain(P, &e2->x.p->arr[2*i+1], &data);
 396         }
 397     }
 398     evalue_free(e2);
 399     if (data.pl_all) {
 400         data.pl_all->sign = options->bernstein_optimize;
 401         remove_redundants(data.pl_all, params, options);
 402     }
 403     return data.pl_all;
 404 }