Quantum++  v0.8
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 inline double rand(double a = 0, double b = 1)
44 {
45  std::uniform_real_distribution<> ud(a, b);
46 
47 #ifndef _NO_THREAD_LOCAL
49 #else
50  return ud(RandomDevices::get_instance()._rng);
51  #endif // _NO_THREAD_LOCAL
52 }
53 
62 inline bigint rand(bigint a = std::numeric_limits<bigint>::min(),
63  bigint b = std::numeric_limits<bigint>::max())
64 {
65  std::uniform_int_distribution<bigint> uid(a, b);
66 
67 #ifndef _NO_THREAD_LOCAL
69 #else
70  return uid(RandomDevices::get_instance()._rng);
71  #endif // _NO_THREAD_LOCAL
72 }
73 
82 inline ubigint rand(ubigint a = std::numeric_limits<ubigint>::min(),
83  ubigint b = std::numeric_limits<ubigint>::max())
84 {
85  std::uniform_int_distribution<ubigint> uid(a, b);
86 
87 #ifndef _NO_THREAD_LOCAL
89 #else
90  return uid(RandomDevices::get_instance()._rng);
91  #endif // _NO_THREAD_LOCAL
92 }
93 
102 inline idx randidx(idx a = std::numeric_limits<idx>::min(),
103  idx b = std::numeric_limits<idx>::max())
104 {
105  std::uniform_int_distribution<idx> uid(a, b);
106 
107 #ifndef _NO_THREAD_LOCAL
108  return uid(RandomDevices::get_thread_local_instance()._rng);
109 #else
110  return uid(RandomDevices::get_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  if (rows == 0 || cols == 0)
161  throw Exception("qpp::rand()", Exception::Type::ZERO_SIZE);
162 
163  return dmat::Zero(rows, cols).unaryExpr(
164  [a, b](double)
165  {
166  return rand(a, b);
167  });
168 }
169 
191 template<>
192 inline cmat rand(idx rows, idx cols, double a, double b)
193 {
194  if (rows == 0 || cols == 0)
195  throw Exception("qpp::rand()", Exception::Type::ZERO_SIZE);
196 
197  return rand<dmat>(rows, cols, a, b).cast<cplx>() +
198  1_i * rand<dmat>(rows, cols, a, b).cast<cplx>();
199 }
200 
212 template<typename Derived>
213 Derived randn(idx rows, idx cols, double mean = 0,
214  double sigma = 1)
215 {
216  // silence -Wunused-parameter in clang++
217  (void) rows;
218  (void) cols;
219  (void) mean;
220  (void) sigma;
221  throw Exception("qpp::randn()", Exception::Type::UNDEFINED_TYPE);
222 }
223 
245 template<>
246 inline dmat randn(idx rows, idx cols,
247  double mean, double sigma)
248 {
249  if (rows == 0 || cols == 0)
250  throw Exception("qpp::randn()", Exception::Type::ZERO_SIZE);
251 
252  std::normal_distribution<> nd(mean, sigma);
253 
254  return dmat::Zero(rows, cols).unaryExpr(
255  [&nd](double)
256  {
257 #ifndef _NO_THREAD_LOCAL
259 #else
260  return nd(RandomDevices::get_instance()._rng);
261  #endif // _NO_THREAD_LOCAL
262  });
263 }
264 
286 template<>
287 inline cmat randn(idx rows, idx cols,
288  double mean, double sigma)
289 {
290  if (rows == 0 || cols == 0)
291  throw Exception("qpp::randn()", Exception::Type::ZERO_SIZE);
292 
293  return randn<dmat>(rows, cols, mean, sigma).cast<cplx>() +
294  1_i * randn<dmat>(rows, cols, mean, sigma).cast<cplx>();
295 }
296 
305 inline double randn(double mean = 0, double sigma = 1)
306 {
307  std::normal_distribution<> nd(mean, sigma);
308 
309 #ifndef _NO_THREAD_LOCAL
311 #else
312  return nd(RandomDevices::get_instance()._rng);
313  #endif // _NO_THREAD_LOCAL
314 }
315 
322 inline cmat randU(idx D)
323 // ~3 times slower than Toby Cubitt's MATLAB corresponding routine,
324 // because Eigen 3 QR algorithm is not parallelized
325 {
326  if (D == 0)
327  throw Exception("qpp::randU()", Exception::Type::DIMS_INVALID);
328 
329  cmat X(D, D);
330 
331  X = 1 / std::sqrt(2.) * randn < cmat > (D, D);
332  Eigen::HouseholderQR<cmat> qr(X);
333 
334  cmat Q = qr.householderQ();
335  // phase correction so that the resultant matrix is
336  // uniformly distributed according to the Haar measure
337 
338  Eigen::VectorXcd phases = (rand<dmat>(D, 1)).cast<cplx>();
339  for (idx i = 0; i < static_cast<idx>(phases.rows()); ++i)
340  phases(i) = std::exp(2 * pi * 1_i * phases(i));
341 
342  Q = Q * phases.asDiagonal();
343 
344  return Q;
345 }
346 
354 inline cmat randV(idx Din, idx Dout)
355 {
356  if (Din == 0 || Dout == 0 || Din > Dout)
357  throw Exception("qpp::randV()", Exception::Type::DIMS_INVALID);
358 
359  return randU(Dout).block(0, 0, Dout, Din);
360 }
361 
372 inline std::vector <cmat> randkraus(idx N, idx D)
373 {
374  if (N == 0)
375  throw Exception("qpp::randkraus()", Exception::Type::OUT_OF_RANGE);
376  if (D == 0)
377  throw Exception("qpp::randkraus()", Exception::Type::DIMS_INVALID);
378 
379  std::vector <cmat> result(N);
380  for (idx i = 0; i < N; ++i)
381  result[i] = cmat::Zero(D, D);
382 
383  cmat Fk(D, D);
384  cmat U = randU(N * D);
385 
386 #pragma omp parallel for collapse(3)
387  for (idx k = 0; k < N; ++k)
388  for (idx a = 0; a < D; ++a)
389  for (idx b = 0; b < D; ++b)
390  result[k](a, b) = U(a * N + k, b * N);
391 
392  return result;
393 }
394 
401 inline cmat randH(idx D)
402 {
403  if (D == 0)
404  throw Exception("qpp::randH()", Exception::Type::DIMS_INVALID);
405 
406  cmat H = 2 * rand<cmat>(D, D) - (1. + 1_i) * cmat::Ones(D, D);
407 
408  return H + adjoint(H);
409 }
410 
417 inline ket randket(idx D)
418 {
419  if (D == 0)
420  throw Exception("qpp::randket()", Exception::Type::DIMS_INVALID);
421 
422  /* slow
423  ket kt = ket::Ones(D);
424  ket result = static_cast<ket>(randU(D) * kt);
425 
426  return result;
427  */
428 
429  ket kt = randn < cmat > (D, 1);
430 
431  return kt / norm(kt);
432 }
433 
440 inline cmat randrho(idx D)
441 {
442  if (D == 0)
443  throw Exception("qpp::randrho()", Exception::Type::DIMS_INVALID);
444 
445  cmat result = 10 * randH(D);
446  result = result * adjoint(result);
447 
448  return result / trace(result);
449 }
450 
460 inline std::vector <idx> randperm(idx n)
461 {
462  if (n == 0)
463  throw Exception("qpp::randperm()", Exception::Type::PERM_INVALID);
464 
465  std::vector <idx> result(n);
466 
467  // fill in increasing order
468  std::iota(std::begin(result), std::end(result), 0);
469  // shuffle
470 #ifndef _NO_THREAD_LOCAL
471  std::shuffle(std::begin(result), std::end(result),
473 #else
474  std::shuffle(std::begin(result), std::end(result),
476  #endif // _NO_THREAD_LOCAL
477 
478  return result;
479 }
480 
481 } /* namespace qpp */
482 
483 #endif /* RANDOM_H_ */
static thread_local RandomDevices & get_thread_local_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:102
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:102
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:72
unsigned long long int ubigint
Non-negative big integer.
Definition: types.h:46
std::vector< cmat > randkraus(idx N, idx D)
Generates a set of random Kraus operators.
Definition: random.h:372
ket randket(idx D)
Generates a random normalized ket (pure state vector)
Definition: random.h:417
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:57
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:460
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:52
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:213
cmat randU(idx D)
Generates a random unitary matrix.
Definition: random.h:322
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:401
double rand(double a=0, double b=1)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:43
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:67
cmat randrho(idx D)
Generates a random density matrix.
Definition: random.h:440
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:354