SHOGUN
v3.2.0
|
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 }