bernstein/src/main.cpp

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3
   4 #include <ginac/ginac.h>
   5 #include <gmp.h>
   6 #include "polylib++.h"
   7
   8 #include <bernstein/bernstein.h>
   9 #include <bernstein/maximize.h>
  10
  11 #define MAXRAYS 1000
  12
  13 using namespace std;
  14 using namespace GiNaC;
  15 using namespace bernstein;
  16
  17 static ex readPolynomial(const exvector& vars, const exvector& params);
  18 static void printCoefficients(lst coeffs);
  19 static int printMaxMinCoefficient(Polyhedron *VD, lst coeffs, const exvector& Params);
  20
  21 /* main function */
  22 int main(void) {
  23         Matrix *a, *b;
  24         Polyhedron *A, *B;              // initial matrices
  25         char **param_name;      // name of the parameters
  26         exvector params, vars;
  27         ex polynomial;
  28
  29         Param_Polyhedron *PP;
  30         Param_Domain   *Q;
  31
  32         unsigned int nb_param, nb_var;
  33
  34         printf("\n===============================================\n");
  35
  36         a = Matrix_Read();
  37         A = Constraints2Polyhedron(a,200);
  38
  39         b = Matrix_Read();
  40         B = Constraints2Polyhedron(b,200);
  41
  42         /* Set the number of parameters */
  43         nb_param        = B->Dimension;
  44         nb_var          = A->Dimension - B->Dimension;
  45
  46         /* Read the name of the parameters */
  47         param_name = Read_ParamNames(stdin, nb_var+nb_param);
  48         vars = constructParameterVector(param_name, nb_var);
  49         params = constructParameterVector(param_name+nb_var, nb_param);
  50
  51         polynomial = readPolynomial(vars, params);
  52
  53         Matrix_Free(a);
  54         Matrix_Free(b);
  55
  56         /* Find the parametrized domains */
  57         PP = Polyhedron2Param_Domain(A,B,MAXRAYS);
  58         for(Q=PP->D;Q;Q=Q->next) {
  59                 Polyhedron *VD;
  60                 lst coeffs;
  61
  62                 printf("\nDomain: \n");
  63                 VD = DomainSimplify(Q->Domain, B, MAXRAYS);
  64                 Print_Domain(stdout, VD, param_name+nb_var);
  65                 matrix VM = domainVertices(PP, Q, params);
  66                 coeffs = bernsteinExpansion(VM, polynomial, vars, params);
  67                 printCoefficients(coeffs);
  68                 printMaxMinCoefficient(Q->Domain, coeffs, params);
  69                 Domain_Free(VD);
  70                 printf("\n\n===============================================\n");
  71         }
  72
  73         Domain_Free(A);
  74         Domain_Free(B);
  75         Param_Polyhedron_Free(PP);
  76         free(param_name);
  77
  78         return 0;
  79 } /* main */
  80
  81
  82 void printCoefficients(lst coeffs)
  83 {
  84         cout << "-----------------------------------------------" << endl;
  85         cout << "Coefficients: " << endl << endl;
  86
  87         for (lst::const_iterator i = coeffs.begin(); i != coeffs.end(); ++i)
  88             cout << "\t" << *i << endl;
  89 }
  90
  91
  92 /* Reads the polynomial matrix, converts it to long long precision and calls ginac functions */
  93 ex readPolynomial(const exvector& vars, const exvector& params)
  94 {
  95         char buffer[1024], *s;
  96         lst allvars;
  97         ex p;
  98
  99         for (int i = 0; i < vars.size(); ++i)
 100                 allvars.append(vars[i]);
 101         for (int i = 0; i < params.size(); ++i)
 102                 allvars.append(params[i]);
 103
 104         do
 105                 s = fgets(buffer, 1024, stdin);
 106         while (s && (s[0] == '#' || s[0] == '\n'));
 107
 108         if (!s)
 109                 return 0;
 110
 111         try {
 112                 p = ex(string(s), allvars);
 113         } catch (exception &p) {
 114                 cerr << p.what() << endl;
 115                 return 0;
 116         }
 117
 118         return p;
 119 }
 120
 121
 122 int printMaxMinCoefficient(Polyhedron *VD, lst coeffs, const exvector& Params)
 123 {
 124         cout << "\tMinimum coefficient(s): " << minimize(VD, coeffs, Params) << endl;
 125         cout << "\tMaximum coefficient(s): " << maximize(VD, coeffs, Params) << endl;
 126         return 0;
 127 }