omega/polyfunc.cc

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