Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
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_DENSE_FEATURE_MATRIX_H__
18 #define __KQP_DENSE_FEATURE_MATRIX_H__
19 
20 #include <numeric>
21 #include <boost/lexical_cast.hpp>
22 
23 #include <kqp/feature_matrix.hpp>
24 #include <kqp/subset.hpp>
25 #include <kqp/intervals.hpp>
26 
27 namespace kqp {
28 
29  template <typename Scalar> class Dense;
30  template <typename Scalar> class DenseSpace;
31 
32 # include <kqp/define_header_logger.hpp>
33  DEFINE_KQP_HLOGGER("kqp.feature-matrix.dense");
34 
35 
40  template <typename Scalar>
41  class Dense : public FeatureMatrixBase<Scalar> {
42  public:
43  typedef Dense<Scalar> Self;
44  KQP_SPACE_TYPEDEFS("dense", Scalar);
45 
46  virtual ~Dense() {}
47 
49  Dense() {}
50 
52  Dense(Index dimension) : m_matrix(dimension, 0) {
53  }
54 
56  Dense(const ScalarMatrix &m) : m_matrix(m) {}
57 
59  Dense(const Self &other) : m_gramMatrix(other.m_gramMatrix), m_matrix(other.m_matrix) {}
60 
61 
62 #ifndef SWIG
63  inline static SelfPtr create(Index dimension) { return SelfPtr(new Self(dimension)); }
64 
66  Dense(ScalarMatrix &&m) : m_matrix(m) {}
67 
69  static FMatrixPtr create(const ScalarMatrix &m) {
70  return FMatrixPtr(new Self(m));
71  }
72 
74  static FMatrixPtr create(const ScalarMatrix &&m) {
75  return FMatrixPtr(new Self(std::move(m)));
76  }
77 #endif
78 
79 
80 
84  template<typename Derived>
85  void add(const Eigen::DenseBase<Derived> &m, const std::vector<bool> *which = NULL) {
86  if (m_matrix.cols() == 0)
87  m_matrix.resize(m.rows(), 0);
88  if (m.rows() != m_matrix.rows())
89  KQP_THROW_EXCEPTION_F(illegal_operation_exception,
90  "Cannot add a vector of dimension %d (dimension is %d)", % m.rows() % m_matrix.rows());
91 
92 
93  Intervals intervals(which, m.cols());
94  Index offset = m_matrix.cols();
95  m_matrix.conservativeResize(m_matrix.rows(), intervals.selected() + m_matrix.cols());
96 
97  for(auto i = intervals.begin(); i != intervals.end(); i++) {
98  Index cols = i->second - i->first + 1;
99  KQP_HLOG_DEBUG_F("Copying cols %d to %d from a %d x %d matrix into (%d-%d) in a %d x %d matrix",
100  %i->first %i->second % m.rows() % m.cols()
101  %offset %(offset+cols-1) % m_matrix.rows() % m_matrix.cols());
102  this->m_matrix.block(0, offset, m_matrix.rows(), cols) = m.block(0, i->first, m_matrix.rows(), cols);
103  offset += cols;
104  }
105 
106  }
107 
109  const ScalarMatrix& getMatrix() const {
110  return this->m_matrix;
111  }
112 
113  const ScalarMatrix &toDense() const {
114  return this->m_matrix;
115  }
116 
117 #ifndef SWIG
118  inline static const Self &cast(const FMatrixBase &m) {
119  return kqp::our_dynamic_cast<const Self&>(m);
120  }
121 #endif
122 
123  void add(const FMatrixBase &other, const std::vector<bool> *which = NULL) override {
124  this->add(cast(other).getMatrix(), which);
125  }
126 
127 
128  virtual Index size() const override {
129  return m_matrix.cols();
130  }
131 
132  Index dimension() const {
133  return m_matrix.rows();
134  }
135 
137  const ScalarMatrix &gramMatrix() const {
138  if (m_matrix.cols() == m_gramMatrix.rows()) return m_gramMatrix;
139 
140  // We lose space here, could be used otherwise???
141  Index current = m_gramMatrix.rows();
142  if (current < m_matrix.cols())
143  m_gramMatrix.conservativeResize(m_matrix.cols(), m_matrix.cols());
144 
145  Index tofill = m_matrix.cols() - current;
146 
147  // Compute the remaining inner products
148  m_gramMatrix.bottomRightCorner(tofill, tofill).noalias() = this->getMatrix().rightCols(tofill).adjoint() * this->getMatrix().rightCols(tofill);
149  m_gramMatrix.topRightCorner(current, tofill).noalias() = this->getMatrix().leftCols(current).adjoint() * this->getMatrix().rightCols(tofill);
150  m_gramMatrix.bottomLeftCorner(tofill, current) = m_gramMatrix.topRightCorner(current, tofill).adjoint().eval();
151 
152  return m_gramMatrix;
153  }
154 
156  template<class DerivedMatrix>
157  void _inner(const Self &other, DerivedMatrix &result) const {
158  result = this->getMatrix().adjoint() * other.getMatrix();
159  }
160 
161 
162  FMatrixBasePtr linearCombination(const ScalarAltMatrix & mA, Scalar alpha, const Self *mY, const ScalarAltMatrix *mB, Scalar beta) const override {
163  ScalarMatrix m(alpha * (getMatrix() * mA));
164  if (mY != 0)
165  m += beta * mY->getMatrix() * *mB;
166 
167  return FMatrixBasePtr(new Self(std::move(m)));
168  }
169 
170 
171  FMatrixBasePtr subset(const std::vector<bool>::const_iterator &begin, const std::vector<bool>::const_iterator &end) const override {
172  ScalarMatrix m;
173  select_columns(begin, end, this->m_matrix, m);
174  return FMatrixBasePtr(new Self(std::move(m)));
175  }
176 
177  virtual FMatrixBasePtr copy() const override {
178  return FMatrixBasePtr(new Self(*this));
179  }
180 
181 
182 
183  Self& operator=(const Self &other) {
184  m_matrix = other.m_matrix;
185  m_gramMatrix = other.m_gramMatrix;
186  return *this;
187  }
188 
189  virtual FMatrixBase& operator=(const FMatrixBase &other) override {
190  return *this = kqp::our_dynamic_cast<const Self&>(other);
191  }
192 
193 
194  const ScalarMatrix * operator->() const {
195  return &m_matrix;
196  }
197  const ScalarMatrix & operator*() const {
198  return m_matrix;
199  }
200 
201  private:
202 
204  mutable ScalarMatrix m_gramMatrix;
205 
207  ScalarMatrix m_matrix;
208 
209  friend class DenseSpace<Scalar>;
210 
211  };
212 
213 
214 
215  template<typename Scalar>
216  std::ostream& operator<<(std::ostream &out, const Dense<Scalar> &f) {
217  return out << "[Dense Matrix with scalar " << KQP_DEMANGLE((Scalar)0) << "]" << std::endl << f.getMatrix();
218  }
219 
220 
224  template<typename Scalar>
225  class DenseSpace : public SpaceBase<Scalar> {
226  public:
227  typedef SpaceBase<Scalar> Self;
228  KQP_SPACE_TYPEDEFS("dense", Scalar);
229 #ifndef SWIG
230  using SpaceBase<Scalar>::k;
231 #endif
232  static FSpacePtr create(Index dimension) { return FSpacePtr(new DenseSpace(dimension)); }
233 
234  DenseSpace(Index dimension) : m_dimension(dimension) {}
235  DenseSpace() : m_dimension(0) {}
236 
237  FSpacePtr copy() const override {
238  return FSpacePtr(new DenseSpace(*this));
239  }
240 
241  inline static const Dense<Scalar>& cast(const FMatrixBase &mX) { return kqp::our_dynamic_cast<const Dense<Scalar> &>(mX); }
242 
243  FMatrixBasePtr newMatrix(const ScalarMatrix &mX) const {
244  return FMatrixBasePtr(new Dense<Scalar>(mX));
245  }
246 
247 
248  // Overriden methods
249 
250  Index dimension() const override { return m_dimension; }
251  void dimension(Index dimension) { m_dimension = dimension; }
252 
253 
254  virtual FMatrixBasePtr newMatrix() const override {
255  return FMatrixBasePtr(new Dense<Scalar>(m_dimension));
256  }
257  virtual FMatrixBasePtr newMatrix(const FMatrixBase &mX) const override {
258  return FMatrixBasePtr(new Dense<Scalar>(cast(mX)));
259  }
260 
261  virtual bool _canLinearlyCombine() const override {
262  return true;
263  }
264 
265  const ScalarMatrix &k(const FeatureMatrixBase<Scalar> &mX) const override {
266  return cast(mX).gramMatrix();
267  }
268 
269  virtual ScalarMatrix k(const FeatureMatrixBase<Scalar> &mX1, const ScalarAltMatrix &mY1, const RealAltVector &mD1,
270  const FeatureMatrixBase<Scalar> &mX2, const ScalarAltMatrix &mY2, const RealAltVector &mD2) const override {
271  return mD1.asDiagonal() * mY1.adjoint() * cast(mX1)->adjoint() * *cast(mX2) * mY2 * mD2.asDiagonal();
272  };
273 
274  virtual FMatrixBasePtr linearCombination(const FMatrixBase &mX, const ScalarAltMatrix &mA, Scalar alpha,
275  const FMatrixBase *mY, const ScalarAltMatrix *mB, Scalar beta) const override {
276  return FMatrixBasePtr(cast(mX).linearCombination(mA, alpha, kqp::our_dynamic_cast<const Dense<Scalar> *>(mY), mB, beta));
277  }
278 
279 
280  virtual void load(const picojson::object &json) override {
281  m_dimension = kqp::getNumeric<Index>("", json, "dimension");
282  }
283 
284  virtual picojson::object save() const override {
285  picojson::object json = SpaceBase<Scalar>::save();
286  json["dimension"] = picojson::value((double)m_dimension);
287  return json;
288  }
289 
290  private:
291  Index m_dimension;
292  };
293 
294 } // end namespace kqp
295 
296 #endif