Quantum++  v0.6
C++11 quantum computing library
random.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 RANDOM_H_
28 #define RANDOM_H_
29 
30 // random matrices/states
31 
32 namespace qpp
33 {
34 
43 double rand(double a = 0, double b = 1)
44 {
45  std::uniform_real_distribution<> ud(a, b);
46 
48 }
49 
58 idx randidx(idx a = std::numeric_limits<idx>::min(),
59  idx b = std::numeric_limits<idx>::max())
60 {
61  std::uniform_int_distribution<idx> uid(a, b);
62 
64 }
65 
77 template<typename Derived>
78 Derived rand(idx rows, idx cols, double a = 0, double b = 1)
79 {
80  throw Exception("qpp::rand()", Exception::Type::UNDEFINED_TYPE);
81 }
82 
104 template<>
105 inline dmat rand(idx rows, idx cols, double a, double b)
106 {
107  if (rows == 0 || cols == 0)
108  throw Exception("qpp::rand()", Exception::Type::ZERO_SIZE);
109 
110  return dmat::Zero(rows, cols).unaryExpr(
111  [a, b](double)
112  {
113  return rand(a, b);
114  });
115 }
116 
138 template<>
139 inline cmat rand(idx rows, idx cols, double a, double b)
140 {
141  if (rows == 0 || cols == 0)
142  throw Exception("qpp::rand()", Exception::Type::ZERO_SIZE);
143 
144  return rand<dmat>(rows, cols, a, b).cast<cplx>() +
145  1_i * rand<dmat>(rows, cols, a, b).cast<cplx>();
146 }
147 
159 template<typename Derived>
160 Derived randn(idx rows, idx cols, double mean = 0,
161  double sigma = 1)
162 {
163  throw Exception("qpp::randn()", Exception::Type::UNDEFINED_TYPE);
164 }
165 
187 template<>
188 inline dmat randn(idx rows, idx cols,
189  double mean, double sigma)
190 {
191  if (rows == 0 || cols == 0)
192  throw Exception("qpp::randn()", Exception::Type::ZERO_SIZE);
193 
194  std::normal_distribution<> nd(mean, sigma);
195 
196  return dmat::Zero(rows, cols).unaryExpr(
197  [&nd](double)
198  {
200  });
201 
202 }
203 
225 template<>
226 inline cmat randn(idx rows, idx cols,
227  double mean, double sigma)
228 {
229  if (rows == 0 || cols == 0)
230  throw Exception("qpp::randn()", Exception::Type::ZERO_SIZE);
231 
232  return randn<dmat>(rows, cols, mean, sigma).cast<cplx>() +
233  1_i * randn<dmat>(rows, cols, mean, sigma).cast<cplx>();
234 }
235 
244 double randn(double mean = 0, double sigma = 1)
245 {
246  std::normal_distribution<> nd(mean, sigma);
247 
249 }
250 
258 // ~3 times slower than Toby Cubitt's MATLAB corresponding routine,
259 // because Eigen 3 QR algorithm is not parallelized
260 {
261  if (D == 0)
262  throw Exception("qpp::randU()", Exception::Type::DIMS_INVALID);
263 
264  cmat X(D, D);
265 
266  X = 1 / std::sqrt(2.) * randn < cmat > (D, D);
267  Eigen::HouseholderQR<cmat> qr(X);
268 
269  cmat Q = qr.householderQ();
270  // phase correction so that the resultant matrix is
271  // uniformly distributed according to the Haar measure
272 
273  Eigen::VectorXcd phases = (rand<dmat>(D, 1)).cast<cplx>();
274  for (idx i = 0; i < static_cast<idx>(phases.rows()); ++i)
275  phases(i) = std::exp(2 * pi * 1_i * phases(i));
276 
277  Q = Q * phases.asDiagonal();
278 
279  return Q;
280 }
281 
289 cmat randV(idx Din, idx Dout)
290 {
291  if (Din == 0 || Dout == 0 || Din > Dout)
292  throw Exception("qpp::randV()", Exception::Type::DIMS_INVALID);
293 
294  return randU(Dout).block(0, 0, Dout, Din);
295 }
296 
307 std::vector<cmat> randkraus(idx N, idx D)
308 {
309  if (N == 0)
310  throw Exception("qpp::randkraus()", Exception::Type::OUT_OF_RANGE);
311  if (D == 0)
312  throw Exception("qpp::randkraus()", Exception::Type::DIMS_INVALID);
313 
314  std::vector<cmat> result(N);
315  for (idx i = 0; i < N; ++i)
316  result[i] = cmat::Zero(D, D);
317 
318  cmat Fk(D, D);
319  cmat U = randU(N * D);
320 
321 #pragma omp parallel for collapse(3)
322  for (idx k = 0; k < N; ++k)
323  for (idx a = 0; a < D; ++a)
324  for (idx b = 0; b < D; ++b)
325  result[k](a, b) = U(a * N + k, b * N);
326 
327  return result;
328 }
329 
337 {
338  if (D == 0)
339  throw Exception("qpp::randH()", Exception::Type::DIMS_INVALID);
340 
341  cmat H = 2 * rand<cmat>(D, D) - (1. + 1_i) * cmat::Ones(D, D);
342 
343  return H + adjoint(H);
344 }
345 
353 {
354  if (D == 0)
355  throw Exception("qpp::randket()", Exception::Type::DIMS_INVALID);
356 
357  /* slow
358  ket kt = ket::Ones(D);
359  ket result = static_cast<ket>(randU(D) * kt);
360 
361  return result;
362  */
363 
364  ket kt = randn < cmat > (D, 1);
365 
366  return kt / norm(kt);
367 }
368 
376 {
377  if (D == 0)
378  throw Exception("qpp::randrho()", Exception::Type::DIMS_INVALID);
379 
380  cmat result = 10 * randH(D);
381  result = result * adjoint(result);
382 
383  return result / trace(result);
384 }
385 
395 std::vector<idx> randperm(idx n)
396 {
397  if (n == 0)
398  throw Exception("qpp::randperm()", Exception::Type::PERM_INVALID);
399 
400  std::vector<idx> result(n);
401 
402  // fill in increasing order
403  std::iota(std::begin(result), std::end(result), 0);
404  // shuffle
405  std::shuffle(std::begin(result), std::end(result),
407 
408  return result;
409 }
410 
411 } /* namespace qpp */
412 
413 #endif /* RANDOM_H_ */
static thread_local RandomDevices & get_thread_local_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:100
idx randidx(idx a=std::numeric_limits< idx >::min(), idx b=std::numeric_limits< idx >::max())
Generates a random index (idx) uniformly distributed in the interval [a, b].
Definition: random.h:58
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:61
std::vector< cmat > randkraus(idx N, idx D)
Generates a set of random Kraus operators.
Definition: random.h:307
ket randket(idx D)
Generates a random normalized ket (pure state vector)
Definition: random.h:352
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:46
Quantum++ main namespace.
Definition: codes.h:30
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:231
std::vector< idx > randperm(idx n)
Generates a random uniformly distributed permutation.
Definition: random.h:395
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:84
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:41
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
Derived randn(idx rows, idx cols, double mean=0, double sigma=1)
Generates a random matrix with entries normally distributed in N(mean, sigma)
Definition: random.h:160
cmat randU(idx D)
Generates a random unitary matrix.
Definition: random.h:257
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:121
constexpr double pi
Definition: constants.h:87
cmat randH(idx D)
Generates a random Hermitian matrix.
Definition: random.h:336
double rand(double a=0, double b=1)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:43
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:56
cmat randrho(idx D)
Generates a random density matrix.
Definition: random.h:375
std::size_t idx
Non-negative integer index.
Definition: types.h:36
cmat randV(idx Din, idx Dout)
Generates a random isometry matrix.
Definition: random.h:289