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