bernstein/src/piecewise_lst.cpp

   1 #include <assert.h>
   2 #include <cln/cln.h>
   3
   4 #include <bernstein/bernstein.h>
   5 #include <bernstein/maximize.h>
   6 #include <bernstein/piecewise_lst.h>
   7
   8 using std::cerr;
   9 using std::endl;
  10
  11 using namespace GiNaC;
  12
  13 namespace bernstein {
  14
  15 static void print_lst(std::ostream & os, int sign, lst l)
  16 {
  17     for (int i = 1; i < l.nops(); ++i) {
  18         if (sign > 0)
  19             os << "max(";
  20         else if (sign < 0)
  21             os << "min(";
  22     }
  23     if (l.nops() > 0)
  24         l.op(0).print(print_csrc(os), 0);
  25     for (int i = 1; i < l.nops(); ++i) {
  26         os << ",";
  27         l.op(i).print(print_csrc(os), 0);
  28         if (sign)
  29             os << ")";
  30     }
  31 }
  32
  33 static void printpoly(std::ostream& o, Polyhedron *D, const exvector& p)
  34 {
  35     for (int i = 0; i < D->NbConstraints; ++i) {
  36         int first = 1;
  37         if (i)
  38             o << " && ";
  39         for (int j = 0; j < D->Dimension; ++j) {
  40             if (value_zero_p(D->Constraint[i][1+j]))
  41                 continue;
  42             if (!first)
  43                 o << " + ";
  44             o << VALUE_TO_INT(D->Constraint[i][1+j]) << "*" << p[j];
  45             first = 0;
  46         }
  47         if (value_notzero_p(D->Constraint[i][1+D->Dimension])) {
  48             if (!first)
  49                 o << " + ";
  50             o << VALUE_TO_INT(D->Constraint[i][1+D->Dimension]);
  51         }
  52         if (value_zero_p(D->Constraint[i][0]))
  53             o << " == 0";
  54         else
  55             o << " >= 0";
  56     }
  57 }
  58
  59 static void printdomain(std::ostream& o, Polyhedron *D, const exvector& p)
  60 {
  61     if (!D->next)
  62         printpoly(o, D, p);
  63     else for (Polyhedron *P = D; P; P = P->next) {
  64         o << "(";
  65         printpoly(o, P, p);
  66         o << ")";
  67         if (P->next)
  68             o << " || ";
  69     }
  70 }
  71
  72 std::ostream & operator<< (std::ostream & os, const piecewise_lst & pl)
  73 {
  74     if (pl.list.size() == 1 && universeQ(pl.list[0].first))
  75         print_lst(os, pl.sign, pl.list[0].second);
  76     else {
  77         for (int i = 0; i < pl.list.size(); ++i) {
  78             os << "(";
  79             printdomain(os, pl.list[i].first, pl.vars);
  80             os << ") ? (";
  81             print_lst(os, pl.sign, pl.list[i].second);
  82             os << ") : ";
  83         }
  84         os << "0";
  85     }
  86     return os;
  87 }
  88
  89 static std::vector<guarded_lst> combine(const std::vector<guarded_lst>& one,
  90                                         const std::vector<guarded_lst>& other,
  91                                         const GiNaC::exvector& vars, int sign)
  92 {
  93     Polyhedron *fd;
  94     std::vector<guarded_lst> newlist;
  95     for (int j = 0; j < other.size(); ++j) {
  96         assert(one.size() >= 1);
  97         fd = DomainDifference(other[j].first, one[0].first, 0);
  98         if (!emptyQ(fd))
  99             for (int i = 1; i < one.size(); ++i) {
 100                 Polyhedron *t = fd;
 101                 fd = DomainDifference(fd, one[i].first, 0);
 102                 Domain_Free(t);
 103                 if (emptyQ(fd))
 104                     break;
 105             }
 106         if (emptyQ(fd)) {
 107             Domain_Free(fd);
 108             continue;
 109         }
 110         newlist.push_back(guarded_lst(fd, other[j].second));
 111     }
 112     for (int i = 0; i < one.size(); ++i) {
 113         fd = one[i].first;
 114         for (int j = 0; j < other.size(); ++j) {
 115             Polyhedron *t, *d;
 116             d = DomainIntersection(other[j].first, one[i].first, 0);
 117             if (emptyQ(d)) {
 118                 Domain_Free(d);
 119                 continue;
 120             }
 121             t = fd;
 122             fd = DomainDifference(fd, other[j].first, 0);
 123             if (t != one[i].first)
 124                 Domain_Free(t);
 125             lst list = remove_redundants(d, one[i].second, other[j].second,
 126                                          vars, sign);
 127             newlist.push_back(guarded_lst(d, list));
 128         }
 129         if (!emptyQ(fd))
 130             newlist.push_back(guarded_lst(fd, one[i].second));
 131         else
 132             Domain_Free(fd);
 133         if (fd != one[i].first)
 134             Domain_Free(one[i].first);
 135     }
 136     return newlist;
 137 }
 138
 139 using std::ostream;
 140
 141 ostream & operator<< (ostream & os, const exvector & v)
 142 {
 143     os << "[";
 144     for (int i = 0; i < v.size(); ++i) {
 145         if (i)
 146             os << ", ";
 147         os << v[i];
 148     }
 149     os << "]";
 150     return os;
 151 }
 152
 153 void piecewise_lst::add_guarded_lst(Polyhedron *D, GiNaC::lst coeffs)
 154 {
 155     coeffs = remove_redundants(D, coeffs, coeffs, vars, sign);
 156     assert(coeffs.nops() > 0);
 157     list.push_back(guarded_lst(D, coeffs));
 158 }
 159
 160 piecewise_lst& piecewise_lst::combine(const piecewise_lst& other)
 161 {
 162     assert(vars == other.vars);
 163     assert(sign == other.sign);
 164     list = bernstein::combine(list, other.list, vars, sign);
 165     return *this;
 166 }
 167
 168
 169 void piecewise_lst::maximize()
 170 {
 171     exvector params;
 172     for (int i = 0; i < list.size(); ++i) {
 173         list[i].second = bernstein::maximize(list[i].first, list[i].second, vars);
 174     }
 175 }
 176
 177 void piecewise_lst::minimize()
 178 {
 179     exvector params;
 180     for (int i = 0; i < list.size(); ++i) {
 181         list[i].second = bernstein::minimize(list[i].first, list[i].second, vars);
 182     }
 183 }
 184
 185 void piecewise_lst::simplify_domains(Polyhedron *ctx, unsigned MaxRays)
 186 {
 187     for (int i = 0; i < list.size(); ++i) {
 188         Polyhedron *D = list[i].first;
 189         list[i].first = DomainSimplify(D, ctx, MaxRays);
 190         Domain_Free(D);
 191     }
 192 }
 193
 194 numeric piecewise_lst::evaluate(const exvector& values)
 195 {
 196     Value *v = new Value[values.size()];
 197     numeric result;
 198
 199     for (int i = 0; i < values.size(); ++i) {
 200         value_init(v[i]);
 201         assert(is_a<numeric>(values[i]));
 202         numeric2value(ex_to<numeric>(values[i]), v[i]);
 203     }
 204     result = evaluate(values, values.size(), v);
 205     for (int i = 0; i < values.size(); ++i)
 206         value_clear(v[i]);
 207     delete [] v;
 208     return result;
 209 }
 210
 211 void piecewise_lst::evaluate(int n, Value *v, Value *num, Value *den)
 212 {
 213     exvector values(n);
 214     numeric result;
 215     for (int i = 0; i < values.size(); ++i)
 216         values[i] = value2numeric(v[i]);
 217     result = evaluate(values, values.size(), v);
 218     numeric2value(result.numer(), *num);
 219     numeric2value(result.denom(), *den);
 220 }
 221
 222 numeric piecewise_lst::evaluate(const exvector& values, int n, Value *v)
 223 {
 224     numeric result = 0;
 225
 226     for (int i = 0; i < list.size(); ++i) {
 227         if (!in_domain(list[i].first, v))
 228             continue;
 229         exmap m;
 230         for (int j = 0; j < n; ++j)
 231             m[vars[j]] = values[j];
 232         ex ex_val = list[i].second.subs(m);
 233         assert(is_a<lst>(ex_val));
 234         lst val = ex_to<lst>(ex_val);;
 235         ex opt = val.op(0);
 236         for (int j = 1; j < val.nops(); ++j) {
 237             assert(sign);
 238             if (sign > 0 && val.op(j) > opt)
 239                 opt = val.op(j);
 240             if (sign < 0 && val.op(j) < opt)
 241                 opt = val.op(j);
 242         }
 243         assert(is_a<numeric>(opt));
 244         result = ex_to<numeric>(opt);
 245         break;
 246     }
 247     return result;
 248 }
 249
 250 void piecewise_lst::add(const GiNaC::ex& poly)
 251 {
 252     for (int i = 0; i < list.size(); ++i) {
 253         lst::const_iterator j;
 254         lst new_coeff;
 255         for (j = list[i].second.begin(); j != list[i].second.end(); ++j)
 256             new_coeff.append(*j + poly);
 257         list[i].second = new_coeff;
 258     }
 259 }
 260
 261 int piecewise_lst::is_equal(const piecewise_lst& other) const
 262 {
 263     if (list.size() != other.list.size())
 264         return 0;
 265
 266     for (int i = 0; i < list.size(); ++i) {
 267         int j;
 268         lst l1, l2;
 269         for (j = 0; j < other.list.size(); ++j) {
 270             if (!PolyhedronIncludes(list[i].first, other.list[j].first))
 271                 continue;
 272             if (!PolyhedronIncludes(other.list[j].first, list[i].first))
 273                 continue;
 274             break;
 275         }
 276         if (j >= other.list.size())
 277             return 0;
 278         l1 = list[i].second;
 279         l2 = other.list[j].second;
 280         l1.sort();
 281         l2.sort();
 282         if (!l1.is_equal(l2))
 283             return 0;
 284     }
 285
 286     return 1;
 287 }
 288
 289 }