Quantum++  v1.0
A modern C++11 quantum computing library
random.h
Go to the documentation of this file.
1 /*
2  * This file is part of Quantum++.
3  *
4  * MIT License
5  *
6  * Copyright (c) 2013 - 2018 Vlad Gheorghiu (vgheorgh@gmail.com)
7  *
8  * Permission is hereby granted, free of charge, to any person obtaining a copy
9  * of this software and associated documentation files (the "Software"), to deal
10  * in the Software without restriction, including without limitation the rights
11  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12  * copies of the Software, and to permit persons to whom the Software is
13  * furnished to do so, subject to the following conditions:
14  *
15  * The above copyright notice and this permission notice shall be included in
16  * all copies or substantial portions of the Software.
17  *
18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24  * SOFTWARE.
25  */
26 
32 #ifndef RANDOM_H_
33 #define RANDOM_H_
34 
35 namespace qpp {
45 inline double rand(double a, double b) {
46  // EXCEPTION CHECKS
47 
48  if (a >= b)
49  throw exception::OutOfRange("qpp::rand()");
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().get_prng());
56 #else
57  return ud(RandomDevices::get_thread_local_instance().get_prng());
58 #endif // NO_THREAD_LOCAL_
59 }
60 
72 inline bigint rand(bigint a, bigint b) {
73  // EXCEPTION CHECKS
74 
75  if (a > b)
76  throw exception::OutOfRange("qpp::rand()");
77  // END EXCEPTION CHECKS
78 
79  std::uniform_int_distribution<bigint> uid(a, b);
80 
81 #ifdef NO_THREAD_LOCAL_
82  return uid(RandomDevices::get_instance().get_prng());
83 #else
84  return uid(RandomDevices::get_thread_local_instance().get_prng());
85 #endif // NO_THREAD_LOCAL_
86 }
87 
96 inline idx randidx(idx a = std::numeric_limits<idx>::min(),
97  idx b = std::numeric_limits<idx>::max()) {
98  // EXCEPTION CHECKS
99 
100  if (a > b)
101  throw exception::OutOfRange("qpp::randidx()");
102  // END EXCEPTION CHECKS
103 
104  std::uniform_int_distribution<idx> uid(a, b);
105 
106 #ifdef NO_THREAD_LOCAL_
107  return uid(RandomDevices::get_instance().get_prng());
108 #else
109  return uid(RandomDevices::get_thread_local_instance().get_prng());
110 #endif // NO_THREAD_LOCAL_
111 }
112 
124 template <typename Derived>
125 Derived rand(idx rows, idx cols, double a = 0, double b = 1) {
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  // EXCEPTION CHECKS
158 
159  if (rows == 0 || cols == 0)
160  throw exception::ZeroSize("qpp::rand()");
161  if (a >= b)
162  throw exception::OutOfRange("qpp::rand()");
163  // END EXCEPTION CHECKS
164 
165  return dmat::Zero(rows, cols).unaryExpr([a, b](double) {
166  return rand(a, b);
167  });
168 }
169 
191 template <>
192 inline cmat rand(idx rows, idx cols, double a, double b) {
193  // EXCEPTION CHECKS
194 
195  if (rows == 0 || cols == 0)
196  throw exception::ZeroSize("qpp::rand()");
197  if (a >= b)
198  throw exception::OutOfRange("qpp::rand()");
199  // END EXCEPTION CHECKS
200 
201  return rand<dmat>(rows, cols, a, b).cast<cplx>() +
202  1_i * rand<dmat>(rows, cols, a, b).cast<cplx>();
203 }
204 
216 template <typename Derived>
217 Derived randn(idx rows, idx cols, double mean = 0, double sigma = 1) {
218  // silence -Wunused-parameter in clang++
219  (void) rows;
220  (void) cols;
221  (void) mean;
222  (void) sigma;
223  throw exception::UndefinedType("qpp::randn()");
224 }
225 
247 template <>
248 inline dmat randn(idx rows, idx cols, double mean, double sigma) {
249  // EXCEPTION CHECKS
250 
251  if (rows == 0 || cols == 0)
252  throw exception::ZeroSize("qpp::randn()");
253  // END EXCEPTION CHECKS
254 
255  std::normal_distribution<> nd(mean, sigma);
256 
257  return dmat::Zero(rows, cols).unaryExpr([&nd](double) {
258 #ifdef NO_THREAD_LOCAL_
259  return nd(RandomDevices::get_instance().get_prng());
260 #else
261  return nd(RandomDevices::get_thread_local_instance().get_prng());
262 #endif // NO_THREAD_LOCAL_
263  });
264 }
265 
287 template <>
288 inline cmat randn(idx rows, idx cols, double mean, double sigma) {
289  // EXCEPTION CHECKS
290 
291  if (rows == 0 || cols == 0)
292  throw exception::ZeroSize("qpp::randn()");
293  // END EXCEPTION CHECKS
294 
295  return randn<dmat>(rows, cols, mean, sigma).cast<cplx>() +
296  1_i * randn<dmat>(rows, cols, mean, sigma).cast<cplx>();
297 }
298 
307 inline double randn(double mean = 0, double sigma = 1) {
308  std::normal_distribution<> nd(mean, sigma);
309 
310 #ifdef NO_THREAD_LOCAL_
311  return nd(RandomDevices::get_instance().get_prng());
312 #else
313  return nd(RandomDevices::get_thread_local_instance().get_prng());
314 #endif // NO_THREAD_LOCAL_
315 }
316 
323 inline cmat randU(idx D = 2)
324 // ~3 times slower than Toby Cubitt's MATLAB corresponding routine,
325 // because Eigen 3 QR algorithm is not parallelized
326 {
327  // EXCEPTION CHECKS
328 
329  if (D == 0)
330  throw exception::DimsInvalid("qpp::randU()");
331  // END EXCEPTION CHECKS
332 
333  cmat X = 1 / std::sqrt(2.) * randn<cmat>(D, D);
334  Eigen::HouseholderQR<cmat> qr(X);
335 
336  cmat Q = qr.householderQ();
337  // phase correction so that the resultant matrix is
338  // uniformly distributed according to the Haar measure
339 
340  Eigen::VectorXcd phases = (rand<dmat>(D, 1)).cast<cplx>();
341  for (idx i = 0; i < static_cast<idx>(phases.rows()); ++i)
342  phases(i) = std::exp(2 * pi * 1_i * phases(i));
343 
344  Q = Q * phases.asDiagonal();
345 
346  return Q;
347 }
348 
356 inline cmat randV(idx Din, idx Dout) {
357  // EXCEPTION CHECKS
358 
359  if (Din == 0 || Dout == 0 || Din > Dout)
360  throw exception::DimsInvalid("qpp::randV()");
361  // END EXCEPTION CHECKS
362 
363  return randU(Dout).block(0, 0, Dout, Din);
364 }
365 
376 inline std::vector<cmat> randkraus(idx N, idx D = 2) {
377  // EXCEPTION CHECKS
378 
379  if (N == 0)
380  throw exception::OutOfRange("qpp::randkraus()");
381  if (D == 0)
382  throw exception::DimsInvalid("qpp::randkraus()");
383  // END EXCEPTION CHECKS
384 
385  std::vector<cmat> result(N);
386  for (idx i = 0; i < N; ++i)
387  result[i] = cmat::Zero(D, D);
388 
389  cmat Fk(D, D);
390  cmat U = randU(N * D);
391 
392 #ifdef WITH_OPENMP_
393 #pragma omp parallel for collapse(3)
394 #endif // WITH_OPENMP_
395  for (idx k = 0; k < N; ++k)
396  for (idx a = 0; a < D; ++a)
397  for (idx b = 0; b < D; ++b)
398  result[k](a, b) = U(a * N + k, b * N);
399 
400  return result;
401 }
402 
409 inline cmat randH(idx D = 2) {
410  // EXCEPTION CHECKS
411 
412  if (D == 0)
413  throw exception::DimsInvalid("qpp::randH()");
414  // END EXCEPTION CHECKS
415 
416  cmat H = 2 * rand<cmat>(D, D) - (1. + 1_i) * cmat::Ones(D, D);
417 
418  return H + adjoint(H);
419 }
420 
427 inline ket randket(idx D = 2) {
428  // EXCEPTION CHECKS
429 
430  if (D == 0)
431  throw exception::DimsInvalid("qpp::randket()");
432  // END EXCEPTION CHECKS
433 
434  /* slow
435  ket kt = ket::Ones(D);
436  ket result = static_cast<ket>(randU(D) * kt);
437 
438  return result;
439  */
440 
441  ket kt = randn<cmat>(D, 1);
442 
443  return kt / norm(kt);
444 }
445 
452 inline cmat randrho(idx D = 2) {
453  // EXCEPTION CHECKS
454 
455  if (D == 0)
456  throw exception::DimsInvalid("qpp::randrho()");
457  // END EXCEPTION CHECKS
458 
459  cmat result = 10 * randH(D);
460  result = result * adjoint(result);
461 
462  return result / trace(result);
463 }
464 
474 inline std::vector<idx> randperm(idx N) {
475  // EXCEPTION CHECKS
476 
477  if (N == 0)
478  throw exception::PermInvalid("qpp::randperm()");
479  // END EXCEPTION CHECKS
480 
481  std::vector<idx> result(N);
482 
483  // fill in increasing order
484  std::iota(std::begin(result), std::end(result), 0);
485 // shuffle
486 #ifdef NO_THREAD_LOCAL_
487  std::shuffle(std::begin(result), std::end(result),
488  RandomDevices::get_instance().get_prng());
489 #else
490  std::shuffle(std::begin(result), std::end(result),
492 
493 #endif // NO_THREAD_LOCAL_
494 
495  return result;
496 }
497 
505 inline std::vector<double> randprob(idx N) {
506  // EXCEPTION CHECKS
507 
508  if (N == 0)
509  throw exception::OutOfRange("qpp::randprob()");
510  // END EXCEPTION CHECKS
511 
512  std::vector<double> result(N);
513 
514  // generate
515  std::exponential_distribution<> ed(1);
516  for (idx i = 0; i < N; ++i) {
517 #ifdef NO_THREAD_LOCAL_
518  result[i] = ed(qpp::RandomDevices::get_instance().get_prng());
519 #else
520  result[i] =
522 #endif // NO_THREAD_LOCAL_
523  }
524 
525  // normalize
526  double sumprob = sum(result);
527  for (idx i = 0; i < N; ++i)
528  result[i] /= sumprob;
529 
530  return result;
531 }
532 
533 } /* namespace qpp */
534 
535 #endif /* RANDOM_H_ */
std::vector< idx > randperm(idx N)
Generates a random uniformly distributed permutation.
Definition: random.h:474
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
Not defined for this type exception.
Definition: exception.h:556
cmat randH(idx D=2)
Generates a random Hermitian matrix.
Definition: random.h:409
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:69
cmat randrho(idx D=2)
Generates a random density matrix.
Definition: random.h:452
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
Quantum++ main namespace.
Definition: codes.h:35
ket randket(idx D=2)
Generates a random normalized ket (pure state vector)
Definition: random.h:427
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:242
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:197
Invalid dimension(s) exception.
Definition: exception.h:269
Invalid permutation exception.
Definition: exception.h:380
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:83
std::vector< double > randprob(idx N)
Generates a random probability vector uniformly distributed over the probability simplex.
Definition: random.h:505
static RandomDevices & get_thread_local_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:103
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:49
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:217
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:45
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:124
constexpr double pi
Definition: constants.h:80
Parameter out of range exception.
Definition: exception.h:515
static RandomDevices & get_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:92
std::vector< cmat > randkraus(idx N, idx D=2)
Generates a set of random Kraus operators.
Definition: random.h:376
long long int bigint
Big integer.
Definition: types.h:44
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
std::size_t idx
Non-negative integer index.
Definition: types.h:39
cmat randV(idx Din, idx Dout)
Generates a random isometry matrix.
Definition: random.h:356
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:201
cmat randU(idx D=2)
Generates a random unitary matrix.
Definition: random.h:323
Object has zero size exception.
Definition: exception.h:134