SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SparseMatrixOperator.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) 2013 Soumyajit De
00008  */
00009 
00010 #include <shogun/lib/config.h>
00011 #include <shogun/lib/SGVector.h>
00012 #include <shogun/lib/SGSparseMatrix.h>
00013 #include <shogun/lib/SGSparseVector.h>
00014 #include <shogun/base/Parameter.h>
00015 #include <shogun/mathematics/linalg/linop/SparseMatrixOperator.h>
00016 #include <shogun/mathematics/eigen3.h>
00017 
00018 namespace shogun
00019 {
00020 
00021 template<class T>
00022 CSparseMatrixOperator<T>::CSparseMatrixOperator()
00023     : CMatrixOperator<T>()
00024     {
00025         init();
00026     }
00027 
00028 template<class T>
00029 CSparseMatrixOperator<T>::CSparseMatrixOperator(SGSparseMatrix<T> op)
00030     : CMatrixOperator<T>(op.num_features),
00031       m_operator(op)
00032     {
00033         init();
00034     }
00035 
00036 template<class T>
00037 CSparseMatrixOperator<T>::CSparseMatrixOperator
00038     (const CSparseMatrixOperator<T>& orig)
00039     : CMatrixOperator<T>(orig.get_dimension())
00040     {
00041         init();
00042 
00043         typedef SGSparseVector<T> vector;
00044         typedef SGSparseVectorEntry<T> entry;
00045 
00046         m_operator=SGSparseMatrix<T>(orig.m_operator.num_vectors, orig.m_operator.num_features);
00047 
00048         vector* rows=SG_MALLOC(vector, m_operator.num_features);
00049         for (index_t i=0; i<m_operator.num_vectors; ++i)
00050         {
00051             entry* features=SG_MALLOC(entry, orig.m_operator[i].num_feat_entries);
00052             for (index_t j=0; j<orig.m_operator[i].num_feat_entries; ++j)
00053             {
00054                 features[j].feat_index=orig.m_operator[i].features[j].feat_index;
00055                 features[j].entry=orig.m_operator[i].features[j].entry;
00056             }
00057             rows[i].features=features;
00058             rows[i].num_feat_entries=m_operator[i].num_feat_entries;
00059         }
00060         m_operator.sparse_matrix=rows;
00061 
00062         SG_SGCDEBUG("%s deep copy created (%p)\n", this->get_name(), this);
00063     }
00064 
00065 template<class T>
00066 void CSparseMatrixOperator<T>::init()
00067     {
00068         CSGObject::set_generic<T>();
00069 
00070         this->m_parameters->add_vector(&m_operator.sparse_matrix,
00071                 &m_operator.num_vectors, "sparse_matrix",
00072                 "The sparse matrix of the linear operator.");
00073         this->m_parameters->add(&m_operator.num_features,
00074                 "m_operator.num_features", "Number of features.");
00075     }
00076 
00077 template<class T>
00078 CSparseMatrixOperator<T>::~CSparseMatrixOperator()
00079     {
00080     }
00081 
00082 template<class T>
00083 SGSparseMatrix<T> CSparseMatrixOperator<T>::get_matrix_operator() const
00084     {
00085         return m_operator;
00086     }
00087 
00088 #ifdef HAVE_EIGEN3
00089 template<class T>
00090 SparsityStructure* CSparseMatrixOperator<T>::get_sparsity_structure(
00091     int64_t power) const
00092     {
00093         REQUIRE(power>0, "matrix-power is non-positive!\n");
00094 
00095         // create casted operator in bool for capturing the sparsity
00096         CSparseMatrixOperator<bool>* sp_str
00097             =static_cast<CSparseMatrixOperator<bool>*>(*this);
00098 
00099         // eigen3 map for this sparse matrix in which we compute current power
00100         Eigen::SparseMatrix<bool> current_power
00101             =EigenSparseUtil<bool>::toEigenSparse(sp_str->get_matrix_operator());
00102 
00103         // final power of the matrix goes into this one
00104         Eigen::SparseMatrix<bool> matrix_power;
00105 
00106         // compute matrix power with O(log(n)) matrix-multiplication!
00107         // traverse the bits of the power and compute the powers of 2 which
00108         // makes up this number. in the process multiply these to get the result
00109         bool lsb=true;
00110         while (power)
00111         {
00112             // if the current bit is on, it should contribute to the final result
00113             if (1 & power)
00114             {
00115                 if (lsb)
00116                 {
00117                     // if seeing a 1 for the first time, then this should be the first
00118                     // power we should consider
00119                     matrix_power=current_power;
00120                     lsb=false;
00121                 }
00122                 else
00123                     matrix_power=matrix_power*current_power;
00124             }
00125             power=power>>1;
00126 
00127             // save unnecessary matrix-multiplication
00128             if (power)
00129                 current_power=current_power*current_power;
00130         }
00131 
00132         // create the sparsity structure using the final power
00133         int32_t* outerIndexPtr=const_cast<int32_t*>(matrix_power.outerIndexPtr());
00134         int32_t* innerIndexPtr=const_cast<int32_t*>(matrix_power.innerIndexPtr());
00135 
00136         SG_UNREF(sp_str);
00137 
00138         return new SparsityStructure(outerIndexPtr, innerIndexPtr,
00139             matrix_power.rows());
00140     }
00141 #else
00142 template<class T>
00143 SparsityStructure* CSparseMatrixOperator<T>::get_sparsity_structure(
00144     int64_t power) const
00145     {
00146         SG_SWARNING("Eigen3 required\n");
00147         return new SparsityStructure();
00148     }
00149 #endif // HAVE_EIGEN3
00150 
00151 template<> \
00152 SparsityStructure* CSparseMatrixOperator<complex128_t>
00153 	::get_sparsity_structure(int64_t power) const
00154   {
00155     SG_SERROR("Not supported for complex128_t\n");
00156     return new SparsityStructure();
00157   }
00158 
00159 template<class T>
00160 SGVector<T> CSparseMatrixOperator<T>::get_diagonal() const
00161     {
00162         REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
00163 
00164         const int32_t diag_size=m_operator.num_vectors>m_operator.num_features ?
00165             m_operator.num_features : m_operator.num_vectors;
00166 
00167         SGVector<T> diag(diag_size);
00168         diag.set_const(static_cast<T>(0));
00169         for (index_t i=0; i<diag_size; ++i)
00170         {
00171             SGSparseVectorEntry<T>* current_row=m_operator[i].features;
00172             for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
00173             {
00174                 if (i==current_row[j].feat_index)
00175                 {
00176                     diag[i]=current_row[j].entry;
00177                     break;
00178                 }
00179             }
00180         }
00181 
00182         return diag;
00183     }
00184 
00185 template<class T>
00186 void CSparseMatrixOperator<T>::set_diagonal(SGVector<T> diag)
00187     {
00188         REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
00189         REQUIRE(diag.vector, "Diagonal not initialized!\n");
00190 
00191         const int32_t diag_size=m_operator.num_vectors>m_operator.num_features ?
00192             m_operator.num_features : m_operator.num_vectors;
00193 
00194         REQUIRE(diag_size==diag.vlen, "Dimension mismatch!\n");
00195 
00196         bool need_sorting=false;
00197         for (index_t i=0; i<diag_size; ++i)
00198         {
00199             SGSparseVectorEntry<T>* current_row=m_operator[i].features;
00200             bool inserted=false;
00201             // we just change the entry if the diagonal element for this row exists
00202             for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
00203             {
00204                 if (i==current_row[j].feat_index)
00205                 {
00206                     current_row[j].entry=diag[i];
00207                     inserted=true;
00208                     break;
00209                 }
00210             }
00211 
00212             // we create a new entry if the diagonal element for this row doesn't exist
00213             if (!inserted)
00214             {
00215                 index_t j=m_operator[i].num_feat_entries;
00216                 m_operator[i].num_feat_entries=j+1;
00217                 m_operator[i].features=SG_REALLOC(SGSparseVectorEntry<T>,
00218                     m_operator[i].features, j, j+1);
00219                 m_operator[i].features[j].feat_index=i;
00220                 m_operator[i].features[j].entry=diag[i];
00221                 need_sorting=true;
00222             }
00223         }
00224 
00225         if (need_sorting)
00226             m_operator.sort_features();
00227     }
00228 
00229 template<class T>
00230 SGVector<T> CSparseMatrixOperator<T>::apply(SGVector<T> b) const
00231     {
00232         REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
00233         REQUIRE(this->get_dimension()==b.vlen,
00234             "Number of rows of vector must be equal to the "
00235             "number of cols of the operator!\n");
00236 
00237         SGVector<T> result(m_operator.num_vectors);
00238         result=m_operator*b;
00239 
00240         return result;
00241     }
00242 
00243 #define UNDEFINED(type) \
00244 template<> \
00245 SGVector<type> CSparseMatrixOperator<type>::apply(SGVector<type> b) const \
00246     {   \
00247         SG_SERROR("Not supported for %s\n", #type);\
00248         return b; \
00249     }
00250 
00251 UNDEFINED(bool)
00252 UNDEFINED(char)
00253 UNDEFINED(int8_t)
00254 UNDEFINED(uint8_t)
00255 UNDEFINED(int16_t)
00256 UNDEFINED(uint16_t)
00257 UNDEFINED(int32_t)
00258 UNDEFINED(uint32_t)
00259 UNDEFINED(int64_t)
00260 UNDEFINED(uint64_t)
00261 UNDEFINED(float32_t)
00262 UNDEFINED(floatmax_t)
00263 #undef UNDEFINED
00264 
00265 template class CSparseMatrixOperator<bool>;
00266 template class CSparseMatrixOperator<char>;
00267 template class CSparseMatrixOperator<int8_t>;
00268 template class CSparseMatrixOperator<uint8_t>;
00269 template class CSparseMatrixOperator<int16_t>;
00270 template class CSparseMatrixOperator<uint16_t>;
00271 template class CSparseMatrixOperator<int32_t>;
00272 template class CSparseMatrixOperator<uint32_t>;
00273 template class CSparseMatrixOperator<int64_t>;
00274 template class CSparseMatrixOperator<uint64_t>;
00275 template class CSparseMatrixOperator<float32_t>;
00276 template class CSparseMatrixOperator<float64_t>;
00277 template class CSparseMatrixOperator<floatmax_t>;
00278 template class CSparseMatrixOperator<complex128_t>;
00279 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation