random.h

Go to the documentation of this file.

/*
 * This file is part of Quantum++.
 *
 * MIT License
 *
 * Copyright (c) 2013 - 2020 Vlad Gheorghiu.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef RANDOM_H_
#define RANDOM_H_

namespace qpp {
inline double rand(double a, double b) {
    // EXCEPTION CHECKS

    if (a >= b)
        throw exception::OutOfRange("qpp::rand()");
    // END EXCEPTION CHECKS

    std::uniform_real_distribution<> ud(a, b);
    auto& gen =
#ifdef NO_THREAD_LOCAL_
        RandomDevices::get_instance().get_prng();
#else
        RandomDevices::get_thread_local_instance().get_prng();
#endif

    return ud(gen);
}

inline bigint rand(bigint a, bigint b) {
    // EXCEPTION CHECKS

    if (a > b)
        throw exception::OutOfRange("qpp::rand()");
    // END EXCEPTION CHECKS

    std::uniform_int_distribution<bigint> uid(a, b);

    auto& gen =
#ifdef NO_THREAD_LOCAL_
        RandomDevices::get_instance().get_prng();
#else
        RandomDevices::get_thread_local_instance().get_prng();
#endif

    return uid(gen);
}

inline idx randidx(idx a = std::numeric_limits<idx>::min(),
                   idx b = std::numeric_limits<idx>::max()) {
    // EXCEPTION CHECKS

    if (a > b)
        throw exception::OutOfRange("qpp::randidx()");
    // END EXCEPTION CHECKS

    std::uniform_int_distribution<idx> uid(a, b);
    auto& gen =
#ifdef NO_THREAD_LOCAL_
        RandomDevices::get_instance().get_prng();
#else
        RandomDevices::get_thread_local_instance().get_prng();
#endif

    return uid(gen);
}

template <typename Derived>
Derived rand(idx rows QPP_UNUSED_, idx cols QPP_UNUSED_,
             double a QPP_UNUSED_ = 0, double b QPP_UNUSED_ = 1) {
    throw exception::UndefinedType("qpp::rand()");
}

template <>
inline dmat rand(idx rows, idx cols, double a, double b) {
    // EXCEPTION CHECKS

    if (rows == 0 || cols == 0)
        throw exception::ZeroSize("qpp::rand()");
    if (a >= b)
        throw exception::OutOfRange("qpp::rand()");
    // END EXCEPTION CHECKS

    return dmat::Zero(rows, cols).unaryExpr([a, b](double) {
        return rand(a, b);
    });
}

template <>
inline cmat rand(idx rows, idx cols, double a, double b) {
    // EXCEPTION CHECKS

    if (rows == 0 || cols == 0)
        throw exception::ZeroSize("qpp::rand()");
    if (a >= b)
        throw exception::OutOfRange("qpp::rand()");
    // END EXCEPTION CHECKS

    return rand<dmat>(rows, cols, a, b).cast<cplx>() +
           1_i * rand<dmat>(rows, cols, a, b).cast<cplx>();
}

template <typename Derived>
Derived randn(idx rows QPP_UNUSED_, idx cols QPP_UNUSED_,
              double mean QPP_UNUSED_ = 0, double sigma QPP_UNUSED_ = 1) {
    throw exception::UndefinedType("qpp::randn()");
}

template <>
inline dmat randn(idx rows, idx cols, double mean, double sigma) {
    // EXCEPTION CHECKS

    if (rows == 0 || cols == 0)
        throw exception::ZeroSize("qpp::randn()");
    // END EXCEPTION CHECKS

    std::normal_distribution<> nd(mean, sigma);
    auto& gen =
#ifdef NO_THREAD_LOCAL_
        RandomDevices::get_instance().get_prng();
#else
        RandomDevices::get_thread_local_instance().get_prng();
#endif

    return dmat::Zero(rows, cols).unaryExpr([&nd, &gen](double) {
        return nd(gen);
    });
}

template <>
inline cmat randn(idx rows, idx cols, double mean, double sigma) {
    // EXCEPTION CHECKS

    if (rows == 0 || cols == 0)
        throw exception::ZeroSize("qpp::randn()");
    // END EXCEPTION CHECKS

    return randn<dmat>(rows, cols, mean, sigma).cast<cplx>() +
           1_i * randn<dmat>(rows, cols, mean, sigma).cast<cplx>();
}

inline double randn(double mean = 0, double sigma = 1) {
    std::normal_distribution<> nd(mean, sigma);
    auto& gen =
#ifdef NO_THREAD_LOCAL_
        RandomDevices::get_instance().get_prng();
#else
        RandomDevices::get_thread_local_instance().get_prng();
#endif

    return nd(gen);
}

inline cmat randU(idx D = 2)
// ~3 times slower than Toby Cubitt's MATLAB corresponding routine,
// because Eigen 3 QR algorithm is not parallelized
{
    // EXCEPTION CHECKS

    if (D == 0)
        throw exception::DimsInvalid("qpp::randU()");
    // END EXCEPTION CHECKS

    cmat X = 1 / std::sqrt(2.) * randn<cmat>(D, D);
    Eigen::HouseholderQR<cmat> qr(X);

    cmat Q = qr.householderQ();
    // phase correction so that the resultant matrix is
    // uniformly distributed according to the Haar measure

    Eigen::VectorXcd phases = (rand<dmat>(D, 1)).cast<cplx>();
    for (idx i = 0; i < static_cast<idx>(phases.rows()); ++i)
        phases(i) = std::exp(2 * pi * 1_i * phases(i));

    Q = Q * phases.asDiagonal();

    return Q;
}

inline cmat randV(idx Din, idx Dout) {
    // EXCEPTION CHECKS

    if (Din == 0 || Dout == 0 || Din > Dout)
        throw exception::DimsInvalid("qpp::randV()");
    // END EXCEPTION CHECKS

    return randU(Dout).block(0, 0, Dout, Din);
}

inline std::vector<cmat> randkraus(idx N, idx D = 2) {
    // EXCEPTION CHECKS

    if (N == 0)
        throw exception::OutOfRange("qpp::randkraus()");
    if (D == 0)
        throw exception::DimsInvalid("qpp::randkraus()");
    // END EXCEPTION CHECKS

    std::vector<cmat> result(N);
    for (idx i = 0; i < N; ++i)
        result[i] = cmat::Zero(D, D);

    cmat Fk(D, D);
    cmat U = randU(N * D);

#ifdef WITH_OPENMP_
#pragma omp parallel for collapse(3)
#endif // WITH_OPENMP_
    for (idx k = 0; k < N; ++k)
        for (idx a = 0; a < D; ++a)
            for (idx b = 0; b < D; ++b)
                result[k](a, b) = U(a * N + k, b * N);

    return result;
}

inline cmat randH(idx D = 2) {
    // EXCEPTION CHECKS

    if (D == 0)
        throw exception::DimsInvalid("qpp::randH()");
    // END EXCEPTION CHECKS

    cmat H = 2 * rand<cmat>(D, D) - (1. + 1_i) * cmat::Ones(D, D);

    return H + H.adjoint();
}

inline ket randket(idx D = 2) {
    // EXCEPTION CHECKS

    if (D == 0)
        throw exception::DimsInvalid("qpp::randket()");
    // END EXCEPTION CHECKS

    /* slow
     ket kt = ket::Ones(D);
     ket result = static_cast<ket>(randU(D) * kt);

     return result;
     */

    ket kt = randn<cmat>(D, 1);

    return kt / kt.norm();
}

inline cmat randrho(idx D = 2) {
    // EXCEPTION CHECKS

    if (D == 0)
        throw exception::DimsInvalid("qpp::randrho()");
    // END EXCEPTION CHECKS

    cmat result = 10 * randH(D);
    result = result * result.adjoint();

    return result / result.trace();
}

inline std::vector<idx> randperm(idx N) {
    // EXCEPTION CHECKS

    if (N == 0)
        throw exception::PermInvalid("qpp::randperm()");
    // END EXCEPTION CHECKS

    std::vector<idx> result(N);

    // fill in increasing order
    std::iota(std::begin(result), std::end(result), 0);

    // shuffle
    auto& gen =
#ifdef NO_THREAD_LOCAL_
        RandomDevices::get_instance().get_prng();
#else
        RandomDevices::get_thread_local_instance().get_prng();
#endif
    std::shuffle(std::begin(result), std::end(result), gen);

    return result;
}

inline std::vector<double> randprob(idx N) {
    // EXCEPTION CHECKS

    if (N == 0)
        throw exception::OutOfRange("qpp::randprob()");
    // END EXCEPTION CHECKS

    std::vector<double> result(N);

    // generate
    std::exponential_distribution<> ed(1);
    auto& gen =
#ifdef NO_THREAD_LOCAL_
        RandomDevices::get_instance().get_prng();
#else
        RandomDevices::get_thread_local_instance().get_prng();
#endif
    for (idx i = 0; i < N; ++i)
        result[i] = ed(gen);

    // normalize
    double sumprob = std::accumulate(std::begin(result), std::end(result), 0.0);
    for (idx i = 0; i < N; ++i)
        result[i] /= sumprob;

    return result;
}

} /* namespace qpp */

#endif /* RANDOM_H_ */