Quantum++  v1.0-rc2
A modern 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 namespace qpp
31 {
41 inline double rand(double a, double b)
42 {
43  // EXCEPTION CHECKS
44 
45  if (a >= b)
46  throw exception::OutOfRange("qpp::rand()");
47  // END EXCEPTION CHECKS
48 
49  std::uniform_real_distribution<> ud(a, b);
50 
51 #ifdef NO_THREAD_LOCAL_
52  return ud(RandomDevices::get_instance().get_prng());
53 #else
54  return ud(RandomDevices::get_thread_local_instance().get_prng());
55 #endif // NO_THREAD_LOCAL_
56 }
57 
69 inline bigint rand(bigint a, bigint b)
70 {
71  // EXCEPTION CHECKS
72 
73  if (a > b)
74  throw exception::OutOfRange("qpp::rand()");
75  // END EXCEPTION CHECKS
76 
77  std::uniform_int_distribution<bigint> uid(a, b);
78 
79 #ifdef NO_THREAD_LOCAL_
80  return uid(RandomDevices::get_instance().get_prng());
81 #else
82  return uid(RandomDevices::get_thread_local_instance().get_prng());
83 #endif // NO_THREAD_LOCAL_
84 }
85 
94 inline idx randidx(idx a = std::numeric_limits<idx>::min(),
95  idx b = std::numeric_limits<idx>::max())
96 {
97  // EXCEPTION CHECKS
98 
99  if (a > b)
100  throw exception::OutOfRange("qpp::randidx()");
101  // END EXCEPTION CHECKS
102 
103  std::uniform_int_distribution<idx> uid(a, b);
104 
105 #ifdef NO_THREAD_LOCAL_
106  return uid(RandomDevices::get_instance().get_prng());
107 #else
108  return uid(RandomDevices::get_thread_local_instance().get_prng());
109 #endif // NO_THREAD_LOCAL_
110 }
111 
123 template<typename Derived>
124 Derived rand(idx rows, idx cols, double a = 0, double b = 1)
125 {
126  // silence -Wunused-parameter in clang++
127  (void) rows;
128  (void) cols;
129  (void) a;
130  (void) b;
131  throw exception::UndefinedType("qpp::rand()");
132 }
133 
155 template<>
156 inline dmat rand(idx rows, idx cols, double a, double b)
157 {
158  // EXCEPTION CHECKS
159 
160  if (rows == 0 || cols == 0)
161  throw exception::ZeroSize("qpp::rand()");
162  if (a >= b)
163  throw exception::OutOfRange("qpp::rand()");
164  // END EXCEPTION CHECKS
165 
166  return dmat::Zero(rows, cols).unaryExpr(
167  [a, b](double)
168  {
169  return rand(a, b);
170  });
171 }
172 
194 template<>
195 inline cmat rand(idx rows, idx cols, double a, double b)
196 {
197  // EXCEPTION CHECKS
198 
199  if (rows == 0 || cols == 0)
200  throw exception::ZeroSize("qpp::rand()");
201  if (a >= b)
202  throw exception::OutOfRange("qpp::rand()");
203  // END EXCEPTION CHECKS
204 
205  return rand<dmat>(rows, cols, a, b).cast<cplx>() +
206  1_i * rand<dmat>(rows, cols, a, b).cast<cplx>();
207 }
208 
220 template<typename Derived>
221 Derived randn(idx rows, idx cols, double mean = 0, double sigma = 1)
222 {
223  // silence -Wunused-parameter in clang++
224  (void) rows;
225  (void) cols;
226  (void) mean;
227  (void) sigma;
228  throw exception::UndefinedType("qpp::randn()");
229 }
230 
252 template<>
253 inline dmat randn(idx rows, idx cols, double mean, double sigma)
254 {
255  // EXCEPTION CHECKS
256 
257  if (rows == 0 || cols == 0)
258  throw exception::ZeroSize("qpp::randn()");
259  // END EXCEPTION CHECKS
260 
261  std::normal_distribution<> nd(mean, sigma);
262 
263  return dmat::Zero(rows, cols).unaryExpr(
264  [&nd](double)
265  {
266 #ifdef NO_THREAD_LOCAL_
267  return nd(RandomDevices::get_instance().get_prng());
268 #else
269  return nd(
271 #endif // NO_THREAD_LOCAL_
272  });
273 }
274 
296 template<>
297 inline cmat randn(idx rows, idx cols, double mean, double sigma)
298 {
299  // EXCEPTION CHECKS
300 
301  if (rows == 0 || cols == 0)
302  throw exception::ZeroSize("qpp::randn()");
303  // END EXCEPTION CHECKS
304 
305  return randn<dmat>(rows, cols, mean, sigma).cast<cplx>() +
306  1_i * randn<dmat>(rows, cols, mean, sigma).cast<cplx>();
307 }
308 
317 inline double randn(double mean = 0, double sigma = 1)
318 {
319  std::normal_distribution<> nd(mean, sigma);
320 
321 #ifdef NO_THREAD_LOCAL_
322  return nd(RandomDevices::get_instance().get_prng());
323 #else
324  return nd(RandomDevices::get_thread_local_instance().get_prng());
325 #endif // NO_THREAD_LOCAL_
326 }
327 
334 inline cmat randU(idx D = 2)
335 // ~3 times slower than Toby Cubitt's MATLAB corresponding routine,
336 // because Eigen 3 QR algorithm is not parallelized
337 {
338  // EXCEPTION CHECKS
339 
340  if (D == 0)
341  throw exception::DimsInvalid("qpp::randU()");
342  // END EXCEPTION CHECKS
343 
344  cmat X = 1 / std::sqrt(2.) * randn<cmat>(D, D);
345  Eigen::HouseholderQR<cmat> qr(X);
346 
347  cmat Q = qr.householderQ();
348  // phase correction so that the resultant matrix is
349  // uniformly distributed according to the Haar measure
350 
351  Eigen::VectorXcd phases = (rand<dmat>(D, 1)).cast<cplx>();
352  for (idx i = 0; i < static_cast<idx>(phases.rows()); ++i)
353  phases(i) = std::exp(2 * pi * 1_i * phases(i));
354 
355  Q = Q * phases.asDiagonal();
356 
357  return Q;
358 }
359 
367 inline cmat randV(idx Din, idx Dout)
368 {
369  // EXCEPTION CHECKS
370 
371  if (Din == 0 || Dout == 0 || Din > Dout)
372  throw exception::DimsInvalid("qpp::randV()");
373  // END EXCEPTION CHECKS
374 
375  return randU(Dout).block(0, 0, Dout, Din);
376 }
377 
388 inline std::vector<cmat> randkraus(idx N, idx D = 2)
389 {
390  // EXCEPTION CHECKS
391 
392  if (N == 0)
393  throw exception::OutOfRange("qpp::randkraus()");
394  if (D == 0)
395  throw exception::DimsInvalid("qpp::randkraus()");
396  // END EXCEPTION CHECKS
397 
398  std::vector<cmat> result(N);
399  for (idx i = 0; i < N; ++i)
400  result[i] = cmat::Zero(D, D);
401 
402  cmat Fk(D, D);
403  cmat U = randU(N * D);
404 
405 #ifdef WITH_OPENMP_
406 #pragma omp parallel for collapse(3)
407 #endif // WITH_OPENMP_
408  for (idx k = 0; k < N; ++k)
409  for (idx a = 0; a < D; ++a)
410  for (idx b = 0; b < D; ++b)
411  result[k](a, b) = U(a * N + k, b * N);
412 
413  return result;
414 }
415 
422 inline cmat randH(idx D = 2)
423 {
424  // EXCEPTION CHECKS
425 
426  if (D == 0)
427  throw exception::DimsInvalid("qpp::randH()");
428  // END EXCEPTION CHECKS
429 
430  cmat H = 2 * rand<cmat>(D, D) - (1. + 1_i) * cmat::Ones(D, D);
431 
432  return H + adjoint(H);
433 }
434 
441 inline ket randket(idx D = 2)
442 {
443  // EXCEPTION CHECKS
444 
445  if (D == 0)
446  throw exception::DimsInvalid("qpp::randket()");
447  // END EXCEPTION CHECKS
448 
449  /* slow
450  ket kt = ket::Ones(D);
451  ket result = static_cast<ket>(randU(D) * kt);
452 
453  return result;
454  */
455 
456  ket kt = randn<cmat>(D, 1);
457 
458  return kt / norm(kt);
459 }
460 
467 inline cmat randrho(idx D = 2)
468 {
469  // EXCEPTION CHECKS
470 
471  if (D == 0)
472  throw exception::DimsInvalid("qpp::randrho()");
473  // END EXCEPTION CHECKS
474 
475  cmat result = 10 * randH(D);
476  result = result * adjoint(result);
477 
478  return result / trace(result);
479 }
480 
490 inline std::vector<idx> randperm(idx N)
491 {
492  // EXCEPTION CHECKS
493 
494  if (N == 0)
495  throw exception::PermInvalid("qpp::randperm()");
496  // END EXCEPTION CHECKS
497 
498  std::vector<idx> result(N);
499 
500  // fill in increasing order
501  std::iota(std::begin(result), std::end(result), 0);
502  // shuffle
503 #ifdef NO_THREAD_LOCAL_
504  std::shuffle(std::begin(result), std::end(result),
505  RandomDevices::get_instance().get_prng());
506 #else
507  std::shuffle(std::begin(result), std::end(result),
509 
510 #endif // NO_THREAD_LOCAL_
511 
512  return result;
513 }
514 
522 inline std::vector<double> randprob(idx N)
523 {
524  // EXCEPTION CHECKS
525 
526  if (N == 0)
527  throw exception::OutOfRange("qpp::randprob()");
528  // END EXCEPTION CHECKS
529 
530  std::vector<double> result(N);
531 
532  // generate
533  std::exponential_distribution<> ed(1);
534  for (idx i = 0; i < N; ++i)
535  {
536 #ifdef NO_THREAD_LOCAL_
537  result[i] = ed(qpp::RandomDevices::get_instance().get_prng());
538 #else
539  result[i] =
541 #endif // NO_THREAD_LOCAL_
542  }
543 
544  // normalize
545  double sumprob = sum(result);
546  for (idx i = 0; i < N; ++i)
547  result[i] /= sumprob;
548 
549  return result;
550 }
551 
552 } /* namespace qpp */
553 
554 #endif /* RANDOM_H_ */
std::vector< idx > randperm(idx N)
Generates a random uniformly distributed permutation.
Definition: random.h:490
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:94
Not defined for this type exception.
Definition: exception.h:619
cmat randH(idx D=2)
Generates a random Hermitian matrix.
Definition: random.h:422
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:65
cmat randrho(idx D=2)
Generates a random density matrix.
Definition: random.h:467
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:50
Quantum++ main namespace.
Definition: codes.h:30
ket randket(idx D=2)
Generates a random normalized ket (pure state vector)
Definition: random.h:441
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:248
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:204
Invalid dimension(s) exception.
Definition: exception.h:287
Invalid permutation exception.
Definition: exception.h:412
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:81
std::vector< double > randprob(idx N)
Generates a random probability vector uniformly distributed over the probability simplex.
Definition: random.h:522
static RandomDevices & get_thread_local_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:101
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:45
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:221
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:41
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:124
constexpr double pi
Definition: constants.h:76
Parameter out of range exception.
Definition: exception.h:567
static RandomDevices & get_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:89
std::vector< cmat > randkraus(idx N, idx D=2)
Generates a set of random Kraus operators.
Definition: random.h:388
long long int bigint
Big integer.
Definition: types.h:40
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:60
std::size_t idx
Non-negative integer index.
Definition: types.h:35
cmat randV(idx Din, idx Dout)
Generates a random isometry matrix.
Definition: random.h:367
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:205
cmat randU(idx D=2)
Generates a random unitary matrix.
Definition: random.h:334
Object has zero size exception.
Definition: exception.h:134