bernstein/src/main.cpp

   1 #include <assert.h>
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <unistd.h>
   5
   6 #include <ginac/ginac.h>
   7 #include <gmp.h>
   8 #include "polylib++.h"
   9
  10 #include <bernstein/bernstein.h>
  11 #include <bernstein/maximize.h>
  12 #include <bernstein/piecewise_lst.h>
  13
  14 #define MAXRAYS 1000
  15
  16 using namespace std;
  17 using namespace GiNaC;
  18 using namespace bernstein;
  19
  20 static ex readPolynomial(const exvector& vars, const exvector& params);
  21 static void printCoefficients(lst coeffs);
  22 static int printMaxMinCoefficient(Polyhedron *VD, lst coeffs, const exvector& Params);
  23
  24 static char *readLine(void)
  25 {
  26         static char buffer[1024];
  27         char *s;
  28
  29         do
  30                 s = fgets(buffer, 1024, stdin);
  31         while (s && (s[0] == '#' || s[0] == '\n'));
  32
  33         return s;
  34 }
  35
  36 static void readExpected(const exvector& vars, const exvector& params,
  37                          piecewise_lst **exp_all,
  38                          piecewise_lst **exp_min,
  39                          piecewise_lst **exp_max)
  40 {
  41         int n = 0;
  42         char *s;
  43         piecewise_lst *pl_all, *pl_min, *pl_max;
  44
  45         pl_all = new piecewise_lst(params);
  46         pl_min = new piecewise_lst(params);
  47         pl_max = new piecewise_lst(params);
  48
  49         s = readLine();
  50         assert(s);
  51
  52         sscanf(s, "%d", &n);
  53
  54         for (int i = 0 ; i < n ; ++i) {
  55                 ex all, min, max;
  56                 Polyhedron *D;
  57                 Matrix *M = Matrix_Read();
  58
  59                 D = Constraints2Polyhedron(M, MAXRAYS);
  60                 Matrix_Free(M);
  61
  62                 all = readPolynomial(vars, params);
  63                 min = readPolynomial(vars, params);
  64                 max = readPolynomial(vars, params);
  65                 assert(is_a<lst>(all));
  66                 assert(is_a<lst>(min));
  67                 assert(is_a<lst>(max));
  68                 pl_all->add_guarded_lst(Polyhedron_Copy(D), ex_to<lst>(all));
  69                 pl_min->add_guarded_lst(Polyhedron_Copy(D), ex_to<lst>(min));
  70                 pl_max->add_guarded_lst(D, ex_to<lst>(max));
  71         }
  72         *exp_all = pl_all;
  73         *exp_min = pl_min;
  74         *exp_max = pl_max;
  75 }
  76
  77 /* main function */
  78 int main(int argc, char *argv[])
  79 {
  80         Matrix *a, *b;
  81         Polyhedron *A, *B;              // initial matrices
  82         const char **param_name;        // name of the parameters
  83         exvector params, vars;
  84         ex polynomial;
  85         piecewise_lst *exp_all, *exp_min, *exp_max;
  86         piecewise_lst *got_all, *got_min, *got_max;
  87
  88         Param_Polyhedron *PP;
  89         Param_Domain   *Q;
  90
  91         unsigned int nb_param, nb_var;
  92
  93         int c;
  94         int verify = 0;
  95
  96         while ((c = getopt(argc, argv, "T")) != -1) {
  97                 switch(c) {
  98                 case 'T':
  99                         verify = 1;
 100                         break;
 101                 }
 102         }
 103
 104         printf("\n===============================================\n");
 105
 106         a = Matrix_Read();
 107         A = Constraints2Polyhedron(a,200);
 108
 109         b = Matrix_Read();
 110         B = Constraints2Polyhedron(b,200);
 111
 112         /* Set the number of parameters */
 113         nb_param        = B->Dimension;
 114         nb_var          = A->Dimension - B->Dimension;
 115
 116         /* Read the name of the parameters */
 117         param_name = Read_ParamNames(stdin, nb_var+nb_param);
 118         vars = constructParameterVector(param_name, nb_var);
 119         params = constructParameterVector(param_name+nb_var, nb_param);
 120
 121         polynomial = readPolynomial(vars, params);
 122
 123         if (verify) {
 124                 readExpected(vars, params, &exp_all, &exp_min, &exp_max);
 125                 got_all = new piecewise_lst(params);
 126                 got_min = new piecewise_lst(params);
 127                 got_max = new piecewise_lst(params);
 128         }
 129
 130         Matrix_Free(a);
 131         Matrix_Free(b);
 132
 133         /* Find the parametrized domains */
 134         PP = Polyhedron2Param_Domain(A,B,MAXRAYS);
 135         for(Q=PP->D;Q;Q=Q->next) {
 136                 Polyhedron *VD;
 137                 lst coeffs;
 138
 139                 printf("\nDomain: \n");
 140                 VD = DomainSimplify(Q->Domain, B, MAXRAYS);
 141                 Print_Domain(stdout, VD, param_name+nb_var);
 142                 matrix VM = domainVertices(PP, Q, params);
 143                 coeffs = bernsteinExpansion(VM, polynomial, vars, params);
 144                 printCoefficients(coeffs);
 145                 printMaxMinCoefficient(Q->Domain, coeffs, params);
 146                 Domain_Free(VD);
 147                 printf("\n\n===============================================\n");
 148
 149                 if (!verify)
 150                         continue;
 151
 152                 got_all->add_guarded_lst(Polyhedron_Copy(Q->Domain), coeffs);
 153                 got_min->add_guarded_lst(Polyhedron_Copy(Q->Domain),
 154                                          minimize(Q->Domain, coeffs, params));
 155                 got_max->add_guarded_lst(Polyhedron_Copy(Q->Domain),
 156                                          maximize(Q->Domain, coeffs, params));
 157         }
 158
 159         Domain_Free(A);
 160         Domain_Free(B);
 161         Param_Polyhedron_Free(PP);
 162         free(param_name);
 163
 164         if (verify) {
 165                 if (!got_all->is_equal(*exp_all)) {
 166                         cerr << "expected:" << endl;
 167                         cerr << *exp_all << endl;
 168                         cerr << "got:" << endl;
 169                         cerr << *got_all << endl;
 170                         return 1;
 171                 }
 172                 if (!got_min->is_equal(*exp_min)) {
 173                         cerr << "expected:" << endl;
 174                         cerr << *exp_min << endl;
 175                         cerr << "got:" << endl;
 176                         cerr << *got_min << endl;
 177                         return 1;
 178                 }
 179                 if (!got_max->is_equal(*exp_max)) {
 180                         cerr << "expected:" << endl;
 181                         cerr << *exp_max << endl;
 182                         cerr << "got:" << endl;
 183                         cerr << *got_max << endl;
 184                         return 1;
 185                 }
 186
 187                 delete exp_all;
 188                 delete exp_max;
 189                 delete exp_min;
 190                 delete got_all;
 191                 delete got_max;
 192                 delete got_min;
 193         }
 194
 195         return 0;
 196 } /* main */
 197
 198
 199 void printCoefficients(lst coeffs)
 200 {
 201         cout << "-----------------------------------------------" << endl;
 202         cout << "Coefficients: " << endl << endl;
 203
 204         for (lst::const_iterator i = coeffs.begin(); i != coeffs.end(); ++i)
 205             cout << "\t" << *i << endl;
 206 }
 207
 208
 209 /* Reads the polynomial matrix, converts it to long long precision and calls ginac functions */
 210 ex readPolynomial(const exvector& vars, const exvector& params)
 211 {
 212         char *s;
 213         lst allvars;
 214         ex p;
 215
 216         for (int i = 0; i < vars.size(); ++i)
 217                 allvars.append(vars[i]);
 218         for (int i = 0; i < params.size(); ++i)
 219                 allvars.append(params[i]);
 220
 221         s = readLine();
 222         if (!s)
 223                 return 0;
 224
 225         try {
 226                 p = ex(string(s), allvars);
 227         } catch (exception &p) {
 228                 cerr << p.what() << endl;
 229                 return 0;
 230         }
 231
 232         return p;
 233 }
 234
 235
 236 int printMaxMinCoefficient(Polyhedron *VD, lst coeffs, const exvector& Params)
 237 {
 238         cout << "\tMinimum coefficient(s): " << minimize(VD, coeffs, Params) << endl;
 239         cout << "\tMaximum coefficient(s): " << maximize(VD, coeffs, Params) << endl;
 240         return 0;
 241 }