Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
sparse_dense.hpp
1 /*
2  This file is part of the Kernel Quantum Probability library (KQP).
3 
4  KQP is free software: you can redistribute it and/or modify
5  it under the terms of the GNU General Public License as published by
6  the Free Software Foundation, either version 3 of the License, or
7  (at your option) any later version.
8 
9  KQP is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with KQP. If not, see <http://www.gnu.org/licenses/>.
16  */
17 #ifndef __KQP_SPARSE_DENSE_FEATURE_MATRIX_H__
18 #define __KQP_SPARSE_DENSE_FEATURE_MATRIX_H__
19 
20 #include <boost/serialization/map.hpp>
21 
22 #include <boost/lexical_cast.hpp>
23 #include <boost/function_output_iterator.hpp>
24 
25 #include <kqp/subset.hpp>
26 #include <kqp/feature_matrix.hpp>
27 #include <Eigen/Sparse>
28 
29 
30 namespace kqp {
31 
42  template <typename Scalar>
43  class SparseDense : public FeatureMatrixBase<Scalar> {
44  public:
45  KQP_SCALAR_TYPEDEFS(Scalar);
46  typedef SparseDense<Scalar> Self;
47 
48  typedef std::map<Index,Index> RowMap;
49 
50  virtual ~SparseDense() {}
51 
52  SparseDense() : m_dimension(0) {}
53  SparseDense(Index dimension) : m_dimension(dimension) {}
54  SparseDense(const Self &other) : m_dimension(other.m_dimension), m_map(other.m_map), m_matrix(other.m_matrix), m_gramMatrix(other.m_gramMatrix) {}
55 
57  static FMatrixPtr create(const ScalarMatrix &m) {
58  return FMatrixPtr(new Self(m));
59  }
60 
62  static FMatrixPtr create(const Eigen::SparseMatrix<Scalar, Eigen::ColMajor> &m, double threshold = epsilon()) {
63  return FMatrixPtr(new Self(m, threshold));
64  }
65 
67  static FMatrixPtr create(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &m, double threshold = epsilon()) {
68  return FMatrixPtr(new Self(m, threshold));
69  }
70 
71 
72 #ifndef SWIG
73  SparseDense(Index dimension, RowMap &&map, ScalarMatrix &&matrix)
74  : m_dimension(dimension), m_map(std::move(map)), m_matrix(std::move(matrix)) {
75 
76  }
77 #endif
78 
80  SparseDense(const Eigen::SparseMatrix<Scalar, Eigen::ColMajor> &mat, double threshold = epsilon()) : m_dimension(mat.rows()) {
81  // --- Compute which rows we need
82 
83  // Computing the column norms
84  RealVector norms(mat.cols());
85  norms.setZero();
86  for (Index k=0; k<mat.cols(); ++k) { // Loop on cols
87  // FIXME: use norm() when Eigen fixed
88  // norms[k] = mat.innerVector(k).norm();
89  for (typename Eigen::SparseMatrix<Scalar, Eigen::ColMajor>::InnerIterator it(mat,k); it; ++it) { // Loop on rows
90  norms[k] += Eigen::internal::abs2(it.value());
91  }
92  }
93 
94  norms = norms.cwiseSqrt();
95 
96  // Computing selected rows map
97  Index rows = 0;
98  for (int k=0; k<mat.outerSize(); ++k) { // Loop on cols
99  for (typename Eigen::SparseMatrix<Scalar, Eigen::ColMajor>::InnerIterator it(mat,k); it; ++it) { // Loop on rows
100  if (std::abs(it.value()) / norms[it.col()] > threshold) {
101  if (m_map.find(it.row()) == m_map.end()) {
102  m_map[it.row()] = rows;
103  rows += 1;
104  }
105  }
106  }
107  }
108 
109  // --- Copying
110  assert(rows == (Index)m_map.size());
111  m_matrix.resize(rows, mat.cols());
112  m_matrix.setZero();
113 
114  if (rows > 0)
115  for (int k=0; k<mat.outerSize(); ++k) { // Loop on cols
116  for (typename Eigen::SparseMatrix<Scalar, Eigen::ColMajor>::InnerIterator it(mat,k); it; ++it) { // Loop on rows
117  if (std::abs(it.value()) / norms[it.col()] > threshold)
118  m_matrix(m_map[it.row()], it.col()) = it.value();
119  }
120  }
121 
122  }
123 
124 
126  SparseDense(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &mat, double threshold = epsilon()) : m_dimension(mat.rows()) {
127  // --- Compute which rows we need
128 
129 
130  // Computing the column norms
131  RealVector norms(mat.cols());
132  norms.setZero();
133  for (int k=0; k<mat.outerSize(); ++k) // Loop on rows
134  for (typename Eigen::SparseMatrix<Scalar, Eigen::RowMajor>::InnerIterator it(mat,k); it; ++it) // Loop on cols
135  norms[it.col()] += Eigen::internal::abs2(it.value());
136 
137  // Computing selected rows
138  norms = norms.cwiseSqrt();
139  Index rows = 0;
140  for (int k=0; k<mat.outerSize(); ++k) {
141  bool any = false;
142  for (typename Eigen::SparseMatrix<Scalar, Eigen::RowMajor>::InnerIterator it(mat,k); it && !any; ++it) {
143  any |= std::abs(it.value()) / norms[it.col()] > threshold;
144  }
145  if (any) {
146  m_map[k] = rows;
147  rows += 1;
148  }
149  }
150 
151  // --- Copying
152  m_matrix.resize(rows, mat.cols());
153  m_matrix.setZero();
154 
155  for(auto i = m_map.begin(); i != m_map.end(); i++) {
156  for (typename Eigen::SparseMatrix<Scalar, Eigen::RowMajor>::InnerIterator it(mat,i->first); it; ++it)
157  m_matrix(i->second, it.col()) = it.value();
158 
159  }
160 
161  }
162 
171  SparseDense(const ScalarMatrix &mat, double threshold = Eigen::NumTraits<Scalar>::epsilon()) : m_dimension(mat.rows()) {
172  // Compute which rows we need
173  Matrix<Real, 1, Dynamic> norms = mat.cwiseAbs2().colwise().sum();
174 
175  size_t rows = 0;
176  for(Index i = 0; i < mat.rows(); i++) {
177  if ((std::abs(mat.row(i).array()) / norms.array() > threshold).any()) {
178  m_map[i] = rows;
179  rows += 1;
180  }
181  }
182 
183  // Copy the selected rows
184  m_matrix.resize(rows, mat.cols());
185  if ((size_t)mat.rows() == m_map.size())
186  m_matrix = mat;
187  else
188  for(auto i = m_map.begin(); i != m_map.end(); i++) {
189  m_matrix.row(i->second) = mat.row(i->first);
190  }
191 
192  }
193 
194 
196  ScalarMatrix toDense() const {
197  ScalarMatrix mat(m_dimension, size());
198  mat.setZero();
199 
200  for(auto i = m_map.begin(); i != m_map.end(); i++)
201  mat.row(i->first) = m_matrix.row(i->second);
202 
203  return mat;
204  }
205 
207  Index denseDimension() const {
208  return m_matrix.rows();
209  }
210 
218  void cleanup(Real threshold = epsilon()) {
219  Matrix<Real, 1, Dynamic> norms = m_matrix.cwiseAbs2().colwise().sum();
220 
221  std::vector<bool> selected(m_matrix.rows(), true);
222  for(auto i = m_map.begin(); i != m_map.end(); ) {
223  auto j = i;
224  i++;
225  if ((std::abs(m_matrix.row(j->second).array()) / norms.array() <= threshold).all()) {
226  selected[j->second] = false;
227  m_map.erase(j);
228  }
229  }
230 
231  Index newSize = 0;
232  std::vector<Index> ix(m_matrix.rows(), 0);
233  for(Index i = 0; i < m_matrix.rows(); i++) {
234  ix[i] = newSize;
235  if (selected[i]) newSize++;
236  }
237  for(auto i = m_map.begin(); i != m_map.end(); i++)
238  i->second = ix[i->second];
239 
240 
241  select_rows(selected, m_matrix, m_matrix);
242  assert(m_matrix.rows() == newSize);
243  }
244 
245  // --- Base methods
246  virtual Index size() const {
247  return m_matrix.cols();
248  }
249 
250  Index dimension() const {
251  return m_dimension;
252  }
253 
254 
255 #ifndef SWIG
256  struct Insert {
257  RowMap &m_map;
258  void operator()(const RowMap::const_reference &x) const {
259  size_t s = m_map.size();
260  m_map[x.first] = s;
261  }
262  };
263 
264  struct KeyComparator {
265  bool operator()(const RowMap::const_reference &a, const RowMap::const_reference b) {
266  return a.first < b.first;
267  }
268  };
269 #endif
270 
271  void add(const FMatrixBase &_other, const std::vector<bool> *which = NULL) override {
272  const Self &other = kqp::our_dynamic_cast<const Self&>(_other);
273 
274  if (m_dimension != other.m_dimension)
275  KQP_THROW_EXCEPTION_F(illegal_argument_exception, "Cannot add vectors of different sizes (%d vs %d)", %m_dimension %other.m_dimension);
276 
277  // Compute which vectors to add
278  std::vector<Index> ix;
279  Index toAdd = 0;
280  if (which) {
281  for(size_t i = 0; i < which->size(); i++)
282  if ((*which)[i]) ix.push_back(i);
283  toAdd = ix.size();
284  } else toAdd = other.size();
285 
286  if (toAdd == 0) return;
287 
288  Index oldRows = m_map.size();
289  std::set_difference(other.m_map.begin(), other.m_map.end(), m_map.begin(), m_map.end(),
290  boost::make_function_output_iterator(Insert({m_map})), KeyComparator());
291 
292  Index newRows = m_map.size() - oldRows;
293 
294  // Add
295  Index offset = m_matrix.cols();
296  m_matrix.conservativeResize(m_map.size(), m_matrix.cols() + toAdd);
297  m_matrix.topRightCorner(oldRows, toAdd).setZero();
298  m_matrix.bottomRows(newRows).setZero();
299 
300  for(auto i = other.m_map.begin(); i != other.m_map.end(); i++) {
301  Index otherRow = i->second;
302  Index selfRow = m_map[i->first];
303 
304  if (which) {
305  for(size_t k = 0; k < ix.size(); k++)
306  m_matrix(selfRow, offset + k) = other.m_matrix(otherRow, ix[k]);
307  } else
308  m_matrix.row(selfRow).tail(toAdd) = other.m_matrix.row(otherRow);
309  }
310 
311  }
312 
313 
314  const ScalarMatrix &gramMatrix() const {
315  if (size() == m_gramMatrix.rows()) return m_gramMatrix;
316 
317  // We lose space here, could be used otherwise???
318  Index current = m_gramMatrix.rows();
319  if (current < size())
320  m_gramMatrix.conservativeResize(size(), size());
321 
322  Index tofill = size() - current;
323 
324  // Compute the remaining inner products
325  m_gramMatrix.bottomRightCorner(tofill, tofill).noalias() = m_matrix.rightCols(tofill).adjoint() * m_matrix.rightCols(tofill);
326  m_gramMatrix.topRightCorner(current, tofill).noalias() = m_matrix.leftCols(current).adjoint() * m_matrix.rightCols(tofill);
327  m_gramMatrix.bottomLeftCorner(tofill, current) = m_gramMatrix.topRightCorner(current, tofill).adjoint().eval();
328 
329  return m_gramMatrix;
330  }
331 
332 
334  template<class DerivedMatrix>
335  void inner(const Self &other, DerivedMatrix &result) const {
336  if (result.rows() != this->size() || result.cols() != other.size())
337  result.resize(this->size(), other.size());
338  result.setZero();
339 
340  struct Collector {
341  const ScalarMatrix &mat1, &mat2;
342  const RowMap &map1, &map2;
343  DerivedMatrix &mat;
344 
345  void operator()(const RowMap::const_reference & x) const {
346  Index ix1 = this->map1.find(x.first)->second;
347  Index ix2 = this->map2.find(x.first)->second;
348  this->mat += this->mat1.row(ix1).adjoint() * this->mat2.row(ix2);
349  }
350  } collector({m_matrix, other.m_matrix, m_map, other.m_map, result});
351 
352  std::set_intersection(m_map.begin(), m_map.end(), other.m_map.begin(), other.m_map.end(),
353  boost::make_function_output_iterator(collector), KeyComparator());
354 
355  }
356 
357 
358  // Computes alpha * X * A + beta * Y * B (X = *this)
359  FMatrixBasePtr linearCombination(const ScalarAltMatrix &mA, Scalar alpha, const Self *mY, const ScalarAltMatrix *mB, Scalar beta) const {
360  // Simple case: we don't have to add another matrix
361  if (!mY)
362  return FMatrixBasePtr(new Self(m_dimension, RowMap(m_map), alpha * m_matrix * mA));
363 
364  // Add the keys
365  RowMap newMap;
366  std::set_union(mY->m_map.begin(), mY->m_map.end(), m_map.begin(), m_map.end(),
367  boost::make_function_output_iterator(Insert({newMap})),
368  KeyComparator());
369 
370 
371  // Perform the linear combination
372  ScalarMatrix mat(newMap.size(), mA.cols());
373  mat.setZero();
374 
375  for(auto i = m_map.begin(), end = m_map.end(); i != end; i++)
376  mat.row(newMap[i->first]) = alpha * m_matrix.row(i->second) * mA;
377 
378  for(auto i = mY->m_map.begin(), end = mY->m_map.end(); i != end; i++)
379  mat.row(newMap[i->first]) += beta * mY->m_matrix.row(i->second) * *mB;
380 
381  // Move and cleanup before returning
382  FMatrixBasePtr sdMat(new Self(m_dimension, std::move(newMap), std::move(mat)));
383  kqp::our_dynamic_cast<Self&>(*sdMat).cleanup(epsilon());
384  return sdMat;
385  }
386 
387 
388 
389  FMatrixBasePtr subset(const std::vector<bool>::const_iterator &begin, const std::vector<bool>::const_iterator &end) const override {
390  boost::shared_ptr<Self> dest(new Self());
391  select_columns(begin, end, m_matrix, dest->m_matrix);
392 
393  dest->m_dimension = m_dimension;
394  dest->m_map = m_map;
395 
396  return dest;
397  }
398 
399  virtual FMatrixBasePtr copy() const override {
400  return FMatrixBasePtr(new Self(*this));
401  }
402 
403  Self &operator=(const Self &other) {
404  m_dimension = other.m_dimension;
405  m_matrix = other.m_matrix;
406  m_map = other.m_map;
407  m_gramMatrix = other.m_gramMatrix;
408  return *this;
409  }
410 
411 
412  virtual FMatrixBase &operator=(const FMatrixBase &_other) override {
413  return *this = kqp::our_dynamic_cast<const Self&>(_other);
414  }
415 
417  const std::map<Index, Index> &map() const {
418  return m_map;
419  }
420 
422  const ScalarMatrix &getMatrix() const {
423  return m_matrix;
424  }
425 
426  template<typename Archive>
427  inline void serialize(Archive & ar, const unsigned int /*file_version*/) {
428  ar & m_dimension;
429  ar & m_map;
430  ar & m_matrix;
431  }
432  private:
434  Index m_dimension;
435 
437  std::map<Index, Index> m_map;
438 
440  ScalarMatrix m_matrix;
441 
443  mutable ScalarMatrix m_gramMatrix;
444 
445  };
446 
447 
448  template<typename Scalar>
449  class SparseDenseSpace : public SpaceBase<Scalar> {
450  public:
452  KQP_SPACE_TYPEDEFS("sparse-dense", Scalar);
453 #ifndef SWIG
454  using SpaceBase<Scalar>::k;
455 #endif
456  static FSpacePtr create(Index dimension) { return FSpacePtr(new SparseDenseSpace(dimension)); }
457 
458  SparseDenseSpace(Index dimension) : m_dimension(dimension) {}
459  SparseDenseSpace() : m_dimension(0) {}
460 
461  inline static const SparseDense<Scalar>& cast(const FeatureMatrixBase<Scalar> &mX) { return kqp::our_dynamic_cast<const SparseDense<Scalar> &>(mX); }
462 
463 
464  Index dimension() const override { return m_dimension; }
465  void dimension(Index dimension) { m_dimension = dimension; }
466 
467  virtual FSpacePtr copy() const override { return FSpacePtr(new SparseDenseSpace(m_dimension)); }
468 
469  virtual FMatrixBasePtr newMatrix() const override {
470  return FMatrixBasePtr(new SparseDense<Scalar>(m_dimension));
471  }
472  virtual FMatrixBasePtr newMatrix(const FMatrixBase &mX) const override {
473  return FMatrixBasePtr(new SparseDense<Scalar>(cast(mX)));
474  }
475 
476  virtual bool _canLinearlyCombine() const override {
477  return true;
478  }
479 
480  const ScalarMatrix &k(const FeatureMatrixBase<Scalar> &mX) const override {
481  return cast(mX).gramMatrix();
482  }
483 
484  virtual ScalarMatrix k(const FeatureMatrixBase<Scalar> &mX1, const ScalarAltMatrix &mY1, const RealAltVector &mD1,
485  const FeatureMatrixBase<Scalar> &mX2, const ScalarAltMatrix &mY2, const RealAltVector &mD2) const override {
486  ScalarMatrix inner;
487  cast(mX1).inner(cast(mX2), inner);
488  return mD1.asDiagonal() * mY1.adjoint() * inner * mY2 * mD2.asDiagonal();
489  };
490 
491  virtual FMatrixBasePtr linearCombination(const FeatureMatrixBase<Scalar> &mX, const ScalarAltMatrix &mA, Scalar alpha,
492  const FeatureMatrixBase<Scalar> *mY, const ScalarAltMatrix *mB, Scalar beta) const override {
493  return cast(mX).linearCombination(mA, alpha, kqp::our_dynamic_cast<const SparseDense<Scalar> *>(mY), mB, beta);
494  }
495 
496  virtual void load(const picojson::object &json) override {
497  m_dimension = kqp::getNumeric<Index>("", json, "dimension");
498  }
499 
500  virtual picojson::object save() const override {
501  picojson::object json = SpaceBase<Scalar>::save();
502  json["dimension"] = picojson::value((double)m_dimension);
503  return json;
504  }
505 
506  private:
507  Index m_dimension;
508 
509  };
510 
511 } // end namespace kqp
512 
513 #endif