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