SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SGMatrix.cpp
Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Written (W) 2011-2013 Heiko Strathmann
00008  * Written (W) 2012 Fernando Jose Iglesias Garcia
00009  * Written (W) 2010,2012 Soeren Sonnenburg
00010  * Copyright (C) 2010 Berlin Institute of Technology
00011  * Copyright (C) 2012 Soeren Sonnenburg
00012  */
00013 
00014 #include <shogun/lib/config.h>
00015 #include <shogun/io/SGIO.h>
00016 #include <shogun/io/File.h>
00017 #include <shogun/lib/SGMatrix.h>
00018 #include <shogun/lib/SGVector.h>
00019 #include <shogun/mathematics/Math.h>
00020 #include <shogun/mathematics/lapack.h>
00021 #include <shogun/lib/SGMatrixList.h>
00022 
00023 namespace shogun {
00024 
00025 template <class T>
00026 SGMatrix<T>::SGMatrix() : SGReferencedData()
00027 {
00028     init_data();
00029 }
00030 
00031 template <class T>
00032 SGMatrix<T>::SGMatrix(bool ref_counting) : SGReferencedData(ref_counting)
00033 {
00034     init_data();
00035 }
00036 
00037 template <class T>
00038 SGMatrix<T>::SGMatrix(T* m, index_t nrows, index_t ncols, bool ref_counting)
00039     : SGReferencedData(ref_counting), matrix(m),
00040     num_rows(nrows), num_cols(ncols) { }
00041 
00042 template <class T>
00043 SGMatrix<T>::SGMatrix(index_t nrows, index_t ncols, bool ref_counting)
00044     : SGReferencedData(ref_counting), num_rows(nrows), num_cols(ncols)
00045 {
00046     matrix=SG_MALLOC(T, ((int64_t) nrows)*ncols);
00047 }
00048 
00049 template <class T>
00050 SGMatrix<T>::SGMatrix(const SGMatrix &orig) : SGReferencedData(orig)
00051 {
00052     copy_data(orig);
00053 }
00054 
00055 template <class T>
00056 SGMatrix<T>::~SGMatrix()
00057 {
00058     unref();
00059 }
00060 
00061 template <class T>
00062 bool SGMatrix<T>::operator==(SGMatrix<T>& other)
00063 {
00064     if (num_rows!=other.num_rows || num_cols!=other.num_cols)
00065         return false;
00066 
00067     if (matrix!=other.matrix)
00068         return false;
00069 
00070     return true;
00071 }
00072 
00073 template <class T>
00074 bool SGMatrix<T>::equals(SGMatrix<T>& other)
00075 {
00076     if (num_rows!=other.num_rows || num_cols!=other.num_cols)
00077         return false;
00078 
00079     for (int64_t i=0; i<int64_t(num_rows)*num_cols; ++i)
00080     {
00081         if (matrix[i]!=other.matrix[i])
00082             return false;
00083     }
00084 
00085     return true;
00086 }
00087 
00088 template <class T>
00089 void SGMatrix<T>::set_const(T const_elem)
00090 {
00091     for (int64_t i=0; i<int64_t(num_rows)*num_cols; i++)
00092         matrix[i]=const_elem ;
00093 }
00094 
00095 template <class T>
00096 void SGMatrix<T>::zero()
00097 {
00098     if (matrix && (int64_t(num_rows)*num_cols))
00099         set_const(0);
00100 }
00101 
00102 template <>
00103 void SGMatrix<complex128_t>::zero()
00104 {
00105     if (matrix && (int64_t(num_rows)*num_cols))
00106         set_const(complex128_t(0.0));
00107 }
00108 
00109 template <class T>
00110 T SGMatrix<T>::max_single()
00111 {
00112     T max=matrix[0];
00113     for (int64_t i=1; i<int64_t(num_rows)*num_cols; ++i)
00114     {
00115         if (matrix[i]>max)
00116             max=matrix[i];
00117     }
00118 
00119     return max;
00120 }
00121 
00122 template <>
00123 complex128_t SGMatrix<complex128_t>::max_single()
00124 {
00125     SG_SERROR("SGMatrix::max_single():: Not supported for complex128_t\n");
00126     return complex128_t(0.0);
00127 }
00128 
00129 template <class T>
00130 SGMatrix<T> SGMatrix<T>::clone()
00131 {
00132     return SGMatrix<T>(clone_matrix(matrix, num_rows, num_cols),
00133             num_rows, num_cols);
00134 }
00135 
00136 template <class T>
00137 T* SGMatrix<T>::clone_matrix(const T* matrix, int32_t nrows, int32_t ncols)
00138 {
00139     T* result = SG_MALLOC(T, int64_t(nrows)*ncols);
00140     for (int64_t i=0; i<int64_t(nrows)*ncols; i++)
00141         result[i]=matrix[i];
00142 
00143     return result;
00144 }
00145 
00146 template <class T>
00147 void SGMatrix<T>::transpose_matrix(
00148     T*& matrix, int32_t& num_feat, int32_t& num_vec)
00149 {
00150     /* this should be done in-place! Heiko */
00151     T* transposed=SG_MALLOC(T, int64_t(num_vec)*num_feat);
00152     for (int64_t i=0; i<num_vec; i++)
00153     {
00154         for (int64_t j=0; j<num_feat; j++)
00155             transposed[i+j*num_vec]=matrix[i*num_feat+j];
00156     }
00157 
00158     SG_FREE(matrix);
00159     matrix=transposed;
00160 
00161     CMath::swap(num_feat, num_vec);
00162 }
00163 
00164 template <class T>
00165 void SGMatrix<T>::create_diagonal_matrix(T* matrix, T* v,int32_t size)
00166 {
00167     for(int64_t i=0;i<size;i++)
00168     {
00169         for(int64_t j=0;j<size;j++)
00170         {
00171             if(i==j)
00172                 matrix[j*size+i]=v[i];
00173             else
00174                 matrix[j*size+i]=0;
00175         }
00176     }
00177 }
00178 
00179 template <class T>
00180 float64_t SGMatrix<T>::trace(
00181     float64_t* mat, int32_t cols, int32_t rows)
00182 {
00183     float64_t trace=0;
00184     for (int64_t i=0; i<rows; i++)
00185         trace+=mat[i*cols+i];
00186     return trace;
00187 }
00188 
00189 template <class T>
00190 T* SGMatrix<T>::get_row_sum(T* matrix, int32_t m, int32_t n)
00191 {
00192     T* rowsums=SG_CALLOC(T, n);
00193 
00194     for (int64_t i=0; i<n; i++)
00195     {
00196         for (int64_t j=0; j<m; j++)
00197             rowsums[i]+=matrix[j+i*m];
00198     }
00199     return rowsums;
00200 }
00201 
00202 template <class T>
00203 T* SGMatrix<T>::get_column_sum(T* matrix, int32_t m, int32_t n)
00204 {
00205     T* colsums=SG_CALLOC(T, m);
00206 
00207     for (int64_t i=0; i<n; i++)
00208     {
00209         for (int64_t j=0; j<m; j++)
00210             colsums[j]+=matrix[j+i*m];
00211     }
00212     return colsums;
00213 }
00214 
00215 template <class T>
00216 void SGMatrix<T>::center()
00217 {
00218     center_matrix(matrix, num_rows, num_cols);
00219 }
00220 
00221 template <class T>
00222 void SGMatrix<T>::center_matrix(T* matrix, int32_t m, int32_t n)
00223 {
00224     float64_t num_data=n;
00225 
00226     T* colsums=get_column_sum(matrix, m,n);
00227     T* rowsums=get_row_sum(matrix, m,n);
00228 
00229     for (int32_t i=0; i<m; i++)
00230         colsums[i]/=num_data;
00231     for (int32_t j=0; j<n; j++)
00232         rowsums[j]/=num_data;
00233 
00234     T s=SGVector<T>::sum(rowsums, n)/num_data;
00235 
00236     for (int64_t i=0; i<n; i++)
00237     {
00238         for (int64_t j=0; j<m; j++)
00239             matrix[i*m+j]+=s-colsums[j]-rowsums[i];
00240     }
00241 
00242     SG_FREE(rowsums);
00243     SG_FREE(colsums);
00244 }
00245 
00246 template <class T>
00247 void SGMatrix<T>::remove_column_mean()
00248 {
00249     /* compute "row" sums (which I would call column sums), i.e. sum of all
00250      * elements in a fixed column */
00251     T* means=get_row_sum(matrix, num_rows, num_cols);
00252 
00253     /* substract column mean from every corresponding entry */
00254     for (int64_t i=0; i<num_cols; ++i)
00255     {
00256         means[i]/=num_rows;
00257         for (int64_t j=0; j<num_rows; ++j)
00258             matrix[i*num_rows+j]-=means[i];
00259     }
00260 
00261     SG_FREE(means);
00262 }
00263 
00264 template<class T> void SGMatrix<T>::display_matrix(const char* name) const
00265 {
00266     display_matrix(matrix, num_rows, num_cols, name);
00267 }
00268 
00269 template <class T>
00270 void SGMatrix<T>::display_matrix(
00271     const SGMatrix<T> matrix, const char* name,
00272     const char* prefix)
00273 {
00274     matrix.display_matrix();
00275 }
00276 
00277 template <>
00278 void SGMatrix<bool>::display_matrix(
00279     const bool* matrix, int32_t rows, int32_t cols, const char* name,
00280     const char* prefix)
00281 {
00282     ASSERT(rows>=0 && cols>=0)
00283     SG_SPRINT("%s%s=[\n", prefix, name)
00284     for (int64_t i=0; i<rows; i++)
00285     {
00286         SG_SPRINT("%s[", prefix)
00287         for (int64_t j=0; j<cols; j++)
00288             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i] ? 1 : 0,
00289                 j==cols-1? "" : ",");
00290         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00291     }
00292     SG_SPRINT("%s]\n", prefix)
00293 }
00294 
00295 template <>
00296 void SGMatrix<char>::display_matrix(
00297     const char* matrix, int32_t rows, int32_t cols, const char* name,
00298     const char* prefix)
00299 {
00300     ASSERT(rows>=0 && cols>=0)
00301     SG_SPRINT("%s%s=[\n", prefix, name)
00302     for (int64_t i=0; i<rows; i++)
00303     {
00304         SG_SPRINT("%s[", prefix)
00305         for (int64_t j=0; j<cols; j++)
00306             SG_SPRINT("%s\t%c%s", prefix, matrix[j*rows+i],
00307                 j==cols-1? "" : ",");
00308         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00309     }
00310     SG_SPRINT("%s]\n", prefix)
00311 }
00312 
00313 template <>
00314 void SGMatrix<int8_t>::display_matrix(
00315     const int8_t* matrix, int32_t rows, int32_t cols, const char* name,
00316     const char* prefix)
00317 {
00318     ASSERT(rows>=0 && cols>=0)
00319     SG_SPRINT("%s%s=[\n", prefix, name)
00320     for (int64_t i=0; i<rows; i++)
00321     {
00322         SG_SPRINT("%s[", prefix)
00323         for (int64_t j=0; j<cols; j++)
00324             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
00325                 j==cols-1? "" : ",");
00326         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00327     }
00328     SG_SPRINT("%s]\n", prefix)
00329 }
00330 
00331 template <>
00332 void SGMatrix<uint8_t>::display_matrix(
00333     const uint8_t* matrix, int32_t rows, int32_t cols, const char* name,
00334     const char* prefix)
00335 {
00336     ASSERT(rows>=0 && cols>=0)
00337     SG_SPRINT("%s%s=[\n", prefix, name)
00338     for (int64_t i=0; i<rows; i++)
00339     {
00340         SG_SPRINT("%s[", prefix)
00341         for (int64_t j=0; j<cols; j++)
00342             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
00343                 j==cols-1? "" : ",");
00344         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00345     }
00346     SG_SPRINT("%s]\n", prefix)
00347 }
00348 
00349 template <>
00350 void SGMatrix<int16_t>::display_matrix(
00351     const int16_t* matrix, int32_t rows, int32_t cols, const char* name,
00352     const char* prefix)
00353 {
00354     ASSERT(rows>=0 && cols>=0)
00355     SG_SPRINT("%s%s=[\n", prefix, name)
00356     for (int64_t i=0; i<rows; i++)
00357     {
00358         SG_SPRINT("%s[", prefix)
00359         for (int64_t j=0; j<cols; j++)
00360             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
00361                 j==cols-1? "" : ",");
00362         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00363     }
00364     SG_SPRINT("%s]\n", prefix)
00365 }
00366 
00367 template <>
00368 void SGMatrix<uint16_t>::display_matrix(
00369     const uint16_t* matrix, int32_t rows, int32_t cols, const char* name,
00370     const char* prefix)
00371 {
00372     ASSERT(rows>=0 && cols>=0)
00373     SG_SPRINT("%s%s=[\n", prefix, name)
00374     for (int64_t i=0; i<rows; i++)
00375     {
00376         SG_SPRINT("%s[", prefix)
00377         for (int64_t j=0; j<cols; j++)
00378             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
00379                 j==cols-1? "" : ",");
00380         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00381     }
00382     SG_SPRINT("%s]\n", prefix)
00383 }
00384 
00385 
00386 template <>
00387 void SGMatrix<int32_t>::display_matrix(
00388     const int32_t* matrix, int32_t rows, int32_t cols, const char* name,
00389     const char* prefix)
00390 {
00391     ASSERT(rows>=0 && cols>=0)
00392     SG_SPRINT("%s%s=[\n", prefix, name)
00393     for (int64_t i=0; i<rows; i++)
00394     {
00395         SG_SPRINT("%s[", prefix)
00396         for (int64_t j=0; j<cols; j++)
00397             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
00398                 j==cols-1? "" : ",");
00399         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00400     }
00401     SG_SPRINT("%s]\n", prefix)
00402 }
00403 
00404 template <>
00405 void SGMatrix<uint32_t>::display_matrix(
00406     const uint32_t* matrix, int32_t rows, int32_t cols, const char* name,
00407     const char* prefix)
00408 {
00409     ASSERT(rows>=0 && cols>=0)
00410     SG_SPRINT("%s%s=[\n", prefix, name)
00411     for (int64_t i=0; i<rows; i++)
00412     {
00413         SG_SPRINT("%s[", prefix)
00414         for (int64_t j=0; j<cols; j++)
00415             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
00416                 j==cols-1? "" : ",");
00417         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00418     }
00419     SG_SPRINT("%s]\n", prefix)
00420 }
00421 template <>
00422 void SGMatrix<int64_t>::display_matrix(
00423     const int64_t* matrix, int32_t rows, int32_t cols, const char* name,
00424     const char* prefix)
00425 {
00426     ASSERT(rows>=0 && cols>=0)
00427     SG_SPRINT("%s%s=[\n", prefix, name)
00428     for (int64_t i=0; i<rows; i++)
00429     {
00430         SG_SPRINT("%s[", prefix)
00431         for (int64_t j=0; j<cols; j++)
00432             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
00433                 j==cols-1? "" : ",");
00434         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00435     }
00436     SG_SPRINT("%s]\n", prefix)
00437 }
00438 
00439 template <>
00440 void SGMatrix<uint64_t>::display_matrix(
00441     const uint64_t* matrix, int32_t rows, int32_t cols, const char* name,
00442     const char* prefix)
00443 {
00444     ASSERT(rows>=0 && cols>=0)
00445     SG_SPRINT("%s%s=[\n", prefix, name)
00446     for (int64_t i=0; i<rows; i++)
00447     {
00448         SG_SPRINT("%s[", prefix)
00449         for (int64_t j=0; j<cols; j++)
00450             SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
00451                 j==cols-1? "" : ",");
00452         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00453     }
00454     SG_SPRINT("%s]\n", prefix)
00455 }
00456 
00457 template <>
00458 void SGMatrix<float32_t>::display_matrix(
00459     const float32_t* matrix, int32_t rows, int32_t cols, const char* name,
00460     const char* prefix)
00461 {
00462     ASSERT(rows>=0 && cols>=0)
00463     SG_SPRINT("%s%s=[\n", prefix, name)
00464     for (int64_t i=0; i<rows; i++)
00465     {
00466         SG_SPRINT("%s[", prefix)
00467         for (int64_t j=0; j<cols; j++)
00468             SG_SPRINT("%s\t%.18g%s", prefix, (float) matrix[j*rows+i],
00469                 j==cols-1? "" : ",");
00470         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00471     }
00472     SG_SPRINT("%s]\n", prefix)
00473 }
00474 
00475 template <>
00476 void SGMatrix<float64_t>::display_matrix(
00477     const float64_t* matrix, int32_t rows, int32_t cols, const char* name,
00478     const char* prefix)
00479 {
00480     ASSERT(rows>=0 && cols>=0)
00481     SG_SPRINT("%s%s=[\n", prefix, name)
00482     for (int64_t i=0; i<rows; i++)
00483     {
00484         SG_SPRINT("%s[", prefix)
00485         for (int64_t j=0; j<cols; j++)
00486             SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i],
00487                 j==cols-1? "" : ",");
00488         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00489     }
00490     SG_SPRINT("%s]\n", prefix)
00491 }
00492 
00493 template <>
00494 void SGMatrix<floatmax_t>::display_matrix(
00495     const floatmax_t* matrix, int32_t rows, int32_t cols, const char* name,
00496     const char* prefix)
00497 {
00498     ASSERT(rows>=0 && cols>=0)
00499     SG_SPRINT("%s%s=[\n", prefix, name)
00500     for (int64_t i=0; i<rows; i++)
00501     {
00502         SG_SPRINT("%s[", prefix)
00503         for (int64_t j=0; j<cols; j++)
00504             SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i],
00505                 j==cols-1? "" : ",");
00506         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00507     }
00508     SG_SPRINT("%s]\n", prefix)
00509 }
00510 
00511 template <>
00512 void SGMatrix<complex128_t>::display_matrix(
00513     const complex128_t* matrix, int32_t rows, int32_t cols, const char* name,
00514     const char* prefix)
00515 {
00516     ASSERT(rows>=0 && cols>=0)
00517     SG_SPRINT("%s%s=[\n", prefix, name)
00518     for (int64_t i=0; i<rows; i++)
00519     {
00520         SG_SPRINT("%s[", prefix)
00521         for (int64_t j=0; j<cols; j++)
00522             SG_SPRINT("%s\t(%.18g+i%.18g)%s", prefix, matrix[j*rows+i].real(),
00523                 matrix[j*rows+i].imag(), j==cols-1? "" : ",");
00524         SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
00525     }
00526     SG_SPRINT("%s]\n", prefix)
00527 }
00528 
00529 template <>
00530 SGMatrix<char> SGMatrix<char>::create_identity_matrix(index_t size, char scale)
00531 {
00532     SG_SNOTIMPLEMENTED
00533     return SGMatrix<char>();
00534 }
00535 
00536 template <>
00537 SGMatrix<int8_t> SGMatrix<int8_t>::create_identity_matrix(index_t size, int8_t scale)
00538 {
00539     SGMatrix<int8_t> I(size, size);
00540     for (index_t i=0; i<size; ++i)
00541     {
00542         for (index_t j=0; j<size; ++j)
00543             I(i,j)=i==j ? scale : 0.0;
00544     }
00545 
00546     return I;
00547 }
00548 
00549 template <>
00550 SGMatrix<uint8_t> SGMatrix<uint8_t>::create_identity_matrix(index_t size, uint8_t scale)
00551 {
00552     SGMatrix<uint8_t> I(size, size);
00553     for (index_t i=0; i<size; ++i)
00554     {
00555         for (index_t j=0; j<size; ++j)
00556             I(i,j)=i==j ? scale : 0.0;
00557     }
00558 
00559     return I;
00560 }
00561 
00562 template <>
00563 SGMatrix<bool> SGMatrix<bool>::create_identity_matrix(index_t size, bool scale)
00564 {
00565     SGMatrix<bool> I(size, size);
00566     for (index_t i=0; i<size; ++i)
00567     {
00568         for (index_t j=0; j<size; ++j)
00569             I(i,j)=i==j ? scale : (!scale);
00570     }
00571 
00572     return I;
00573 }
00574 
00575 template <>
00576 SGMatrix<int16_t> SGMatrix<int16_t>::create_identity_matrix(index_t size, int16_t scale)
00577 {
00578     SGMatrix<int16_t> I(size, size);
00579     for (index_t i=0; i<size; ++i)
00580     {
00581         for (index_t j=0; j<size; ++j)
00582             I(i,j)=i==j ? scale : 0.0;
00583     }
00584 
00585     return I;
00586 }
00587 
00588 template <>
00589 SGMatrix<uint16_t> SGMatrix<uint16_t>::create_identity_matrix(index_t size, uint16_t scale)
00590 {
00591     SGMatrix<uint16_t> I(size, size);
00592     for (index_t i=0; i<size; ++i)
00593     {
00594         for (index_t j=0; j<size; ++j)
00595             I(i,j)=i==j ? scale : 0.0;
00596     }
00597 
00598     return I;
00599 }
00600 
00601 template <>
00602 SGMatrix<int32_t> SGMatrix<int32_t>::create_identity_matrix(index_t size, int32_t scale)
00603 {
00604     SGMatrix<int32_t> I(size, size);
00605     for (index_t i=0; i<size; ++i)
00606     {
00607         for (index_t j=0; j<size; ++j)
00608             I(i,j)=i==j ? scale : 0.0;
00609     }
00610 
00611     return I;
00612 }
00613 
00614 template <>
00615 SGMatrix<uint32_t> SGMatrix<uint32_t>::create_identity_matrix(index_t size, uint32_t scale)
00616 {
00617     SGMatrix<uint32_t> I(size, size);
00618     for (index_t i=0; i<size; ++i)
00619     {
00620         for (index_t j=0; j<size; ++j)
00621             I(i,j)=i==j ? scale : 0.0;
00622     }
00623 
00624     return I;
00625 }
00626 
00627 template <>
00628 SGMatrix<int64_t> SGMatrix<int64_t>::create_identity_matrix(index_t size, int64_t scale)
00629 {
00630     SGMatrix<int64_t> I(size, size);
00631     for (index_t i=0; i<size; ++i)
00632     {
00633         for (index_t j=0; j<size; ++j)
00634             I(i,j)=i==j ? scale : 0.0;
00635     }
00636 
00637     return I;
00638 }
00639 
00640 template <>
00641 SGMatrix<uint64_t> SGMatrix<uint64_t>::create_identity_matrix(index_t size, uint64_t scale)
00642 {
00643     SGMatrix<uint64_t> I(size, size);
00644     for (index_t i=0; i<size; ++i)
00645     {
00646         for (index_t j=0; j<size; ++j)
00647             I(i,j)=i==j ? scale : 0.0;
00648     }
00649 
00650     return I;
00651 }
00652 
00653 template <>
00654 SGMatrix<float32_t> SGMatrix<float32_t>::create_identity_matrix(index_t size, float32_t scale)
00655 {
00656     SGMatrix<float32_t> I(size, size);
00657     for (index_t i=0; i<size; ++i)
00658     {
00659         for (index_t j=0; j<size; ++j)
00660             I(i,j)=i==j ? scale : 0.0;
00661     }
00662 
00663     return I;
00664 }
00665 
00666 template <>
00667 SGMatrix<float64_t> SGMatrix<float64_t>::create_identity_matrix(index_t size, float64_t scale)
00668 {
00669     SGMatrix<float64_t> I(size, size);
00670     for (index_t i=0; i<size; ++i)
00671     {
00672         for (index_t j=0; j<size; ++j)
00673             I(i,j)=i==j ? scale : 0.0;
00674     }
00675 
00676     return I;
00677 }
00678 
00679 template <>
00680 SGMatrix<floatmax_t> SGMatrix<floatmax_t>::create_identity_matrix(index_t size, floatmax_t scale)
00681 {
00682     SGMatrix<floatmax_t> I(size, size);
00683     for (index_t i=0; i<size; ++i)
00684     {
00685         for (index_t j=0; j<size; ++j)
00686             I(i,j)=i==j ? scale : 0.0;
00687     }
00688 
00689     return I;
00690 }
00691 
00692 template <>
00693 SGMatrix<complex128_t> SGMatrix<complex128_t>::create_identity_matrix(index_t size, complex128_t scale)
00694 {
00695     SGMatrix<complex128_t> I(size, size);
00696     for (index_t i=0; i<size; ++i)
00697     {
00698         for (index_t j=0; j<size; ++j)
00699             I(i,j)=i==j ? scale : complex128_t(0.0);
00700     }
00701 
00702     return I;
00703 }
00704 
00705 
00706 template <class T>
00707 SGMatrix<float64_t> SGMatrix<T>::create_centering_matrix(index_t size)
00708 {
00709     SGMatrix<float64_t> H=SGMatrix<float64_t>::create_identity_matrix(size, 1.0);
00710 
00711     float64_t subtract=1.0/size;
00712     for (index_t i=0; i<size; ++i)
00713     {
00714         for (index_t j=0; j<0; ++j)
00715             H(i,j)-=subtract;
00716     }
00717 
00718     return H;
00719 }
00720 
00721 //Howto construct the pseudo inverse (from "The Matrix Cookbook")
00722 //
00723 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
00724 //
00725 //The matrix A+ known as the pseudo inverse is unique and does always exist.
00726 //
00727 //The pseudo inverse A+ can be constructed from the singular value
00728 //decomposition A = UDV^T , by  A^+ = V(D+)U^T.
00729 
00730 #ifdef HAVE_LAPACK
00731 template <class T>
00732 float64_t* SGMatrix<T>::pinv(
00733         float64_t* matrix, int32_t rows, int32_t cols, float64_t* target)
00734 {
00735     if (!target)
00736         target=SG_MALLOC(float64_t, rows*cols);
00737 
00738     char jobu='A';
00739     char jobvt='A';
00740     int m=rows; /* for calling external lib */
00741     int n=cols; /* for calling external lib */
00742     int lda=m; /* for calling external lib */
00743     int ldu=m; /* for calling external lib */
00744     int ldvt=n; /* for calling external lib */
00745     int info=-1; /* for calling external lib */
00746     int32_t lsize=CMath::min((int32_t) m, (int32_t) n);
00747     double* s=SG_MALLOC(double, lsize);
00748     double* u=SG_MALLOC(double, m*m);
00749     double* vt=SG_MALLOC(double, n*n);
00750 
00751     wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
00752     ASSERT(info==0)
00753 
00754     for (int64_t i=0; i<n; i++)
00755     {
00756         for (int64_t j=0; j<lsize; j++)
00757             vt[i*n+j]=vt[i*n+j]/s[j];
00758     }
00759 
00760     cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
00761 
00762     SG_FREE(u);
00763     SG_FREE(vt);
00764     SG_FREE(s);
00765 
00766     return target;
00767 }
00768 
00770 template <class T>
00771 void SGMatrix<T>::inverse(SGMatrix<float64_t> matrix)
00772 {
00773     ASSERT(matrix.num_cols==matrix.num_rows)
00774     int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols);
00775     clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
00776     clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
00777     SG_FREE(ipiv);
00778 }
00779 
00780 template <class T>
00781 SGVector<float64_t> SGMatrix<T>::compute_eigenvectors(SGMatrix<float64_t> matrix)
00782 {
00783     if (matrix.num_rows!=matrix.num_rows)
00784     {
00785         SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix"
00786                 " rows and columns are not equal!\n");
00787     }
00788 
00789     /* use reference counting for SGVector */
00790     SGVector<float64_t> result(NULL, 0, true);
00791     result.vlen=matrix.num_rows;
00792     result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows,
00793             matrix.num_rows);
00794     return result;
00795 }
00796 
00797 template <class T>
00798 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m)
00799 {
00800     ASSERT(n == m)
00801 
00802     char V='V';
00803     char U='U';
00804     int info;
00805     int ord=n;
00806     int lda=n;
00807     double* eigenvalues=SG_CALLOC(float64_t, n+1);
00808 
00809     // lapack sym matrix eigenvalues+vectors
00810     wrap_dsyev(V, U,  ord, matrix, lda,
00811             eigenvalues, &info);
00812 
00813     if (info!=0)
00814         SG_SERROR("DSYEV failed with code %d\n", info)
00815 
00816     return eigenvalues;
00817 }
00818 
00819 template <class T>
00820 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors,
00821                                            int n, int il, int iu)
00822 {
00823     eigenvalues = SG_MALLOC(double, n);
00824     eigenvectors = SG_MALLOC(double, (iu-il+1)*n);
00825     int status = 0;
00826     wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status);
00827     ASSERT(status==0)
00828 }
00829 
00830 #endif //HAVE_LAPACK
00831 
00832 template <class T>
00833 SGMatrix<float64_t> SGMatrix<T>::matrix_multiply(
00834         SGMatrix<float64_t> A, SGMatrix<float64_t> B,
00835         bool transpose_A, bool transpose_B, float64_t scale)
00836 {
00837     /* these variables store size of transposed matrices*/
00838     index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
00839     index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
00840     index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
00841     index_t cols_B=transpose_B ? B.num_rows : B.num_cols;
00842 
00843     /* do a dimension check */
00844     if (cols_A!=rows_B)
00845     {
00846         SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: "
00847                 "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
00848     }
00849 
00850     /* allocate result matrix */
00851     SGMatrix<float64_t> C(rows_A, cols_B);
00852     C.zero();
00853 #ifdef HAVE_LAPACK
00854     /* multiply */
00855     cblas_dgemm(CblasColMajor,
00856             transpose_A ? CblasTrans : CblasNoTrans,
00857             transpose_B ? CblasTrans : CblasNoTrans,
00858             rows_A, cols_B, cols_A, scale,
00859             A.matrix, A.num_rows, B.matrix, B.num_rows,
00860             0.0, C.matrix, C.num_rows);
00861 #else
00862     for (int32_t i=0; i<rows_A; i++)
00863     {
00864         for (int32_t j=0; j<cols_B; j++)
00865         {
00866             for (int32_t k=0; k<cols_A; k++)
00867                 C(i,j) += A(i,k)*B(k,j);
00868         }
00869     }
00870 #endif //HAVE_LAPACK
00871 
00872     return C;
00873 }
00874 
00875 template<class T>
00876 SGMatrix<T> SGMatrix<T>::get_allocated_matrix(index_t num_rows,
00877         index_t num_cols, SGMatrix<T> pre_allocated)
00878 {
00879     SGMatrix<T> result;
00880 
00881     /* evtl use pre-allocated space */
00882     if (pre_allocated.matrix)
00883     {
00884         result=pre_allocated;
00885 
00886         /* check dimension consistency */
00887         if (pre_allocated.num_rows!=num_rows ||
00888                 pre_allocated.num_cols!=num_cols)
00889         {
00890             SG_SERROR("SGMatrix<T>::get_allocated_matrix(). Provided target"
00891                     "matrix dimensions (%dx%d) do not match passed data "
00892                     "dimensions (%dx%d)!\n", pre_allocated.num_rows,
00893                     pre_allocated.num_cols, num_rows, num_cols);
00894         }
00895     }
00896     else
00897     {
00898         /* otherwise, allocate space */
00899         result=SGMatrix<T>(num_rows, num_cols);
00900     }
00901 
00902     return result;
00903 }
00904 
00905 template<class T>
00906 void SGMatrix<T>::copy_data(const SGReferencedData &orig)
00907 {
00908     matrix=((SGMatrix*)(&orig))->matrix;
00909     num_rows=((SGMatrix*)(&orig))->num_rows;
00910     num_cols=((SGMatrix*)(&orig))->num_cols;
00911 }
00912 
00913 template<class T>
00914 void SGMatrix<T>::init_data()
00915 {
00916     matrix=NULL;
00917     num_rows=0;
00918     num_cols=0;
00919 }
00920 
00921 template<class T>
00922 void SGMatrix<T>::free_data()
00923 {
00924     SG_FREE(matrix);
00925     matrix=NULL;
00926     num_rows=0;
00927     num_cols=0;
00928 }
00929 
00930 template<class T>
00931 void SGMatrix<T>::load(CFile* loader)
00932 {
00933     ASSERT(loader)
00934     unref();
00935 
00936     SG_SET_LOCALE_C;
00937     SGMatrix<T> mat;
00938     loader->get_matrix(mat.matrix, mat.num_rows, mat.num_cols);
00939     copy_data(mat);
00940     copy_refcount(mat);
00941     ref();
00942     SG_RESET_LOCALE;
00943 }
00944 
00945 template<>
00946 void SGMatrix<complex128_t>::load(CFile* loader)
00947 {
00948     SG_SERROR("SGMatrix::load():: Not supported for complex128_t\n");
00949 }
00950 
00951 template<class T>
00952 void SGMatrix<T>::save(CFile* writer)
00953 {
00954     ASSERT(writer)
00955     SG_SET_LOCALE_C;
00956     writer->set_matrix(matrix, num_rows, num_cols);
00957     SG_RESET_LOCALE;
00958 }
00959 
00960 template<>
00961 void SGMatrix<complex128_t>::save(CFile* saver)
00962 {
00963     SG_SERROR("SGMatrix::save():: Not supported for complex128_t\n");
00964 }
00965 
00966 template<class T>
00967 SGVector<T> SGMatrix<T>::get_row_vector(index_t row) const
00968 {
00969     SGVector<T> rowv(num_cols);
00970     for (int64_t i = 0; i < num_cols; i++)
00971     {
00972         rowv[i] = matrix[i*num_rows+row];
00973     }
00974     return rowv;
00975 }
00976 
00977 template<class T>
00978 SGVector<T> SGMatrix<T>::get_diagonal_vector() const
00979 {
00980     index_t diag_vlen=CMath::min(num_cols, num_rows);
00981     SGVector<T> diag(diag_vlen);
00982 
00983     for (int64_t i=0; i<diag_vlen; i++)
00984     {
00985         diag[i]=matrix[i*num_rows+i];
00986     }
00987 
00988     return diag;
00989 }
00990 
00991 template class SGMatrix<bool>;
00992 template class SGMatrix<char>;
00993 template class SGMatrix<int8_t>;
00994 template class SGMatrix<uint8_t>;
00995 template class SGMatrix<int16_t>;
00996 template class SGMatrix<uint16_t>;
00997 template class SGMatrix<int32_t>;
00998 template class SGMatrix<uint32_t>;
00999 template class SGMatrix<int64_t>;
01000 template class SGMatrix<uint64_t>;
01001 template class SGMatrix<float32_t>;
01002 template class SGMatrix<float64_t>;
01003 template class SGMatrix<floatmax_t>;
01004 template class SGMatrix<complex128_t>;
01005 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation