The Quantum Exact Simulation Toolkit v4.2.0
Loading...
Searching...
No Matches
test_calculations.cpp
1/** @file
2 * Ported tests of the deprecated QuEST v3 interface,
3 * unit testing the "calculations" module.
4 *
5 * This file should be excluded from doxygen parsing so
6 * as not to conflict with the doc of the v4 unit tests.
7 *
8 * @author Tyson Jones
9 * @author Oliver Thomson Brown (ported to Catch2 v3)
10 * @author Ali Rezaei (tested porting to QuEST v4)
11 */
12
13#include <catch2/catch_test_macros.hpp>
14#include <catch2/catch_approx.hpp>
15#include <catch2/matchers/catch_matchers_string.hpp>
16#include <catch2/generators/catch_generators_range.hpp>
17
18#include "quest.h"
19#include "test_utilities.hpp"
20
21/* allows concise use of ContainsSubstring in catch's REQUIRE_THROWS_WITH */
22using Catch::Matchers::ContainsSubstring;
23using Catch::Approx;
24
25
26/** @sa calcDensityInnerProduct
27 * @ingroup deprecatedtests
28 * @author Tyson Jones
29 */
30TEST_CASE( "calcDensityInnerProduct", "[calculations]" ) {
31
32 Qureg mat1 = createForcedDensityQureg(NUM_QUBITS);
33 Qureg mat2 = createForcedDensityQureg(NUM_QUBITS);
34
35 SECTION( "correctness" ) {
36
37 // repeat these random tests 10 times
38 GENERATE( range(0,10) );
39
40 SECTION( "density-matrix" ) {
41
42 SECTION( "pure" ) {
43
44 // mat1 = |r1><r1|
45 QVector r1 = getRandomStateVector(NUM_QUBITS);
46 toQureg(mat1, getKetBra(r1,r1));
47
48 // mat2 = |r2><r2|
49 QVector r2 = getRandomStateVector(NUM_QUBITS);
50 toQureg(mat2, getKetBra(r2,r2));
51
52 // prod( |r1><r1|, |r2><r2| ) = |<r1|r2>|^2
53 qcomp prod = 0;
54 for (size_t i=0; i<r1.size(); i++)
55 prod += conj(r1[i]) * r2[i];
56 qreal densProd = pow(abs(prod),2);
57
58 REQUIRE( real(calcDensityInnerProduct(mat1,mat2)) == Approx(densProd).margin(100 * REAL_EPS) );
59 }
60 SECTION( "mixed" ) {
61
62 QMatrix ref1 = getRandomDensityMatrix(NUM_QUBITS);
63 QMatrix ref2 = getRandomDensityMatrix(NUM_QUBITS);
64 toQureg(mat1, ref1);
65 toQureg(mat2, ref2);
66
67 // prod(mat1, mat2) = sum_{ij} conj(mat1_{ij}) * mat2_{ij}
68 qcomp refProd = 0;
69 for (size_t i=0; i<ref1.size(); i++)
70 for (size_t j=0; j<ref1.size(); j++)
71 refProd += conj(ref1[i][j]) * ref2[i][j];
72 REQUIRE( imag(refProd) == Approx(0).margin(REAL_EPS) );
73
74 REQUIRE( real(calcDensityInnerProduct(mat1,mat2)) == Approx(real(refProd)).margin(100 * REAL_EPS) );
75
76 // should be invariant under ordering
77 REQUIRE( real(calcDensityInnerProduct(mat1,mat2)) == Approx(real(calcDensityInnerProduct(mat2,mat1))).margin(100 * REAL_EPS) );
78 }
79 SECTION( "unnormalised" ) {
80
81 // set both to random (non-Hermitian) complex matrices
82 QMatrix ref1 = getRandomQMatrix(1<<NUM_QUBITS);
83 QMatrix ref2 = getRandomQMatrix(1<<NUM_QUBITS);
84 toQureg(mat1, ref1);
85 toQureg(mat2, ref2);
86
87 // prod(mat1, mat2) = real(sum_{ij} conj(mat1_{ij}) * mat2_{ij})
88 qcomp refProd = 0;
89 for (size_t i=0; i<ref1.size(); i++)
90 for (size_t j=0; j<ref1.size(); j++)
91 refProd += conj(ref1[i][j]) * ref2[i][j];
92
93 REQUIRE( real(calcDensityInnerProduct(mat1,mat2)) == Approx(real(refProd)).margin(100 * REAL_EPS) );
94 }
95 }
96 }
97 SECTION( "input validation" ) {
98
99 SECTION( "dimensions" ) {
100
101 Qureg mat3 = createDensityQureg(NUM_QUBITS + 1);
102 REQUIRE_THROWS_WITH( calcDensityInnerProduct(mat1,mat3), ContainsSubstring("differing numbers of qubits") );
103 destroyQureg(mat3);
104 }
105
106 // in v4, calcDensityInnerProduct() redirects to calcInnerProduct() which is
107 // valid for both statevectors and density matrices (and combinations thereof!)
108
109 // SECTION( "state-vectors" ) {
110
111 // Qureg vec = createForcedQureg(NUM_QUBITS);
112
113 // REQUIRE_THROWS_WITH( calcDensityInnerProduct(mat1,vec), ContainsSubstring("valid only for density matrices") );
114 // REQUIRE_THROWS_WITH( calcDensityInnerProduct(vec,mat1), ContainsSubstring("valid only for density matrices") );
115 // REQUIRE_THROWS_WITH( calcDensityInnerProduct(vec,vec), ContainsSubstring("valid only for density matrices") );
116
117 // destroyQureg(vec);
118 // }
119 }
120 destroyQureg(mat1);
121 destroyQureg(mat2);
122}
123
124
125
126/** @sa calcExpecDiagonalOp
127 * @ingroup deprecatedtests
128 * @author Tyson Jones
129 */
130TEST_CASE( "calcExpecDiagonalOp", "[calculations]" ) {
131
132 Qureg vec = createForcedQureg(NUM_QUBITS);
133 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
134 initDebugState(vec);
135 initDebugState(mat);
136 QVector vecRef = toQVector(vec);
137 QMatrix matRef = toQMatrix(mat);
138
139 SECTION( "correctness" ) {
140
141 // try 10 random operators
142 GENERATE( range(0,10) );
143
144 // make a totally random (non-Hermitian) diagonal oeprator
145 DiagonalOp op = createDiagonalOp(NUM_QUBITS, getQuESTEnv());
146 for (long long int i=0; i<op.numElemsPerNode; i++)
147 op.cpuElems[i] = getRandomComplex();
148 syncDiagonalOp(op);
149
150 SECTION( "state-vector" ) {
151
152 /* calcExpecDiagOp calculates <qureg|diag|qureg> */
153
154 QVector sumRef = toQMatrix(op) * vecRef;
155 qcomp prod = 0;
156 for (size_t i=0; i<vecRef.size(); i++)
157 prod += conj(vecRef[i]) * sumRef[i];
158
159 qcomp res = calcExpecDiagonalOp(vec, op);
160 REQUIRE( real(res) == Approx(real(prod)).margin(REAL_EPS) );
161 REQUIRE( imag(res) == Approx(imag(prod)).margin(REAL_EPS) );
162 }
163 SECTION( "density-matrix" ) {
164
165 /* calcExpecDiagOp calculates Trace( diag * qureg ) */
166 matRef = toQMatrix(op) * matRef;
167 qcomp tr = 0;
168 for (size_t i=0; i<matRef.size(); i++)
169 tr += matRef[i][i];
170
171 qcomp res = calcExpecDiagonalOp(mat, op);
172 REQUIRE( real(res) == Approx(real(tr)).margin(100*REAL_EPS) );
173 REQUIRE( imag(res) == Approx(imag(tr)).margin(100*REAL_EPS) );
174 }
175
176 destroyDiagonalOp(op, getQuESTEnv());
177 }
178 SECTION( "input validation" ) {
179
180 SECTION( "mismatching size" ) {
181
182 DiagonalOp op = createDiagonalOp(NUM_QUBITS + 1, getQuESTEnv());
183
184 REQUIRE_THROWS_WITH( calcExpecDiagonalOp(vec, op), ContainsSubstring("different number of qubits"));
185 REQUIRE_THROWS_WITH( calcExpecDiagonalOp(mat, op), ContainsSubstring("different number of qubits"));
186
187 destroyDiagonalOp(op, getQuESTEnv());
188 }
189 }
190 destroyQureg(vec);
191 destroyQureg(mat);
192}
193
194
195
196// calcExpecPauliHamil removed because PauliHamil is deprecated,
197// and replacement PauliStrSum presently has no compatible constructor
198
199 // /** @sa calcExpecPauliHamil
200 // * @ingroup deprecatedtests
201 // * @author Tyson Jones
202 // */
203 // TEST_CASE( "calcExpecPauliHamil", "[calculations]" ) {
204
205 // Qureg vec = createForcedQureg(NUM_QUBITS);
206 // Qureg mat = createForcedDensityQureg(NUM_QUBITS);
207 // initDebugState(vec);
208 // initDebugState(mat);
209 // QVector vecRef = toQVector(vec);
210 // QMatrix matRef = toQMatrix(mat);
211
212 // Qureg vecWork = createForcedQureg(NUM_QUBITS);
213 // Qureg matWork = createForcedDensityQureg(NUM_QUBITS);
214
215 // SECTION( "correctness" ) {
216
217 // /* it's too expensive to try every possible Pauli configuration, so
218 // * we'll try 10 random codes, and for each, random coefficients
219 // */
220 // GENERATE( range(0,10) );
221
222 // int numTerms = GENERATE( 1, 2, 10, 15 );
223 // PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
224 // setRandomPauliSum(hamil);
225 // QMatrix refHamil = toQMatrix(hamil);
226
227 // SECTION( "state-vector" ) {
228
229 // /* calcExpecPauliHamil calculates <qureg|pauliHum|qureg> */
230
231 // QVector sumRef = refHamil * vecRef;
232 // qcomp prod = 0;
233 // for (size_t i=0; i<vecRef.size(); i++)
234 // prod += conj(vecRef[i]) * sumRef[i];
235 // REQUIRE( imag(prod) == Approx(0).margin(10*REAL_EPS) );
236
237 // qreal res = calcExpecPauliHamil(vec, hamil, vecWork);
238 // REQUIRE( res == Approx(real(prod)).margin(10*REAL_EPS) );
239 // }
240 // SECTION( "density-matrix" ) {
241
242 // /* calcExpecPauliHamil calculates Trace( pauliHamil * qureg ) */
243 // matRef = refHamil * matRef;
244 // qreal tr = 0;
245 // for (size_t i=0; i<matRef.size(); i++)
246 // tr += real(matRef[i][i]);
247 // // (get real, since we start in a non-Hermitian state, hence diagonal isn't real)
248
249 // qreal res = calcExpecPauliHamil(mat, hamil, matWork);
250 // REQUIRE( res == Approx(tr).margin(1E2*REAL_EPS) );
251 // }
252
253 // destroyPauliHamil(hamil);
254 // }
255 // SECTION( "input validation" ) {
256
257 // SECTION( "pauli codes" ) {
258
259 // int numTerms = 3;
260 // PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
261
262 // // make one pauli code wrong
263 // hamil.pauliCodes[GENERATE_COPY( range(0,numTerms*NUM_QUBITS) )] = (pauliOpType) GENERATE( -1, 4 );
264 // REQUIRE_THROWS_WITH( calcExpecPauliHamil(vec, hamil, vecWork), ContainsSubstring("Invalid Pauli code") );
265
266 // destroyPauliHamil(hamil);
267 // }
268 // SECTION( "workspace type" ) {
269
270 // int numTerms = 1;
271 // PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
272
273 // REQUIRE_THROWS_WITH( calcExpecPauliHamil(vec, hamil, mat), ContainsSubstring("Registers must both be state-vectors or both be density matrices") );
274 // REQUIRE_THROWS_WITH( calcExpecPauliHamil(mat, hamil, vec), ContainsSubstring("Registers must both be state-vectors or both be density matrices") );
275
276 // destroyPauliHamil(hamil);
277 // }
278 // SECTION( "workspace dimensions" ) {
279
280 // int numTerms = 1;
281 // PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
282
283 // Qureg vec2 = createQureg(NUM_QUBITS + 1);
284 // REQUIRE_THROWS_WITH( calcExpecPauliHamil(vec, hamil, vec2), ContainsSubstring("Dimensions") && ContainsSubstring("don't match") );
285 // destroyQureg(vec2);
286
287 // Qureg mat2 = createDensityQureg(NUM_QUBITS + 1);
288 // REQUIRE_THROWS_WITH( calcExpecPauliHamil(mat, hamil, mat2), ContainsSubstring("Dimensions") && ContainsSubstring("don't match") );
289 // destroyQureg(mat2);
290
291 // destroyPauliHamil(hamil);
292 // }
293 // SECTION( "matching hamiltonian qubits" ) {
294
295 // int numTerms = 1;
296 // PauliHamil hamil = createPauliHamil(NUM_QUBITS + 1, numTerms);
297
298 // REQUIRE_THROWS_WITH( calcExpecPauliHamil(vec, hamil, vecWork), ContainsSubstring("same number of qubits") );
299 // REQUIRE_THROWS_WITH( calcExpecPauliHamil(mat, hamil, matWork), ContainsSubstring("same number of qubits") );
300
301 // destroyPauliHamil(hamil);
302 // }
303 // }
304 // destroyQureg(vec);
305 // destroyQureg(mat);
306 // destroyQureg(vecWork);
307 // destroyQureg(matWork);
308 // }
309
310
311
312/** @sa calcExpecPauliProd
313 * @ingroup deprecatedtests
314 * @author Tyson Jones
315 */
316TEST_CASE( "calcExpecPauliProd", "[calculations]" ) {
317
318 QuESTEnv env = getQuESTEnv();
319 Qureg vec = createCustomQureg(NUM_QUBITS, 0, env.isDistributed, env.isGpuAccelerated, env.isMultithreaded);
320 Qureg mat = createCustomQureg(NUM_QUBITS, 1, env.isDistributed, env.isGpuAccelerated, env.isMultithreaded);
321
322 initDebugState(vec);
323 initDebugState(mat);
324 QVector vecRef = toQVector(vec);
325 QMatrix matRef = toQMatrix(mat);
326
327 Qureg vecWork = createForcedQureg(NUM_QUBITS);
328 Qureg matWork = createForcedDensityQureg(NUM_QUBITS);
329
330 SECTION( "correctness" ) {
331
332 int numTargs = GENERATE( range(1,NUM_QUBITS+1) );
333 int* targs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numTargs) );
334
335 /* it's too expensive to try ALL Pauli sequences, via
336 * pauliOpType* paulis = GENERATE_COPY( pauliseqs(numTargs) );.
337 * Furthermore, take(10, pauliseqs(numTargs)) will try the same pauli codes.
338 * Hence, we instead opt to repeatedlyrandomly generate pauliseqs
339 */
340 GENERATE( range(0,10) ); // gen 10 random pauli-codes for every targs
341 vector<pauliOpType> paulis(numTargs);
342 for (int i=0; i<numTargs; i++)
343 paulis[i] = (pauliOpType) getRandomInt(0,4);
344
345 // produce a numTargs-big matrix 'pauliProd' by pauli-matrix tensoring
346 QMatrix iMatr{{1,0},{0,1}};
347 QMatrix xMatr{{0,1},{1,0}};
348 QMatrix yMatr{{0,-qcomp(0,1)},{qcomp(0,1),0}};
349 QMatrix zMatr{{1,0},{0,-1}};
350 QMatrix pauliProd{{1}};
351 for (int i=0; i<numTargs; i++) {
352 QMatrix fac;
353 if (paulis[i] == PAULI_I) fac = iMatr;
354 if (paulis[i] == PAULI_X) fac = xMatr;
355 if (paulis[i] == PAULI_Y) fac = yMatr;
356 if (paulis[i] == PAULI_Z) fac = zMatr;
357 pauliProd = getKroneckerProduct(fac, pauliProd);
358 }
359
360 SECTION( "state-vector" ) {
361
362 /* calcExpecPauliProd calculates <qureg|pauliProd|qureg> */
363
364 QVector prodRef = vecRef;
365 applyReferenceOp(prodRef, targs, numTargs, pauliProd);
366 qcomp prod = 0;
367 for (size_t i=0; i<vecRef.size(); i++)
368 prod += conj(vecRef[i]) * prodRef[i];
369 REQUIRE( imag(prod) == Approx(0).margin(REAL_EPS) );
370
371 qreal res = calcExpecPauliProd(vec, targs, paulis.data(), numTargs, vecWork);
372 REQUIRE( res == Approx(real(prod)).margin(REAL_EPS) );
373 }
374 SECTION( "density-matrix" ) {
375
376 /* calcExpecPauliProd calculates Trace( pauliProd * qureg ) */
377
378 // produce (pauliProd * mat)
379 QMatrix fullOp = getFullOperatorMatrix(NULL, 0, targs, numTargs, pauliProd, NUM_QUBITS);
380 matRef = fullOp * matRef;
381
382 // compute real(trace(pauliProd * mat))
383 qreal tr = 0;
384 for (size_t i=0; i<matRef.size(); i++)
385 tr += real(matRef[i][i]);
386 // (get real, since we start in a non-Hermitian state, hence diagonal isn't real)
387
388 // disable validation during call, because result is non-real and will upset post-check
389 setValidationOff();
390 qreal res = calcExpecPauliProd(mat, targs, paulis.data(), numTargs, matWork);
391 setValidationOn();
392
393 REQUIRE( res == Approx(tr).margin(10*REAL_EPS) );
394 }
395 }
396 SECTION( "input validation" ) {
397
398 pauliOpType pauliOpsToAvoidDeprecSegFault[100];
399 for (int i=0; i<100; i++)
400 pauliOpsToAvoidDeprecSegFault[i] = PAULI_X;
401
402 SECTION( "number of targets" ) {
403
404 int targs[NUM_QUBITS+1];
405 for (int i=0; i<NUM_QUBITS+1; i++)
406 targs[i] = i;
407
408 // too few
409 int numTargs = GENERATE( -1, 0 );
410 REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, pauliOpsToAvoidDeprecSegFault, numTargs, vecWork), ContainsSubstring("Invalid number of Paulis") );
411
412 // too many
413 numTargs = NUM_QUBITS + 1;
414 REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, pauliOpsToAvoidDeprecSegFault, numTargs, vecWork), ContainsSubstring("highest-index non-identity Pauli operator") && ContainsSubstring("exceeds the maximum target") );
415 }
416 SECTION( "target indices" ) {
417
418 int numTargs = 3;
419 int targs[3] = {0, 1, 2};
420
421 int badInd = GENERATE( range(0,3) );
422
423 // make one index too small
424 targs[badInd] = -1;
425 REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, pauliOpsToAvoidDeprecSegFault, numTargs, vecWork), ContainsSubstring("Pauli indices must be non-negative") );
426
427 // make one index too big
428 targs[badInd] = NUM_QUBITS;
429 REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, pauliOpsToAvoidDeprecSegFault, numTargs, vecWork), ContainsSubstring("highest-index non-identity Pauli operator") );
430
431 // make one index WAY too big
432 targs[badInd] = 65;
433 REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, pauliOpsToAvoidDeprecSegFault, numTargs, vecWork), ContainsSubstring("exceed the maximum number of representable Pauli operators") );
434 }
435 SECTION( "repetition in targets" ) {
436
437 int numTargs = 3;
438 int targs[3] = {0, 1, 1};
439 REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, pauliOpsToAvoidDeprecSegFault, numTargs, vecWork), ContainsSubstring("Indices must be unique") );
440 }
441 SECTION( "pauli codes" ) {
442
443 int numTargs = 3;
444 int targs[3] = {0, 1, 2};
445 pauliOpType codes[3] = {PAULI_X, PAULI_Y, PAULI_Z};
446
447 // make one pauli wrong
448 codes[GENERATE( range(0,3) )] = (pauliOpType) GENERATE( -1, 4 );
449 REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, codes, numTargs, vecWork), ContainsSubstring("invalid Pauli code") );
450 }
451
452 // workspace is no longer needed; argument is ignored
453
454 // SECTION( "workspace type" ) {
455
456 // int numTargs = 1;
457 // int targs[1] = {0};
458 // pauliOpType codes[1] = {PAULI_I};
459
460 // REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, codes, numTargs, matWork), ContainsSubstring("Registers must both be state-vectors or both be density matrices") );
461 // REQUIRE_THROWS_WITH( calcExpecPauliProd(mat, targs, codes, numTargs, vecWork), ContainsSubstring("Registers must both be state-vectors or both be density matrices") );
462 // }
463 // SECTION( "workspace dimensions" ) {
464
465 // int numTargs = 1;
466 // int targs[1] = {0};
467 // pauliOpType codes[1] = {PAULI_I};
468
469 // Qureg vec2 = createQureg(NUM_QUBITS + 1);
470 // REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, codes, numTargs, vec2), ContainsSubstring("Dimensions") && ContainsSubstring("don't match") );
471 // destroyQureg(vec2);
472
473 // Qureg mat2 = createDensityQureg(NUM_QUBITS + 1);
474 // REQUIRE_THROWS_WITH( calcExpecPauliProd(mat, targs, codes, numTargs, mat2), ContainsSubstring("Dimensions") && ContainsSubstring("don't match") );
475 // destroyQureg(mat2);
476 // }
477 }
478 destroyQureg(vec);
479 destroyQureg(mat);
480 destroyQureg(vecWork);
481 destroyQureg(matWork);
482}
483
484
485
486/** @sa calcExpecPauliSum
487 * @ingroup deprecatedtests
488 * @author Tyson Jones
489 */
490TEST_CASE( "calcExpecPauliSum", "[calculations]" ) {
491
492 Qureg vec = createForcedQureg(NUM_QUBITS);
493 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
494
495 QVector vecRef = getRandomStateVector(NUM_QUBITS);
496 QMatrix matRef = getRandomDensityMatrix(NUM_QUBITS);
497 toQureg(vec, vecRef);
498 toQureg(mat, matRef);
499
500 // accepted by v3 deprecated API but discarded by v4
501 Qureg vecWork = createForcedQureg(NUM_QUBITS);
502 Qureg matWork = createForcedDensityQureg(NUM_QUBITS);
503
504 SECTION( "correctness" ) {
505
506 int numSumTerms = GENERATE( 1, 2, 10, 15 );
507
508 /* it's too expensive to try every possible Pauli configuration, so
509 * we'll try 10 random codes, and for each, random coefficients
510 */
511 GENERATE( range(0,10) );
512 int totNumCodes = numSumTerms * NUM_QUBITS;
513 vector<pauliOpType> paulis(totNumCodes);
514 vector<qreal> coeffs(numSumTerms);
515 setRandomPauliSum(coeffs.data(), paulis.data(), NUM_QUBITS, numSumTerms);
516
517 // produce a numTargs-big matrix 'pauliSum' by pauli-matrix tensoring and summing
518 QMatrix pauliSum = toQMatrix(coeffs.data(), paulis.data(), NUM_QUBITS, numSumTerms);
519
520 SECTION( "state-vector" ) {
521
522 /* calcExpecPauliSum calculates <qureg|pauliSum|qureg> */
523
524 QVector sumRef = pauliSum * vecRef;
525 qcomp prod = 0;
526 for (size_t i=0; i<vecRef.size(); i++)
527 prod += conj(vecRef[i]) * sumRef[i];
528 REQUIRE( imag(prod) == Approx(0).margin(10*REAL_EPS) );
529
530 qreal res = calcExpecPauliSum(vec, paulis.data(), coeffs.data(), numSumTerms, vecWork);
531 REQUIRE( res == Approx(real(prod)).margin(10*REAL_EPS) );
532 }
533 SECTION( "density-matrix" ) {
534
535 /* calcExpecPauliSum calculates Trace( pauliSum * qureg ) */
536 matRef = pauliSum * matRef;
537 qreal tr = 0;
538 for (size_t i=0; i<matRef.size(); i++)
539 tr += real(matRef[i][i]);
540 // (get real, since we start in a non-Hermitian state, hence diagonal isn't real)
541
542 qreal res = calcExpecPauliSum(mat, paulis.data(), coeffs.data(), numSumTerms, matWork);
543 REQUIRE( res == Approx(tr).margin(1E2*REAL_EPS) );
544 }
545 }
546 SECTION( "input validation" ) {
547
548 // cannot be validated; deprecated API copies before validating numSumTerms, causing segfault
549
550 // SECTION( "number of sum terms" ) {
551
552 // int numSumTerms = GENERATE( -1, 0 );
553 // REQUIRE_THROWS_WITH( calcExpecPauliSum(vec, NULL, NULL, numSumTerms, vecWork), ContainsSubstring("The number of terms must be a positive integer") );
554 // }
555
556 SECTION( "pauli codes" ) {
557
558 // make valid params
559 int numSumTerms = 3;
560 vector<qreal> coeffs(numSumTerms);
561 vector<pauliOpType> codes(numSumTerms*NUM_QUBITS);
562 for (int i=0; i<numSumTerms*NUM_QUBITS; i++)
563 codes[i] = PAULI_I;
564
565 // make one pauli wrong
566 codes[GENERATE_COPY( range(0,numSumTerms*NUM_QUBITS) )] = (pauliOpType) GENERATE( -1, 4 );
567 REQUIRE_THROWS_WITH( calcExpecPauliSum(vec, codes.data(), coeffs.data(), numSumTerms, vecWork), ContainsSubstring("invalid Pauli code") );
568 }
569
570 // the v4 API does not use a workspace, so it is discarded by the v3 deprecation layer
571
572 // SECTION( "workspace type" ) {
573
574 // // make valid params
575 // int numSumTerms = 1;
576 // qreal coeffs[1] = {0};
577 // pauliOpType codes[NUM_QUBITS];
578 // for (int i=0; i<NUM_QUBITS; i++)
579 // codes[i] = PAULI_I;
580
581 // REQUIRE_THROWS_WITH( calcExpecPauliSum(vec, codes, coeffs, numSumTerms, mat), ContainsSubstring("Registers must both be state-vectors or both be density matrices") );
582 // REQUIRE_THROWS_WITH( calcExpecPauliSum(mat, codes, coeffs, numSumTerms, vec), ContainsSubstring("Registers must both be state-vectors or both be density matrices") );
583 // }
584
585 // SECTION( "workspace dimensions" ) {
586
587 // // make valid params
588 // int numSumTerms = 1;
589 // qreal coeffs[1] = {0};
590 // pauliOpType codes[NUM_QUBITS];
591 // for (int i=0; i<NUM_QUBITS; i++)
592 // codes[i] = PAULI_I;
593
594 // Qureg vec2 = createQureg(NUM_QUBITS + 1);
595 // REQUIRE_THROWS_WITH( calcExpecPauliSum(vec, codes, coeffs, numSumTerms, vec2), ContainsSubstring("Dimensions") && ContainsSubstring("don't match") );
596 // destroyQureg(vec2);
597
598 // Qureg mat2 = createDensityQureg(NUM_QUBITS + 1);
599 // REQUIRE_THROWS_WITH( calcExpecPauliSum(mat, codes, coeffs, numSumTerms, mat2), ContainsSubstring("Dimensions") && ContainsSubstring("don't match") );
600 // destroyQureg(mat2);
601 // }
602 }
603 destroyQureg(vec);
604 destroyQureg(mat);
605 destroyQureg(vecWork);
606 destroyQureg(matWork);
607}
608
609
610
611/** @sa calcFidelity
612 * @ingroup deprecatedtests
613 * @author Tyson Jones
614 */
615TEST_CASE( "calcFidelity", "[calculations]" ) {
616
617 Qureg vec = createForcedQureg(NUM_QUBITS);
618 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
619 Qureg pure = createForcedQureg(NUM_QUBITS);
620
621 SECTION( "correctness" ) {
622
623 // repeat the below random tests 10 times
624 GENERATE( range(0,10) );
625
626 SECTION( "state-vector" ) {
627
628 /* calcFidelity computes |<vec|pure>|^2 */
629
630 SECTION( "normalised" ) {
631
632 // random L2 vectors
633 QVector vecRef = getRandomStateVector(NUM_QUBITS);
634 QVector pureRef = getRandomStateVector(NUM_QUBITS);
635 toQureg(vec, vecRef);
636 toQureg(pure, pureRef);
637
638 // |<vec|vec>|^2 = |1|^2 = 1
639 REQUIRE( calcFidelity(vec,vec) == Approx(1) );
640
641 // |<vec|pure>|^2 = |sum_j conj(vec_j) * pure_j|^2
642 qcomp dotProd = 0;
643 for (size_t i=0; i<vecRef.size(); i++)
644 dotProd += conj(vecRef[i]) * pureRef[i];
645 qreal refFid = pow(abs(dotProd), 2);
646
647 REQUIRE( calcFidelity(vec,pure) == Approx(refFid) );
648 }
649
650 // unnormalised test is no longer supported, since v4 calcFidelity
651 // validates thet the fidelity is correctly approximately real
652
653 // SECTION( "unnormalised" ) {
654
655 // // random unnormalised vectors
656 // QVector vecRef = getRandomQVector(1<<NUM_QUBITS);
657 // QVector pureRef = getRandomQVector(1<<NUM_QUBITS);
658 // toQureg(vec, vecRef);
659 // toQureg(pure, pureRef);
660
661 // // Let nv be magnitude of vec, hence |unit-vec> = 1/sqrt(nv)|vec>
662 // qreal nv = 0;
663 // for (size_t i=0; i<vecRef.size(); i++)
664 // nv += pow(abs(vecRef[i]), 2);
665 // // then <vec|vec> = sqrt(nv)*sqrt(nv) <unit-vec|unit-vec> = nv,
666 // // hence |<vec|vec>|^2 = nv*nv
667 // REQUIRE( calcFidelity(vec,vec) == Approx( nv*nv ) );
668
669 // qcomp dotProd = 0;
670 // for (size_t i=0; i<vecRef.size(); i++)
671 // dotProd += conj(vecRef[i]) * pureRef[i];
672 // qreal refFid = pow(abs(dotProd), 2);
673
674 // REQUIRE( calcFidelity(vec,pure) == Approx(refFid) );
675 // }
676 }
677 SECTION( "density-matrix" ) {
678
679 /* calcFidelity computes <pure|mat|pure> */
680
681 SECTION( "pure" ) {
682
683 QVector pureRef = getRandomStateVector(NUM_QUBITS);
684 toQureg(pure, pureRef);
685
686 // test when density matrix is the same pure state
687 QMatrix matRef = getKetBra(pureRef, pureRef);
688 toQureg(mat, matRef);
689 REQUIRE( calcFidelity(mat,pure) == Approx(1) );
690
691 // test when density matrix is a random pure state
692 QVector r1 = getRandomStateVector(NUM_QUBITS);
693 matRef = getKetBra(r1, r1); // actually pure |r1><r1|
694 toQureg(mat, matRef);
695
696 // <pure|r1><r1|pure> = |<r1|pure>|^2 = |sum_j conj(r1_j) * pure_j|^2
697 qcomp dotProd = 0;
698 for (size_t i=0; i<r1.size(); i++)
699 dotProd += conj(r1[i]) * pureRef[i];
700 qreal refFid = pow(abs(dotProd), 2);
701
702 REQUIRE( calcFidelity(mat,pure) == Approx(refFid).margin(100 * REAL_EPS) );
703 }
704 SECTION( "mixed" ) {
705
706 QVector pureRef = getRandomStateVector(NUM_QUBITS);
707 toQureg(pure, pureRef);
708
709 // test when density matrix is mixed
710 QMatrix matRef = getRandomDensityMatrix(NUM_QUBITS);
711 toQureg(mat, matRef);
712
713 // <pure|mat|pure> = <pure| (Mat|pure>)
714 QVector rhs = matRef * pureRef;
715 qcomp dotProd = 0;
716 for (size_t i=0; i<rhs.size(); i++)
717 dotProd += conj(pureRef[i]) * rhs[i];
718
719 REQUIRE( imag(dotProd) == Approx(0).margin(REAL_EPS) );
720 REQUIRE( calcFidelity(mat,pure) == Approx(real(dotProd)).margin(100 * REAL_EPS) );
721 }
722
723 // unnormalised test is no longer supported, since v4 calcFidelity
724 // validates thet the fidelity is correctly approximately real
725
726 // SECTION( "unnormalised" ) {
727
728 // // test when both density matrix and pure state are unnormalised
729 // QVector pureRef = getRandomQVector(1<<NUM_QUBITS);
730 // QMatrix matRef = getRandomQMatrix(1<<NUM_QUBITS);
731 // toQureg(pure, pureRef);
732 // toQureg(mat, matRef);
733
734 // // real[ <pure|mat|pure> ] = real[ <pure| (Mat|pure>) ]
735 // QVector rhs = matRef * pureRef;
736 // qcomp dotProd = 0;
737 // for (size_t i=0; i<rhs.size(); i++)
738 // dotProd += conj(pureRef[i]) * rhs[i];
739
740 // REQUIRE( calcFidelity(mat,pure) == Approx(real(dotProd)) );
741 // }
742 }
743 }
744 SECTION( "input validation" ) {
745
746 SECTION( "dimensions" ) {
747
748 // two state-vectors
749 Qureg vec2 = createQureg(vec.numQubits + 1);
750 REQUIRE_THROWS_WITH( calcFidelity(vec2,vec), ContainsSubstring("differing numbers of qubits") );
751 destroyQureg(vec2);
752
753 // density-matrix and state-vector
754 Qureg mat2 = createDensityQureg(vec.numQubits + 1);
755 REQUIRE_THROWS_WITH( calcFidelity(mat2,vec), ContainsSubstring("differing numbers of qubits") );
756 destroyQureg(mat2);
757 }
758 SECTION( "density-matrices" ) {
759
760 // two mixed statess
761 REQUIRE_THROWS_WITH( calcFidelity(mat,mat), ContainsSubstring("Quregs cannot both be density matrices") );
762 }
763 }
764 destroyQureg(vec);
765 destroyQureg(mat);
766 destroyQureg(pure);
767}
768
769
770
771/** @sa calcHilbertSchmidtDistance
772 * @ingroup deprecatedtests
773 * @author Tyson Jones
774 */
775TEST_CASE( "calcHilbertSchmidtDistance", "[calculations]" ) {
776
777 Qureg mat1 = createForcedDensityQureg(NUM_QUBITS);
778 Qureg mat2 = createForcedDensityQureg(NUM_QUBITS);
779
780 SECTION( "correctness" ) {
781
782 // perform these random tests 10 times
783 GENERATE( range(0,10) );
784
785 SECTION( "density-matrix" ) {
786
787 SECTION( "pure" ) {
788
789 // create random |r1><r1| and |r2><r2| states
790 QVector r1 = getRandomStateVector(NUM_QUBITS);
791 QMatrix m1 = getKetBra(r1,r1);
792 toQureg(mat1, m1);
793 QVector r2 = getRandomStateVector(NUM_QUBITS);
794 QMatrix m2 = getKetBra(r2,r2);
795 toQureg(mat2, m2);
796
797 // Tr{ (a-b)(a-b)^dagger } = sum_{ij} |a_{ij} - b_{ij}|^2
798 qreal tr = 0;
799 for (size_t i=0; i<m1.size(); i++)
800 for (size_t j=0; j<m1.size(); j++)
801 tr += pow(abs(m1[i][j] - m2[i][j]), 2);
802
803 qreal res = calcHilbertSchmidtDistance(mat1, mat2);
804 REQUIRE( res == Approx(sqrt(tr)) );
805
806 }
807 SECTION( "normalised" ) {
808
809 QMatrix ref1 = getRandomDensityMatrix(NUM_QUBITS);
810 QMatrix ref2 = getRandomDensityMatrix(NUM_QUBITS);
811 toQureg(mat1, ref1);
812 toQureg(mat2, ref2);
813
814 // Tr{ (a-b)(a-b)^dagger } = sum_{ij} |a_{ij} - b_{ij}|^2
815 qreal tr = 0;
816 for (size_t i=0; i<ref1.size(); i++)
817 for (size_t j=0; j<ref1.size(); j++)
818 tr += pow(abs(ref1[i][j] - ref2[i][j]), 2);
819
820 qreal res = calcHilbertSchmidtDistance(mat1, mat2);
821 REQUIRE( res == Approx(sqrt(tr)) );
822 }
823 SECTION( "unnormalised" ) {
824
825 // mat1 and mat2 are both random matrices
826 QMatrix ref1 = getRandomQMatrix(1<<NUM_QUBITS);
827 QMatrix ref2 = getRandomQMatrix(1<<NUM_QUBITS);
828 toQureg(mat1, ref1);
829 toQureg(mat2, ref2);
830
831 // Tr{ (a-b)(a-b)^dagger } = sum_{ij} |a_{ij} - b_{ij}|^2
832 qreal tr = 0;
833 for (size_t i=0; i<ref1.size(); i++)
834 for (size_t j=0; j<ref1.size(); j++)
835 tr += pow(abs(ref1[i][j] - ref2[i][j]), 2);
836
837 qreal res = calcHilbertSchmidtDistance(mat1, mat2);
838 REQUIRE( res == Approx(sqrt(tr)) );
839 }
840 }
841 }
842 SECTION( "input validation") {
843
844 SECTION( "dimensions" ) {
845
846 Qureg mat3 = createDensityQureg(NUM_QUBITS + 1);
847 REQUIRE_THROWS_WITH( calcHilbertSchmidtDistance(mat1,mat3), ContainsSubstring("differing numbers of qubits") );
848 destroyQureg(mat3);
849 }
850
851 // in v4, this function calls calcDistance() which has state-vector overloads
852
853 // SECTION( "state-vector" ) {
854
855 // Qureg vec = createForcedQureg(NUM_QUBITS);
856
857 // REQUIRE_THROWS_WITH( calcHilbertSchmidtDistance(vec,mat1), ContainsSubstring("valid only for density matrices") );
858 // REQUIRE_THROWS_WITH( calcHilbertSchmidtDistance(mat1,vec), ContainsSubstring("valid only for density matrices") );
859 // REQUIRE_THROWS_WITH( calcHilbertSchmidtDistance(vec,vec), ContainsSubstring("valid only for density matrices") );
860
861 // destroyQureg(vec);
862 // }
863 }
864 destroyQureg(mat1);
865 destroyQureg(mat2);
866}
867
868
869
870/** @sa calcInnerProduct
871 * @ingroup deprecatedtests
872 * @author Tyson Jones
873 */
874TEST_CASE( "calcInnerProduct", "[calculations]" ) {
875
876 Qureg vec1 = createForcedQureg(NUM_QUBITS);
877 Qureg vec2 = createForcedQureg(NUM_QUBITS);
878
879 SECTION( "correctness" ) {
880
881 // perform these random tests 10 times
882 GENERATE( range(0,10) );
883
884 SECTION( "state-vector" ) {
885
886 SECTION( "normalised" ) {
887
888 // <r1|r2> = sum_j conj(r1_j) * r2_j
889 QVector r1 = getRandomStateVector(NUM_QUBITS);
890 QVector r2 = getRandomStateVector(NUM_QUBITS);
891 qcomp prod = 0;
892 for (size_t i=0; i<r1.size(); i++)
893 prod += conj(r1[i]) * r2[i];
894
895 toQureg(vec1, r1);
896 toQureg(vec2, r2);
897 qcomp res = calcInnerProduct(vec1,vec2);
898
899 REQUIRE( real(res) == Approx(real(prod)) );
900 REQUIRE( imag(res) == Approx(imag(prod)) );
901 }
902 SECTION( "unnormalised" ) {
903
904 // <r1|r2> = sum_j conj(r1_j) * r2_j
905 QVector r1 = getRandomQVector(1<<NUM_QUBITS);
906 QVector r2 = getRandomQVector(1<<NUM_QUBITS);
907 qcomp prod = 0;
908 for (size_t i=0; i<r1.size(); i++)
909 prod += conj(r1[i]) * r2[i];
910
911 toQureg(vec1, r1);
912 toQureg(vec2, r2);
913 qcomp res = calcInnerProduct(vec1,vec2);
914
915 REQUIRE( real(res) == Approx(real(prod)).margin(REAL_EPS) );
916 REQUIRE( imag(res) == Approx(imag(prod)).margin(REAL_EPS) );
917 }
918 }
919 }
920 SECTION( "input validation" ) {
921
922 SECTION( "dimensions" ) {
923
924 Qureg vec3 = createQureg(NUM_QUBITS + 1);
925 REQUIRE_THROWS_WITH( calcInnerProduct(vec1,vec3), ContainsSubstring("differing numbers of qubits") );
926 destroyQureg(vec3);
927 }
928
929 // density-matrix arguments are permitted in v4
930
931 // SECTION( "density-matrix" ) {
932
933 // Qureg mat = createForcedDensityQureg(NUM_QUBITS);
934
935 // REQUIRE_THROWS_WITH( calcInnerProduct(vec1,mat), ContainsSubstring("valid only for state-vectors") );
936 // REQUIRE_THROWS_WITH( calcInnerProduct(mat,vec1), ContainsSubstring("valid only for state-vectors") );
937 // REQUIRE_THROWS_WITH( calcInnerProduct(mat,mat), ContainsSubstring("valid only for state-vectors") );
938
939 // destroyQureg(mat);
940 // }
941 }
942 destroyQureg(vec1);
943 destroyQureg(vec2);
944}
945
946
947
948/** @sa calcProbOfAllOutcomes
949 * @ingroup deprecatedtests
950 * @author Tyson Jones
951 */
952TEST_CASE( "calcProbOfAllOutcomes", "[calculations]" ) {
953
954 Qureg vec = createForcedQureg(NUM_QUBITS);
955 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
956
957 SECTION( "correctness" ) {
958
959 // generate all possible qubit arrangements
960 int numQubits = GENERATE_COPY( range(1,NUM_QUBITS+1) );
961 int* qubits = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numQubits) );
962
963 int numOutcomes = 1<<numQubits;
964 vector<qreal> probs(numOutcomes);
965 QVector refProbs = QVector(numOutcomes);
966
967 SECTION( "state-vector" ) {
968
969 SECTION( "normalised" ) {
970
971 QVector ref = getRandomStateVector(NUM_QUBITS);
972 toQureg(vec, ref);
973
974 // prob is sum of |amp|^2 of basis states which encode outcome
975 for (size_t i=0; i<ref.size(); i++) {
976 int outcome = 0;
977 for (int q=0; q<numQubits; q++) {
978 int bit = (i >> qubits[q]) & 1;
979 outcome += bit * (1 << q);
980 }
981 refProbs[outcome] += pow(abs(ref[i]), 2);
982 }
983
984 calcProbOfAllOutcomes(probs.data(), vec, qubits, numQubits);
985 REQUIRE( areEqual(refProbs, probs.data()) );
986 }
987 SECTION( "unnormalised" ) {
988
989 QVector ref = getRandomQVector(1<<NUM_QUBITS);
990 toQureg(vec, ref);
991
992 // prob is sum of |amp|^2 of basis states which encode outcome
993 for (size_t i=0; i<ref.size(); i++) {
994 int outcome = 0;
995 for (int q=0; q<numQubits; q++) {
996 int bit = (i >> qubits[q]) & 1;
997 outcome += bit * (1 << q);
998 }
999 refProbs[outcome] += pow(abs(ref[i]), 2);
1000 }
1001
1002 calcProbOfAllOutcomes(probs.data(), vec, qubits, numQubits);
1003 REQUIRE( areEqual(refProbs, probs.data()) );
1004 }
1005 }
1006 SECTION( "density-matrix" ) {
1007
1008 SECTION( "normalised" ) {
1009
1010 QMatrix ref = getRandomDensityMatrix(NUM_QUBITS);
1011 toQureg(mat, ref);
1012
1013 // prob is sum of diagonals which encode outcome
1014 for (size_t i=0; i<ref.size(); i++) {
1015 int outcome = 0;
1016 for (int q=0; q<numQubits; q++) {
1017 int bit = (i >> qubits[q]) & 1;
1018 outcome += bit * (1 << q);
1019 }
1020 refProbs[outcome] += real(ref[i][i]);
1021 }
1022
1023 calcProbOfAllOutcomes(probs.data(), mat, qubits, numQubits);
1024 REQUIRE( areEqual(refProbs, probs.data()) );
1025 }
1026 SECTION( "unnormalised" ) {
1027
1028 QMatrix ref = getRandomQMatrix(1<<NUM_QUBITS);
1029 toQureg(mat, ref);
1030
1031 // prob is sum of diagonals which encode outcome
1032 for (size_t i=0; i<ref.size(); i++) {
1033 int outcome = 0;
1034 for (int q=0; q<numQubits; q++) {
1035 int bit = (i >> qubits[q]) & 1;
1036 outcome += bit * (1 << q);
1037 }
1038 refProbs[outcome] += real(ref[i][i]);
1039 }
1040
1041 calcProbOfAllOutcomes(probs.data(), mat, qubits, numQubits);
1042 REQUIRE( areEqual(refProbs, probs.data()) );
1043 }
1044 }
1045 }
1046 SECTION( "input validation" ) {
1047
1048 int numQubits = 3;
1049 int qubits[] = {0, 1, 2};
1050 qreal probs[8];
1051
1052 SECTION( "number of qubits" ) {
1053
1054 // too small
1055 numQubits = GENERATE( -1, 0 );
1056 REQUIRE_THROWS_WITH( calcProbOfAllOutcomes(probs, mat, qubits, numQubits), ContainsSubstring("specified number of target qubits") && ContainsSubstring("invalid.") );
1057
1058 // too big
1059 numQubits = NUM_QUBITS + 1;
1060 REQUIRE_THROWS_WITH( calcProbOfAllOutcomes(probs, mat, qubits, numQubits), ContainsSubstring("target qubits") && ContainsSubstring("exceeds the number of qubits in the Qureg") );
1061 }
1062 SECTION( "qubit indices" ) {
1063
1064 qubits[GENERATE_COPY(range(0,numQubits))] = GENERATE( -1, NUM_QUBITS );
1065 REQUIRE_THROWS_WITH( calcProbOfAllOutcomes(probs, mat, qubits, numQubits), ContainsSubstring("Invalid target qubit") );
1066 }
1067 SECTION( "repetition of qubits" ) {
1068
1069 qubits[GENERATE_COPY(1,2)] = qubits[0];
1070 REQUIRE_THROWS_WITH( calcProbOfAllOutcomes(probs, mat, qubits, numQubits), ContainsSubstring("qubits must be unique") );
1071 }
1072 }
1073 destroyQureg(vec);
1074 destroyQureg(mat);
1075}
1076
1077
1078
1079
1080/** @sa calcProbOfOutcome
1081 * @ingroup deprecatedtests
1082 * @author Tyson Jones
1083 */
1084TEST_CASE( "calcProbOfOutcome", "[calculations]" ) {
1085
1086 Qureg vec = createForcedQureg(NUM_QUBITS);
1087 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1088
1089 SECTION( "correctness" ) {
1090
1091 int target = GENERATE( range(0,NUM_QUBITS) );
1092 int outcome = GENERATE( 0, 1 );
1093
1094 SECTION( "state-vector" ) {
1095
1096 SECTION( "normalised" ) {
1097
1098 QVector ref = getRandomStateVector(NUM_QUBITS);
1099 toQureg(vec, ref);
1100
1101 // prob is sum of |amp|^2 of amplitudes where target bit is outcome
1102 qreal prob = 0;
1103 for (size_t ind=0; ind<ref.size(); ind++) {
1104 int bit = (ind >> target) & 1; // target-th bit
1105 if (bit == outcome)
1106 prob += pow(abs(ref[ind]), 2);
1107 }
1108
1109 REQUIRE( calcProbOfOutcome(vec, target, outcome) == Approx(prob) );
1110 }
1111 SECTION( "unnormalised" ) {
1112
1113 QVector ref = getRandomQVector(1<<NUM_QUBITS);
1114 toQureg(vec, ref);
1115
1116 // prob is sum of |amp|^2 of amplitudes where target bit is outcome
1117 qreal prob = 0;
1118 for (size_t ind=0; ind<ref.size(); ind++) {
1119 int bit = (ind >> target) & 1; // target-th bit
1120 if (bit == outcome)
1121 prob += pow(abs(ref[ind]), 2);
1122 }
1123
1124 REQUIRE( calcProbOfOutcome(vec, target, outcome) == Approx(prob) );
1125 }
1126 }
1127 SECTION( "density-matrix" ) {
1128
1129 SECTION( "pure" ) {
1130
1131 // set mat to a random |r><r|
1132 QVector ref = getRandomStateVector(NUM_QUBITS);
1133 toQureg(mat, getKetBra(ref, ref));
1134
1135 // calc prob of the state-vector
1136 qreal prob = 0;
1137 for (size_t ind=0; ind<ref.size(); ind++) {
1138 int bit = (ind >> target) & 1; // target-th bit
1139 if (bit == outcome)
1140 prob += pow(abs(ref[ind]), 2);
1141 }
1142
1143 REQUIRE( calcProbOfOutcome(mat, target, outcome) == Approx(prob) );
1144 }
1145 SECTION( "mixed" ) {
1146
1147 QMatrix ref = getRandomDensityMatrix(NUM_QUBITS);
1148 toQureg(mat, ref);
1149
1150 // prob is sum of diagonal amps (should be real) where target bit is outcome
1151 qcomp tr = 0;
1152 for (size_t ind=0; ind<ref.size(); ind++) {
1153 int bit = (ind >> target) & 1; // target-th bit
1154 if (bit == outcome)
1155 tr += ref[ind][ind];
1156 }
1157
1158 REQUIRE( imag(tr) == Approx(0).margin(REAL_EPS) );
1159
1160 REQUIRE( calcProbOfOutcome(mat, target, outcome) == Approx(real(tr)) );
1161 }
1162 SECTION( "unnormalised" ) {
1163
1164 QMatrix ref = getRandomQMatrix(1<<NUM_QUBITS);
1165 toQureg(mat, ref);
1166
1167 // prob is (sum of real of diagonal amps where target bit is outcome)
1168 qreal tr = 0;
1169 for (size_t ind=0; ind<ref.size(); ind++) {
1170 int bit = (ind >> target) & 1; // target-th bit
1171 if (bit == outcome)
1172 tr += real(ref[ind][ind]);
1173 }
1174
1175 REQUIRE( calcProbOfOutcome(mat, target, outcome) == Approx(tr) );
1176 }
1177 }
1178 }
1179 SECTION( "input validation" ) {
1180
1181 SECTION( "qubit indices" ) {
1182
1183 int target = GENERATE( -1, NUM_QUBITS );
1184 REQUIRE_THROWS_WITH( calcProbOfOutcome(vec, target, 0), ContainsSubstring("Invalid target qubit") );
1185 }
1186 SECTION( "outcome value" ) {
1187
1188 int outcome = GENERATE( -1, 2 );
1189 REQUIRE_THROWS_WITH( calcProbOfOutcome(vec, 0, outcome), ContainsSubstring("measurement outcome") && ContainsSubstring("invalid") );
1190 }
1191 }
1192 destroyQureg(vec);
1193 destroyQureg(mat);
1194}
1195
1196
1197
1198/** @sa calcPurity
1199 * @ingroup deprecatedtests
1200 * @author Tyson Jones
1201 */
1202TEST_CASE( "calcPurity", "[calculations]" ) {
1203
1204 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1205
1206 SECTION( "correctness" ) {
1207
1208 // perform the following random tests 10 times
1209 GENERATE( range(1,10) );
1210
1211 SECTION( "density-matrix" ) {
1212
1213 SECTION( "pure" ) {
1214
1215 // pure states have unity purity
1216 initZeroState(mat);
1217 REQUIRE( calcPurity(mat) == 1 );
1218
1219 // (try also a pure random L2-vector)
1220 QVector r1 = getRandomStateVector(NUM_QUBITS); // |r>
1221 QMatrix m1 = getKetBra(r1, r1); // |r><r|
1222 toQureg(mat, m1);
1223 REQUIRE( calcPurity(mat) == Approx(1) );
1224
1225 }
1226 SECTION( "mixed" ) {
1227
1228 // mixed states have 1/2^N < purity < 1
1229 QMatrix ref = getRandomDensityMatrix(NUM_QUBITS);
1230 toQureg(mat, ref);
1231 qreal purity = calcPurity(mat);
1232 REQUIRE( purity < 1 );
1233 REQUIRE( purity >= 1/pow(2.,NUM_QUBITS) );
1234
1235 // compare to Tr(rho^2)
1236 QMatrix prod = ref*ref;
1237 qreal tr = 0;
1238 for (size_t i=0; i<prod.size(); i++)
1239 tr += real(prod[i][i]);
1240 REQUIRE( purity == Approx(tr) );
1241 }
1242 SECTION( "unnormalised" ) {
1243
1244 // unphysical states give sum_{ij} |rho_ij|^2
1245 QMatrix ref = getRandomQMatrix(1<<NUM_QUBITS);
1246 qreal tot = 0;
1247 for (size_t i=0; i<ref.size(); i++)
1248 for (size_t j=0; j<ref.size(); j++)
1249 tot += pow(abs(ref[i][j]), 2);
1250
1251 toQureg(mat, ref);
1252 REQUIRE( calcPurity(mat) == Approx(tot) );
1253 }
1254 }
1255 }
1256 SECTION( "input validation" ) {
1257
1258 // in v4, this accepts state-vectors
1259
1260 // SECTION( "state-vector" ) {
1261
1262 // Qureg vec = createForcedQureg(NUM_QUBITS);
1263 // REQUIRE_THROWS_WITH( calcPurity(vec), ContainsSubstring("valid only for density matrices") );
1264 // destroyQureg(vec);
1265 // }
1266 }
1267 destroyQureg(mat);
1268}
1269
1270
1271
1272/** @sa calcTotalProb
1273 * @ingroup deprecatedtests
1274 * @author Tyson Jones
1275 */
1276TEST_CASE( "calcTotalProb", "[calculations]" ) {
1277
1278 Qureg vec = createForcedQureg(NUM_QUBITS);
1279 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1280
1281 SECTION( "correctness" ) {
1282
1283 SECTION( "state-vector" ) {
1284
1285 // normalised: prob(vec) = 1
1286 initPlusState(vec);
1287 REQUIRE( calcTotalProb(vec) == Approx(1) );
1288
1289 // zero norm: prob(vec) = 0
1290 initBlankState(vec);
1291 REQUIRE( calcTotalProb(vec) == 0 );
1292
1293 // random L2 state: prob(vec) = 1
1294 toQureg(vec, getRandomStateVector(NUM_QUBITS));
1295 REQUIRE( calcTotalProb(vec) == Approx(1) );
1296
1297 // unnormalised: prob(vec) = sum_i |vec_i|^2
1298 initDebugState(vec);
1299 QVector ref = toQVector(vec);
1300 qreal refProb = 0;
1301 for (size_t i=0; i<ref.size(); i++)
1302 refProb += pow(abs(ref[i]), 2);
1303 REQUIRE( calcTotalProb(vec) == Approx(refProb) );
1304 }
1305 SECTION( "density-matrix" ) {
1306
1307 // normalised: prob(mat) = 1
1308 initPlusState(mat);
1309 REQUIRE( calcTotalProb(mat) == Approx(1) );
1310
1311 // zero norm: prob(mat) = 0
1312 initBlankState(mat);
1313 REQUIRE( calcTotalProb(mat) == 0 );
1314
1315 // random density matrix: prob(mat) = 1
1316 toQureg(mat, getRandomDensityMatrix(NUM_QUBITS));
1317 REQUIRE( calcTotalProb(mat) == Approx(1) );
1318
1319 // unnormalised: prob(mat) = sum_i real(mat_{ii})
1320 initDebugState(mat);
1321 QMatrix ref = toQMatrix(mat);
1322 qreal refProb = 0;
1323 for (size_t i=0; i<ref.size(); i++)
1324 refProb += real(ref[i][i]);
1325 REQUIRE( calcTotalProb(mat) == Approx(refProb) );
1326 }
1327 }
1328 SECTION( "input validation" ) {
1329
1330 // no validation
1331 SUCCEED();
1332 }
1333 destroyQureg(vec);
1334 destroyQureg(mat);
1335}
1336
1337
1338
1339/** @sa getAmp
1340 * @ingroup deprecatedtests
1341 * @author Tyson Jones
1342 */
1343TEST_CASE( "getAmp", "[calculations]" ) {
1344
1345 Qureg vec = createForcedQureg(NUM_QUBITS);
1346
1347 SECTION( "correctness" ) {
1348
1349 SECTION( "state-vector" ) {
1350
1351 initDebugState(vec);
1352 QVector ref = toQVector(vec);
1353
1354 int ind = GENERATE( range(0,1<<NUM_QUBITS) );
1355 qcomp amp = getAmp(vec,ind);
1356 REQUIRE( amp == ref[ind] );
1357 }
1358 }
1359 SECTION( "input validation" ) {
1360
1361 SECTION( "state index" ) {
1362
1363 int ind = GENERATE( -1, 1<<NUM_QUBITS );
1364 REQUIRE_THROWS_WITH( getAmp(vec,ind), ContainsSubstring("Basis state index") && ContainsSubstring("invalid") );
1365 }
1366 SECTION( "density-matrix" ) {
1367
1368 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1369 REQUIRE_THROWS_WITH( getAmp(mat,0), ContainsSubstring("Expected a statevector Qureg but received a density matrix") );
1370 destroyQureg(mat);
1371 }
1372 }
1373 destroyQureg(vec);
1374}
1375
1376
1377
1378/** @sa getDensityAmp
1379 * @ingroup deprecatedtests
1380 * @author Tyson Jones
1381 */
1382TEST_CASE( "getDensityAmp", "[calculations]" ) {
1383
1384 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1385
1386 SECTION( "correctness" ) {
1387
1388 SECTION( "density-matrix" ) {
1389
1390 initDebugState(mat);
1391 QMatrix ref = toQMatrix(mat);
1392
1393 int row = GENERATE( range(0,1<<NUM_QUBITS) );
1394 int col = GENERATE( range(0,1<<NUM_QUBITS) );
1395
1396 qcomp amp = getDensityAmp(mat,row,col);
1397 REQUIRE( amp == ref[row][col] );
1398 }
1399 }
1400 SECTION( "input validation" ) {
1401
1402 SECTION( "state index" ) {
1403
1404 int ind = GENERATE( -1, 1<<NUM_QUBITS );
1405 REQUIRE_THROWS_WITH( getDensityAmp(mat,ind,0), ContainsSubstring("The row and column indices") && ContainsSubstring("invalid") );
1406 REQUIRE_THROWS_WITH( getDensityAmp(mat,0,ind), ContainsSubstring("The row and column indices") && ContainsSubstring("invalid") );
1407
1408 }
1409 SECTION( "state-vector" ) {
1410
1411 Qureg vec = createForcedQureg(NUM_QUBITS);
1412 REQUIRE_THROWS_WITH( getDensityAmp(vec,0,0), ContainsSubstring("Expected a density matrix Qureg but received a statevector") );
1413 destroyQureg(vec);
1414 }
1415 }
1416 destroyQureg(mat);
1417}
1418
1419
1420
1421/** @sa getImagAmp
1422 * @ingroup deprecatedtests
1423 * @author Tyson Jones
1424 */
1425TEST_CASE( "getImagAmp", "[calculations]" ) {
1426
1427 Qureg vec = createForcedQureg(NUM_QUBITS);
1428
1429 SECTION( "correctness" ) {
1430
1431 SECTION( "state-vector" ) {
1432
1433 initDebugState(vec);
1434 QVector ref = toQVector(vec);
1435
1436 int ind = GENERATE( range(0,1<<NUM_QUBITS) );
1437 REQUIRE( getImagAmp(vec,ind) == imag(ref[ind]) );
1438 }
1439 }
1440 SECTION( "input validation" ) {
1441
1442 SECTION( "state index" ) {
1443
1444 int ind = GENERATE( -1, 1<<NUM_QUBITS );
1445 REQUIRE_THROWS_WITH( getImagAmp(vec,ind), ContainsSubstring("Basis state index") && ContainsSubstring("invalid") );
1446 }
1447 SECTION( "density-matrix" ) {
1448
1449 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1450 REQUIRE_THROWS_WITH( getImagAmp(mat,0), ContainsSubstring("Expected a statevector Qureg but received a density matrix") );
1451 destroyQureg(mat);
1452 }
1453 }
1454 destroyQureg(vec);
1455}
1456
1457
1458
1459/** @sa getNumAmps
1460 * @ingroup deprecatedtests
1461 * @author Tyson Jones
1462 */
1463TEST_CASE( "getNumAmps", "[calculations]" ) {
1464
1465 SECTION( "correctness" ) {
1466
1467 // test >= NUM_QUBITS so as not to limit distribution size
1468 int numQb = GENERATE( range(NUM_QUBITS, NUM_QUBITS+10) );
1469
1470 SECTION( "state-vector" ) {
1471
1472 Qureg vec = createForcedQureg(numQb);
1473 REQUIRE( getNumAmps(vec) == (1<<numQb) );
1474 destroyQureg(vec);
1475 }
1476 }
1477
1478 // in v4, argument can be a density matrix (return total num amps)
1479
1480 // SECTION( "input validation" ) {
1481
1482 // SECTION( "density-matrix" ) {
1483 // Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1484 // REQUIRE_THROWS_WITH( getNumAmps(mat), ContainsSubstring("Expected a statevector Qureg but received a density matrix") );
1485 // destroyQureg(mat);
1486 // }
1487 // }
1488}
1489
1490
1491
1492/** @sa getNumQubits
1493 * @ingroup deprecatedtests
1494 * @author Tyson Jones
1495 */
1496TEST_CASE( "getNumQubits", "[calculations]" ) {
1497
1498 SECTION( "correctness" ) {
1499
1500 // test >= NUM_QUBITS so as not to limit distribution size
1501 int numQb = GENERATE( range(NUM_QUBITS, NUM_QUBITS+10) );
1502
1503 SECTION( "state-vector" ) {
1504
1505 Qureg vec = createForcedQureg(numQb);
1506 REQUIRE( getNumQubits(vec) == numQb );
1507 destroyQureg(vec);
1508 }
1509 SECTION( "density-matrix" ) {
1510
1511 // density matrices use square as much memory; we must be careful not to seg-fault!
1512 if (2*numQb > 25)
1513 numQb = 13; // max size
1514
1515 Qureg mat = createForcedDensityQureg(numQb);
1516 REQUIRE( getNumQubits(mat) == numQb );
1517 destroyQureg(mat);
1518 }
1519 }
1520 SECTION( "input validation" ) {
1521
1522 // no validation
1523 SUCCEED();
1524 }
1525}
1526
1527
1528
1529/** @sa getProbAmp
1530 * @ingroup deprecatedtests
1531 * @author Tyson Jones
1532 */
1533TEST_CASE( "getProbAmp", "[calculations]" ) {
1534
1535 Qureg vec = createForcedQureg(NUM_QUBITS);
1536
1537 SECTION( "correctness" ) {
1538
1539 SECTION( "state-vector" ) {
1540
1541 initDebugState(vec);
1542 QVector ref = toQVector(vec);
1543
1544 int ind = GENERATE( range(0,1<<NUM_QUBITS) );
1545 qreal refCalc = pow(abs(ref[ind]), 2);
1546 REQUIRE( getProbAmp(vec,ind) == Approx(refCalc) );
1547 }
1548 }
1549 SECTION( "input validation" ) {
1550
1551 SECTION( "state index" ) {
1552
1553 int ind = GENERATE( -1, 1<<NUM_QUBITS );
1554 REQUIRE_THROWS_WITH( getProbAmp(vec,ind), ContainsSubstring("Basis state index") && ContainsSubstring("invalid") );
1555 }
1556
1557 // in v4, this redirects to calcProbOfBasisState() which accepts
1558 // both statevectors and density matrices
1559
1560 // SECTION( "density-matrix" ) {
1561
1562 // Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1563 // REQUIRE_THROWS_WITH( getProbAmp(mat,0), ContainsSubstring("Expected a statevector Qureg but received a density matrix") );
1564 // destroyQureg(mat);
1565 // }
1566 }
1567 destroyQureg(vec);
1568}
1569
1570
1571
1572/** @sa getRealAmp
1573 * @ingroup deprecatedtests
1574 * @author Tyson Jones
1575 */
1576TEST_CASE( "getRealAmp", "[calculations]" ) {
1577
1578 Qureg vec = createForcedQureg(NUM_QUBITS);
1579
1580 SECTION( "correctness" ) {
1581
1582 SECTION( "state-vector" ) {
1583
1584 initDebugState(vec);
1585 QVector ref = toQVector(vec);
1586
1587 int ind = GENERATE( range(0,1<<NUM_QUBITS) );
1588 REQUIRE( getRealAmp(vec,ind) == real(ref[ind]) );
1589 }
1590 }
1591 SECTION( "input validation" ) {
1592
1593 SECTION( "state index" ) {
1594
1595 int ind = GENERATE( -1, 1<<NUM_QUBITS );
1596 REQUIRE_THROWS_WITH( getRealAmp(vec,ind), ContainsSubstring("Basis state index") && ContainsSubstring("invalid") );
1597 }
1598 SECTION( "density-matrix" ) {
1599
1600 Qureg mat = createForcedDensityQureg(NUM_QUBITS);
1601 REQUIRE_THROWS_WITH( getRealAmp(mat,0), ContainsSubstring("Expected a statevector Qureg but received a density matrix") );
1602 destroyQureg(mat);
1603 }
1604 }
1605 destroyQureg(vec);
1606}
1607
1608
calcFidelity
qreal calcFidelity(Qureg qureg, Qureg other)
Definition calculations.cpp:303
calcInnerProduct
qcomp calcInnerProduct(Qureg qureg, Qureg other)
Definition calculations.cpp:49
calcPurity
qreal calcPurity(Qureg qureg)
Definition calculations.cpp:293
calcTotalProb
qreal calcTotalProb(Qureg qureg)
Definition calculations.cpp:278
setValidationOff
void setValidationOff()
Definition debug.cpp:86
setValidationOn
void setValidationOn()
Definition debug.cpp:80
TEST_CASE
TEST_CASE("calcDensityInnerProduct", "[calculations]")
Definition test_calculations.cpp:30
setRandomPauliSum
void setRandomPauliSum(qreal *coeffs, pauliOpType *codes, int numQubits, int numTerms)
Definition test_utilities.cpp:1370
getFullOperatorMatrix
QMatrix getFullOperatorMatrix(int *ctrls, int numCtrls, int *targs, int numTargs, QMatrix op, int numQubits)
Definition test_utilities.cpp:417
getRandomQVector
QVector getRandomQVector(int dim)
Definition test_utilities.cpp:549
getKetBra
QMatrix getKetBra(QVector ket, QVector bra)
Definition test_utilities.cpp:274
areEqual
bool areEqual(QVector a, QVector b)
Definition test_utilities.cpp:512
QMatrix
vector< vector< qcomp > > QMatrix
Definition test_utilities.hpp:60
applyReferenceOp
void applyReferenceOp(QVector &state, int *ctrls, int numCtrls, int *targs, int numTargs, QMatrix op)
Definition test_utilities.cpp:869
getRandomQMatrix
QMatrix getRandomQMatrix(int dim)
Definition test_utilities.cpp:492
toQMatrix
QMatrix toQMatrix(CompMatr1 src)
Definition test_utilities.cpp:1199
QVector
vector< qcomp > QVector
Definition test_utilities.hpp:71
toQVector
QVector toQVector(Qureg qureg)
Definition test_utilities.cpp:1250
toQureg
void toQureg(Qureg qureg, QVector vec)
Definition test_utilities.cpp:1327
getQuESTEnv
QuESTEnv getQuESTEnv()
Definition environment.cpp:416
initPlusState
void initPlusState(Qureg qureg)
Definition initialisations.cpp:59
initZeroState
void initZeroState(Qureg qureg)
Definition initialisations.cpp:50
initDebugState
void initDebugState(Qureg qureg)
Definition initialisations.cpp:92
initBlankState
void initBlankState(Qureg qureg)
Definition initialisations.cpp:41
createDensityQureg
Qureg createDensityQureg(int numQubits)
Definition qureg.cpp:291
createForcedQureg
Qureg createForcedQureg(int numQubits)
Definition qureg.cpp:299
createForcedDensityQureg
Qureg createForcedDensityQureg(int numQubits)
Definition qureg.cpp:309
createCustomQureg
Qureg createCustomQureg(int numQubits, int isDensMatr, int useDistrib, int useGpuAccel, int useMultithread)
Definition qureg.cpp:277
createQureg
Qureg createQureg(int numQubits)
Definition qureg.cpp:283
destroyQureg
void destroyQureg(Qureg qureg)
Definition qureg.cpp:334
getKroneckerProduct
qmatrix getKroneckerProduct(qmatrix a, qmatrix b)
Definition linalg.cpp:523
getRandomStateVector
qvector getRandomStateVector(int numQb)
Definition random.cpp:296
getRandomComplex
qcomp getRandomComplex()
Definition random.cpp:107
getRandomDensityMatrix
qmatrix getRandomDensityMatrix(int numQb)
Definition random.cpp:308
getRandomInt
int getRandomInt(int min, int maxExcl)
Definition random.cpp:90
Qureg
Definition qureg.h:49