bernstein/src/maximize.cpp

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