Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
eigen_identity.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 
18 #ifndef _KQP_EIGEN_IDENTITY_H_
19 #define _KQP_EIGEN_IDENTITY_H_
20 
21 #include <kqp/kqp.hpp>
22 
23 namespace kqp {
24  template<typename Derived> struct DiagonalBlockWrapper;
25  template<typename Derived> class AltMatrixBase;
26  struct IdentityStorage {};
27 }
28 
29 namespace Eigen {
30 
31  template<typename Derived> class SparseMatrixBase;
32 
34  template<typename Scalar>
35  class Identity {
36  public:
37 #ifndef SWIG
38  typedef typename MatrixXd::Index Index;
39  typedef Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<Scalar>, Eigen::Matrix<Scalar,Dynamic,1> > VectorType;
40 
41  class RowWise {
42  public:
43  RowWise(Index rows, Index cols) : m_rows(rows), m_cols(cols) {}
44 
45  auto squaredNorm() const -> decltype(Eigen::Matrix<Scalar, Eigen::Dynamic, 1>::Ones(-1)) {
46  eigen_assert(m_rows == m_cols); // otherwise, it is a bit more complex! We need a constant expression
47  return Eigen::Matrix<Scalar, Eigen::Dynamic, 1>::Ones(m_rows);
48  }
49  private:
50  Index m_rows;
51  Index m_cols;
52  };
53 
54  RowWise rowwise() const { return RowWise(m_rows, m_cols); }
55 #endif
56 
57 
58  Identity() : m_rows(0), m_cols(0) {}
59  Identity(Index size) : m_rows(size), m_cols(size) {}
60 
61 
62  void swap(Identity &other) {
63  std::swap(m_rows, other.m_rows);
64  std::swap(m_cols, other.m_cols);
65  }
66 
67  Index rows() const { return this->m_rows; }
68  Index cols() const { return this->m_cols; }
69 
70 
71 
72  Scalar trace() const { return std::min(m_rows, m_cols); }
73  Scalar sum() const { return std::min(m_rows, m_cols); }
74  Scalar squaredNorm() const { return std::min(m_rows, m_cols); }
75 
76  void resize(Index rows, Index cols) {
77  m_rows = rows;
78  m_cols = cols;
79  }
80 
81  void conservativeResize(Index rows, Index cols) {
82  this->resize(rows, cols);
83  }
84 
85  Scalar operator()(Index i, Index j) const {
86  if (i == j) return 1;
87  return 0;
88  }
89 
90 #ifndef SWIG
91 
92  Identity(Index rows, Index cols) : m_rows(rows), m_cols(cols) { eigen_assert(rows == cols); }
93 
94  kqp::DiagonalBlockWrapper<VectorType> block(Index startCol, Index startRow, Index blockRows, Index blockCols) const {
95  return kqp::DiagonalBlockWrapper<VectorType>(VectorType(std::min(m_rows,m_cols),1,1), startCol, startRow, blockRows, blockCols);
96  }
97 
98  template<typename CwiseUnaryOp>
99  Identity unaryExpr(const CwiseUnaryOp &) {
101  //KQP_THROW_EXCEPTION(kqp::illegal_argument_exception, "Cannot apply a unary expression on the Identity matrix");
102  }
103 
104  auto getVectorIdentity() const -> decltype(Eigen::Matrix<Scalar,Dynamic,1>::Ones(0).asDiagonal()) {
105  return Eigen::Matrix<Scalar,Dynamic,1>::Ones(this->rows()).asDiagonal();
106  }
107 
108  auto getIdentityMatrix() const -> decltype(Eigen::Matrix<Scalar,Dynamic,Dynamic>::Identity(0,0)) {
109  return Eigen::Matrix<Scalar,Dynamic,Dynamic>::Identity(this->rows(),this->rows());
110  }
111 #endif
112 
113 
114  template<class Archive>
115  inline void serialize(Archive & ar, const unsigned int /*file_version*/) {
116  ar & m_rows;
117  ar & m_cols;
118  }
119 
120 
121  private:
122  Index m_rows;
123  Index m_cols;
124  };
125 
126 
127 #ifndef SWIG
128  template <typename Derived, typename Scalar> \
129  const DiagonalWrapper<Derived> operator* (const Identity<Scalar> &KQP_M_DEBUG(lhs), const DiagonalWrapper<Derived> &rhs) { \
130  eigen_assert(lhs.cols() == rhs.rows() \
131  && "invalid matrix product" \
132  && "if you wanted a coeff-wise or a dot product use the respective explicit functions"); \
133  return rhs.derived(); \
134  }
135 
136  template <typename Derived, typename Scalar> \
137  const DiagonalWrapper<Derived> operator* (const DiagonalWrapper<Derived> &lhs, const Identity<Scalar> &KQP_M_DEBUG(rhs)) { \
138  eigen_assert(lhs.cols() == rhs.rows() \
139  && "invalid matrix product" \
140  && "if you wanted a coeff-wise or a dot product use the respective explicit functions"); \
141  return lhs.derived();\
142  }
143 
144 
145 
147 # define KQP_IDENTITY_PRE_MULT(matrixtype) \
148  template <typename Derived, typename Scalar> \
149  const typename Eigen::internal::ref_selector<Derived>::type operator* (const Identity<Scalar> &KQP_M_DEBUG(lhs), const matrixtype &rhs) { \
150  eigen_assert(lhs.cols() == rhs.rows() \
151  && "invalid matrix product" \
152  && "if you wanted a coeff-wise or a dot product use the respective explicit functions"); \
153  return rhs.derived(); \
154  }
155 
157 # define KQP_IDENTITY_POST_MULT(matrixtype) \
158  template <typename Derived, typename Scalar> \
159  const typename Eigen::internal::ref_selector<Derived>::type operator* (const matrixtype &lhs, const Identity<Scalar> &KQP_M_DEBUG(rhs)) { \
160  eigen_assert(lhs.cols() == rhs.rows() \
161  && "invalid matrix product" \
162  && "if you wanted a coeff-wise or a dot product use the respective explicit functions"); \
163  return lhs.derived();\
164  }
165 
167 # define KQP_IDENTITY_MULT(type) KQP_IDENTITY_PRE_MULT(type) KQP_IDENTITY_POST_MULT(type)
168 
169 
170  KQP_IDENTITY_MULT(kqp::AltMatrixBase<Derived>)
171  KQP_IDENTITY_MULT(MatrixBase<Derived>)
172  KQP_IDENTITY_MULT(SparseMatrixBase<Derived>)
173 
174 
176 # define KQP_IDENTITY_ADD(op, matrixtype) \
177  template <typename Derived, typename Scalar> \
178  auto operator op (const Identity<Scalar> &lhs, const matrixtype &rhs) -> decltype(lhs.getIdentityMatrix() op rhs.derived()) { \
179  eigen_assert(lhs.cols() == rhs.cols()); \
180  eigen_assert(lhs.rows() == rhs.rows()); \
181  return lhs.getIdentityMatrix() op rhs.derived(); \
182  } \
183  template <typename Derived, typename Scalar> \
184  auto operator op (const matrixtype &lhs, const Identity<Scalar> &rhs) -> decltype(lhs.derived() op rhs.getIdentityMatrix()) { \
185  eigen_assert(lhs.cols() == rhs.cols()); \
186  eigen_assert(lhs.rows() == rhs.rows()); \
187  return lhs.derived() op rhs.getIdentityMatrix();\
188  }
189 
190 # define KQP_IDENTITY_OP(mtype) KQP_IDENTITY_ADD(+, mtype) KQP_IDENTITY_ADD(-, mtype)
191 
192  KQP_IDENTITY_OP(MatrixBase<Derived>);
193  KQP_IDENTITY_OP(SparseMatrixBase<Derived>);
194  KQP_IDENTITY_OP(kqp::AltMatrixBase<Derived>)
195 
196  namespace internal {
197  template<typename _Scalar>
198  struct traits<Identity<_Scalar>> {
200  typedef typename MatrixXd::Index Index;
201  typedef _Scalar Scalar;
202  };
203  }
204 #endif // SWIG
205 
206 }
207 
208 #endif