Quantum++  v0.8.6
C++11 quantum computing library
util.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2016 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 all sizes match
119 template<typename T1, typename T2>
120 bool _check_matching_sizes(const T1& lhs, const T2& rhs) noexcept
121 {
122  return lhs.size() == rhs.size();
123 }
124 
125 // check that dims is a valid dimension vector
126 inline bool _check_dims(const std::vector<idx>& dims)
127 {
128  if (dims.size() == 0)
129  return false;
130 
131  return std::find_if(std::begin(dims), std::end(dims),
132  [dims](idx i) -> bool
133  {
134  if (i == 0) return true;
135  else return false;
136  }) == std::end(dims);
137 }
138 
139 // check that valid dims match the dimensions
140 // of valid (non-zero sized) square matrix
141 template<typename Derived>
142 bool _check_dims_match_mat(const std::vector<idx>& dims,
143  const Eigen::MatrixBase<Derived>& A)
144 {
145  const dyn_mat<typename Derived::Scalar>& rA = A;
146 
147  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
148  static_cast<idx>(1), std::multiplies<idx>());
149 
150  return proddim == static_cast<idx>(rA.rows());
151 }
152 
153 // check that valid dims match the dimensions of valid column vector
154 template<typename Derived>
155 bool _check_dims_match_cvect(const std::vector<idx>& dims,
156  const Eigen::MatrixBase<Derived>& V)
157 {
158  const dyn_mat<typename Derived::Scalar>& rV = V;
159 
160  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
161  static_cast<idx>(1), std::multiplies<idx>());
162 
163  return proddim == static_cast<idx>(rV.rows());
164 }
165 
166 // check that valid dims match the dimensions of valid row vector
167 template<typename Derived>
168 bool _check_dims_match_rvect(const std::vector<idx>& dims,
169  const Eigen::MatrixBase<Derived>& V)
170 {
171  const dyn_mat<typename Derived::Scalar>& rV = V;
172 
173  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
174  static_cast<idx>(1), std::multiplies<idx>());;
175 
176  return proddim == static_cast<idx>(rV.cols());
177 }
178 
179 // check that all elements in valid dims equal to dim
180 inline bool _check_eq_dims(const std::vector<idx>& dims, idx dim) noexcept
181 {
182  for (idx i : dims)
183  if (i != dim)
184  return false;
185 
186  return true;
187 }
188 
189 // check that subsys is valid with respect to valid dims
190 inline bool _check_subsys_match_dims(const std::vector<idx>& subsys,
191  const std::vector<idx>& dims)
192 {
193  // // check non-zero sized subsystems
194  // if (subsys.size() == 0)
195  // return false;
196 
197  // check valid number of subsystems
198  if (subsys.size() > dims.size())
199  return false;
200 
201  // sort the subsystems
202  std::vector<idx> subsyssort = subsys;
203  std::sort(std::begin(subsyssort), std::end(subsyssort));
204 
205  // check duplicates
206  if (std::unique(std::begin(subsyssort), std::end(subsyssort))
207  != std::end(subsyssort))
208  return false;
209 
210  // check range of subsystems
211  return std::find_if(std::begin(subsyssort), std::end(subsyssort),
212  [dims](idx i) -> bool
213  {
214  return i > dims.size() - 1;
215  }) == std::end(subsyssort);
216 }
217 
218 // check matrix is 2 x 2
219 template<typename Derived>
220 bool _check_qubit_matrix(const Eigen::MatrixBase<Derived>& A) noexcept
221 {
222  const dyn_mat<typename Derived::Scalar>& rA = A;
223 
224  return rA.rows() == 2 && rA.cols() == 2;
225 }
226 
227 // check column vector is 2 x 1
228 template<typename Derived>
229 bool _check_qubit_cvector(const Eigen::MatrixBase<Derived>& V) noexcept
230 {
231  const dyn_mat<typename Derived::Scalar>& rV = V;
232 
233  return rV.rows() == 2 && rV.cols() == 1;
234 }
235 
236 // check row vector is 1 x 2
237 template<typename Derived>
238 bool _check_qubit_rvector(const Eigen::MatrixBase<Derived>& V) noexcept
239 {
240  const dyn_mat<typename Derived::Scalar>& rV = V;
241 
242  return rV.rows() == 1 && rV.cols() == 2;
243 }
244 
245 // check row vector is 1 x 2 or 2 x 1
246 template<typename Derived>
247 bool _check_qubit_vector(const Eigen::MatrixBase<Derived>& V) noexcept
248 {
249  const dyn_mat<typename Derived::Scalar>& rV = V;
250 
251  return (rV.rows() == 1 && rV.cols() == 2) ||
252  (rV.rows() == 2 && rV.cols() == 1);
253 }
254 
255 
256 // check valid permutation
257 inline bool _check_perm(const std::vector<idx>& perm)
258 {
259  if (perm.size() == 0)
260  return false;
261 
262  std::vector<idx> ordered(perm.size());
263  std::iota(std::begin(ordered), std::end(ordered), 0);
264 
265  return std::is_permutation(std::begin(ordered), std::end(ordered),
266  std::begin(perm));
267 }
268 
269 // Kronecker product of 2 matrices, preserve return type
270 // internal function for the variadic template function wrapper kron()
271 template<typename Derived1, typename Derived2>
272 dyn_mat<typename Derived1::Scalar> _kron2(const Eigen::MatrixBase<Derived1>& A,
273  const Eigen::MatrixBase<Derived2>& B)
274 {
277 
278  // EXCEPTION CHECKS
279 
280  // check types
281  if (!std::is_same<typename Derived1::Scalar,
282  typename Derived2::Scalar>::value)
283  throw Exception("qpp::kron()", Exception::Type::TYPE_MISMATCH);
284 
285  // check zero-size
287  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
288 
289  // check zero-size
291  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
292  // END EXCEPTION CHECKS
293 
294  idx Acols = static_cast<idx>(rA.cols());
295  idx Arows = static_cast<idx>(rA.rows());
296  idx Bcols = static_cast<idx>(rB.cols());
297  idx Brows = static_cast<idx>(rB.rows());
298 
300  result.resize(Arows * Brows, Acols * Bcols);
301 
302 #pragma omp parallel for collapse(2)
303  for (idx j = 0; j < Acols; ++j) // column major order for speed
304  for (idx i = 0; i < Arows; ++i)
305  result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
306 
307  return result;
308 }
309 
310 // Direct sum of 2 matrices, preserve return type
311 // internal function for the variadic template function wrapper dirsum()
312 template<typename Derived1, typename Derived2>
314  const Eigen::MatrixBase<Derived1>& A,
315  const Eigen::MatrixBase<Derived2>& B)
316 {
319 
320  // EXCEPTION CHECKS
321 
322  // check types
323  if (!std::is_same<typename Derived1::Scalar,
324  typename Derived2::Scalar>::value)
325  throw Exception("qpp::dirsum()", Exception::Type::TYPE_MISMATCH);
326 
327  // check zero-size
329  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
330 
331  // check zero-size
333  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
334  // END EXCEPTION CHECKS
335 
336  idx Acols = static_cast<idx>(rA.cols());
337  idx Arows = static_cast<idx>(rA.rows());
338  idx Bcols = static_cast<idx>(rB.cols());
339  idx Brows = static_cast<idx>(rB.rows());
340 
343  Acols + Bcols);
344 
345  result.block(0, 0, Arows, Acols) = rA;
346  result.block(Arows, Acols, Brows, Bcols) = rB;
347 
348  return result;
349 }
350 
351 // may be useful, extracts variadic template argument pack into a std::vector
352 template<typename T>
353 // ends the recursion
354 void variadic_vector_emplace(std::vector<T>&)
355 {
356 }
357 
358 // may be useful, extracts variadic template argument pack into a std::vector
359 template<typename T, typename First, typename ... Args>
360 void variadic_vector_emplace(std::vector<T>& v, First&& first, Args&& ... args)
361 {
362  v.emplace_back(std::forward<First>(first));
363  variadic_vector_emplace(v, std::forward<Args>(args)...);
364 }
365 
366 } /* namespace internal */
367 } /* namespace qpp */
368 
369 #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:229
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:80
bool _check_matching_sizes(const T1 &lhs, const T2 &rhs) noexcept
Definition: util.h:120
bool _check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:190
bool _check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:142
bool _check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:155
dyn_mat< typename Derived1::Scalar > _kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:272
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:83
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:168
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:220
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:257
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:354
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:313
bool _check_qubit_vector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:247
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:126
bool _check_eq_dims(const std::vector< idx > &dims, idx dim) noexcept
Definition: util.h:180
bool _check_qubit_rvector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:238