polylib: testCompressParms.c Source File

00001 /** 
00002  * $Id: testCompressParms.c,v 1.4 2006/09/18 03:09:03 meister Exp $
00003  * 
00004  * Test routines for kernel/compress_parms.c functions
00005  * @author B. Meister, 3/2006
00006  * 
00007  */
00008 
00009 #include <polylib/polylib.h>
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012 #include <string.h>
00013 
00014 #ifdef dbg
00015 #undef dbg
00016 #endif
00017 #define dbg 1
00018 
00019 #define TEST(a) if (isOk = a) { \
00020                   printf(#a" tested ok.\n"); \
00021                 } \
00022                 else { \
00023                   printf(#a" NOT OK\n"); \
00024                 } 
00025 
00026 #define maxRays 200
00027 
00028 const char *origNames[] =
00029         {"n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z"};
00030 
00031 int main(int argc, char ** argv) {
00032   int isOk = 0;
00033   Matrix * A, * B;
00034   if (argc>1) {
00035     printf("Warning: No arguments taken into account: testing"
00036            "remove_parm_eqs().\n");
00037   }
00038 
00039   A = Matrix_Read();
00040   B = Matrix_Read();
00041   TEST( test_Constraints_Remove_parm_eqs(A, B) )
00042   TEST( test_Polyhedron_Remove_parm_eqs(A, B) )
00043   TEST( test_Constraints_fullDimensionize(A, B, 4) )
00044   Matrix_Free(A);
00045   Matrix_Free(B);
00046   return (1-isOk);
00047 }
00048 
00049 
00050 /** extracts the equalities involving the parameters only, try to introduce
00051     them back and compare the two polyhedra.
00052     Reads a polyhedron and a context.
00053  */
00054 int test_Constraints_Remove_parm_eqs(Matrix * A, Matrix * B) {
00055   int isOk = 1;
00056   Matrix * M, *C, *Cp, * Eqs, *M1, *C1;
00057   Polyhedron *Pm, *Pc, *Pcp, *Peqs, *Pint;
00058   unsigned int * elimParms;
00059   printf("----- test_Constraints_Remove_parm_eqs() -----\n");
00060   M1 = Matrix_Copy(A);
00061   C1 = Matrix_Copy(B);
00062 
00063   M = Matrix_Copy(M1);
00064   C = Matrix_Copy(C1);
00065 
00066    /* compute the combined polyhedron */
00067   Pm = Constraints2Polyhedron(M, maxRays);
00068   Pc = Constraints2Polyhedron(C, maxRays);
00069   Pcp = align_context(Pc, Pm->Dimension, maxRays);
00070   Polyhedron_Free(Pc);
00071   Pc = DomainIntersection(Pm, Pcp, maxRays);
00072   Polyhedron_Free(Pm);
00073   Polyhedron_Free(Pcp);
00074   Matrix_Free(M);
00075   Matrix_Free(C);
00076 
00077   /* extract the parm-equalities, expressed in the combined space */
00078   Eqs = Constraints_Remove_parm_eqs(&M1, &C1, 1, &elimParms);
00079 
00080   printf("Removed equalities: \n");
00081   show_matrix(Eqs); 
00082   printf("Polyhedron without equalities involving only parameters: \n");
00083   show_matrix(M1);  
00084   printf("Context without equalities: \n");
00085   show_matrix(C1);  
00086   
00087   /* compute the supposedly-same polyhedron, using the extracted equalities */
00088   Pm = Constraints2Polyhedron(M1, maxRays);
00089   Pcp = Constraints2Polyhedron(C1, maxRays);
00090   Peqs = align_context(Pcp, Pm->Dimension, maxRays);
00091   Polyhedron_Free(Pcp);
00092   Pcp = DomainIntersection(Pm, Peqs, maxRays);
00093   Polyhedron_Free(Peqs);
00094   Polyhedron_Free(Pm);
00095   Peqs = Constraints2Polyhedron(Eqs, maxRays);
00096   Matrix_Free(Eqs);
00097   Matrix_Free(M1);
00098   Matrix_Free(C1);
00099   Pint = DomainIntersection(Pcp, Peqs, maxRays);
00100   Polyhedron_Free(Pcp);
00101   Polyhedron_Free(Peqs);
00102 
00103   /* test their equality */
00104   if (!PolyhedronIncludes(Pint, Pc)) {
00105     isOk = 0;
00106   }
00107   else {
00108     if (!PolyhedronIncludes(Pc, Pint)) {
00109       isOk = 0;
00110     }
00111   }
00112   Polyhedron_Free(Pc);
00113   Polyhedron_Free(Pint);
00114   return isOk;
00115 } /* test_Constraints_Remove_parm_eqs() */
00116 
00117 
00118 /** extracts the equalities holding on the parameters only, try to introduce
00119     them back and compare the two polyhedra.
00120     Reads a polyhedron and a context.
00121  */
00122 int test_Polyhedron_Remove_parm_eqs(Matrix * A, Matrix * B) {
00123   int isOk = 1;
00124   Matrix * M, *C;
00125   Polyhedron *Pm, *Pc, *Pcp, *Peqs, *Pint, *Pint1;
00126   unsigned int * elimParms;
00127   printf("----- test_Polyhedron_Remove_parm_eqs() -----\n");
00128 
00129   M = Matrix_Copy(A);
00130   C = Matrix_Copy(B);
00131 
00132    /* compute the combined polyhedron */
00133   Pm = Constraints2Polyhedron(M, maxRays);
00134   Pc = Constraints2Polyhedron(C, maxRays);
00135   Pcp = align_context(Pc, Pm->Dimension, maxRays);
00136   Polyhedron_Free(Pc);
00137   Pint1 = DomainIntersection(Pm, Pcp, maxRays);
00138   Polyhedron_Free(Pm);
00139   Polyhedron_Free(Pcp);
00140   Matrix_Free(M);
00141   Matrix_Free(C);
00142 
00143   M = Matrix_Copy(A);
00144   C = Matrix_Copy(B);
00145   /* extract the parm-equalities, expressed in the combined space */
00146   Pm = Constraints2Polyhedron(M, maxRays);
00147   Pc = Constraints2Polyhedron(C, maxRays);
00148   Matrix_Free(M);
00149   Matrix_Free(C);
00150   Peqs = Polyhedron_Remove_parm_eqs(&Pm, &Pc, 1, &elimParms, 200);
00151   
00152   /* compute the supposedly-same polyhedron, using the extracted equalities */
00153   Pcp = align_context(Pc, Pm->Dimension, maxRays);
00154   Polyhedron_Free(Pc);
00155   Pc = DomainIntersection(Pm, Pcp, maxRays);
00156   Polyhedron_Free(Pm);
00157   Polyhedron_Free(Pcp);
00158  
00159   Pint = DomainIntersection(Pc, Peqs, maxRays);
00160   Polyhedron_Free(Pc);
00161   Polyhedron_Free(Peqs);
00162 
00163   /* test their equality */
00164   if (!PolyhedronIncludes(Pint, Pint1)) {
00165     isOk = 0;
00166   }
00167   else {
00168     if (!PolyhedronIncludes(Pint1, Pint)) {
00169       isOk = 0;
00170     }
00171   }
00172   Polyhedron_Free(Pint1);
00173   Polyhedron_Free(Pint);
00174   return isOk;
00175 } /* test_Polyhedron_remove_parm_eqs() */
00176 
00177 
00178 /** 
00179  * Eliminates certain parameters from a vector of values for parameters
00180  * @param origParms the initial vector of values of parameters
00181  * @param elimParms the list of parameters to be eliminated in the vector
00182  * @param newParms the vector of values without the eliminated ones.
00183  */
00184 void valuesWithoutElim(Matrix * origParms, unsigned int * elimParms, 
00185                        Matrix ** newParms) {
00186   unsigned int i, j=0;
00187   if (*newParms==NULL) {
00188     *newParms = Matrix_Alloc(1, origParms->NbColumns-elimParms[0]);
00189   } /* else assume enough space is allocated */
00190   if (elimParms[0] ==0) {
00191     for (i=0; i< origParms->NbColumns; i++) {
00192       value_assign((*newParms)->p[0][i], origParms->p[0][i]);
00193     }
00194   }
00195   for (i=0; i< origParms->NbColumns; i++) {
00196     if (i!=elimParms[j+1]) {
00197       value_assign((*newParms)->p[0][i-j], origParms->p[0][i]);
00198     }
00199     else {
00200       j++;
00201     }
00202   }
00203 }/* valuesWithoutElim */
00204 
00205 
00206 /**
00207  * takes a list of parameter names, a list ofparameters to eliminate, and
00208  * returns the list of parameters without the eliminated ones.
00209  * @param parms the original parameter names
00210  * @param nbParms the number of original parmaeters
00211  * @param elimParms the array-list of indices of parameters to eliminate (first
00212  * element set to the number of its elements)
00213  * @param newParms the returned list of parm names. Allocated if set to NULL,
00214  * reused if not.
00215  * @return the number of names in the returned list.
00216  */
00217 unsigned int namesWithoutElim(const char **parms, unsigned nbParms,
00218                               unsigned int * elimParms,
00219                               const char ***newParms)
00220 {
00221   unsigned int i, j=0;
00222   unsigned int newSize = nbParms -elimParms[0];
00223   if (dbg) {
00224     printf("Size of the new parm vector: %d\n", newSize);
00225   }
00226   if (*newParms==NULL) {
00227     *newParms = malloc(newSize*sizeof(char *));
00228   }
00229   if (elimParms[0]==0) {
00230     for (i=0; i< nbParms; i++) {
00231       (*newParms)[i] = strdup(parms[i]);
00232     }
00233     return newSize;
00234   }
00235   for (i=0; i< nbParms; i++) {
00236     if (i!=elimParms[j+1]) {
00237       (*newParms)[i-j] = strdup(parms[i]);
00238     }
00239     else {
00240       j++;
00241     }
00242   }
00243   return newSize;
00244 }
00245 
00246 
00247 /**
00248  * Tests Constraints_fullDimensionize by comparing the Ehrhart polynomials 
00249  * @param A the input set of constraints
00250  * @param B the corresponding context
00251  * @param the number of samples to generate for the test
00252  * @return 1 if the Ehrhart polynomial had the same value for the
00253  * full-dimensional and non-full-dimensional sets of constraints, for their
00254  * corresponding sample parameters values.
00255  */
00256 int test_Constraints_fullDimensionize(Matrix * A, Matrix * B, 
00257                                       unsigned int nbSamples) {
00258   Matrix * Eqs= NULL, *ParmEqs=NULL, *VL=NULL;
00259   unsigned int * elimVars=NULL, * elimParms=NULL;
00260   Matrix * sample, * smallerSample=NULL;
00261   Matrix * transfSample=NULL;
00262   Matrix * parmVL=NULL;
00263   unsigned int i, j, r, nbOrigParms, nbParms;
00264   Value div, mod, *origVal=NULL, *fullVal=NULL;
00265   Matrix * VLInv;
00266   Polyhedron * P, *PC;
00267   Matrix * M, *C;
00268   Enumeration * origEP, * fullEP=NULL;
00269   const char **fullNames = NULL;
00270   int isOk = 1; /* holds the result */
00271 
00272   /* compute the origial Ehrhart polynomial */
00273   M = Matrix_Copy(A);
00274   C = Matrix_Copy(B);
00275   P = Constraints2Polyhedron(M, maxRays);
00276   PC = Constraints2Polyhedron(C, maxRays);
00277   origEP = Polyhedron_Enumerate(P, PC, maxRays, origNames);
00278   Matrix_Free(M);
00279   Matrix_Free(C);
00280   Polyhedron_Free(P);
00281   Polyhedron_Free(PC);
00282 
00283   /* compute the full-dimensional polyhedron corresponding to A and its Ehrhart
00284      polynomial */
00285   M = Matrix_Copy(A);
00286   C = Matrix_Copy(B);
00287   nbOrigParms = B->NbColumns-2;
00288   Constraints_fullDimensionize(&M, &C, &VL, &Eqs, &ParmEqs, 
00289                                &elimVars, &elimParms, maxRays);
00290   if ((Eqs->NbRows==0) && (ParmEqs->NbRows==0)) {
00291     Matrix_Free(M);
00292     Matrix_Free(C);
00293     Matrix_Free(Eqs);
00294     Matrix_Free(ParmEqs);
00295     free(elimVars);
00296     free(elimParms);
00297     return 1;
00298   }
00299   nbParms = C->NbColumns-2;
00300   P = Constraints2Polyhedron(M, maxRays);
00301   PC = Constraints2Polyhedron(C, maxRays);
00302   namesWithoutElim(origNames, nbOrigParms, elimParms, &fullNames);
00303   fullEP = Polyhedron_Enumerate(P, PC, maxRays, fullNames);
00304   Matrix_Free(M);
00305   Matrix_Free(C);
00306   Polyhedron_Free(P);
00307   Polyhedron_Free(PC);
00308   
00309   /* make a set of sample parameter values and compare the corresponding
00310      Ehrhart polnomials */
00311   sample = Matrix_Alloc(1,nbOrigParms);
00312   transfSample = Matrix_Alloc(1, nbParms);
00313   Lattice_extractSubLattice(VL, nbParms, &parmVL);
00314   VLInv = Matrix_Alloc(parmVL->NbRows, parmVL->NbRows+1);
00315   MatInverse(parmVL, VLInv);
00316   if (dbg) {
00317     show_matrix(parmVL);
00318     show_matrix(VLInv);
00319   }
00320   srand(nbSamples);
00321   value_init(mod);
00322   value_init(div);
00323   for (i = 0; i< nbSamples; i++) {
00324     /* create a random sample */
00325     for (j=0; j< nbOrigParms; j++) {
00326       value_set_si(sample->p[0][j], rand()%100);
00327     }
00328     /* compute the corresponding value for the full-dimensional
00329        constraints */
00330     valuesWithoutElim(sample, elimParms, &smallerSample); 
00331     /* (N' i' 1)^T = VLinv.(N i 1)^T*/
00332     for (r = 0; r < nbParms; r++) {
00333       Inner_Product(&(VLInv->p[r][0]), smallerSample->p[0], nbParms,
00334                     &(transfSample->p[0][r]));
00335       /* add the constant part */
00336       value_addto(transfSample->p[0][r], transfSample->p[0][r], 
00337                                          VLInv->p[r][VLInv->NbColumns-2]);
00338       value_pdivision(div, transfSample->p[0][r], 
00339                          VLInv->p[r][VLInv->NbColumns-1]);
00340       value_subtract(mod, transfSample->p[0][r], div);
00341       /* if the parameters value does not belong to the validity lattice, the
00342          Ehrhart polynomial is zero. */
00343       if (!value_zero_p(mod)) {
00344         fullEP = Enumeration_zero(nbParms, maxRays);
00345         break;
00346       }
00347     }
00348     /* compare the two forms of the Ehrhart polynomial.*/
00349     if (origEP ==NULL) break; /* NULL has loose semantics for EPs */
00350     origVal = compute_poly(origEP, sample->p[0]);
00351     fullVal = compute_poly(fullEP, transfSample->p[0]);
00352     if (!value_eq(*origVal, *fullVal)) {
00353       isOk = 0;
00354       printf("EPs don't match. \n Original value = ");
00355       value_print(stdout, VALUE_FMT, *origVal);
00356       printf("\n Original sample = [");
00357       for (j=0; j<sample->NbColumns; j++) {
00358         value_print(stdout, VALUE_FMT, sample->p[0][j]);
00359         printf(" ");
00360       }
00361       printf("] \n EP = ");
00362       if(origEP!=NULL) {
00363         print_evalue(stdout, &(origEP->EP), origNames);
00364       }
00365       else {
00366         printf("NULL");
00367       }
00368       printf(" \n Full-dimensional value = ");
00369       value_print(stdout, P_VALUE_FMT, *fullVal);
00370       printf("\n full-dimensional sample = [");
00371       for (j=0; j<sample->NbColumns; j++) {
00372         value_print(stdout, VALUE_FMT, transfSample->p[0][j]);
00373         printf(" ");
00374       }
00375       printf("] \n EP = ");
00376       if(origEP!=NULL) {
00377         print_evalue(stdout, &(origEP->EP), fullNames);
00378       }
00379       else {
00380         printf("NULL");
00381       }
00382     }
00383     if (dbg) {
00384       printf("\nOriginal value = ");
00385       value_print(stdout, VALUE_FMT, *origVal);
00386       printf("\nFull-dimensional value = ");
00387       value_print(stdout, P_VALUE_FMT, *fullVal);
00388       printf("\n");
00389     }
00390     value_clear(*origVal);
00391     value_clear(*fullVal);
00392   }
00393   value_clear(mod);
00394   value_clear(div);
00395   Matrix_Free(sample);
00396   Matrix_Free(smallerSample);
00397   Matrix_Free(transfSample);
00398   Enumeration_Free(origEP);
00399   Enumeration_Free(fullEP);
00400   return isOk;
00401 } /* test_Constraints_fullDimensionize */