SHOGUN  v3.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SparseMatrixOperator.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Soumyajit De
8  */
9 
10 #include <shogun/lib/config.h>
11 #include <shogun/lib/SGVector.h>
14 #include <shogun/base/Parameter.h>
17 
18 namespace shogun
19 {
20 
21 template<class T>
23  : CMatrixOperator<T>()
24  {
25  init();
26  }
27 
28 template<class T>
30  : CMatrixOperator<T>(op.num_features),
31  m_operator(op)
32  {
33  init();
34  }
35 
36 template<class T>
40  {
41  init();
42 
43  typedef SGSparseVector<T> vector;
44  typedef SGSparseVectorEntry<T> entry;
45 
46  m_operator=SGSparseMatrix<T>(orig.m_operator.num_vectors, orig.m_operator.num_features);
47 
48  vector* rows=SG_MALLOC(vector, m_operator.num_features);
49  for (index_t i=0; i<m_operator.num_vectors; ++i)
50  {
51  entry* features=SG_MALLOC(entry, orig.m_operator[i].num_feat_entries);
52  for (index_t j=0; j<orig.m_operator[i].num_feat_entries; ++j)
53  {
54  features[j].feat_index=orig.m_operator[i].features[j].feat_index;
55  features[j].entry=orig.m_operator[i].features[j].entry;
56  }
57  rows[i].features=features;
58  rows[i].num_feat_entries=m_operator[i].num_feat_entries;
59  }
60  m_operator.sparse_matrix=rows;
61 
62  SG_SGCDEBUG("%s deep copy created (%p)\n", this->get_name(), this);
63  }
64 
65 template<class T>
67  {
68  CSGObject::set_generic<T>();
69 
70  this->m_parameters->add_vector(&m_operator.sparse_matrix,
71  &m_operator.num_vectors, "sparse_matrix",
72  "The sparse matrix of the linear operator.");
73  this->m_parameters->add(&m_operator.num_features,
74  "m_operator.num_features", "Number of features.");
75  }
76 
77 template<class T>
79  {
80  }
81 
82 template<class T>
84  {
85  return m_operator;
86  }
87 
88 #ifdef HAVE_EIGEN3
89 template<class T>
91  int64_t power) const
92  {
93  REQUIRE(power>0, "matrix-power is non-positive!\n");
94 
95  // create casted operator in bool for capturing the sparsity
97  =static_cast<CSparseMatrixOperator<bool>*>(*this);
98 
99  // eigen3 map for this sparse matrix in which we compute current power
100  Eigen::SparseMatrix<bool> current_power
102 
103  // final power of the matrix goes into this one
104  Eigen::SparseMatrix<bool> matrix_power;
105 
106  // compute matrix power with O(log(n)) matrix-multiplication!
107  // traverse the bits of the power and compute the powers of 2 which
108  // makes up this number. in the process multiply these to get the result
109  bool lsb=true;
110  while (power)
111  {
112  // if the current bit is on, it should contribute to the final result
113  if (1 & power)
114  {
115  if (lsb)
116  {
117  // if seeing a 1 for the first time, then this should be the first
118  // power we should consider
119  matrix_power=current_power;
120  lsb=false;
121  }
122  else
123  matrix_power=matrix_power*current_power;
124  }
125  power=power>>1;
126 
127  // save unnecessary matrix-multiplication
128  if (power)
129  current_power=current_power*current_power;
130  }
131 
132  // create the sparsity structure using the final power
133  int32_t* outerIndexPtr=const_cast<int32_t*>(matrix_power.outerIndexPtr());
134  int32_t* innerIndexPtr=const_cast<int32_t*>(matrix_power.innerIndexPtr());
135 
136  SG_UNREF(sp_str);
137 
138  return new SparsityStructure(outerIndexPtr, innerIndexPtr,
139  matrix_power.rows());
140  }
141 #else
142 template<class T>
144  int64_t power) const
145  {
146  SG_SWARNING("Eigen3 required\n");
147  return new SparsityStructure();
148  }
149 #endif // HAVE_EIGEN3
150 
151 template<> \
152 SparsityStructure* CSparseMatrixOperator<complex128_t>
153  ::get_sparsity_structure(int64_t power) const
154  {
155  SG_SERROR("Not supported for complex128_t\n");
156  return new SparsityStructure();
157  }
158 
159 template<class T>
161  {
162  REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
163 
164  const int32_t diag_size=m_operator.num_vectors>m_operator.num_features ?
165  m_operator.num_features : m_operator.num_vectors;
166 
167  SGVector<T> diag(diag_size);
168  diag.set_const(static_cast<T>(0));
169  for (index_t i=0; i<diag_size; ++i)
170  {
171  SGSparseVectorEntry<T>* current_row=m_operator[i].features;
172  for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
173  {
174  if (i==current_row[j].feat_index)
175  {
176  diag[i]=current_row[j].entry;
177  break;
178  }
179  }
180  }
181 
182  return diag;
183  }
184 
185 template<class T>
187  {
188  REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
189  REQUIRE(diag.vector, "Diagonal not initialized!\n");
190 
191  const int32_t diag_size=m_operator.num_vectors>m_operator.num_features ?
192  m_operator.num_features : m_operator.num_vectors;
193 
194  REQUIRE(diag_size==diag.vlen, "Dimension mismatch!\n");
195 
196  bool need_sorting=false;
197  for (index_t i=0; i<diag_size; ++i)
198  {
199  SGSparseVectorEntry<T>* current_row=m_operator[i].features;
200  bool inserted=false;
201  // we just change the entry if the diagonal element for this row exists
202  for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
203  {
204  if (i==current_row[j].feat_index)
205  {
206  current_row[j].entry=diag[i];
207  inserted=true;
208  break;
209  }
210  }
211 
212  // we create a new entry if the diagonal element for this row doesn't exist
213  if (!inserted)
214  {
215  index_t j=m_operator[i].num_feat_entries;
216  m_operator[i].num_feat_entries=j+1;
217  m_operator[i].features=SG_REALLOC(SGSparseVectorEntry<T>,
218  m_operator[i].features, j, j+1);
219  m_operator[i].features[j].feat_index=i;
220  m_operator[i].features[j].entry=diag[i];
221  need_sorting=true;
222  }
223  }
224 
225  if (need_sorting)
226  m_operator.sort_features();
227  }
228 
229 template<class T>
231  {
232  REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
233  REQUIRE(this->get_dimension()==b.vlen,
234  "Number of rows of vector must be equal to the "
235  "number of cols of the operator!\n");
236 
237  SGVector<T> result(m_operator.num_vectors);
238  result=m_operator*b;
239 
240  return result;
241  }
242 
243 #define UNDEFINED(type) \
244 template<> \
245 SGVector<type> CSparseMatrixOperator<type>::apply(SGVector<type> b) const \
246  { \
247  SG_SERROR("Not supported for %s\n", #type);\
248  return b; \
249  }
250 
251 UNDEFINED(bool)
252 UNDEFINED(char)
253 UNDEFINED(int8_t)
254 UNDEFINED(uint8_t)
255 UNDEFINED(int16_t)
256 UNDEFINED(uint16_t)
257 UNDEFINED(int32_t)
258 UNDEFINED(uint32_t)
259 UNDEFINED(int64_t)
260 UNDEFINED(uint64_t)
263 #undef UNDEFINED
264 
265 template class CSparseMatrixOperator<bool>;
266 template class CSparseMatrixOperator<char>;
267 template class CSparseMatrixOperator<int8_t>;
268 template class CSparseMatrixOperator<uint8_t>;
269 template class CSparseMatrixOperator<int16_t>;
270 template class CSparseMatrixOperator<uint16_t>;
271 template class CSparseMatrixOperator<int32_t>;
272 template class CSparseMatrixOperator<uint32_t>;
273 template class CSparseMatrixOperator<int64_t>;
274 template class CSparseMatrixOperator<uint64_t>;
275 template class CSparseMatrixOperator<float32_t>;
276 template class CSparseMatrixOperator<float64_t>;
277 template class CSparseMatrixOperator<floatmax_t>;
278 template class CSparseMatrixOperator<complex128_t>;
279 }

SHOGUN Machine Learning Toolbox - Documentation