Quantum++  v0.8
C++11 quantum computing library
util.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2015 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef INTERNAL_UTIL_H_
28 #define INTERNAL_UTIL_H_
29 
30 namespace qpp
31 {
36 namespace internal
37 {
38 
39 // integer index to multi-index, use C-style array for speed
40 // standard lexicographical order, e.g. 00, 01, 10, 11
41 inline void _n2multiidx(idx n, idx numdims, const idx* dims, idx* result)
42 noexcept
43 {
44  // no error checks to improve speed
45  for (idx i = 0; i < numdims; ++i)
46  {
47  result[numdims - i - 1] = n % (dims[numdims - i - 1]);
48  n /= (dims[numdims - i - 1]);
49  }
50 }
51 
52 // multi-index to integer index, use C-style array for speed,
53 // standard lexicographical order, e.g. 00->0, 01->1, 10->2, 11->3
54 inline idx _multiidx2n(const idx* midx, idx numdims, const idx* dims)
55 noexcept
56 {
57  // no error checks to improve speed
58 
59  // Static allocation for speed!
60  // double the size for matrices reshaped as vectors
61  idx part_prod[2 * maxn];
62 
63  idx result = 0;
64  part_prod[numdims - 1] = 1;
65  for (idx i = 1; i < numdims; ++i)
66  {
67  part_prod[numdims - i - 1] =
68  part_prod[numdims - i] * dims[numdims - i];
69  result += midx[numdims - i - 1] * part_prod[numdims - i - 1];
70  }
71 
72  return result + midx[numdims - 1];
73 }
74 
75 // check square matrix
76 template<typename Derived>
77 bool _check_square_mat(const Eigen::MatrixBase<Derived>& A)
78 {
80 
81  return rA.rows() == rA.cols();
82 }
83 
84 // check whether input is a vector or not
85 template<typename Derived>
86 bool _check_vector(const Eigen::MatrixBase<Derived>& A)
87 {
89 
90  return rA.rows() == 1 || rA.cols() == 1;
91 }
92 
93 // check whether input is a row vector or not
94 template<typename Derived>
95 bool _check_rvector(const Eigen::MatrixBase<Derived>& A)
96 {
98 
99  return rA.rows() == 1;
100 }
101 
102 // check whether input is a column vector or not
103 template<typename Derived>
104 bool _check_cvector(const Eigen::MatrixBase<Derived>& A)
105 {
106  const dyn_mat<typename Derived::Scalar>& rA = A;
107 
108  return rA.cols() == 1;
109 }
110 
111 // check non-zero size of object that supports size() function
112 template<typename T>
113 bool _check_nonzero_size(const T& x) noexcept
114 {
115  return x.size() != 0;
116 }
117 
118 // check that dims is a valid dimension vector
119 inline bool _check_dims(const std::vector <idx>& dims)
120 {
121  if (dims.size() == 0)
122  return false;
123 
124  return std::find_if(std::begin(dims), std::end(dims),
125  [dims](idx i) -> bool
126  {
127  if (i == 0) return true;
128  else return false;
129  }) == std::end(dims);
130 }
131 
132 // check that valid dims match the dimensions
133 // of valid (non-zero sized) quare matrix
134 template<typename Derived>
135 bool _check_dims_match_mat(const std::vector <idx>& dims,
136  const Eigen::MatrixBase<Derived>& A)
137 {
138  const dyn_mat<typename Derived::Scalar>& rA = A;
139 
140  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
141  static_cast<idx>(1), std::multiplies<idx>());
142 
143  return proddim == static_cast<idx>(rA.rows());
144 }
145 
146 // check that valid dims match the dimensions of valid column vector
147 template<typename Derived>
148 bool _check_dims_match_cvect(const std::vector <idx>& dims,
149  const Eigen::MatrixBase<Derived>& V)
150 {
151  const dyn_mat<typename Derived::Scalar>& rV = V;
152 
153  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
154  static_cast<idx>(1), std::multiplies<idx>());
155 
156  return proddim == static_cast<idx>(rV.rows());
157 }
158 
159 // check that valid dims match the dimensions of valid row vector
160 template<typename Derived>
161 bool _check_dims_match_rvect(const std::vector <idx>& dims,
162  const Eigen::MatrixBase<Derived>& V)
163 {
164  const dyn_mat<typename Derived::Scalar>& rV = V;
165 
166  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
167  static_cast<idx>(1), std::multiplies<idx>());;
168 
169  return proddim == static_cast<idx>(rV.cols());
170 }
171 
172 // check that all elements in valid dims equal to dim
173 inline bool _check_eq_dims(const std::vector <idx>& dims, idx dim) noexcept
174 {
175  for (idx i : dims)
176  if (i != dim)
177  return false;
178 
179  return true;
180 }
181 
182 // check that subsys is valid with respect to valid dims
183 inline bool _check_subsys_match_dims(const std::vector <idx>& subsys,
184  const std::vector <idx>& dims)
185 {
186  // // check non-zero sized subsystems
187  // if (subsys.size() == 0)
188  // return false;
189 
190  // check valid number of subsystems
191  if (subsys.size() > dims.size())
192  return false;
193 
194  // sort the subsystems
195  std::vector <idx> subsyssort = subsys;
196  std::sort(std::begin(subsyssort), std::end(subsyssort));
197 
198  // check duplicates
199  if (std::unique(std::begin(subsyssort), std::end(subsyssort))
200  != std::end(subsyssort))
201  return false;
202 
203  // check range of subsystems
204  return std::find_if(std::begin(subsyssort), std::end(subsyssort),
205  [dims](idx i) -> bool
206  {
207  return i > dims.size() - 1;
208  }) == std::end(subsyssort);
209 }
210 
211 // check matrix is 2 x 2
212 template<typename Derived>
213 bool _check_qubit_matrix(const Eigen::MatrixBase<Derived>& A) noexcept
214 {
215  const dyn_mat<typename Derived::Scalar>& rA = A;
216 
217  return rA.rows() == 2 && rA.cols() == 2;
218 }
219 
220 // check column vector is 2 x 1
221 template<typename Derived>
222 bool _check_qubit_cvector(const Eigen::MatrixBase<Derived>& V) noexcept
223 {
224  const dyn_mat<typename Derived::Scalar>& rV = V;
225 
226  return rV.rows() == 2 && rV.cols() == 1;
227 }
228 
229 // check row vector is 1 x 2
230 template<typename Derived>
231 bool _check_qubit_rvector(const Eigen::MatrixBase<Derived>& V) noexcept
232 {
233  const dyn_mat<typename Derived::Scalar>& rV = V;
234 
235  return rV.rows() == 1 && rV.cols() == 2;
236 }
237 
238 // check row vector is 1 x 2 or 2 x 1
239 template<typename Derived>
240 bool _check_qubit_vector(const Eigen::MatrixBase<Derived>& V) noexcept
241 {
242  const dyn_mat<typename Derived::Scalar>& rV = V;
243 
244  return (rV.rows() == 1 && rV.cols() == 2) ||
245  (rV.rows() == 2 && rV.cols() == 1);
246 }
247 
248 
249 // check valid permutation
250 inline bool _check_perm(const std::vector <idx>& perm)
251 {
252  if (perm.size() == 0)
253  return false;
254 
255  std::vector <idx> ordered(perm.size());
256  std::iota(std::begin(ordered), std::end(ordered), 0);
257 
258  return std::is_permutation(std::begin(ordered), std::end(ordered),
259  std::begin(perm));
260 }
261 
262 // Kronecker product of 2 matrices, preserve return type
263 // internal function for the variadic template function wrapper kron()
264 template<typename Derived1, typename Derived2>
265 dyn_mat<typename Derived1::Scalar> _kron2(const Eigen::MatrixBase<Derived1>& A,
266  const Eigen::MatrixBase<Derived2>& B)
267 {
270 
271  // EXCEPTION CHECKS
272 
273  // check types
274  if (!std::is_same<typename Derived1::Scalar,
275  typename Derived2::Scalar>::value)
276  throw Exception("qpp::kron()", Exception::Type::TYPE_MISMATCH);
277 
278  // check zero-size
280  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
281 
282  // check zero-size
284  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
285 
286  idx Acols = static_cast<idx>(rA.cols());
287  idx Arows = static_cast<idx>(rA.rows());
288  idx Bcols = static_cast<idx>(rB.cols());
289  idx Brows = static_cast<idx>(rB.rows());
290 
292  result.resize(Arows * Brows, Acols * Bcols);
293 
294 #pragma omp parallel for collapse(2)
295  for (idx j = 0; j < Acols; ++j) // column major order for speed
296  for (idx i = 0; i < Arows; ++i)
297  result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
298 
299  return result;
300 }
301 
302 // Direct sum of 2 matrices, preserve return type
303 // internal function for the variadic template function wrapper dirsum()
304 template<typename Derived1, typename Derived2>
306  const Eigen::MatrixBase<Derived1>& A,
307  const Eigen::MatrixBase<Derived2>& B)
308 {
311 
312  // EXCEPTION CHECKS
313 
314  // check types
315  if (!std::is_same<typename Derived1::Scalar,
316  typename Derived2::Scalar>::value)
317  throw Exception("qpp::dirsum()", Exception::Type::TYPE_MISMATCH);
318 
319  // check zero-size
321  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
322 
323  // check zero-size
325  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
326 
327  idx Acols = static_cast<idx>(rA.cols());
328  idx Arows = static_cast<idx>(rA.rows());
329  idx Bcols = static_cast<idx>(rB.cols());
330  idx Brows = static_cast<idx>(rB.rows());
331 
333  dyn_mat < typename Derived1::Scalar > ::Zero(Arows + Brows,
334  Acols + Bcols);
335 
336  result.block(0, 0, Arows, Acols) = rA;
337  result.block(Arows, Acols, Brows, Bcols) = rB;
338 
339  return result;
340 }
341 
342 // may be useful, extracts variadic template argument pack into a std::vector
343 template<typename T>
344 // ends the recursion
345 void variadic_vector_emplace(std::vector <T>&)
346 {
347 }
348 
349 // may be useful, extracts variadic template argument pack into a std::vector
350 template<typename T, typename First, typename ... Args>
351 void variadic_vector_emplace(std::vector <T>& v, First&& first, Args&& ... args)
352 {
353  v.emplace_back(std::forward<First>(first));
354  variadic_vector_emplace(v, std::forward<Args>(args)...);
355 }
356 
357 } /* namespace internal */
358 } /* namespace qpp */
359 
360 #endif /* INTERNAL_UTIL_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:104
bool _check_qubit_cvector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:222
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:82
bool _check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:183
bool _check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:135
bool _check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:148
dyn_mat< typename Derived1::Scalar > _kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:265
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:84
bool _check_vector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:86
Quantum++ main namespace.
Definition: codes.h:30
bool _check_dims_match_rvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:161
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:77
bool _check_rvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:95
bool _check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:213
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:250
idx _multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:54
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
void variadic_vector_emplace(std::vector< T > &)
Definition: util.h:345
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:113
dyn_mat< typename Derived1::Scalar > _dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:305
bool _check_qubit_vector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:240
void _n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:41
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:119
bool _check_eq_dims(const std::vector< idx > &dims, idx dim) noexcept
Definition: util.h:173
bool _check_qubit_rvector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:231