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