Robotsystem of the Robot Group Leipzig: matrix.tests.hpp Source File

00001 /***************************************************************************
00002                           matrix.tests.cpp  -  description
00003                              -------------------
00004     email                : georg.martius@web.de
00005 ***************************************************************************/
00006 // 
00007 // $Log: matrix.tests.hpp,v $
00008 // Revision 1.5  2008/06/18 13:46:20  martius
00009 // beside and toBeside added
00010 // addColumns and addRows is now based on toAbove and toBeside
00011 // bug fix in removeColumns
00012 //
00013 // Revision 1.4  2007/07/17 07:28:06  martius
00014 // store and restore test
00015 //
00016 // Revision 1.3  2006/07/20 17:14:36  martius
00017 // removed std namespace from matrix.h
00018 // storable interface
00019 // abstract model and invertablemodel as superclasses for networks
00020 //
00021 // Revision 1.2  2006/07/14 12:24:01  martius
00022 // selforg becomes HEAD
00023 //
00024 // Revision 1.1.2.1  2006/07/10 12:01:02  martius
00025 // Matrixlib moved to selforg
00026 //
00027 // Revision 1.15  2005/10/21 11:58:25  martius
00028 // map2 (similar to map but for 2 matrices)
00029 // changed naming of functions to be more consistent.
00030 //  Functions with "to" in front indicate the change of this. (Still not consistent with add, mult ...)
00031 //
00032 // Revision 1.14  2005/10/06 17:10:06  martius
00033 // convertToList
00034 // above and toAbove
00035 //
00036 // Revision 1.13  2005/08/06 20:47:36  martius
00037 // Commented
00038 //
00039 // Revision 1.12  2005/06/17 15:18:09  martius
00040 // 2x2 Inversion tested
00041 //
00042 // Revision 1.11  2005/06/10 08:21:50  martius
00043 // mult???wise are copy operations now!
00044 // toMult???wise are inplace instead
00045 //
00046 // Revision 1.10  2005/06/09 11:52:03  martius
00047 // multMT (M * M^T) and multTM (M^T * M)
00048 //
00049 // Revision 1.9  2005/06/02 22:49:33  martius
00050 // speed tests extended
00051 //
00052 // Revision 1.8  2005/06/02 08:49:26  martius
00053 // mult_row/column_wise
00054 //
00055 // Revision 1.7  2005/05/30 22:40:56  martius
00056 // map becomes toMap and the new map returns a new matrix
00057 // exp becomes toExp
00058 //
00059 // Revision 1.6  2005/05/30 21:46:54  martius
00060 // *** empty log message ***
00061 //
00062 // Revision 1.5  2005/05/30 16:43:02  martius
00063 // map function included (component-wise function application)
00064 //
00065 // Revision 1.4  2005/05/30 10:15:47  martius
00066 // proper log entry in header
00067 //
00068 /***************************************************************************/
00069 
00070 // //////////////////////////////////////////////////////////////////////////////
00071 // ////////////////// UNIT TESTS for matrix lib /////////////////////////////////////////////////
00072 
00073 #ifdef UNITTEST
00074 #include "unit_test.hpp"
00075 using namespace matrix;
00076 using namespace std;
00077 
00078 const D EPS=1e-9;
00079 bool comparetoidentity(const Matrix& m)  {
00080   int worstdiagonal = 0;
00081   D maxunitydeviation = 0.0;
00082   D currentunitydeviation;
00083   for ( int i = 0; i < m.getM(); i++ )  {
00084     currentunitydeviation = m.val(i,i) - 1.;
00085     if ( currentunitydeviation < 0.0) currentunitydeviation *= -1.;
00086     if ( currentunitydeviation > maxunitydeviation )  {
00087       maxunitydeviation = currentunitydeviation;
00088       worstdiagonal = i;
00089     }
00090   }
00091   int worstoffdiagonalrow = 0;
00092   int worstoffdiagonalcolumn = 0;
00093   D maxzerodeviation = 0.0;
00094   D currentzerodeviation ;
00095   for ( int i = 0; i < m.getM(); i++ )  {
00096     for ( int j = 0; j < m.getN(); j++ )  {
00097       if ( i == j ) continue;  // we look only at non-diagonal terms
00098       currentzerodeviation = m.val(i,j);
00099       if ( currentzerodeviation < 0.0) currentzerodeviation *= -1.0;
00100       if ( currentzerodeviation > maxzerodeviation )  {
00101         maxzerodeviation = currentzerodeviation;
00102         worstoffdiagonalrow = i;
00103         worstoffdiagonalcolumn = j;
00104       }
00105 
00106     }
00107   }
00108 //   cout << "\tWorst diagonal value deviation from unity: " 
00109 //        << maxunitydeviation << " at row/column " << worstdiagonal << endl;
00110 //   cout << "\tWorst off-diagonal value deviation from zero: " 
00111 //        << maxzerodeviation << " at row = " << worstoffdiagonalrow 
00112 //        << ", column = " << worstoffdiagonalcolumn << endl;
00113   return (maxunitydeviation < EPS && maxzerodeviation < EPS); 
00114 }
00115 
00116 
00117 UNIT_TEST_DEFINES
00118 
00119 DEFINE_TEST( check_creation ) {  
00120   cout << "\n -[ Creation and Access ]-\n";
00121   Matrix M1(3,3);  
00122   D testdata[9]={1,0,0, 0,1,0, 0,0,1};
00123   Matrix M2(3,3);
00124   M2.set(testdata);
00125   M1.toId();
00126   unit_assert( "id=identity_set", M1 == M2 );
00127   D testdata2[12]={1,2,3, 4,5,6, 7,8,9, 10,11,12};
00128   Matrix M3(4,3, testdata2);
00129   unit_assert( "dimension of matrix", M3.getM() == 4 &&  M3.getN() == 3 );
00130   unit_assert( "field check", 
00131                M3.val(0,2) == 3 &&  M3.val(2,1) == 8 &&  M3.val(3,0) == 10 );
00132   Matrix M4(M3);  
00133   unit_assert( "copy constructor", M3 == M4 );
00134   Matrix M5(1,3,testdata2+3);  
00135   unit_assert( "row", M3.row(1) == M5 );
00136 
00137   unit_pass();
00138 }
00139 
00140 DEFINE_TEST( check_vector_operation ) {  
00141   cout << "\n -[ Vector Operations ]-\n";
00142   D testdata[3]={-1,3,2};
00143   const Matrix V1(1,3, testdata);
00144   const Matrix V2(3,1, testdata);  
00145   Matrix V3(V1);  
00146 
00147   V3.toTranspose();
00148   unit_assert( "transpose", V3 == V2 );
00149   V3.toTranspose();
00150   unit_assert( "double transpose", V3 == V1 );
00151   D testdata2[3]={-2,6,4};
00152   Matrix V4(1,3,testdata2);    
00153   V3.add(V1,V1);
00154   unit_assert( "add", V3 == V4 );
00155   D testdata3[3]={1,-3,-2};
00156   Matrix V5(1,3,testdata3);
00157   V4.sub(V1,V3);
00158   unit_assert( "sub", V4 == V5 );
00159 
00160   D testdata4[3]={3,-9,-6};
00161   V5.set(testdata4);
00162   V4.copy(V1);
00163   V4.toMult(-3.0);
00164   unit_assert( "mult with scalar I", V4 == V5 );
00165   V4.mult(V1,-3.0);
00166   unit_assert( "mult with scalar II", V4 == V5 );
00167 
00168   double f = 14;
00169   Matrix V6(1,1,&f);
00170   V3.copy(V1);  
00171   V3.toTranspose();
00172   V4.mult(V1,V3);
00173   unit_assert( "scalarproduct", V4 == V6 );
00174   V4.copy(V1);
00175   unit_assert( "scalarproduct with exp(2)", V4.multMT() == V6 );
00176   V4.mult(V3,V1);
00177   unit_assert( "vector^T*vector=matrix", V4.getM() == 3);
00178   unit_assert( "vector^T*vector=exp(2)", V3.multMT().getM() == 3);
00179     
00180   unit_pass();
00181 }
00182 
00183 DEFINE_TEST( check_matrix_operation ) {  
00184   cout << "\n -[ Matrix Operations ]-\n";
00185   D testdata[6]={1,2,3, 4,5,6 };
00186   const Matrix M1(2,3, testdata);
00187   D testdata2[6]={1,4, 2,5, 3,6 };
00188   const Matrix M2(3,2, testdata2);  
00189   Matrix M3(M1);  
00190 
00191   M3.toTranspose();
00192   unit_assert( "transpose", M3 == M2 );
00193   D testdata3[6]={2,4,6, 8,10,12};
00194   Matrix M4(2,3,testdata3);    
00195   M3.add(M1,M1);
00196   unit_assert( "add", M3 == M4 );
00197 
00198   D testdata4[6]={-0.5, -1, -1.5,  -2, -2.5, -3};
00199   M3.set(testdata4);
00200   M4.copy(M1);
00201   M4.toMult(-0.5);
00202   unit_assert( "mult with scalar", M4 == M3 );
00203 
00204   D testdata5[6] = {2, -1,  0, 3,  2, -2};
00205   D testdata6[4] = {8, -1,  20, -1};
00206   Matrix M5 (3,2, testdata5);
00207   Matrix M6(2,2,testdata6);
00208   M4.mult(M1,M5);
00209   unit_assert( "mult(matrix, matrix)", M4 == M6 );
00210 
00211   M3.copy(M1);  
00212   M4.copy(M1);  
00213   M4.toExp(1);
00214   unit_assert( "exp(1)", M3 == M4 );
00215   M3.toTranspose();
00216   M4.toExp(T);
00217   unit_assert( "exp(T)=transpose", M3 == M4 );
00218   M3.toId();
00219   M4.toExp(0);
00220   unit_assert( "exp(0)=id", M3 == M4 );
00221 
00222   D testdata7[16] = {1,2,3,4, -4,2,1,3, 0.3,-0.9, 4, -3, 1,0.5,0.3,5.0};
00223   Matrix M7(4,4,testdata7);
00224   Matrix M8(M7);  
00225   M7.toExp(-1);  
00226   M4.mult(M7,M8);
00227   unit_assert( "exp(-1)*exp(1)=id",   comparetoidentity(M4) );
00228 
00229   D testdata9[6] = {sin(1.0),sin(2.0),sin(3.0), sin(4.0),sin(5.0),sin(6.0) };
00230   Matrix M9(2,3,testdata9);  
00231   M4.copy(M1);
00232   M4.toMap(sin);
00233   unit_assert( "map(sin)",   M4 == M9 );
00234 
00235   D testdata10[6] = {2,4,6, -0.4,-0.5,-0.6 };
00236   D testdata11[2] = {2,-0.1};
00237   Matrix M10(2,3,testdata10);  
00238   Matrix M11(2,1,testdata11);  
00239   M4.copy(M1);
00240   M4.toMultrowwise(M11);
00241   unit_assert( "multrowwise()",   M4 == M10 );
00242   D testdata12[6] = {2,1,0, 8, 2.5, 0 };
00243   D testdata13[3] = {2, 0.5, 0};
00244   M10.set(2,3,testdata12);  
00245   Matrix M12(3,1,testdata13);  
00246   M4.copy(M1);
00247   M4.toMultcolwise(M12);
00248   unit_assert( "multrowwise()",   M4 == M10 );
00249 
00250   M3.copy(M1);  
00251   M4.copy(M1);
00252   M4.toTranspose();
00253   M5 = M3.multMT();
00254   M6.mult(M1,M4);
00255   unit_assert( "multMT() ",   M5 == M6 );
00256   M5 = M3.multTM();
00257   M6.mult(M4,M1);
00258   unit_assert( "multTM() ",   M5 == M6 );
00259 
00260   D testdata20[12]={1,2,3, 4,5,6, 1,2,3, 4,5,6 };
00261   const Matrix M20(4,3, testdata20);
00262   const Matrix M21 = M1.above(M1);
00263   unit_assert( "above() ",   M20 == M21 );
00264   D testdata22[8]={1,2,3, 7, 4,5,6, 8};
00265   D testdata23[2]={7, 8};
00266   const Matrix M22(2,4, testdata22);
00267   const Matrix M23(2,1, testdata23);
00268   const Matrix M24 = M1.beside(M23);
00269   unit_assert( "beside() ",   M24 == M22 );
00270 
00271   Matrix M30 = M24;
00272   const Matrix M31 = M30.removeColumns(1);
00273   unit_assert( "removeColumns() ",   M31 == M1 );
00274   Matrix M32 = M20;  
00275   const Matrix M33 = M32.removeRows(2);
00276   unit_assert( "removeRows() ",  M33 == M1 );
00277    
00278   unit_pass();
00279 }
00280 
00281 DEFINE_TEST( check_matrix_operators ) {  
00282   cout << "\n -[ Matrix Operators (+ - * ^)]-\n";
00283   D testdata[6]={1,2,3, 4,5,6 };
00284   const Matrix M1(2,3, testdata);
00285   D testdata2[6]={1,4, 2,5, 3,6 };
00286   const Matrix M2(3,2, testdata2);  
00287   unit_assert( "^T ",  (M1^T) == M2 );
00288   D testdata3[6]={2,4,6, 8,10,12};
00289   Matrix M4(2,3,testdata3);    
00290   unit_assert( "+  ", M1+M1 == M4 );
00291   unit_assert( "-  ", M1+M1-M1 == M1 );
00292 
00293   D testdata4[6]={-0.5, -1, -1.5,  -2, -2.5, -3};
00294   Matrix M3(2,3, testdata4);
00295   unit_assert( "* scalar", M1*(-0.5) == M3 );
00296 
00297   D testdata5[6] = {2, -1,  0, 3,  2, -2};
00298   D testdata6[4] = {8, -1,  20, -1};
00299   Matrix M5 (3,2, testdata5);
00300   Matrix M6(2,2,testdata6);
00301   unit_assert( "*  ", M1*M5 == M6 );
00302 
00303   unit_assert( "^1 ", (M1^1) == M1 );
00304   M3.toId();
00305   unit_assert( "^0=id ", (M1^0) == M3 );
00306 
00307   D testdata7[16] = {1,2,3,4, -4,2,1,3, 0.3,-0.9, 4, -3, 1,0.5,0.3,5.0};
00308   Matrix M7(4,4,testdata7);
00309   unit_assert( "^1 * ^-1=id ",   comparetoidentity(M7*(M7^-1)) );
00310   unit_pass();
00311 }
00312 
00313 DEFINE_TEST( speed ) {  
00314   cout << "\n -[ Speed: Inverion]-\n";  
00315 #ifndef NDEBUG
00316   cout << "   DEBUG MODE! use -DNDEBUG -O3 (not -g) to get full performance\n";
00317 #endif
00318   Matrix M1;
00319   srand(time(0));
00320   D testdata0[9] = {1,2, -4,2}; 
00321   Matrix M2(2,2,testdata0);
00322   UNIT_MEASURE_START("2x2 Matrix inversion", 100000)
00323     M1 = (M2^-1);
00324   UNIT_MEASURE_STOP("");
00325   unit_assert( "validation", comparetoidentity(M1*M2));
00326   /* LU version:  555648/s */
00327   /* Explicit  : 1428775/s */
00328   D testdata1[9] = {1,2,3, -4,2,1, 0.3,-0.9}; 
00329   Matrix M3(3,3,testdata1);
00330   UNIT_MEASURE_START("3x3 Matrix inversion", 100000)
00331     M1 = (M3^-1);
00332   UNIT_MEASURE_STOP("");
00333   unit_assert( "validation", comparetoidentity(M1*M3));
00334 
00335   D testdata2[16] = {1,2,3,4, -4,2,1,3, 0.3,-0.9, 4, -3, 1,0.5,0.3,5.0};
00336   Matrix M4(4,4,testdata2);
00337   UNIT_MEASURE_START("4x4 Matrix inversion", 100000)
00338     M1 = (M4^-1);
00339   UNIT_MEASURE_STOP("");
00340   unit_assert( "validation", comparetoidentity(M1*M4));
00341 
00342   Matrix M20(20,20); 
00343   for (int i=0; i < M20.getM(); i++)  // define random values for initial matrix
00344     for (int j=0; j < M20.getN(); j++) {
00345       M20.val(i,j) = -22+(100. * rand())/RAND_MAX;
00346     }
00347   UNIT_MEASURE_START("20x20 Matrix inversion",1000)
00348     M1 = (M20^-1);  
00349   UNIT_MEASURE_STOP("");
00350   unit_assert( "validation", comparetoidentity(M1*M20));
00351 
00352   Matrix M200(200,200); 
00353   rand();  // eliminates the first (= zero) call
00354   for (int i=0; i < M200.getM(); i++)  // define random values for initial matrix
00355     for (int j=0; j < M200.getN(); j++) {
00356       M200.val(i,j) = -22+(100. * rand())/RAND_MAX;
00357     }
00358   UNIT_MEASURE_START("200x200 Matrix inversion",2)
00359     M1 = (M200^-1);  
00360   UNIT_MEASURE_STOP("");
00361   unit_assert( "validation", comparetoidentity(M1*M200));
00362 
00363   cout << "\n -[ Speed: Other Operations]-\n";    
00364   UNIT_MEASURE_START("20x20 Matrix multiplication with assignment",5000)
00365     M1 = M20*M20;  
00366   UNIT_MEASURE_STOP("");
00367   UNIT_MEASURE_START("20x20 Matrix addition with assignment",100000)
00368     M1= (M1 + M20);  
00369   UNIT_MEASURE_STOP("");
00370   UNIT_MEASURE_START("20x20 Matrix inplace addition",100000)
00371     M1 += M20;  
00372   UNIT_MEASURE_STOP("");
00373   UNIT_MEASURE_START("20x20 Matrix transposition",100000)
00374     M1 += M20;  
00375   UNIT_MEASURE_STOP("");
00376 
00377   unit_pass();  
00378 }
00379 
00380 
00381 DEFINE_TEST( store_restore ) {  
00382   cout << "\n -[ Store and Restore]-\n";  
00383   Matrix M1(32,1);
00384   for(int i =0; i<32; i++){
00385     M1.val(0,0) = (double)rand()/RAND_MAX;
00386   }
00387   Matrix M2(32,2);
00388   for(int i =0; i<64; i++){
00389     M2.val(i%32,i/32) = (double)rand()/RAND_MAX;
00390   }
00391   FILE* f;
00392   f=fopen("test.dat","wb");
00393   M1.store(f);
00394   M2.store(f);
00395   fclose(f);
00396   f=fopen("test.dat","rb");  
00397   Matrix M3,M4;
00398   M3.restore(f); 
00399   M4.restore(f); 
00400   fclose(f);  
00401   unit_assert( "validation", (M1-M3).map(fabs).elementSum()==0);
00402   unit_assert( "validation", (M2-M4).map(fabs).elementSum()==0);
00403 
00404   unit_pass();  
00405 }
00406 
00407 
00408 
00409 UNIT_TEST_RUN( "Matrix Tests" )
00410   ADD_TEST( check_creation )
00411   ADD_TEST( check_vector_operation )
00412   ADD_TEST( check_matrix_operation )
00413   ADD_TEST( check_matrix_operators )
00414   ADD_TEST( speed )
00415   ADD_TEST( store_restore )
00416 
00417   UNIT_TEST_END
00418 
00419 #endif // UNITTEST