Quantum++  v1.0.0-beta3
C++11 quantum computing library
random.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2017 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 
44 inline double rand(double a, double b)
45 {
46  // EXCEPTION CHECKS
47 
48  if (a >= b)
50  // END EXCEPTION CHECKS
51 
52  std::uniform_real_distribution<> ud(a, b);
53 
54 #ifdef NO_THREAD_LOCAL_
55  return ud(RandomDevices::get_instance().rng_);
56 #else
58 #endif // NO_THREAD_LOCAL_
59 }
60 
70 inline bigint rand(bigint a, bigint b)
71 {
72  // EXCEPTION CHECKS
73 
74  if (a > b)
76  // END EXCEPTION CHECKS
77 
78  std::uniform_int_distribution<bigint> uid(a, b);
79 
80 #ifdef NO_THREAD_LOCAL_
81  return uid(RandomDevices::get_instance().rng_);
82 #else
84 #endif // NO_THREAD_LOCAL_
85 }
86 
96 inline idx randidx(idx a = std::numeric_limits<idx>::min(),
97  idx b = std::numeric_limits<idx>::max())
98 {
99  // EXCEPTION CHECKS
100 
101  if (a > b)
103  // END EXCEPTION CHECKS
104 
105  std::uniform_int_distribution<idx> uid(a, b);
106 
107 #ifdef NO_THREAD_LOCAL_
108  return uid(RandomDevices::get_instance().rng_);
109 #else
110  return uid(RandomDevices::get_thread_local_instance().rng_);
111 #endif // NO_THREAD_LOCAL_
112 }
113 
125 template<typename Derived>
126 Derived rand(idx rows, idx cols, double a = 0, double b = 1)
127 {
128  // silence -Wunused-parameter in clang++
129  (void) rows;
130  (void) cols;
131  (void) a;
132  (void) b;
133  throw Exception("qpp::rand()", Exception::Type::UNDEFINED_TYPE);
134 }
135 
157 template<>
158 inline dmat rand(idx rows, idx cols, double a, double b)
159 {
160  // EXCEPTION CHECKS
161 
162  if (rows == 0 || cols == 0)
163  throw Exception("qpp::rand()", Exception::Type::ZERO_SIZE);
164  if (a >= b)
166  // END EXCEPTION CHECKS
167 
168  return dmat::Zero(rows, cols).unaryExpr(
169  [a, b](double)
170  {
171  return rand(a, b);
172  });
173 }
174 
196 template<>
197 inline cmat rand(idx rows, idx cols, double a, double b)
198 {
199  // EXCEPTION CHECKS
200 
201  if (rows == 0 || cols == 0)
202  throw Exception("qpp::rand()", Exception::Type::ZERO_SIZE);
203  if (a >= b)
205  // END EXCEPTION CHECKS
206 
207  return rand<dmat>(rows, cols, a, b).cast<cplx>() +
208  1_i * rand<dmat>(rows, cols, a, b).cast<cplx>();
209 }
210 
222 template<typename Derived>
223 Derived randn(idx rows, idx cols, double mean = 0,
224  double sigma = 1)
225 {
226  // silence -Wunused-parameter in clang++
227  (void) rows;
228  (void) cols;
229  (void) mean;
230  (void) sigma;
231  throw Exception("qpp::randn()", Exception::Type::UNDEFINED_TYPE);
232 }
233 
255 template<>
256 inline dmat randn(idx rows, idx cols,
257  double mean, double sigma)
258 {
259  // EXCEPTION CHECKS
260 
261  if (rows == 0 || cols == 0)
262  throw Exception("qpp::randn()", Exception::Type::ZERO_SIZE);
263  // END EXCEPTION CHECKS
264 
265  std::normal_distribution<> nd(mean, sigma);
266 
267  return dmat::Zero(rows, cols).unaryExpr(
268  [&nd](double)
269  {
270 #ifdef NO_THREAD_LOCAL_
271  return nd(RandomDevices::get_instance().rng_);
272 #else
274 #endif // NO_THREAD_LOCAL_
275  });
276 }
277 
299 template<>
300 inline cmat randn(idx rows, idx cols, double mean, double sigma)
301 {
302  // EXCEPTION CHECKS
303 
304  if (rows == 0 || cols == 0)
305  throw Exception("qpp::randn()", Exception::Type::ZERO_SIZE);
306  // END EXCEPTION CHECKS
307 
308  return randn<dmat>(rows, cols, mean, sigma).cast<cplx>() +
309  1_i * randn<dmat>(rows, cols, mean, sigma).cast<cplx>();
310 }
311 
320 inline double randn(double mean = 0, double sigma = 1)
321 {
322  std::normal_distribution<> nd(mean, sigma);
323 
324 #ifdef NO_THREAD_LOCAL_
325  return nd(RandomDevices::get_instance().rng_);
326 #else
328 #endif // NO_THREAD_LOCAL_
329 }
330 
337 inline cmat randU(idx D)
338 // ~3 times slower than Toby Cubitt's MATLAB corresponding routine,
339 // because Eigen 3 QR algorithm is not parallelized
340 {
341  // EXCEPTION CHECKS
342 
343  if (D == 0)
344  throw Exception("qpp::randU()", Exception::Type::DIMS_INVALID);
345  // END EXCEPTION CHECKS
346 
347  cmat X = 1 / std::sqrt(2.) * randn<cmat>(D, D);
348  Eigen::HouseholderQR<cmat> qr(X);
349 
350  cmat Q = qr.householderQ();
351  // phase correction so that the resultant matrix is
352  // uniformly distributed according to the Haar measure
353 
354  Eigen::VectorXcd phases = (rand<dmat>(D, 1)).cast<cplx>();
355  for (idx i = 0; i < static_cast<idx>(phases.rows()); ++i)
356  phases(i) = std::exp(2 * pi * 1_i * phases(i));
357 
358  Q = Q * phases.asDiagonal();
359 
360  return Q;
361 }
362 
370 inline cmat randV(idx Din, idx Dout)
371 {
372  // EXCEPTION CHECKS
373 
374  if (Din == 0 || Dout == 0 || Din > Dout)
375  throw Exception("qpp::randV()", Exception::Type::DIMS_INVALID);
376  // END EXCEPTION CHECKS
377 
378  return randU(Dout).block(0, 0, Dout, Din);
379 }
380 
391 inline std::vector<cmat> randkraus(idx N, idx D)
392 {
393  // EXCEPTION CHECKS
394 
395  if (N == 0)
396  throw Exception("qpp::randkraus()", Exception::Type::OUT_OF_RANGE);
397  if (D == 0)
398  throw Exception("qpp::randkraus()", Exception::Type::DIMS_INVALID);
399  // END EXCEPTION CHECKS
400 
401  std::vector<cmat> result(N);
402  for (idx i = 0; i < N; ++i)
403  result[i] = cmat::Zero(D, D);
404 
405  cmat Fk(D, D);
406  cmat U = randU(N * D);
407 
408 #ifdef WITH_OPENMP_
409 #pragma omp parallel for collapse(3)
410 #endif // WITH_OPENMP_
411  for (idx k = 0; k < N; ++k)
412  for (idx a = 0; a < D; ++a)
413  for (idx b = 0; b < D; ++b)
414  result[k](a, b) = U(a * N + k, b * N);
415 
416  return result;
417 }
418 
425 inline cmat randH(idx D)
426 {
427  // EXCEPTION CHECKS
428 
429  if (D == 0)
430  throw Exception("qpp::randH()", Exception::Type::DIMS_INVALID);
431  // END EXCEPTION CHECKS
432 
433  cmat H = 2 * rand<cmat>(D, D) - (1. + 1_i) * cmat::Ones(D, D);
434 
435  return H + adjoint(H);
436 }
437 
444 inline ket randket(idx D)
445 {
446  // EXCEPTION CHECKS
447 
448  if (D == 0)
449  throw Exception("qpp::randket()", Exception::Type::DIMS_INVALID);
450  // END EXCEPTION CHECKS
451 
452  /* slow
453  ket kt = ket::Ones(D);
454  ket result = static_cast<ket>(randU(D) * kt);
455 
456  return result;
457  */
458 
459  ket kt = randn<cmat>(D, 1);
460 
461  return kt / norm(kt);
462 }
463 
470 inline cmat randrho(idx D)
471 {
472  // EXCEPTION CHECKS
473 
474  if (D == 0)
475  throw Exception("qpp::randrho()", Exception::Type::DIMS_INVALID);
476  // END EXCEPTION CHECKS
477 
478  cmat result = 10 * randH(D);
479  result = result * adjoint(result);
480 
481  return result / trace(result);
482 }
483 
493 inline std::vector<idx> randperm(idx N)
494 {
495  // EXCEPTION CHECKS
496 
497  if (N == 0)
498  throw Exception("qpp::randperm()", Exception::Type::PERM_INVALID);
499  // END EXCEPTION CHECKS
500 
501  std::vector<idx> result(N);
502 
503  // fill in increasing order
504  std::iota(std::begin(result), std::end(result), 0);
505  // shuffle
506 #ifdef NO_THREAD_LOCAL_
507  std::shuffle(std::begin(result), std::end(result),
509 #else
510  std::shuffle(std::begin(result), std::end(result),
512 
513 #endif // NO_THREAD_LOCAL_
514 
515  return result;
516 }
517 
525 inline std::vector<double> randprob(idx N)
526 {
527  // EXCEPTION CHECKS
528 
529  if (N == 0)
530  throw Exception("qpp::randprob()", Exception::Type::PERM_INVALID);
531  // END EXCEPTION CHECKS
532 
533  std::vector<double> result(N);
534 
535  // generate
536  std::exponential_distribution<> ed(1);
537  for (idx i = 0; i < N; ++i)
538  {
539 #ifdef NO_THREAD_LOCAL_
540  result[i] = ed(qpp::RandomDevices::get_instance().rng_);
541 #else
542  result[i] = ed(qpp::RandomDevices::get_thread_local_instance().rng_);
543 #endif // NO_THREAD_LOCAL_
544  }
545 
546  // normalize
547  double sumprob = sum(result);
548  for (idx i = 0; i < N; ++i)
549  result[i] /= sumprob;
550 
551  return result;
552 }
553 
554 } /* namespace qpp */
555 
556 #endif /* RANDOM_H_ */
std::vector< idx > randperm(idx N)
Generates a random uniformly distributed permutation.
Definition: random.h:493
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:96
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:66
std::vector< cmat > randkraus(idx N, idx D)
Generates a set of random Kraus operators.
Definition: random.h:391
ket randket(idx D)
Generates a random normalized ket (pure state vector)
Definition: random.h:444
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:51
Quantum++ main namespace.
Definition: codes.h:30
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:252
double sigma(const std::vector< double > &prob, const Container &X, typename std::enable_if< is_iterable< Container >::value >::type *=nullptr)
Standard deviation.
Definition: statistics.h:207
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:84
std::vector< double > randprob(idx N)
Generates a random probability vector uniformly distributed over the probability simplex.
Definition: random.h:525
static RandomDevices & get_thread_local_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:102
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:46
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:223
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:44
cmat randU(idx D)
Generates a random unitary matrix.
Definition: random.h:337
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:127
constexpr double pi
Definition: constants.h:79
cmat randH(idx D)
Generates a random Hermitian matrix.
Definition: random.h:425
static RandomDevices & get_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:90
long long int bigint
Big integer.
Definition: types.h:41
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:61
cmat randrho(idx D)
Generates a random density matrix.
Definition: random.h:470
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:370
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:209