bernstein/src/bernstein.cpp

   1 /*
   2  *      Bernstein Expansion
   3  *
   4  *      - ginac functions
   5  */
   6
   7
   8 #include <iostream>
   9 #include <cln/cln.h>
  10 #include <string>
  11
  12 #include <bernstein/bernstein.h>
  13 #include "polynomial.h"
  14
  15 using namespace std;
  16 using namespace GiNaC;
  17
  18 #define DIGITS 10
  19
  20 namespace bernstein {
  21
  22 static std::string int2String(int n);
  23 static unsigned int findMaxDegree(ex polynomial, const exvector& Vars);
  24 static matrix getAiMatrix(unsigned int nbVert);
  25 static ex getBasis(unsigned int nbVert, matrix &A);
  26 static ex replaceVariablesInPolynomial(const ex &poly, const exvector& V,
  27                                        ex &variables);
  28 static ex powerMonomials(polynomial &poly, matrix &A, unsigned int nbVert,
  29                          unsigned int maxDegree, ex &basis);
  30 static lst getCoefficients(ex hom, unsigned d, const matrix& A);
  31
  32 exvector constructParameterVector(const char * const *param_names, unsigned nbParams)
  33 {
  34         exvector P(nbParams);
  35         for (int i = 0; i < nbParams; ++i) {
  36                 P[i] = symbol(param_names[i]);
  37 #ifdef DEBUG
  38                 cout << "P: " << P[i] << endl;
  39 #endif
  40         }
  41         return P;
  42 }
  43
  44
  45 exvector constructVariableVector(unsigned nbVariables, const char *prefix)
  46 {
  47         exvector V(nbVariables);
  48         for (int i = 0; i < nbVariables; ++i) {
  49                 V[i] = symbol(prefix + int2String(i));
  50 #ifdef DEBUG
  51                 cout << "V: " << V[i] << endl;
  52 #endif
  53         }
  54         return V;
  55 }
  56
  57
  58 numeric value2numeric(const Value v)
  59 {
  60     int sa = v[0]._mp_size;
  61     if (!sa)
  62         return 0;
  63     int abs_sa = sa < 0 ? -sa : sa;
  64     cln::cl_I res = 0;
  65     for (int i = abs_sa-1; i >= 0; --i) {
  66         res = res << GMP_LIMB_BITS;
  67         res = res + v[0]._mp_d[i];
  68     }
  69     return numeric(sa < 0 ? -res : res);
  70 }
  71
  72
  73 matrix domainVertices(Param_Polyhedron *PP, Param_Domain *Q, const exvector& params)
  74 {
  75         Param_Vertices *V;
  76         unsigned nbVertices = 0;
  77         unsigned nbRows, nbColumns;
  78         int v;
  79
  80         assert(PP->nbV > 0);
  81         nbRows = PP->V->Vertex->NbRows;
  82         nbColumns = PP->V->Vertex->NbColumns;
  83
  84         FORALL_PVertex_in_ParamPolyhedron(V, Q, PP)
  85                 ++nbVertices;
  86         END_FORALL_PVertex_in_ParamPolyhedron;
  87
  88         matrix VM(nbVertices, nbRows);          // vertices matrix
  89
  90         v = 0;
  91         FORALL_PVertex_in_ParamPolyhedron(V, Q, PP)
  92                 for (unsigned i = 0; i < nbRows; i++) {
  93                         ex t;
  94                         for (unsigned j = 0; j < nbColumns-2; j++)
  95                                 t += value2numeric(V->Vertex->p[i][j]) * params[j];
  96                         t += value2numeric(V->Vertex->p[i][nbColumns-2]);
  97                         t /= value2numeric(V->Vertex->p[i][nbColumns-1]);
  98 #ifdef DEBUG
  99                         cout << "T: " << t << endl;
 100 #endif
 101                         VM(v, i) = t;
 102                 }
 103                 ++v;
 104         END_FORALL_PVertex_in_ParamPolyhedron;
 105
 106         return VM;
 107 }
 108
 109 /*
 110  * Do the Bernstein Expansion
 111  *
 112  *      vert: vertices matrix
 113  *      poly: polynomial expression
 114  *      vars: vector of variables
 115  *      Params: vector of parameter
 116  */
 117 lst bernsteinExpansion(const matrix& vert, const ex& poly, const exvector& vars,
 118                        const exvector& Params)
 119 {
 120         if (is_exactly_a<lst>(poly)) {
 121                 lst polylst = ex_to<lst>(poly);
 122                 lst::const_iterator j = polylst.begin();
 123
 124                 lst coeff = bernsteinExpansion(vert, *j, vars, Params);
 125
 126                 for (++j; j != polylst.end(); ++j) {
 127                         lst::const_iterator k;
 128                         lst new_coeff = bernsteinExpansion(vert, *j, vars, Params);
 129                         for (k = new_coeff.begin(); k != new_coeff.end(); ++k)
 130                             coeff.append(*k);
 131                         coeff.sort().unique();
 132                 }
 133                 return coeff;
 134         }
 135
 136         unsigned maxDegree = findMaxDegree(poly, vars);
 137         matrix A = getAiMatrix(vert.rows());
 138
 139 #ifdef DEBUG
 140         cout << endl << "Polynomial: " << poly << endl << endl;
 141         cout << vert << endl << endl;
 142         cout << A << endl << endl;
 143 #endif
 144
 145         // obtain the variables value on the basis and replace
 146         ex substitutions = evalm(A * vert);
 147         ex polynom = replaceVariablesInPolynomial(poly, vars, substitutions);
 148
 149 #ifdef DEBUG
 150         cout << variables << endl << endl;
 151         cout << "Preliminar Expansion: " << polynom << endl << endl;
 152 #endif
 153
 154         ex basis = getBasis(vert.rows(), A);
 155
 156 #ifdef DEBUG
 157         cout << "Basis: " << basis << endl<< endl;
 158 #endif
 159
 160         // monomials to n degree
 161         polynomial p(polynom);
 162         ex maxDegreePolynomial = powerMonomials(p, A, vert.rows(), maxDegree, basis);
 163
 164         return getCoefficients(maxDegreePolynomial, maxDegree, A);
 165 }
 166
 167 /*
 168  * Construct A_i matrix
 169  *
 170  *      nbVert: number of vertices
 171  */
 172 matrix getAiMatrix(unsigned int nbVert)
 173 {
 174         matrix A(1, nbVert);            // a_i matrix
 175         for(unsigned int i = 0; i < nbVert; i++) {
 176                 A(0,i) = symbol("a" + int2String(i));
 177 #ifdef DEBUG
 178                 cout << "A: " << A(0,i) << endl;
 179 #endif
 180         }
 181         return A;
 182 }
 183
 184
 185 /*
 186  * Construct the basis
 187  *
 188  *      A: a_i matrix
 189  *      nbVert: number of vertices
 190  */
 191 ex getBasis(unsigned int nbVert, matrix &A)
 192 {
 193         ex basis;
 194         for(unsigned int i = 0; i < nbVert; i++) {
 195                 basis += A(0,i);
 196         }
 197         return basis;
 198 }
 199
 200
 201 /*
 202  * Finds the less than maxDegree monomials and multiply them
 203  * by the basis
 204  *
 205  *      polynomial: original polynomial
 206  *      A: a_i matrix
 207  *      nbVert: number of vertices
 208  *      maxDegree: maximal polynomial multidegree
 209  *      basis: basis of the polytope
 210  */
 211 ex powerMonomials(polynomial &poly, matrix &A, unsigned int nbVert
 212                   , unsigned int maxDegree, ex &basis)
 213 {
 214         ex maxDegreePolynomial;
 215 #ifdef DEBUG
 216         cout << "- Degree --------------------------------------" << endl;
 217 #endif
 218         for (size_t i = 0; i != poly.nbTerms(); ++i) {
 219                 unsigned int degree = 0;                // degree of the monomial
 220
 221                 for(unsigned int j = 0; j < nbVert; j++) {
 222                         degree += poly.term(i).degree(A(0,j));
 223                 }
 224 #ifdef DEBUG
 225                 cout << poly.term(i) << " Degree: " << degree;
 226 #endif
 227                 if(degree < maxDegree) {
 228                         ex degreeUp = poly.term(i) * pow(basis, maxDegree - degree);
 229 #ifdef DEBUG
 230                         cout << "   --> New Term: " <<  degreeUp.expand();
 231 #endif
 232                         maxDegreePolynomial += degreeUp.expand();
 233                 } else {
 234                         maxDegreePolynomial += poly.term(i);
 235                 }
 236 #ifdef DEBUG
 237                 cout << endl << "-----------------------------------------------" << endl;
 238 #endif
 239
 240         }
 241 #ifdef DEBUG
 242         cout << endl << "Final Expansion: " << maxDegreePolynomial << endl << endl;
 243 #endif
 244         return maxDegreePolynomial;
 245 }
 246
 247
 248
 249 /*
 250  * Finds the coefficients of the polynomial in terms of the Bernstein basis
 251  *
 252  *      hom: homogeneous polynomial of degree d
 253  *      A: a_i matrix
 254  *
 255  *      For each monomial of multi-degree (k[0], k[1], ..., k[n-1])
 256  *      we divide the corresponding coefficient by the multinomial
 257  *      coefficient d!/k[0]! k[1]! ... k[n-1]!
 258  *
 259  *      The code looks a bit complicated because it's an iterative
 260  *      implementation of a recursive procedure.
 261  *      For each variable from 0 to n-1, we loop over the possible
 262  *      powers: 0..d for the first; 0..d-k[0]=left[0] for the second; ...
 263  *      Just for fun, we loop through these in the opposite direction
 264  *      only for the first variable.
 265  *
 266  *      c[i] contains the coefficient of the selected powers of the first i+1 vars
 267  *      multinom[i] contains the partial multinomial coefficient.
 268  */
 269 lst getCoefficients(ex hom, unsigned d, const matrix& A)
 270 {
 271         lst coeff;
 272         int n = A.cols();
 273
 274         /* we should probably notice sooner that there is just one vertex */
 275         if (n == 1) {
 276             coeff.append(hom.coeff(A(0, 0), d));
 277             return coeff;
 278         }
 279
 280         ex c[n];
 281         int left[n];
 282         int k[n];
 283         numeric multinom[n];
 284         assert(n >= 2);
 285
 286         multinom[0] = 1;
 287         for (k[0] = d; k[0] >= 0; --k[0]) {
 288                 c[0] = hom.coeff(A(0, 0), k[0]);
 289                 left[0] = d - k[0];
 290                 int i = 1;
 291                 k[i] = -1;
 292                 multinom[i] = multinom[0];
 293                 while (i > 0) {
 294                         if (i == n-1) {
 295                                 for (int j = 2; j <= left[i-1]; ++j)
 296                                         multinom[i] /= j;
 297                                 coeff.append(c[i-1].coeff(A(0, i), left[i-1]) /
 298                                                 multinom[i]);
 299                                 --i;
 300                                 continue;
 301                         }
 302                         if (k[i] >= left[i-1]) {
 303                                 --i;
 304                                 continue;
 305                         }
 306                         ++k[i];
 307                         if (k[i])
 308                                 multinom[i] /= k[i];
 309                         c[i] = c[i-1].coeff(A(0, i), k[i]);
 310                         left[i] = left[i-1] - k[i];
 311                         k[i+1] = -1;
 312                         multinom[i+1] = multinom[i];
 313                         ++i;
 314                 }
 315                 multinom[0] *= k[0];
 316         }
 317         return coeff.sort().unique();
 318 }
 319
 320
 321 // replace the variables in the polynomial
 322 ex replaceVariablesInPolynomial(const ex &poly, const exvector &vars,
 323                                 ex &substitutions)
 324 {
 325         lst replace;
 326
 327         for(unsigned int i = 0; i < vars.size(); i++) {
 328 #ifdef DEBUG
 329                 cout << "Replacing: " << vars[i] << " by " << substitutions[i] << endl;
 330 #endif
 331                 replace.append(vars[i] == substitutions[i]);
 332         }
 333         ex polyRepl = poly.subs(replace);
 334
 335         return(polyRepl.expand());
 336 }
 337
 338
 339 /* Converts int n to string */
 340 string int2String(int n)
 341 {
 342         char numeroVariable[DIGITS];
 343         snprintf(numeroVariable, DIGITS, "%d", n);
 344         string nroV(numeroVariable);
 345
 346         return nroV;
 347 }
 348
 349
 350 /*
 351  * Find the maximum multi-degree of the polinomial
 352  *
 353  *      polynomial: polynomial
 354  *      Vars: variables matrix
 355  *      nbVar: number of variables
 356  */
 357 static unsigned int findMaxDegree(ex polynomial, const exvector& Vars, int pos)
 358 {
 359         if (pos >= Vars.size())
 360                 return 0;
 361
 362         unsigned max = 0;
 363         for (int i = 0; i <= polynomial.degree(Vars[pos]); ++i) {
 364                 unsigned degree = i + findMaxDegree(polynomial.coeff(Vars[pos], i),
 365                                                     Vars, pos+1);
 366                 if (degree > max)
 367                         max = degree;
 368         }
 369         return max;
 370 }
 371
 372 unsigned int findMaxDegree(ex polynomial, const exvector& Vars)
 373 {
 374     return findMaxDegree(polynomial, Vars, 0);
 375 }
 376
 377 }