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