applications/testCompressParms.c

   1 /**
   2  * $Id: testCompressParms.c,v 1.4 2006/09/18 03:09:03 meister Exp $
   3  *
   4  * Test routines for kernel/compress_parms.c functions
   5  * @author B. Meister, 3/2006
   6  *
   7  */
   8 /*
   9     This file is part of PolyLib.
  10
  11     PolyLib is free software: you can redistribute it and/or modify
  12     it under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     PolyLib is distributed in the hope that it will be useful,
  17     but WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19     GNU General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with PolyLib.  If not, see <http://www.gnu.org/licenses/>.
  23 */
  24
  25
  26 #include <polylib/polylib.h>
  27 #include <stdio.h>
  28 #include <stdlib.h>
  29 #include <string.h>
  30
  31 #ifdef dbg
  32 #undef dbg
  33 #endif
  34 #define dbg 1
  35
  36 #define TEST(a) if (isOk = a) { \
  37                   printf(#a" tested ok.\n"); \
  38                 } \
  39                 else { \
  40                   printf(#a" NOT OK\n"); \
  41                 }
  42
  43 #define maxRays 200
  44
  45 const char *origNames[] =
  46         {"n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z"};
  47
  48 int main(int argc, char ** argv) {
  49   int isOk = 0;
  50   Matrix * A, * B;
  51   if (argc>1) {
  52     printf("Warning: No arguments taken into account: testing"
  53            "remove_parm_eqs().\n");
  54   }
  55
  56   A = Matrix_Read();
  57   B = Matrix_Read();
  58   TEST( test_Constraints_Remove_parm_eqs(A, B) )
  59   TEST( test_Polyhedron_Remove_parm_eqs(A, B) )
  60   TEST( test_Constraints_fullDimensionize(A, B, 4) )
  61   Matrix_Free(A);
  62   Matrix_Free(B);
  63   return (1-isOk);
  64 }
  65
  66
  67 /** extracts the equalities involving the parameters only, try to introduce
  68     them back and compare the two polyhedra.
  69     Reads a polyhedron and a context.
  70  */
  71 int test_Constraints_Remove_parm_eqs(Matrix * A, Matrix * B) {
  72   int isOk = 1;
  73   Matrix * M, *C, *Cp, * Eqs, *M1, *C1;
  74   Polyhedron *Pm, *Pc, *Pcp, *Peqs, *Pint;
  75   unsigned int * elimParms;
  76   printf("----- test_Constraints_Remove_parm_eqs() -----\n");
  77   M1 = Matrix_Copy(A);
  78   C1 = Matrix_Copy(B);
  79
  80   M = Matrix_Copy(M1);
  81   C = Matrix_Copy(C1);
  82
  83    /* compute the combined polyhedron */
  84   Pm = Constraints2Polyhedron(M, maxRays);
  85   Pc = Constraints2Polyhedron(C, maxRays);
  86   Pcp = align_context(Pc, Pm->Dimension, maxRays);
  87   Polyhedron_Free(Pc);
  88   Pc = DomainIntersection(Pm, Pcp, maxRays);
  89   Polyhedron_Free(Pm);
  90   Polyhedron_Free(Pcp);
  91   Matrix_Free(M);
  92   Matrix_Free(C);
  93
  94   /* extract the parm-equalities, expressed in the combined space */
  95   Eqs = Constraints_Remove_parm_eqs(&M1, &C1, 1, &elimParms);
  96
  97   printf("Removed equalities: \n");
  98   show_matrix(Eqs);
  99   printf("Polyhedron without equalities involving only parameters: \n");
 100   show_matrix(M1);
 101   printf("Context without equalities: \n");
 102   show_matrix(C1);
 103
 104   /* compute the supposedly-same polyhedron, using the extracted equalities */
 105   Pm = Constraints2Polyhedron(M1, maxRays);
 106   Pcp = Constraints2Polyhedron(C1, maxRays);
 107   Peqs = align_context(Pcp, Pm->Dimension, maxRays);
 108   Polyhedron_Free(Pcp);
 109   Pcp = DomainIntersection(Pm, Peqs, maxRays);
 110   Polyhedron_Free(Peqs);
 111   Polyhedron_Free(Pm);
 112   Peqs = Constraints2Polyhedron(Eqs, maxRays);
 113   Matrix_Free(Eqs);
 114   Matrix_Free(M1);
 115   Matrix_Free(C1);
 116   Pint = DomainIntersection(Pcp, Peqs, maxRays);
 117   Polyhedron_Free(Pcp);
 118   Polyhedron_Free(Peqs);
 119
 120   /* test their equality */
 121   if (!PolyhedronIncludes(Pint, Pc)) {
 122     isOk = 0;
 123   }
 124   else {
 125     if (!PolyhedronIncludes(Pc, Pint)) {
 126       isOk = 0;
 127     }
 128   }
 129   Polyhedron_Free(Pc);
 130   Polyhedron_Free(Pint);
 131   return isOk;
 132 } /* test_Constraints_Remove_parm_eqs() */
 133
 134
 135 /** extracts the equalities holding on the parameters only, try to introduce
 136     them back and compare the two polyhedra.
 137     Reads a polyhedron and a context.
 138  */
 139 int test_Polyhedron_Remove_parm_eqs(Matrix * A, Matrix * B) {
 140   int isOk = 1;
 141   Matrix * M, *C;
 142   Polyhedron *Pm, *Pc, *Pcp, *Peqs, *Pint, *Pint1;
 143   unsigned int * elimParms;
 144   printf("----- test_Polyhedron_Remove_parm_eqs() -----\n");
 145
 146   M = Matrix_Copy(A);
 147   C = Matrix_Copy(B);
 148
 149    /* compute the combined polyhedron */
 150   Pm = Constraints2Polyhedron(M, maxRays);
 151   Pc = Constraints2Polyhedron(C, maxRays);
 152   Pcp = align_context(Pc, Pm->Dimension, maxRays);
 153   Polyhedron_Free(Pc);
 154   Pint1 = DomainIntersection(Pm, Pcp, maxRays);
 155   Polyhedron_Free(Pm);
 156   Polyhedron_Free(Pcp);
 157   Matrix_Free(M);
 158   Matrix_Free(C);
 159
 160   M = Matrix_Copy(A);
 161   C = Matrix_Copy(B);
 162   /* extract the parm-equalities, expressed in the combined space */
 163   Pm = Constraints2Polyhedron(M, maxRays);
 164   Pc = Constraints2Polyhedron(C, maxRays);
 165   Matrix_Free(M);
 166   Matrix_Free(C);
 167   Peqs = Polyhedron_Remove_parm_eqs(&Pm, &Pc, 1, &elimParms, 200);
 168
 169   /* compute the supposedly-same polyhedron, using the extracted equalities */
 170   Pcp = align_context(Pc, Pm->Dimension, maxRays);
 171   Polyhedron_Free(Pc);
 172   Pc = DomainIntersection(Pm, Pcp, maxRays);
 173   Polyhedron_Free(Pm);
 174   Polyhedron_Free(Pcp);
 175
 176   Pint = DomainIntersection(Pc, Peqs, maxRays);
 177   Polyhedron_Free(Pc);
 178   Polyhedron_Free(Peqs);
 179
 180   /* test their equality */
 181   if (!PolyhedronIncludes(Pint, Pint1)) {
 182     isOk = 0;
 183   }
 184   else {
 185     if (!PolyhedronIncludes(Pint1, Pint)) {
 186       isOk = 0;
 187     }
 188   }
 189   Polyhedron_Free(Pint1);
 190   Polyhedron_Free(Pint);
 191   return isOk;
 192 } /* test_Polyhedron_remove_parm_eqs() */
 193
 194
 195 /**
 196  * Eliminates certain parameters from a vector of values for parameters
 197  * @param origParms the initial vector of values of parameters
 198  * @param elimParms the list of parameters to be eliminated in the vector
 199  * @param newParms the vector of values without the eliminated ones.
 200  */
 201 void valuesWithoutElim(Matrix * origParms, unsigned int * elimParms,
 202                        Matrix ** newParms) {
 203   unsigned int i, j=0;
 204   if (*newParms==NULL) {
 205     *newParms = Matrix_Alloc(1, origParms->NbColumns-elimParms[0]);
 206   } /* else assume enough space is allocated */
 207   if (elimParms[0] ==0) {
 208     for (i=0; i< origParms->NbColumns; i++) {
 209       value_assign((*newParms)->p[0][i], origParms->p[0][i]);
 210     }
 211   }
 212   for (i=0; i< origParms->NbColumns; i++) {
 213     if (i!=elimParms[j+1]) {
 214       value_assign((*newParms)->p[0][i-j], origParms->p[0][i]);
 215     }
 216     else {
 217       j++;
 218     }
 219   }
 220 }/* valuesWithoutElim */
 221
 222
 223 /**
 224  * takes a list of parameter names, a list ofparameters to eliminate, and
 225  * returns the list of parameters without the eliminated ones.
 226  * @param parms the original parameter names
 227  * @param nbParms the number of original parmaeters
 228  * @param elimParms the array-list of indices of parameters to eliminate (first
 229  * element set to the number of its elements)
 230  * @param newParms the returned list of parm names. Allocated if set to NULL,
 231  * reused if not.
 232  * @return the number of names in the returned list.
 233  */
 234 unsigned int namesWithoutElim(const char **parms, unsigned nbParms,
 235                               unsigned int * elimParms,
 236                               const char ***newParms)
 237 {
 238   unsigned int i, j=0;
 239   unsigned int newSize = nbParms -elimParms[0];
 240   if (dbg) {
 241     printf("Size of the new parm vector: %d\n", newSize);
 242   }
 243   if (*newParms==NULL) {
 244     *newParms = malloc(newSize*sizeof(char *));
 245   }
 246   if (elimParms[0]==0) {
 247     for (i=0; i< nbParms; i++) {
 248       (*newParms)[i] = strdup(parms[i]);
 249     }
 250     return newSize;
 251   }
 252   for (i=0; i< nbParms; i++) {
 253     if (i!=elimParms[j+1]) {
 254       (*newParms)[i-j] = strdup(parms[i]);
 255     }
 256     else {
 257       j++;
 258     }
 259   }
 260   return newSize;
 261 }
 262
 263
 264 /**
 265  * Tests Constraints_fullDimensionize by comparing the Ehrhart polynomials
 266  * @param A the input set of constraints
 267  * @param B the corresponding context
 268  * @param the number of samples to generate for the test
 269  * @return 1 if the Ehrhart polynomial had the same value for the
 270  * full-dimensional and non-full-dimensional sets of constraints, for their
 271  * corresponding sample parameters values.
 272  */
 273 int test_Constraints_fullDimensionize(Matrix * A, Matrix * B,
 274                                       unsigned int nbSamples) {
 275   Matrix * Eqs= NULL, *ParmEqs=NULL, *VL=NULL;
 276   unsigned int * elimVars=NULL, * elimParms=NULL;
 277   Matrix * sample, * smallerSample=NULL;
 278   Matrix * transfSample=NULL;
 279   Matrix * parmVL=NULL;
 280   unsigned int i, j, r, nbOrigParms, nbParms;
 281   Value div, mod, *origVal=NULL, *fullVal=NULL;
 282   Matrix * VLInv;
 283   Polyhedron * P, *PC;
 284   Matrix * M, *C;
 285   Enumeration * origEP, * fullEP=NULL;
 286   const char **fullNames = NULL;
 287   int isOk = 1; /* holds the result */
 288
 289   /* compute the origial Ehrhart polynomial */
 290   M = Matrix_Copy(A);
 291   C = Matrix_Copy(B);
 292   P = Constraints2Polyhedron(M, maxRays);
 293   PC = Constraints2Polyhedron(C, maxRays);
 294   origEP = Polyhedron_Enumerate(P, PC, maxRays, origNames);
 295   Matrix_Free(M);
 296   Matrix_Free(C);
 297   Polyhedron_Free(P);
 298   Polyhedron_Free(PC);
 299
 300   /* compute the full-dimensional polyhedron corresponding to A and its Ehrhart
 301      polynomial */
 302   M = Matrix_Copy(A);
 303   C = Matrix_Copy(B);
 304   nbOrigParms = B->NbColumns-2;
 305   Constraints_fullDimensionize(&M, &C, &VL, &Eqs, &ParmEqs,
 306                                &elimVars, &elimParms, maxRays);
 307   if ((Eqs->NbRows==0) && (ParmEqs->NbRows==0)) {
 308     Matrix_Free(M);
 309     Matrix_Free(C);
 310     Matrix_Free(Eqs);
 311     Matrix_Free(ParmEqs);
 312     free(elimVars);
 313     free(elimParms);
 314     return 1;
 315   }
 316   nbParms = C->NbColumns-2;
 317   P = Constraints2Polyhedron(M, maxRays);
 318   PC = Constraints2Polyhedron(C, maxRays);
 319   namesWithoutElim(origNames, nbOrigParms, elimParms, &fullNames);
 320   fullEP = Polyhedron_Enumerate(P, PC, maxRays, fullNames);
 321   Matrix_Free(M);
 322   Matrix_Free(C);
 323   Polyhedron_Free(P);
 324   Polyhedron_Free(PC);
 325
 326   /* make a set of sample parameter values and compare the corresponding
 327      Ehrhart polnomials */
 328   sample = Matrix_Alloc(1,nbOrigParms);
 329   transfSample = Matrix_Alloc(1, nbParms);
 330   Lattice_extractSubLattice(VL, nbParms, &parmVL);
 331   VLInv = Matrix_Alloc(parmVL->NbRows, parmVL->NbRows+1);
 332   MatInverse(parmVL, VLInv);
 333   if (dbg) {
 334     show_matrix(parmVL);
 335     show_matrix(VLInv);
 336   }
 337   srand(nbSamples);
 338   value_init(mod);
 339   value_init(div);
 340   for (i = 0; i< nbSamples; i++) {
 341     /* create a random sample */
 342     for (j=0; j< nbOrigParms; j++) {
 343       value_set_si(sample->p[0][j], rand()%100);
 344     }
 345     /* compute the corresponding value for the full-dimensional
 346        constraints */
 347     valuesWithoutElim(sample, elimParms, &smallerSample);
 348     /* (N' i' 1)^T = VLinv.(N i 1)^T*/
 349     for (r = 0; r < nbParms; r++) {
 350       Inner_Product(&(VLInv->p[r][0]), smallerSample->p[0], nbParms,
 351                     &(transfSample->p[0][r]));
 352       /* add the constant part */
 353       value_addto(transfSample->p[0][r], transfSample->p[0][r],
 354                                          VLInv->p[r][VLInv->NbColumns-2]);
 355       value_pdivision(div, transfSample->p[0][r],
 356                          VLInv->p[r][VLInv->NbColumns-1]);
 357       value_subtract(mod, transfSample->p[0][r], div);
 358       /* if the parameters value does not belong to the validity lattice, the
 359          Ehrhart polynomial is zero. */
 360       if (!value_zero_p(mod)) {
 361         fullEP = Enumeration_zero(nbParms, maxRays);
 362         break;
 363       }
 364     }
 365     /* compare the two forms of the Ehrhart polynomial.*/
 366     if (origEP ==NULL) break; /* NULL has loose semantics for EPs */
 367     origVal = compute_poly(origEP, sample->p[0]);
 368     fullVal = compute_poly(fullEP, transfSample->p[0]);
 369     if (!value_eq(*origVal, *fullVal)) {
 370       isOk = 0;
 371       printf("EPs don't match. \n Original value = ");
 372       value_print(stdout, VALUE_FMT, *origVal);
 373       printf("\n Original sample = [");
 374       for (j=0; j<sample->NbColumns; j++) {
 375         value_print(stdout, VALUE_FMT, sample->p[0][j]);
 376         printf(" ");
 377       }
 378       printf("] \n EP = ");
 379       if(origEP!=NULL) {
 380         print_evalue(stdout, &(origEP->EP), origNames);
 381       }
 382       else {
 383         printf("NULL");
 384       }
 385       printf(" \n Full-dimensional value = ");
 386       value_print(stdout, P_VALUE_FMT, *fullVal);
 387       printf("\n full-dimensional sample = [");
 388       for (j=0; j<sample->NbColumns; j++) {
 389         value_print(stdout, VALUE_FMT, transfSample->p[0][j]);
 390         printf(" ");
 391       }
 392       printf("] \n EP = ");
 393       if(origEP!=NULL) {
 394         print_evalue(stdout, &(origEP->EP), fullNames);
 395       }
 396       else {
 397         printf("NULL");
 398       }
 399     }
 400     if (dbg) {
 401       printf("\nOriginal value = ");
 402       value_print(stdout, VALUE_FMT, *origVal);
 403       printf("\nFull-dimensional value = ");
 404       value_print(stdout, P_VALUE_FMT, *fullVal);
 405       printf("\n");
 406     }
 407     value_clear(*origVal);
 408     value_clear(*fullVal);
 409   }
 410   value_clear(mod);
 411   value_clear(div);
 412   Matrix_Free(sample);
 413   Matrix_Free(smallerSample);
 414   Matrix_Free(transfSample);
 415   Enumeration_Free(origEP);
 416   Enumeration_Free(fullEP);
 417   return isOk;
 418 } /* test_Constraints_fullDimensionize */