SHOGUN  5.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
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 template<class T>
90  int64_t power) const
91  {
92  REQUIRE(power>0, "matrix-power is non-positive!\n");
93 
94  // create casted operator in bool for capturing the sparsity
96  =static_cast<CSparseMatrixOperator<bool>*>(*this);
97 
98  // eigen3 map for this sparse matrix in which we compute current power
99  Eigen::SparseMatrix<bool> current_power
101 
102  // final power of the matrix goes into this one
103  Eigen::SparseMatrix<bool> matrix_power;
104 
105  // compute matrix power with O(log(n)) matrix-multiplication!
106  // traverse the bits of the power and compute the powers of 2 which
107  // makes up this number. in the process multiply these to get the result
108  bool lsb=true;
109  while (power)
110  {
111  // if the current bit is on, it should contribute to the final result
112  if (1 & power)
113  {
114  if (lsb)
115  {
116  // if seeing a 1 for the first time, then this should be the first
117  // power we should consider
118  matrix_power=current_power;
119  lsb=false;
120  }
121  else
122  matrix_power=matrix_power*current_power;
123  }
124  power=power>>1;
125 
126  // save unnecessary matrix-multiplication
127  if (power)
128  current_power=current_power*current_power;
129  }
130 
131  // create the sparsity structure using the final power
132  int32_t* outerIndexPtr=const_cast<int32_t*>(matrix_power.outerIndexPtr());
133  int32_t* innerIndexPtr=const_cast<int32_t*>(matrix_power.innerIndexPtr());
134 
135  SG_UNREF(sp_str);
136 
137  return new SparsityStructure(outerIndexPtr, innerIndexPtr,
138  matrix_power.rows());
139  }
140 
141 template<> \
142 SparsityStructure* CSparseMatrixOperator<complex128_t>
143  ::get_sparsity_structure(int64_t power) const
144  {
145  SG_SERROR("Not supported for complex128_t\n");
146  return new SparsityStructure();
147  }
148 
149 template<class T>
151  {
152  REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
153 
154  const int32_t diag_size=m_operator.num_vectors>m_operator.num_features ?
155  m_operator.num_features : m_operator.num_vectors;
156 
157  SGVector<T> diag(diag_size);
158  diag.set_const(static_cast<T>(0));
159  for (index_t i=0; i<diag_size; ++i)
160  {
161  SGSparseVectorEntry<T>* current_row=m_operator[i].features;
162  for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
163  {
164  if (i==current_row[j].feat_index)
165  {
166  diag[i]=current_row[j].entry;
167  break;
168  }
169  }
170  }
171 
172  return diag;
173  }
174 
175 template<class T>
177  {
178  REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
179  REQUIRE(diag.vector, "Diagonal not initialized!\n");
180 
181  const int32_t diag_size=m_operator.num_vectors>m_operator.num_features ?
182  m_operator.num_features : m_operator.num_vectors;
183 
184  REQUIRE(diag_size==diag.vlen, "Dimension mismatch!\n");
185 
186  bool need_sorting=false;
187  for (index_t i=0; i<diag_size; ++i)
188  {
189  SGSparseVectorEntry<T>* current_row=m_operator[i].features;
190  bool inserted=false;
191  // we just change the entry if the diagonal element for this row exists
192  for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
193  {
194  if (i==current_row[j].feat_index)
195  {
196  current_row[j].entry=diag[i];
197  inserted=true;
198  break;
199  }
200  }
201 
202  // we create a new entry if the diagonal element for this row doesn't exist
203  if (!inserted)
204  {
205  index_t j=m_operator[i].num_feat_entries;
206  m_operator[i].num_feat_entries=j+1;
207  m_operator[i].features=SG_REALLOC(SGSparseVectorEntry<T>,
208  m_operator[i].features, j, j+1);
209  m_operator[i].features[j].feat_index=i;
210  m_operator[i].features[j].entry=diag[i];
211  need_sorting=true;
212  }
213  }
214 
215  if (need_sorting)
216  m_operator.sort_features();
217  }
218 
219 template<class T>
221  {
222  REQUIRE(m_operator.sparse_matrix, "Operator not initialized!\n");
223  REQUIRE(this->get_dimension()==b.vlen,
224  "Number of rows of vector must be equal to the "
225  "number of cols of the operator!\n");
226 
227  SGVector<T> result(m_operator.num_vectors);
228  result=m_operator*b;
229 
230  return result;
231  }
232 
233 #define UNDEFINED(type) \
234 template<> \
235 SGVector<type> CSparseMatrixOperator<type>::apply(SGVector<type> b) const \
236  { \
237  SG_SERROR("Not supported for %s\n", #type);\
238  return b; \
239  }
240 
241 UNDEFINED(bool)
242 UNDEFINED(char)
243 UNDEFINED(int8_t)
244 UNDEFINED(uint8_t)
245 UNDEFINED(int16_t)
246 UNDEFINED(uint16_t)
247 UNDEFINED(int32_t)
248 UNDEFINED(uint32_t)
249 UNDEFINED(int64_t)
250 UNDEFINED(uint64_t)
253 #undef UNDEFINED
254 
255 template class CSparseMatrixOperator<bool>;
256 template class CSparseMatrixOperator<char>;
257 template class CSparseMatrixOperator<int8_t>;
258 template class CSparseMatrixOperator<uint8_t>;
259 template class CSparseMatrixOperator<int16_t>;
260 template class CSparseMatrixOperator<uint16_t>;
261 template class CSparseMatrixOperator<int32_t>;
262 template class CSparseMatrixOperator<uint32_t>;
263 template class CSparseMatrixOperator<int64_t>;
264 template class CSparseMatrixOperator<uint64_t>;
265 template class CSparseMatrixOperator<float32_t>;
266 template class CSparseMatrixOperator<float64_t>;
267 template class CSparseMatrixOperator<floatmax_t>;
268 template class CSparseMatrixOperator<complex128_t>;
269 }
SparsityStructure * get_sparsity_structure(int64_t power=1) const
template class SGSparseMatrix
const index_t get_dimension() const
SGSparseMatrix< T > get_matrix_operator() const
int32_t index_t
Definition: common.h:62
#define REQUIRE(x,...)
Definition: SGIO.h:206
virtual SGVector< T > get_diagonal() const
index_t vlen
Definition: SGVector.h:494
virtual SGVector< T > apply(SGVector< T > b) const
shogun vector
static Eigen::SparseMatrix< T > toEigenSparse(SGSparseMatrix< T > sg_matrix)
Definition: eigen3.cpp:22
long double floatmax_t
Definition: common.h:51
SGSparseVector< T > * sparse_matrix
array of sparse vectors of size num_vectors
virtual void set_diagonal(SGVector< T > diag)
float float32_t
Definition: common.h:49
#define SG_UNREF(x)
Definition: SGObject.h:55
Struct that represents the sparsity structure of the Sparse Matrix in CRS. Implementation has been ad...
#define UNDEFINED(type)
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
#define SG_SGCDEBUG(...)
Definition: SGIO.h:163
Class that represents a sparse-matrix linear operator. It computes matrix-vector product in its appl...
template class SGSparseVectorEntry
Definition: File.h:23
#define SG_SERROR(...)
Definition: SGIO.h:179
template class SGSparseVector The assumtion is that the stored SGSparseVectorEntry* vector is orde...
Abstract base class that represents a matrix linear operator. It provides an interface to computes ma...
void set_const(T const_elem)
Definition: SGVector.cpp:150

SHOGUN Machine Learning Toolbox - Documentation