random.cpp

/** @file
 * Testing utilities which generate random objects
 * independently of QuEST's internal generators. 
 * It is important that this file uses only its
 * internal RNG, and not C-library functions like
 * rand(), so that it cannot be interfered with by
 * external files calling e.g. srand().
 *
 * @author Tyson Jones
 */
 
#include "qvector.hpp"
#include "qmatrix.hpp"
#include "macros.hpp"
#include "linalg.hpp"
#include "lists.hpp"
#include "random.hpp"
#include "quest/include/quest.h"
 
#include <cmath>
#include <vector>
#include <tuple>
#include <random>
#include <algorithm>
 
using std::vector;
using std::tuple;
 
 
 
/*
 * RNG
 */
 
 
static std::mt19937 RNG;
 
 
void setRandomTestStateSeeds() {
    DEMAND( isQuESTEnvInit() );
 
    // generate a random seed from hardware rng
    std::random_device cspnrg;
    unsigned seed = cspnrg();
    
    // seed QuEST, which uses only the root node's seed
    setSeeds(&seed, 1);
 
    // broadcast root node seed to all nodes
    getSeeds(&seed);
 
    // seed RNG
    RNG.seed(seed);
}
 
 
 
/*
 * SCALAR
 */
 
 
qreal getRandomReal(qreal min, qreal maxExcl) {
    DEMAND( min < maxExcl );
 
    // advance RNG on every node, identically
    std::uniform_real_distribution<qreal> dist(min,maxExcl);
    qreal out = dist(RNG);
 
    // note despite the doc asserting maxExcl is exclusive, 
    // uniform_real_distribution() can indeed return it! In that
    // case, we substract machine-eps for caller integrity
    if (out >= maxExcl)
        out = std::nextafter(maxExcl, min);
 
    DEMAND( out >= min );
    DEMAND( out < maxExcl );
    return out;
}
 
 
qreal getRandomPhase() {
 
    // accuracy of PI does not matter here
    qreal pi = 3.14159265358979323846;
    return getRandomReal(-2*pi, 2*pi);
}
 
 
int getRandomInt(int min, int maxExcl) {
    DEMAND( maxExcl >= min );
 
    // permit this out of convenience
    // for some test generators
    if (min == maxExcl)
        return min;
 
    qreal r = std::floor(getRandomReal(min, maxExcl));
    int out = static_cast<int>(r);
 
    DEMAND( out >= min );
    DEMAND( out < maxExcl );
    return out;
}
 
 
qcomp getRandomComplex() {
    qreal re = getRandomReal(-1,1);
    qreal im = getRandomReal(-1,1);
    return qcomp(re, im);
}
 
 
qcomp getRandomUnitComplex() {
    return std::exp(1_i * getRandomPhase());
}
 
 
 
/*
 * LIST
 */
 
 
vector<int> getRandomInts(int min, int maxExcl, int len) {
    DEMAND( len >= 0 ); // permit empty
    DEMAND( min < maxExcl );
 
    vector<int> out(len);
 
    for (auto& x : out)
        x = getRandomInt(min, maxExcl);
 
    return out;
}
 
 
vector<int> getRandomOutcomes(int len) {
 
    int min = 0;
    int max = 1;
    return getRandomInts(min, max+1, len);
}
 
 
vector<int> getRandomSubRange(int start, int endExcl, int numElems) {
    DEMAND( endExcl >= start );
    DEMAND( numElems >= 1 );
    DEMAND( numElems <= endExcl - start );
 
    // shuffle entire range (advances RNG on every node, identically)
    vector<int> range = getRange(start, endExcl);
    std::shuffle(range.begin(), range.end(), RNG);
    
    // return first subrange
    return vector<int>(range.begin(), range.begin() + numElems);
}
 
 
vector<qreal> getRandomProbabilities(int numProbs) {
    
    vector<qreal> probs(numProbs, 0);
 
    // generate random unnormalised scalars
    for (auto& p : probs)
        p = getRandomReal(0, 1);
 
    // normalise
    qreal total = 0;
    for (auto& p : probs)
        total += p;
 
    for (auto& p : probs)
        p /= total;
 
    return probs;
}
 
 
listpair getRandomFixedNumCtrlsTargs(int numQubits, int numCtrls, int numTargs) {
 
    vector<int> targsCtrls = getRandomSubRange(0, numQubits, numTargs + numCtrls);
    vector<int> targs = getSublist(targsCtrls, 0, numTargs);
    vector<int> ctrls = getSublist(targsCtrls, numTargs, numCtrls);
 
    return tuple{ctrls,targs};
}
 
 
listtrio getRandomVariNumCtrlsStatesTargs(int numQubits, int minNumTargs, int maxNumTargsIncl) {
    DEMAND( minNumTargs <= maxNumTargsIncl );
    DEMAND( maxNumTargsIncl <= numQubits );
 
    int numTargs = getRandomInt(minNumTargs, maxNumTargsIncl+1);
 
    // numCtrls in [0, remainingNumQb]
    int minNumCtrls = 0;
    int maxNumCtrls = numQubits - numTargs;
    int numCtrls = getRandomInt(minNumCtrls, maxNumCtrls+1);
 
    // distribute qubits randomly
    auto [ctrls,targs] = getRandomFixedNumCtrlsTargs(numQubits, numCtrls, numTargs);
    vector<int> states = getRandomInts(0, 2, numCtrls);
 
    return tuple{ctrls,states,targs};
}
 
 
 
/*
 * VECTOR
 */
 
 
qvector getRandomVector(size_t dim) { 
 
    qvector vec = getZeroVector(dim);
 
    for (auto& elem : vec)
        elem = getRandomComplex();
        
    return vec;
}
 
 
vector<qvector> getRandomOrthonormalVectors(size_t dim, int numVecs) {
    DEMAND( dim >= 1 );
    DEMAND( numVecs >= 1);
    
    vector<qvector> vecs(numVecs);
 
    // produce each vector in-turn
    for (int n=0; n<numVecs; n++) {
 
        // from a random vector
        vecs[n] = getRandomVector(dim);
 
        // orthogonalise by substracting projections of existing vectors
        for (int m=0; m<n; m++)
            vecs[n] -= vecs[m] * getInnerProduct(vecs[m], vecs[n]);
 
        // then re-normalise
        vecs[n] = getNormalised(vecs[n]);
    }
 
    return vecs;
}
 
 
 
/*
 * MATRIX
 */
 
 
qmatrix getRandomNonSquareMatrix(size_t numRows, size_t numCols) {
 
    // this function is DANGEROUS; it produces a
    // non-square matrix, whereas most test utilities
    // assume qmatrix is square. It should ergo be
    // used very cautiously!
 
    qmatrix out = qmatrix(numRows, qvector(numCols));
 
    for (auto& row : out)
        for (auto& elem : row)
            elem = getRandomComplex();
 
    return out;
}
 
 
qmatrix getRandomMatrix(size_t dim) {
 
    return getRandomNonSquareMatrix(dim, dim);
}
 
 
qmatrix getRandomDiagonalMatrix(size_t dim) {
 
    qmatrix out = getZeroMatrix(dim);
 
    for (size_t i=0; i<dim; i++)
        out[i][i] = getRandomComplex();
 
    return out;
}
 
 
 
/*
 * STATES
 */
 
 
qvector getRandomStateVector(int numQb) {
 
    return getNormalised(getRandomVector(getPow2(numQb)));
}
 
 
vector<qvector> getRandomOrthonormalStateVectors(int numQb, int numStates) {
 
    return getRandomOrthonormalVectors(getPow2(numQb), numStates);
}
 
 
qmatrix getRandomDensityMatrix(int numQb) {
    DEMAND( numQb > 0 );
    
    // generate random probabilities to weight random pure states
    int dim = getPow2(numQb);
    vector<qreal> probs = getRandomProbabilities(dim);
    
    // add random pure states
    qmatrix dens = getZeroMatrix(dim);
    for (int i=0; i<dim; i++) {
        qvector pure = getRandomStateVector(numQb);
        dens += probs[i] * getOuterProduct(pure, pure);
    }
    
    return dens;
}
 
 
qmatrix getRandomPureDensityMatrix(int numQb) {
 
    qvector vec = getRandomStateVector(numQb);
    qmatrix mat = getOuterProduct(vec, vec);
    return mat;
}
 
 
void setToRandomState(qvector& state) {
    state = getRandomStateVector(getLog2(state.size()));
}
void setToRandomState(qmatrix& state) {
    state = getRandomDensityMatrix(getLog2(state.size()));
}
 
 
 
/*
 * OPERATORS
 */
 
 
qmatrix getRandomUnitary(int numQb) {
    DEMAND( numQb >= 1 );
 
    // create Z ~ random complex matrix (distribution not too important)
    size_t dim = getPow2(numQb);
    qmatrix matrZ = getRandomMatrix(dim);
    qmatrix matrZT = getTranspose(matrZ);
 
    // create Z = Q R (via QR decomposition) ...
    qmatrix matrQT = getOrthonormalisedRows(matrZ);
    qmatrix matrQ = getTranspose(matrQT);
    qmatrix matrR = getZeroMatrix(dim);
 
    // ... where R_rc = (columm c of Z) . (column r of Q) = (row c of ZT) . (row r of QT) = <r|c>
    for (size_t r=0; r<dim; r++)
        for (size_t c=r; c<dim; c++)
            matrR[r][c] = getInnerProduct(matrQT[r], matrZT[c]);
 
    // create D = normalised diagonal of R
    qmatrix matrD = getZeroMatrix(dim);
    for (size_t i=0; i<dim; i++)
        matrD[i][i] = matrR[i][i] / std::abs(matrR[i][i]);
 
    // create U = Q D
    qmatrix matrU = matrQ * matrD;
 
    DEMAND( isApproxUnitary(matrU) );
    return matrU;
}
 
 
qmatrix getRandomDiagonalUnitary(int numQb) {
    DEMAND( numQb >= 1 );
 
    qmatrix matr = getZeroMatrix(getPow2(numQb));
 
    // unitary diagonals have unit complex scalars
    for (size_t i=0; i<matr.size(); i++)
        matr[i][i] = getRandomUnitComplex();
 
    return matr;
}
 
 
qmatrix getRandomDiagonalHermitian(int numQb) {
    DEMAND( numQb >= 1 );
 
    qmatrix out = getZeroMatrix(getPow2(numQb));
 
    // Hermitian diagonals are real
    for (size_t i=0; i<out.size(); i++)
        out[i][i] = getRandomReal(-10, 10);
 
    return out;
}
 
 
vector<qmatrix> getRandomKrausMap(int numQb, int numOps) {
    DEMAND( numOps >= 1 );
 
    // generate random unitaries
    vector<qmatrix> ops(numOps);
    for (auto& u : ops)
        u = getRandomUnitary(numQb);
 
    // generate random weights
    vector<qreal> weights(numOps);
    for (auto& w : weights)
        w = getRandomReal(0, 1);
        
    // normalise random weights
    qreal sum = 0;
    for (auto& w : weights)
        sum += w;
    for (auto& w : weights)
        w = std::sqrt(w/sum);
        
    // normalise unitaries according to weights
    for (int i=0; i<numOps; i++)
        ops[i] *= weights[i];
 
    DEMAND( isApproxCPTP(ops) );
    return ops;
}
 
 
 
/*
 * PAULIS
 */
 
 
PauliStr getRandomPauliStr(int numQubits) {
 
    std::string paulis = "";
    for (int i=0; i<numQubits; i++)
        paulis += "IXYZ"[getRandomInt(0,4)];
 
    return getPauliStr(paulis);
}
 
 
PauliStr getRandomPauliStr(vector<int> targs) {
 
    std::string paulis = "";
    for (size_t i=0; i<targs.size(); i++)
        paulis += "IXYZ"[getRandomInt(0,4)];
 
    return getPauliStr(paulis, targs);
}
 
 
PauliStr getRandomDiagPauliStr(int numQubits) {
 
    std::string paulis = "";
    for (int i=0; i<numQubits; i++)
        paulis += "IZ"[getRandomInt(0,2)];
 
    return getPauliStr(paulis);
}
 
 
PauliStrSum createRandomNonHermitianPauliStrSum(int numQubits, int numTerms) {
 
    vector<PauliStr> strings(numTerms);
    for (auto& str : strings)
        str = getRandomPauliStr(numQubits);
 
    vector<qcomp> coeffs(numTerms);
    for (auto& c : coeffs)
        c = getRandomComplex();
 
    return createPauliStrSum(strings, coeffs);
}
 
 
PauliStrSum createRandomPauliStrSum(int numQubits, int numTerms) {
 
    PauliStrSum out = createRandomNonHermitianPauliStrSum(numQubits, numTerms);
 
    for (qindex i=0; i<numTerms; i++)
        out.coeffs[i] = std::real(out.coeffs[i]);
 
    return out;
}