bernstein/src/maximize.cpp

   1 #include <assert.h>
   2 #include <algorithm>
   3 #include <bernstein/bernstein.h>
   4 #include <bernstein/maximize.h>
   5
   6 using namespace GiNaC;
   7
   8 using std::cerr;
   9 using std::endl;
  10
  11 namespace bernstein {
  12
  13 static void domainVerticesAndRays(Polyhedron *P, GiNaC::matrix& VM, GiNaC::matrix& RM)
  14 {
  15     POL_ENSURE_VERTICES(P);
  16     int nv = 0, nr = 0;
  17     for (int i = 0; i < P->NbRays; ++i) {
  18         if (value_zero_p(P->Ray[i][0]))
  19             nr += 2;
  20         else if (value_zero_p(P->Ray[i][1+P->Dimension]))
  21             ++nr;
  22         else
  23             ++nv;
  24     }
  25     VM = GiNaC::matrix(nv, P->Dimension);
  26     RM = GiNaC::matrix(nr, P->Dimension);
  27     nv = nr = 0;
  28     for (int i = 0; i < P->NbRays; ++i) {
  29         if (value_zero_p(P->Ray[i][1+P->Dimension])) {
  30             for (int j = 0; j < P->Dimension; ++j)
  31                 RM(nr,j) = value2numeric(P->Ray[i][1+j]);
  32             ++nr;
  33             if (value_notzero_p(P->Ray[i][0]))
  34                 continue;
  35             for (int j = 0; j < P->Dimension; ++j)
  36                 RM(nr,j) = -RM(nr-1,j);
  37             ++nr;
  38             continue;
  39         }
  40         for (int j = 0; j < P->Dimension; ++j)
  41             VM(nv,j) = value2numeric(P->Ray[i][1+j]) /
  42                                         value2numeric(P->Ray[i][1+P->Dimension]);
  43         ++nv;
  44     }
  45 }
  46
  47 static void find_lst_minmax(const lst& l, ex& m, ex& M)
  48 {
  49     lst::const_iterator i = l.begin();
  50     m = M = *i;
  51     for (++i; i != l.end(); ++i) {
  52         if (*i > M)
  53             M = *i;
  54         else if (*i < m)
  55             m = *i;
  56     }
  57 }
  58
  59 static bool is_linear_constraint(ex constraint, const exvector& vars, matrix& linear)
  60 {
  61     linear = matrix(vars.size()+1, 1);
  62     ex cst = constraint;
  63     for (int i = 0; i < vars.size(); ++i) {
  64         if (cst.degree(vars[i]) > 1)
  65             return false;
  66         ex c = cst.coeff(vars[i], 1);
  67         if (!is_a<numeric>(c))
  68             return false;
  69         linear(i,0) = c;
  70         cst = cst.coeff(vars[i], 0);
  71     }
  72     if (!is_a<numeric>(cst))
  73         return false;
  74     linear(vars.size(),0) = cst;
  75     return true;
  76 }
  77
  78 static void find_linear_minmax(Polyhedron *D, matrix constraint, ex& m, ex& M)
  79 {
  80     // we only care about the sign anyway
  81     M = -1; m = 1;
  82     for (Polyhedron *P = D; P; P = P->next) {
  83         POL_ENSURE_VERTICES(P);
  84         for (int i = 0; i < P->NbRays; ++i) {
  85             ex sum = 0;
  86             assert(value_one_p(P->Ray[i][0]));
  87             for (int j = 0; j < P->Dimension+1; ++j)
  88                 sum += value2numeric(P->Ray[i][1+j]) * constraint(j,0);
  89             if (sum > M)
  90                 M = sum;
  91             if (sum < m)
  92                 m = sum;
  93         }
  94     }
  95 }
  96
  97 enum poly_sign {
  98     zero,
  99     positive,
 100     negative,
 101     unknown,
 102 };
 103
 104 static poly_sign polynomial_sign(ex poly, Polyhedron *D, const GiNaC::matrix& VM,
 105                                  const GiNaC::matrix& RM, const exvector& vars,
 106                                  bool expensive_tests);
 107
 108 /* bernstein requires a bounded domain, so if the domain is unbounded,
 109  * we need a different trick.
 110  * If a parameter N has a lowerbound l, then we write the polynomial as
 111  * q(N,M) * (N-l) + r(M)
 112  * if the signs of q and r can be determined and if they are equal,
 113  * then the whole polynomial has the same sign.
 114  */
 115 static poly_sign polynomial_sign_with_rays(ex poly, Polyhedron *D,
 116                                  const GiNaC::matrix& VM,
 117                                  const GiNaC::matrix& RM, const exvector& vars,
 118                                  bool expensive_tests)
 119 {
 120     int i;
 121     for (i = 0; i < vars.size(); ++i)
 122         if (poly.degree(vars[i]) > 0)
 123             break;
 124     assert(i < vars.size());
 125     ex min = VM(0,i);
 126     for (int j = 1; j < VM.rows(); ++j)
 127         if (VM(j,i) < min)
 128             min = VM(j,i);
 129     if (min < 0)
 130         return unknown;
 131     for (int j = 0; j < RM.rows(); ++j)
 132         if (RM(j,i) < 0)
 133             return unknown;
 134     int d = poly.degree(vars[i]);
 135     ex cum = poly.coeff(vars[i],d);
 136     ex q = cum;
 137
 138     for (--d; d >= 1; --d) {
 139         cum *= min;
 140         cum += poly.coeff(vars[i],d);
 141         q *= vars[i];
 142         q += cum;
 143     }
 144     ex r = cum * min + poly.coeff(vars[i], 0);
 145     poly_sign s1 = polynomial_sign(q.expand(), D, VM, RM, vars, expensive_tests);
 146     if (s1 == unknown)
 147         return unknown;
 148     poly_sign s2 = polynomial_sign(r, D, VM, RM, vars, expensive_tests);
 149     if (s2 == unknown)
 150         return unknown;
 151     if (s1 == zero)
 152         return s2;
 153     if (s2 == zero)
 154         return s1;
 155     if (s1 != s2)
 156         return unknown;
 157     return s1;
 158 }
 159
 160 static ex substitute_equalities(ex poly, Polyhedron *D, const exvector& vars)
 161 {
 162     if (D->NbEq == 0)
 163         return poly;
 164     lst replacements;
 165     for (int i = 0; i < D->NbEq; ++i) {
 166         int pos = First_Non_Zero(D->Constraint[i]+1, D->Dimension);
 167         assert(pos != -1);
 168         ex den = -value2numeric(D->Constraint[i][1+pos]);
 169         ex s = value2numeric(D->Constraint[i][1+D->Dimension]) / den;
 170         for (int j = pos+1; j < D->Dimension; ++j) {
 171             if (value_zero_p(D->Constraint[i][1+j]))
 172                 continue;
 173             s += value2numeric(D->Constraint[i][1+j]) * vars[j] / den;
 174         }
 175         replacements.append(vars[pos] == s);
 176     }
 177     poly = poly.subs(replacements).expand();
 178     return poly;
 179 }
 180
 181 static poly_sign polynomial_sign(ex poly, Polyhedron *D, const GiNaC::matrix& VM,
 182                                  const GiNaC::matrix& RM, const exvector& vars,
 183                                  bool expensive_tests)
 184 {
 185     exvector params;
 186     ex minc, maxc;
 187     matrix linear;
 188     poly = substitute_equalities(poly, D, vars);
 189     if (is_a<numeric>(poly))
 190         minc = maxc = poly;
 191     else if (is_linear_constraint(poly, vars, linear))
 192         find_linear_minmax(D, linear, minc, maxc);
 193     else if (expensive_tests && RM.rows() == 0) {
 194         lst coeffs = bernsteinExpansion(VM, poly, vars, params);
 195         find_lst_minmax(coeffs, minc, maxc);
 196     } else
 197         return polynomial_sign_with_rays(poly, D, VM, RM, vars, expensive_tests);
 198     if (maxc <= 0 && minc >= 0)
 199         return zero;
 200     if (maxc <= 0)
 201         return negative;
 202     if (minc >= 0)
 203         return positive;
 204     return unknown;
 205 }
 206
 207 /* Combine list1 and list2 into a single list with the redundant elements
 208  * removed according to sign.  If sign = 1, obviously smaller elements are
 209  * removed; if sign = -1, obviously bigger elements are removed; if sign = 0,
 210  * no elements are removed.
 211  * If list1 and list2 are the same, the whole list is checked for redundant
 212  * elements.  If the lists are different, the individual lists are assumed
 213  * to have no redundant elements.
 214  */
 215 GiNaC::lst remove_redundants(Polyhedron *domain, GiNaC::lst list1, GiNaC::lst list2,
 216                              const GiNaC::exvector& vars, int sign)
 217 {
 218     bool same_list = list1 == list2;
 219
 220     lst::const_iterator j, k;
 221     if (sign == 0) {
 222         if (same_list)
 223             return list1;
 224         lst list = list1;
 225         for (k = list2.begin(); k != list2.end(); ++k)
 226             list.append(*k);
 227         return list.sort().unique();
 228     }
 229
 230     int sign_better = sign > 0 ? positive : negative;
 231     int sign_worse = sign > 0 ? negative : positive;
 232
 233     GiNaC::matrix VM, RM;
 234     domainVerticesAndRays(domain, VM, RM);
 235     lst newlist;
 236     lst removed;
 237     for (j = list1.begin(); j != list1.end(); ++j) {
 238         if (same_list && find(removed.begin(), removed.end(), *j) != removed.end())
 239             continue;
 240         bool needed = true;
 241         if (same_list) {
 242             k = j; ++k;
 243         } else
 244             k = list2.begin();
 245         for (; k != (same_list ? list1.end() : list2.end()); ++k) {
 246             if (find(removed.begin(), removed.end(), *k) != removed.end())
 247                 continue;
 248             ex diff = *j - *k;
 249             poly_sign s = polynomial_sign(diff, domain, VM, RM, vars, false);
 250             if (s == zero || s == sign_worse)
 251                 needed = false;
 252             else if (s == sign_better)
 253                 removed.append(*k);
 254             if (!needed)
 255                 break;
 256         }
 257         if (needed)
 258             newlist.append(*j);
 259     }
 260     if (!same_list)
 261         for (k = list2.begin(); k != list2.end(); ++k) {
 262             if (find(removed.begin(), removed.end(), *k) != removed.end())
 263                 continue;
 264             newlist.append(*k);
 265         }
 266     return newlist;
 267 }
 268
 269 GiNaC::lst maximize(Polyhedron *domain, GiNaC::lst coeffs,
 270                     const GiNaC::exvector& vars)
 271 {
 272     GiNaC::matrix VM, RM;
 273     domainVerticesAndRays(domain, VM, RM);
 274     lst::const_iterator j, k;
 275     lst newlist;
 276     lst removed;
 277     for (j = coeffs.begin(); j != coeffs.end(); ++j) {
 278         if (find(removed.begin(), removed.end(), *j) != removed.end())
 279             continue;
 280         bool needed = true;
 281         k = j; ++k;
 282         for (; k != coeffs.end(); ++k) {
 283             if (find(removed.begin(), removed.end(), *k) != removed.end())
 284                 continue;
 285             ex diff = *j - *k;
 286             poly_sign s = polynomial_sign(diff, domain, VM, RM, vars, true);
 287             if (s == zero || s == negative)
 288                 needed = false;
 289             else if (s == positive)
 290                 removed.append(*k);
 291             if (!needed)
 292                 break;
 293         }
 294         if (needed)
 295             newlist.append(*j);
 296     }
 297     return newlist;
 298 }
 299
 300 GiNaC::lst minimize(Polyhedron *domain, GiNaC::lst coeffs,
 301                     const GiNaC::exvector& vars)
 302 {
 303     lst negcoeffs, negresult, result;
 304     lst::const_iterator j;
 305     for (j = coeffs.begin(); j != coeffs.end(); ++j)
 306         negcoeffs.append(-*j);
 307     negresult = maximize(domain, negcoeffs, vars);
 308     for (j = negresult.begin(); j != negresult.end(); ++j)
 309         result.append(-*j);
 310     return result;
 311 }
 312
 313 }