omega/polyfunc.cc

   1 #include <bernstein/bernstein.h>
   2 #include <bernstein/piecewise_lst.h>
   3 #include <barvinok/bernstein.h>
   4 #include "polyfunc.h"
   5 #include "omega/convert.h"
   6
   7 using namespace GiNaC;
   8 using namespace bernstein;
   9 using namespace std;
  10
  11 void maximize(PolyFunc *polyfunc, Map<Variable_Ref *, GiNaC::ex>& variableMap)
  12 {
  13     varvector vv;
  14     varvector params;
  15     Param_Polyhedron *PP;
  16     exvector exvars;
  17     exvector exparams;
  18
  19     cout << "maximize " << polyfunc->poly << " over ";
  20     polyfunc->domain.simplify();
  21     polyfunc->domain.print(stdout);
  22
  23     Polyhedron *D = relation2Domain(polyfunc->domain, vv, params);
  24     assert(!D->next);
  25     assert(polyfunc->domain.is_set());
  26     int dim = polyfunc->domain.n_set();
  27     assert(D->Dimension == dim + params.size());
  28     Polyhedron *ctx = Universe_Polyhedron(params.size());
  29
  30     for (int i = 0; i < dim; ++i) {
  31         ex var = 0;
  32         foreach_map(vr,Variable_Ref *,s,ex,variableMap, {
  33             if (vr->vid == vv[i]) {
  34                 var = s;
  35                 break;
  36             }
  37         });
  38         assert(var != 0);
  39         exvars.push_back(var);
  40     }
  41     for (int i = dim; i < vv.size(); ++i) {
  42         Global_Var_ID global = vv[i]->get_global_var();
  43         ex var = 0;
  44         foreach_map(vr,Variable_Ref *,s,ex,variableMap, {
  45             if (vr->g == global) {
  46                 var = s;
  47                 break;
  48             }
  49         });
  50         if (var != 0)
  51             exparams.push_back(var);
  52         else
  53             exparams.push_back(symbol(vv[i]->char_name()));
  54     }
  55
  56     PP = Polyhedron2Param_Domain(D, ctx, 0);
  57     piecewise_lst *pl = new piecewise_lst(exparams);
  58     for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
  59         GiNaC::matrix VM = domainVertices(PP, Q, exparams);
  60         lst coeffs = bernsteinExpansion(VM, polyfunc->poly, exvars, exparams);
  61         pl->list.push_back(guarded_lst(Polyhedron_Copy(Q->Domain), coeffs));
  62     }
  63     cout << "coefficients: " << *pl << endl;
  64     pl->maximize();
  65     cout << "maximum: " << *pl << endl;
  66     delete pl;
  67     Param_Polyhedron_Free(PP);
  68
  69     Polyhedron_Free(ctx);
  70     Domain_Free(D);
  71 }