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 "ex_convert.h"
   7 #include "polyfunc.h"
   8 #include "param_util.h"
   9 #include "omega/convert.h"
  10
  11 using namespace GiNaC;
  12 using namespace bernstein;
  13 using namespace std;
  14
  15 #define ALLOC(type) (type*)malloc(sizeof(type))
  16
  17 static void extract_vars(Map<Variable_Ref *, GiNaC::ex>& variableMap,
  18                          const varvector &vv, exvector &exvars)
  19 {
  20     for (int i = 0; i < vv.size(); ++i) {
  21         ex var = 0;
  22         foreach_map(vr,Variable_Ref *,s,ex,variableMap, {
  23             if (vr->vid == vv[i]) {
  24                 var = s;
  25                 break;
  26             }
  27         });
  28         assert(var != 0);
  29         exvars.push_back(var);
  30     }
  31 }
  32
  33 static void extract_params(Map<Variable_Ref *, GiNaC::ex>& variableMap,
  34                              const varvector &params, exvector &exparams)
  35 {
  36     for (int i = 0; i < params.size(); ++i) {
  37         Global_Var_ID global = params[i]->get_global_var();
  38         ex var = 0;
  39         foreach_map(vr,Variable_Ref *,s,ex,variableMap, {
  40             if (vr->g == global) {
  41                 var = s;
  42                 break;
  43             }
  44         });
  45         if (var != 0)
  46             exparams.push_back(var);
  47         else
  48             exparams.push_back(symbol(params[i]->char_name()));
  49     }
  50 }
  51
  52 void maximize(PolyFunc *polyfunc, Map<Variable_Ref *, GiNaC::ex>& variableMap)
  53 {
  54     varvector vv;
  55     varvector params;
  56     Param_Polyhedron *PP;
  57     exvector exvars;
  58     exvector exparams;
  59     struct barvinok_options *options = barvinok_options_new_with_defaults();
  60
  61     cout << "maximize " << polyfunc->poly << " over ";
  62     polyfunc->domain.simplify();
  63     polyfunc->domain.print(stdout);
  64
  65     Polyhedron *D = relation2Domain(polyfunc->domain, vv, params,
  66                                         options->MaxRays);
  67     assert(!D->next);
  68     assert(polyfunc->domain.is_set());
  69     int dim = polyfunc->domain.n_set();
  70     assert(D->Dimension == dim + params.size());
  71     Polyhedron *ctx = Universe_Polyhedron(params.size());
  72
  73     extract_vars(variableMap, vv, exvars);
  74     extract_params(variableMap, params, exparams);
  75
  76     PP = Polyhedron2Param_Polyhedron(D, ctx, options);
  77     piecewise_lst *pl = new piecewise_lst(exparams);
  78     for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
  79         GiNaC::matrix VM = domainVertices(PP, Q, exparams);
  80         lst coeffs = bernsteinExpansion(VM, polyfunc->poly, exvars, exparams);
  81         pl->list.push_back(guarded_lst(Polyhedron_Copy(Q->Domain), coeffs));
  82     }
  83     cout << "coefficients: " << *pl << endl;
  84     pl->maximize();
  85     cout << "maximum: " << *pl << endl;
  86     delete pl;
  87     Param_Polyhedron_Free(PP);
  88
  89     Polyhedron_Free(ctx);
  90     Domain_Free(D);
  91     barvinok_options_free(options);
  92 }
  93
  94 evalue *summate(PolyFunc *polyfunc, Map<Variable_Ref *, GiNaC::ex>& variableMap)
  95 {
  96     varvector vv;
  97     varvector params;
  98     exvector exvars;
  99     exvector exparams;
 100     struct barvinok_options *options = barvinok_options_new_with_defaults();
 101
 102     cout << "sum " << polyfunc->poly << " over ";
 103     polyfunc->domain.simplify();
 104     polyfunc->domain.print(stdout);
 105
 106     Polyhedron *D = relation2Domain(polyfunc->domain, vv, params,
 107                                         options->MaxRays);
 108     assert(!D->next);
 109     assert(polyfunc->domain.is_set());
 110     int dim = polyfunc->domain.n_set();
 111     assert(D->Dimension == dim + params.size());
 112
 113     extract_vars(variableMap, vv, exvars);
 114     extract_params(variableMap, params, exparams);
 115     for (int i = 0; i < exparams.size(); ++i)
 116         exvars.push_back(exparams[i]);
 117
 118     evalue *e = ALLOC(evalue);
 119     value_init(e->d);
 120     value_set_si(e->d, 0);
 121     e->x.p = new_enode(::partition, 2, D->Dimension);
 122     EVALUE_SET_DOMAIN(e->x.p->arr[0], D);
 123     value_clear(e->x.p->arr[1].d);
 124     e->x.p->arr[1] = *ex2evalue(polyfunc->poly, exvars);
 125
 126     evalue *sum = barvinok_summate(e, dim, options);
 127     evalue_free(e);
 128
 129     barvinok_options_free(options);
 130     return sum;
 131 }