Quantum++  v1.2
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 - 2019 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  auto& gen =
54 #ifdef NO_THREAD_LOCAL_
56 #else
58 #endif
59 
60  return ud(gen);
61 }
62 
74 inline bigint rand(bigint a, bigint b) {
75  // EXCEPTION CHECKS
76 
77  if (a > b)
78  throw exception::OutOfRange("qpp::rand()");
79  // END EXCEPTION CHECKS
80 
81  std::uniform_int_distribution<bigint> uid(a, b);
82 
83  auto& gen =
84 #ifdef NO_THREAD_LOCAL_
86 #else
88 #endif
89 
90  return uid(gen);
91 }
92 
101 inline idx randidx(idx a = std::numeric_limits<idx>::min(),
102  idx b = std::numeric_limits<idx>::max()) {
103  // EXCEPTION CHECKS
104 
105  if (a > b)
106  throw exception::OutOfRange("qpp::randidx()");
107  // END EXCEPTION CHECKS
108 
109  std::uniform_int_distribution<idx> uid(a, b);
110  auto& gen =
111 #ifdef NO_THREAD_LOCAL_
113 #else
115 #endif
116 
117  return uid(gen);
118 }
119 
131 template <typename Derived>
132 Derived rand(idx rows QPP_UNUSED_, idx cols QPP_UNUSED_,
133  double a QPP_UNUSED_ = 0, double b QPP_UNUSED_ = 1) {
134  throw exception::UndefinedType("qpp::rand()");
135 }
136 
158 template <>
159 inline dmat rand(idx rows, idx cols, double a, double b) {
160  // EXCEPTION CHECKS
161 
162  if (rows == 0 || cols == 0)
163  throw exception::ZeroSize("qpp::rand()");
164  if (a >= b)
165  throw exception::OutOfRange("qpp::rand()");
166  // END EXCEPTION CHECKS
167 
168  return dmat::Zero(rows, cols).unaryExpr([a, b](double) {
169  return rand(a, b);
170  });
171 }
172 
194 template <>
195 inline cmat rand(idx rows, idx cols, double a, double b) {
196  // EXCEPTION CHECKS
197 
198  if (rows == 0 || cols == 0)
199  throw exception::ZeroSize("qpp::rand()");
200  if (a >= b)
201  throw exception::OutOfRange("qpp::rand()");
202  // END EXCEPTION CHECKS
203 
204  return rand<dmat>(rows, cols, a, b).cast<cplx>() +
205  1_i * rand<dmat>(rows, cols, a, b).cast<cplx>();
206 }
207 
219 template <typename Derived>
220 Derived randn(idx rows QPP_UNUSED_, idx cols QPP_UNUSED_,
221  double mean QPP_UNUSED_ = 0, double sigma QPP_UNUSED_ = 1) {
222  throw exception::UndefinedType("qpp::randn()");
223 }
224 
246 template <>
247 inline dmat randn(idx rows, idx cols, double mean, double sigma) {
248  // EXCEPTION CHECKS
249 
250  if (rows == 0 || cols == 0)
251  throw exception::ZeroSize("qpp::randn()");
252  // END EXCEPTION CHECKS
253 
254  std::normal_distribution<> nd(mean, sigma);
255  auto& gen =
256 #ifdef NO_THREAD_LOCAL_
258 #else
260 #endif
261 
262  return dmat::Zero(rows, cols).unaryExpr([&nd, &gen](double) {
263  return nd(gen);
264  });
265 }
266 
288 template <>
289 inline cmat randn(idx rows, idx cols, double mean, double sigma) {
290  // EXCEPTION CHECKS
291 
292  if (rows == 0 || cols == 0)
293  throw exception::ZeroSize("qpp::randn()");
294  // END EXCEPTION CHECKS
295 
296  return randn<dmat>(rows, cols, mean, sigma).cast<cplx>() +
297  1_i * randn<dmat>(rows, cols, mean, sigma).cast<cplx>();
298 }
299 
308 inline double randn(double mean = 0, double sigma = 1) {
309  std::normal_distribution<> nd(mean, sigma);
310  auto& gen =
311 #ifdef NO_THREAD_LOCAL_
313 #else
315 #endif
316 
317  return nd(gen);
318 }
319 
326 inline cmat randU(idx D = 2)
327 // ~3 times slower than Toby Cubitt's MATLAB corresponding routine,
328 // because Eigen 3 QR algorithm is not parallelized
329 {
330  // EXCEPTION CHECKS
331 
332  if (D == 0)
333  throw exception::DimsInvalid("qpp::randU()");
334  // END EXCEPTION CHECKS
335 
336  cmat X = 1 / std::sqrt(2.) * randn<cmat>(D, D);
337  Eigen::HouseholderQR<cmat> qr(X);
338 
339  cmat Q = qr.householderQ();
340  // phase correction so that the resultant matrix is
341  // uniformly distributed according to the Haar measure
342 
343  Eigen::VectorXcd phases = (rand<dmat>(D, 1)).cast<cplx>();
344  for (idx i = 0; i < static_cast<idx>(phases.rows()); ++i)
345  phases(i) = std::exp(2 * pi * 1_i * phases(i));
346 
347  Q = Q * phases.asDiagonal();
348 
349  return Q;
350 }
351 
359 inline cmat randV(idx Din, idx Dout) {
360  // EXCEPTION CHECKS
361 
362  if (Din == 0 || Dout == 0 || Din > Dout)
363  throw exception::DimsInvalid("qpp::randV()");
364  // END EXCEPTION CHECKS
365 
366  return randU(Dout).block(0, 0, Dout, Din);
367 }
368 
379 inline std::vector<cmat> randkraus(idx N, idx D = 2) {
380  // EXCEPTION CHECKS
381 
382  if (N == 0)
383  throw exception::OutOfRange("qpp::randkraus()");
384  if (D == 0)
385  throw exception::DimsInvalid("qpp::randkraus()");
386  // END EXCEPTION CHECKS
387 
388  std::vector<cmat> result(N);
389  for (idx i = 0; i < N; ++i)
390  result[i] = cmat::Zero(D, D);
391 
392  cmat Fk(D, D);
393  cmat U = randU(N * D);
394 
395 #ifdef WITH_OPENMP_
396 #pragma omp parallel for collapse(3)
397 #endif // WITH_OPENMP_
398  for (idx k = 0; k < N; ++k)
399  for (idx a = 0; a < D; ++a)
400  for (idx b = 0; b < D; ++b)
401  result[k](a, b) = U(a * N + k, b * N);
402 
403  return result;
404 }
405 
412 inline cmat randH(idx D = 2) {
413  // EXCEPTION CHECKS
414 
415  if (D == 0)
416  throw exception::DimsInvalid("qpp::randH()");
417  // END EXCEPTION CHECKS
418 
419  cmat H = 2 * rand<cmat>(D, D) - (1. + 1_i) * cmat::Ones(D, D);
420 
421  return H + H.adjoint();
422 }
423 
430 inline ket randket(idx D = 2) {
431  // EXCEPTION CHECKS
432 
433  if (D == 0)
434  throw exception::DimsInvalid("qpp::randket()");
435  // END EXCEPTION CHECKS
436 
437  /* slow
438  ket kt = ket::Ones(D);
439  ket result = static_cast<ket>(randU(D) * kt);
440 
441  return result;
442  */
443 
444  ket kt = randn<cmat>(D, 1);
445 
446  return kt / kt.norm();
447 }
448 
455 inline cmat randrho(idx D = 2) {
456  // EXCEPTION CHECKS
457 
458  if (D == 0)
459  throw exception::DimsInvalid("qpp::randrho()");
460  // END EXCEPTION CHECKS
461 
462  cmat result = 10 * randH(D);
463  result = result * result.adjoint();
464 
465  return result / result.trace();
466 }
467 
477 inline std::vector<idx> randperm(idx N) {
478  // EXCEPTION CHECKS
479 
480  if (N == 0)
481  throw exception::PermInvalid("qpp::randperm()");
482  // END EXCEPTION CHECKS
483 
484  std::vector<idx> result(N);
485 
486  // fill in increasing order
487  std::iota(std::begin(result), std::end(result), 0);
488 
489  // shuffle
490  auto& gen =
491 #ifdef NO_THREAD_LOCAL_
493 #else
495 #endif
496  std::shuffle(std::begin(result), std::end(result), gen);
497 
498  return result;
499 }
500 
508 inline std::vector<double> randprob(idx N) {
509  // EXCEPTION CHECKS
510 
511  if (N == 0)
512  throw exception::OutOfRange("qpp::randprob()");
513  // END EXCEPTION CHECKS
514 
515  std::vector<double> result(N);
516 
517  // generate
518  std::exponential_distribution<> ed(1);
519  auto& gen =
520 #ifdef NO_THREAD_LOCAL_
522 #else
524 #endif
525  for (idx i = 0; i < N; ++i)
526  result[i] = ed(gen);
527 
528  // normalize
529  double sumprob = std::accumulate(std::begin(result), std::end(result), 0.0);
530  for (idx i = 0; i < N; ++i)
531  result[i] /= sumprob;
532 
533  return result;
534 }
535 
536 } /* namespace qpp */
537 
538 #endif /* RANDOM_H_ */
std::vector< idx > randperm(idx N)
Generates a random uniformly distributed permutation.
Definition: random.h:477
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:101
Not defined for this type exception.
Definition: exception.h:556
cmat randH(idx D=2)
Generates a random Hermitian matrix.
Definition: random.h:412
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:69
std::mt19937 & get_prng()
Returns a reference to the internal PRNG object.
Definition: random_devices.h:61
cmat randrho(idx D=2)
Generates a random density matrix.
Definition: random.h:455
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
Quantum++ main namespace.
Definition: circuits.h:35
ket randket(idx D=2)
Generates a random normalized ket (pure state vector)
Definition: random.h:430
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
std::vector< double > randprob(idx N)
Generates a random probability vector uniformly distributed over the probability simplex.
Definition: random.h:508
static RandomDevices & get_thread_local_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:103
Derived randn(idx rows QPP_UNUSED_, idx cols QPP_UNUSED_, double mean QPP_UNUSED_=0, double sigma QPP_UNUSED_=1)
Generates a random matrix with entries normally distributed in N(mean, sigma)
Definition: random.h:220
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:49
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:45
constexpr double pi
Definition: constants.h:72
Argument 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:379
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, make sure you use an unsigned type.
Definition: types.h:39
cmat randV(idx Din, idx Dout)
Generates a random isometry matrix.
Definition: random.h:359
cmat randU(idx D=2)
Generates a random unitary matrix.
Definition: random.h:326
Object has zero size exception.
Definition: exception.h:134