test_calculations.cpp
Go to the documentation of this file.
1 
2 #include "catch.hpp"
3 #include "QuEST.h"
4 #include "utilities.hpp"
5 
6 /* allows concise use of Contains in catch's REQUIRE_THROWS_WITH */
7 using Catch::Matchers::Contains;
8 
9 
10 
11 
12 // DEBUG
13 #include <sys/time.h>
14 
15 
16 
21 TEST_CASE( "calcDensityInnerProduct", "[calculations]" ) {
22 
23  Qureg mat1 = createDensityQureg(NUM_QUBITS, QUEST_ENV);
24  Qureg mat2 = createDensityQureg(NUM_QUBITS, QUEST_ENV);
25 
26  SECTION( "correctness" ) {
27 
28  // repeat these random tests 10 times
29  GENERATE( range(0,10) );
30 
31  SECTION( "density-matrix" ) {
32 
33  SECTION( "pure" ) {
34 
35  // mat1 = |r1><r1|
36  QVector r1 = getRandomStateVector(NUM_QUBITS);
37  toQureg(mat1, getKetBra(r1,r1));
38 
39  // mat2 = |r2><r2|
40  QVector r2 = getRandomStateVector(NUM_QUBITS);
41  toQureg(mat2, getKetBra(r2,r2));
42 
43  // prod( |r1><r1|, |r2><r2| ) = |<r1|r2>|^2
44  qcomp prod = 0;
45  for (size_t i=0; i<r1.size(); i++)
46  prod += conj(r1[i]) * r2[i];
47  qreal densProd = pow(abs(prod),2);
48 
49  REQUIRE( calcDensityInnerProduct(mat1,mat2) == Approx(densProd) );
50  }
51  SECTION( "mixed" ) {
52 
53  QMatrix ref1 = getRandomDensityMatrix(NUM_QUBITS);
54  QMatrix ref2 = getRandomDensityMatrix(NUM_QUBITS);
55  toQureg(mat1, ref1);
56  toQureg(mat2, ref2);
57 
58  // prod(mat1, mat2) = sum_{ij} conj(mat1_{ij}) * mat2_{ij}
59  qcomp refProd = 0;
60  for (size_t i=0; i<ref1.size(); i++)
61  for (size_t j=0; j<ref1.size(); j++)
62  refProd += conj(ref1[i][j]) * ref2[i][j];
63  REQUIRE( imag(refProd) == Approx(0).margin(REAL_EPS) );
64 
65  REQUIRE( calcDensityInnerProduct(mat1,mat2) == Approx(real(refProd)) );
66 
67  // should be invariant under ordering
68  REQUIRE( calcDensityInnerProduct(mat1,mat2) == Approx(calcDensityInnerProduct(mat2,mat1)) );
69  }
70  SECTION( "unnormalised" ) {
71 
72  // set both to random (non-Hermitian) complex matrices
73  QMatrix ref1 = getRandomQMatrix(1<<NUM_QUBITS);
74  QMatrix ref2 = getRandomQMatrix(1<<NUM_QUBITS);
75  toQureg(mat1, ref1);
76  toQureg(mat2, ref2);
77 
78  // prod(mat1, mat2) = real(sum_{ij} conj(mat1_{ij}) * mat2_{ij})
79  qcomp refProd = 0;
80  for (size_t i=0; i<ref1.size(); i++)
81  for (size_t j=0; j<ref1.size(); j++)
82  refProd += conj(ref1[i][j]) * ref2[i][j];
83 
84  REQUIRE( calcDensityInnerProduct(mat1,mat2) == Approx(real(refProd)) );
85  }
86  }
87  }
88  SECTION( "input validation" ) {
89 
90  SECTION( "dimensions" ) {
91 
92  Qureg mat3 = createDensityQureg(NUM_QUBITS + 1, QUEST_ENV);
93  REQUIRE_THROWS_WITH( calcDensityInnerProduct(mat1,mat3), Contains("Dimensions") && Contains("don't match") );
94  destroyQureg(mat3, QUEST_ENV);
95  }
96  SECTION( "state-vectors" ) {
97 
98  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
99 
100  REQUIRE_THROWS_WITH( calcDensityInnerProduct(mat1,vec), Contains("valid only for density matrices") );
101  REQUIRE_THROWS_WITH( calcDensityInnerProduct(vec,mat1), Contains("valid only for density matrices") );
102  REQUIRE_THROWS_WITH( calcDensityInnerProduct(vec,vec), Contains("valid only for density matrices") );
103 
104  destroyQureg(vec, QUEST_ENV);
105  }
106  }
107  destroyQureg(mat1, QUEST_ENV);
108  destroyQureg(mat2, QUEST_ENV);
109 }
110 
111 
112 
117 TEST_CASE( "calcExpecDiagonalOp", "[calculations]" ) {
118 
119  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
120  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
121  initDebugState(vec);
122  initDebugState(mat);
123  QVector vecRef = toQVector(vec);
124  QMatrix matRef = toQMatrix(mat);
125 
126  SECTION( "correctness" ) {
127 
128  // try 10 random operators
129  GENERATE( range(0,10) );
130 
131  // make a totally random (non-Hermitian) diagonal oeprator
132  DiagonalOp op = createDiagonalOp(NUM_QUBITS, QUEST_ENV);
133  for (long long int i=0; i<op.numElemsPerChunk; i++) {
134  op.real[i] = getRandomReal(-5, 5);
135  op.imag[i] = getRandomReal(-5, 5);
136  }
137  syncDiagonalOp(op);
138 
139  SECTION( "state-vector" ) {
140 
141  /* calcExpecDiagOp calculates <qureg|diag|qureg> */
142 
143  QVector sumRef = toQMatrix(op) * vecRef;
144  qcomp prod = 0;
145  for (size_t i=0; i<vecRef.size(); i++)
146  prod += conj(vecRef[i]) * sumRef[i];
147 
148  Complex res = calcExpecDiagonalOp(vec, op);
149  REQUIRE( res.real == Approx(real(prod)).margin(REAL_EPS) );
150  REQUIRE( res.imag == Approx(imag(prod)).margin(REAL_EPS) );
151  }
152  SECTION( "density-matrix" ) {
153 
154  /* calcExpecDiagOp calculates Trace( diag * qureg ) */
155  matRef = toQMatrix(op) * matRef;
156  qcomp tr = 0;
157  for (size_t i=0; i<matRef.size(); i++)
158  tr += matRef[i][i];
159 
160  Complex res = calcExpecDiagonalOp(mat, op);
161  REQUIRE( res.real == Approx(real(tr)).margin(100*REAL_EPS) );
162  REQUIRE( res.imag == Approx(imag(tr)).margin(100*REAL_EPS) );
163  }
164 
165  destroyDiagonalOp(op, QUEST_ENV);
166  }
167  SECTION( "input validation" ) {
168 
169  SECTION( "mismatching size" ) {
170 
171  DiagonalOp op = createDiagonalOp(NUM_QUBITS + 1, QUEST_ENV);
172 
173  REQUIRE_THROWS_WITH( calcExpecDiagonalOp(vec, op), Contains("equal number of qubits"));
174  REQUIRE_THROWS_WITH( calcExpecDiagonalOp(mat, op), Contains("equal number of qubits"));
175 
176  destroyDiagonalOp(op, QUEST_ENV);
177  }
178  }
179  destroyQureg(vec, QUEST_ENV);
180  destroyQureg(mat, QUEST_ENV);
181 }
182 
183 
184 
189 TEST_CASE( "calcExpecPauliHamil", "[calculations]" ) {
190 
191  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
192  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
193  initDebugState(vec);
194  initDebugState(mat);
195  QVector vecRef = toQVector(vec);
196  QMatrix matRef = toQMatrix(mat);
197 
198  Qureg vecWork = createQureg(NUM_QUBITS, QUEST_ENV);
199  Qureg matWork = createDensityQureg(NUM_QUBITS, QUEST_ENV);
200 
201  SECTION( "correctness" ) {
202 
203  /* it's too expensive to try every possible Pauli configuration, so
204  * we'll try 10 random codes, and for each, random coefficients
205  */
206  GENERATE( range(0,10) );
207 
208  int numTerms = GENERATE( 1, 2, 10, 15 );
209  PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
210  setRandomPauliSum(hamil);
211  QMatrix refHamil = toQMatrix(hamil);
212 
213  SECTION( "state-vector" ) {
214 
215  /* calcExpecPauliHamil calculates <qureg|pauliHum|qureg> */
216 
217  QVector sumRef = refHamil * vecRef;
218  qcomp prod = 0;
219  for (size_t i=0; i<vecRef.size(); i++)
220  prod += conj(vecRef[i]) * sumRef[i];
221  REQUIRE( imag(prod) == Approx(0).margin(10*REAL_EPS) );
222 
223  qreal res = calcExpecPauliHamil(vec, hamil, vecWork);
224  REQUIRE( res == Approx(real(prod)).margin(10*REAL_EPS) );
225  }
226  SECTION( "density-matrix" ) {
227 
228  /* calcExpecPauliHamil calculates Trace( pauliHamil * qureg ) */
229  matRef = refHamil * matRef;
230  qreal tr = 0;
231  for (size_t i=0; i<matRef.size(); i++)
232  tr += real(matRef[i][i]);
233  // (get real, since we start in a non-Hermitian state, hence diagonal isn't real)
234 
235  qreal res = calcExpecPauliHamil(mat, hamil, matWork);
236  REQUIRE( res == Approx(tr).margin(1E2*REAL_EPS) );
237  }
238 
239  destroyPauliHamil(hamil);
240  }
241  SECTION( "input validation" ) {
242 
243  SECTION( "pauli codes" ) {
244 
245  int numTerms = 3;
246  PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
247 
248  // make one pauli code wrong
249  hamil.pauliCodes[GENERATE_COPY( range(0,numTerms*NUM_QUBITS) )] = (pauliOpType) GENERATE( -1, 4 );
250  REQUIRE_THROWS_WITH( calcExpecPauliHamil(vec, hamil, vecWork), Contains("Invalid Pauli code") );
251 
252  destroyPauliHamil(hamil);
253  }
254  SECTION( "workspace type" ) {
255 
256  int numTerms = 1;
257  PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
258 
259  REQUIRE_THROWS_WITH( calcExpecPauliHamil(vec, hamil, mat), Contains("Registers must both be state-vectors or both be density matrices") );
260  REQUIRE_THROWS_WITH( calcExpecPauliHamil(mat, hamil, vec), Contains("Registers must both be state-vectors or both be density matrices") );
261 
262  destroyPauliHamil(hamil);
263  }
264  SECTION( "workspace dimensions" ) {
265 
266  int numTerms = 1;
267  PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
268 
269  Qureg vec2 = createQureg(NUM_QUBITS + 1, QUEST_ENV);
270  REQUIRE_THROWS_WITH( calcExpecPauliHamil(vec, hamil, vec2), Contains("Dimensions") && Contains("don't match") );
271  destroyQureg(vec2, QUEST_ENV);
272 
273  Qureg mat2 = createDensityQureg(NUM_QUBITS + 1, QUEST_ENV);
274  REQUIRE_THROWS_WITH( calcExpecPauliHamil(mat, hamil, mat2), Contains("Dimensions") && Contains("don't match") );
275  destroyQureg(mat2, QUEST_ENV);
276 
277  destroyPauliHamil(hamil);
278  }
279  SECTION( "matching hamiltonian qubits" ) {
280 
281  int numTerms = 1;
282  PauliHamil hamil = createPauliHamil(NUM_QUBITS + 1, numTerms);
283 
284  REQUIRE_THROWS_WITH( calcExpecPauliHamil(vec, hamil, vecWork), Contains("same number of qubits") );
285  REQUIRE_THROWS_WITH( calcExpecPauliHamil(mat, hamil, matWork), Contains("same number of qubits") );
286 
287  destroyPauliHamil(hamil);
288  }
289  }
290  destroyQureg(vec, QUEST_ENV);
291  destroyQureg(mat, QUEST_ENV);
292  destroyQureg(vecWork, QUEST_ENV);
293  destroyQureg(matWork, QUEST_ENV);
294 }
295 
296 
297 
302 TEST_CASE( "calcExpecPauliProd", "[calculations]" ) {
303 
304  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
305  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
306  initDebugState(vec);
307  initDebugState(mat);
308  QVector vecRef = toQVector(vec);
309  QMatrix matRef = toQMatrix(mat);
310 
311  Qureg vecWork = createQureg(NUM_QUBITS, QUEST_ENV);
312  Qureg matWork = createDensityQureg(NUM_QUBITS, QUEST_ENV);
313 
314  SECTION( "correctness" ) {
315 
316  int numTargs = GENERATE( range(1,NUM_QUBITS+1) );
317  int* targs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numTargs) );
318 
319  /* it's too expensive to try ALL Pauli sequences, via
320  * pauliOpType* paulis = GENERATE_COPY( pauliseqs(numTargs) );.
321  * Furthermore, take(10, pauliseqs(numTargs)) will try the same pauli codes.
322  * Hence, we instead opt to repeatedlyrandomly generate pauliseqs
323  */
324  GENERATE( range(0,10) ); // gen 10 random pauli-codes for every targs
325  pauliOpType paulis[numTargs];
326  for (int i=0; i<numTargs; i++)
327  paulis[i] = (pauliOpType) getRandomInt(0,4);
328 
329  // produce a numTargs-big matrix 'pauliProd' by pauli-matrix tensoring
330  QMatrix iMatr{{1,0},{0,1}};
331  QMatrix xMatr{{0,1},{1,0}};
332  QMatrix yMatr{{0,-qcomp(0,1)},{qcomp(0,1),0}};
333  QMatrix zMatr{{1,0},{0,-1}};
334  QMatrix pauliProd{{1}};
335  for (int i=0; i<numTargs; i++) {
336  QMatrix fac;
337  if (paulis[i] == PAULI_I) fac = iMatr;
338  if (paulis[i] == PAULI_X) fac = xMatr;
339  if (paulis[i] == PAULI_Y) fac = yMatr;
340  if (paulis[i] == PAULI_Z) fac = zMatr;
341  pauliProd = getKroneckerProduct(fac, pauliProd);
342  }
343 
344  SECTION( "state-vector" ) {
345 
346  /* calcExpecPauliProd calculates <qureg|pauliProd|qureg> */
347 
348  QVector prodRef = vecRef;
349  applyReferenceOp(prodRef, targs, numTargs, pauliProd);
350  qcomp prod = 0;
351  for (size_t i=0; i<vecRef.size(); i++)
352  prod += conj(vecRef[i]) * prodRef[i];
353  REQUIRE( imag(prod) == Approx(0).margin(REAL_EPS) );
354 
355  qreal res = calcExpecPauliProd(vec, targs, paulis, numTargs, vecWork);
356  REQUIRE( res == Approx(real(prod)).margin(REAL_EPS) );
357  }
358  SECTION( "density-matrix" ) {
359 
360  /* calcExpecPauliProd calculates Trace( pauliProd * qureg ) */
361 
362  // produce (pauliProd * mat)
363  QMatrix fullOp = getFullOperatorMatrix(NULL, 0, targs, numTargs, pauliProd, NUM_QUBITS);
364  matRef = fullOp * matRef;
365 
366  // compute real(trace(pauliProd * mat))
367  qreal tr = 0;
368  for (size_t i=0; i<matRef.size(); i++)
369  tr += real(matRef[i][i]);
370  // (get real, since we start in a non-Hermitian state, hence diagonal isn't real)
371 
372  qreal res = calcExpecPauliProd(mat, targs, paulis, numTargs, matWork);
373  REQUIRE( res == Approx(tr).margin(10*REAL_EPS) );
374  }
375  }
376  SECTION( "input validation" ) {
377 
378  SECTION( "number of targets" ) {
379 
380  int numTargs = GENERATE( -1, 0, NUM_QUBITS+1 );
381  REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, NULL, NULL, numTargs, vecWork), Contains("Invalid number of target") );
382  }
383  SECTION( "target indices" ) {
384 
385  int numTargs = 3;
386  int targs[3] = {0, 1, 2};
387 
388  // make one index wrong
389  targs[GENERATE( range(0,3) )] = GENERATE( -1, NUM_QUBITS );
390  REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, NULL, numTargs, vecWork), Contains("Invalid target qubit") );
391  }
392  SECTION( "repetition in targets" ) {
393 
394  int numTargs = 3;
395  int targs[3] = {0, 1, 1};
396  REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, NULL, numTargs, vecWork), Contains("target qubits must be unique") );
397  }
398  SECTION( "pauli codes" ) {
399 
400  int numTargs = 3;
401  int targs[3] = {0, 1, 2};
402  pauliOpType codes[3] = {PAULI_X, PAULI_Y, PAULI_Z};
403 
404  // make one pauli wrong
405  codes[GENERATE( range(0,3) )] = (pauliOpType) GENERATE( -1, 4 );
406  REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, codes, numTargs, vecWork), Contains("Invalid Pauli code") );
407  }
408  SECTION( "workspace type" ) {
409 
410  int numTargs = 1;
411  int targs[1] = {0};
412  pauliOpType codes[1] = {PAULI_I};
413 
414  REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, codes, numTargs, matWork), Contains("Registers must both be state-vectors or both be density matrices") );
415  REQUIRE_THROWS_WITH( calcExpecPauliProd(mat, targs, codes, numTargs, vecWork), Contains("Registers must both be state-vectors or both be density matrices") );
416  }
417  SECTION( "workspace dimensions" ) {
418 
419  int numTargs = 1;
420  int targs[1] = {0};
421  pauliOpType codes[1] = {PAULI_I};
422 
423  Qureg vec2 = createQureg(NUM_QUBITS + 1, QUEST_ENV);
424  REQUIRE_THROWS_WITH( calcExpecPauliProd(vec, targs, codes, numTargs, vec2), Contains("Dimensions") && Contains("don't match") );
425  destroyQureg(vec2, QUEST_ENV);
426 
427  Qureg mat2 = createDensityQureg(NUM_QUBITS + 1, QUEST_ENV);
428  REQUIRE_THROWS_WITH( calcExpecPauliProd(mat, targs, codes, numTargs, mat2), Contains("Dimensions") && Contains("don't match") );
429  destroyQureg(mat2, QUEST_ENV);
430  }
431  }
432  destroyQureg(vec, QUEST_ENV);
433  destroyQureg(mat, QUEST_ENV);
434  destroyQureg(vecWork, QUEST_ENV);
435  destroyQureg(matWork, QUEST_ENV);
436 }
437 
438 
439 
444 TEST_CASE( "calcExpecPauliSum", "[calculations]" ) {
445 
446  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
447  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
448  initDebugState(vec);
449  initDebugState(mat);
450  QVector vecRef = toQVector(vec);
451  QMatrix matRef = toQMatrix(mat);
452 
453  Qureg vecWork = createQureg(NUM_QUBITS, QUEST_ENV);
454  Qureg matWork = createDensityQureg(NUM_QUBITS, QUEST_ENV);
455 
456  SECTION( "correctness" ) {
457 
458  int numSumTerms = GENERATE( 1, 2, 10, 15 );
459 
460  /* it's too expensive to try every possible Pauli configuration, so
461  * we'll try 10 random codes, and for each, random coefficients
462  */
463  GENERATE( range(0,10) );
464  int totNumCodes = numSumTerms*NUM_QUBITS;
465  pauliOpType paulis[totNumCodes];
466  qreal coeffs[numSumTerms];
467  setRandomPauliSum(coeffs, paulis, NUM_QUBITS, numSumTerms);
468 
469  // produce a numTargs-big matrix 'pauliSum' by pauli-matrix tensoring and summing
470  QMatrix pauliSum = toQMatrix(coeffs, paulis, NUM_QUBITS, numSumTerms);
471 
472  SECTION( "state-vector" ) {
473 
474  /* calcExpecPauliSum calculates <qureg|pauliSum|qureg> */
475 
476  QVector sumRef = pauliSum * vecRef;
477  qcomp prod = 0;
478  for (size_t i=0; i<vecRef.size(); i++)
479  prod += conj(vecRef[i]) * sumRef[i];
480  REQUIRE( imag(prod) == Approx(0).margin(10*REAL_EPS) );
481 
482  qreal res = calcExpecPauliSum(vec, paulis, coeffs, numSumTerms, vecWork);
483  REQUIRE( res == Approx(real(prod)).margin(10*REAL_EPS) );
484  }
485  SECTION( "density-matrix" ) {
486 
487  /* calcExpecPauliSum calculates Trace( pauliSum * qureg ) */
488  matRef = pauliSum * matRef;
489  qreal tr = 0;
490  for (size_t i=0; i<matRef.size(); i++)
491  tr += real(matRef[i][i]);
492  // (get real, since we start in a non-Hermitian state, hence diagonal isn't real)
493 
494  qreal res = calcExpecPauliSum(mat, paulis, coeffs, numSumTerms, matWork);
495  REQUIRE( res == Approx(tr).margin(1E2*REAL_EPS) );
496  }
497  }
498  SECTION( "input validation" ) {
499 
500  SECTION( "number of sum terms" ) {
501 
502  int numSumTerms = GENERATE( -1, 0 );
503  REQUIRE_THROWS_WITH( calcExpecPauliSum(vec, NULL, NULL, numSumTerms, vecWork), Contains("Invalid number of terms in the Pauli sum") );
504  }
505  SECTION( "pauli codes" ) {
506 
507  // make valid params
508  int numSumTerms = 3;
509  qreal coeffs[numSumTerms];
510  pauliOpType codes[numSumTerms*NUM_QUBITS];
511  for (int i=0; i<numSumTerms*NUM_QUBITS; i++)
512  codes[i] = PAULI_I;
513 
514  // make one pauli wrong
515  codes[GENERATE_COPY( range(0,numSumTerms*NUM_QUBITS) )] = (pauliOpType) GENERATE( -1, 4 );
516  REQUIRE_THROWS_WITH( calcExpecPauliSum(vec, codes, coeffs, numSumTerms, vecWork), Contains("Invalid Pauli code") );
517  }
518  SECTION( "workspace type" ) {
519 
520  // make valid params
521  int numSumTerms = 1;
522  qreal coeffs[1] = {0};
523  pauliOpType codes[NUM_QUBITS];
524  for (int i=0; i<NUM_QUBITS; i++)
525  codes[i] = PAULI_I;
526 
527  REQUIRE_THROWS_WITH( calcExpecPauliSum(vec, codes, coeffs, numSumTerms, mat), Contains("Registers must both be state-vectors or both be density matrices") );
528  REQUIRE_THROWS_WITH( calcExpecPauliSum(mat, codes, coeffs, numSumTerms, vec), Contains("Registers must both be state-vectors or both be density matrices") );
529  }
530  SECTION( "workspace dimensions" ) {
531 
532  // make valid params
533  int numSumTerms = 1;
534  qreal coeffs[1] = {0};
535  pauliOpType codes[NUM_QUBITS];
536  for (int i=0; i<NUM_QUBITS; i++)
537  codes[i] = PAULI_I;
538 
539  Qureg vec2 = createQureg(NUM_QUBITS + 1, QUEST_ENV);
540  REQUIRE_THROWS_WITH( calcExpecPauliSum(vec, codes, coeffs, numSumTerms, vec2), Contains("Dimensions") && Contains("don't match") );
541  destroyQureg(vec2, QUEST_ENV);
542 
543  Qureg mat2 = createDensityQureg(NUM_QUBITS + 1, QUEST_ENV);
544  REQUIRE_THROWS_WITH( calcExpecPauliSum(mat, codes, coeffs, numSumTerms, mat2), Contains("Dimensions") && Contains("don't match") );
545  destroyQureg(mat2, QUEST_ENV);
546  }
547  }
548  destroyQureg(vec, QUEST_ENV);
549  destroyQureg(mat, QUEST_ENV);
550  destroyQureg(vecWork, QUEST_ENV);
551  destroyQureg(matWork, QUEST_ENV);
552 }
553 
554 
555 
560 TEST_CASE( "calcFidelity", "[calculations]" ) {
561 
562  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
563  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
564  Qureg pure = createQureg(NUM_QUBITS, QUEST_ENV);
565 
566  SECTION( "correctness" ) {
567 
568  // repeat the below random tests 10 times
569  GENERATE( range(0,10) );
570 
571  SECTION( "state-vector" ) {
572 
573  /* calcFidelity computes |<vec|pure>|^2 */
574 
575  SECTION( "normalised" ) {
576 
577  // random L2 vectors
578  QVector vecRef = getRandomStateVector(NUM_QUBITS);
579  QVector pureRef = getRandomStateVector(NUM_QUBITS);
580  toQureg(vec, vecRef);
581  toQureg(pure, pureRef);
582 
583  // |<vec|vec>|^2 = |1|^2 = 1
584  REQUIRE( calcFidelity(vec,vec) == Approx(1) );
585 
586  // |<vec|pure>|^2 = |sum_j conj(vec_j) * pure_j|^2
587  qcomp dotProd = 0;
588  for (size_t i=0; i<vecRef.size(); i++)
589  dotProd += conj(vecRef[i]) * pureRef[i];
590  qreal refFid = pow(abs(dotProd), 2);
591 
592  REQUIRE( calcFidelity(vec,pure) == Approx(refFid) );
593  }
594  SECTION( "unnormalised" ) {
595 
596  // random unnormalised vectors
597  QVector vecRef = getRandomQVector(1<<NUM_QUBITS);
598  QVector pureRef = getRandomQVector(1<<NUM_QUBITS);
599  toQureg(vec, vecRef);
600  toQureg(pure, pureRef);
601 
602  // Let nv be magnitude of vec, hence |unit-vec> = 1/sqrt(nv)|vec>
603  qreal nv = 0;
604  for (size_t i=0; i<vecRef.size(); i++)
605  nv += pow(abs(vecRef[i]), 2);
606  // then <vec|vec> = sqrt(nv)*sqrt(nv) <unit-vec|unit-vec> = nv,
607  // hence |<vec|vec>|^2 = nv*nv
608  REQUIRE( calcFidelity(vec,vec) == Approx( nv*nv ) );
609 
610  qcomp dotProd = 0;
611  for (size_t i=0; i<vecRef.size(); i++)
612  dotProd += conj(vecRef[i]) * pureRef[i];
613  qreal refFid = pow(abs(dotProd), 2);
614 
615  REQUIRE( calcFidelity(vec,pure) == Approx(refFid) );
616  }
617  }
618  SECTION( "density-matrix" ) {
619 
620  /* calcFidelity computes <pure|mat|pure> */
621 
622  SECTION( "pure" ) {
623 
624  QVector pureRef = getRandomStateVector(NUM_QUBITS);
625  toQureg(pure, pureRef);
626 
627  // test when density matrix is the same pure state
628  QMatrix matRef = getKetBra(pureRef, pureRef);
629  toQureg(mat, matRef);
630  REQUIRE( calcFidelity(mat,pure) == Approx(1) );
631 
632  // test when density matrix is a random pure state
633  QVector r1 = getRandomStateVector(NUM_QUBITS);
634  matRef = getKetBra(r1, r1); // actually pure |r1><r1|
635  toQureg(mat, matRef);
636 
637  // <pure|r1><r1|pure> = |<r1|pure>|^2 = |sum_j conj(r1_j) * pure_j|^2
638  qcomp dotProd = 0;
639  for (size_t i=0; i<r1.size(); i++)
640  dotProd += conj(r1[i]) * pureRef[i];
641  qreal refFid = pow(abs(dotProd), 2);
642 
643  REQUIRE( calcFidelity(mat,pure) == Approx(refFid) );
644  }
645  SECTION( "mixed" ) {
646 
647  QVector pureRef = getRandomStateVector(NUM_QUBITS);
648  toQureg(pure, pureRef);
649 
650  // test when density matrix is mixed
651  QMatrix matRef = getRandomDensityMatrix(NUM_QUBITS);
652  toQureg(mat, matRef);
653 
654  // <pure|mat|pure> = <pure| (Mat|pure>)
655  QVector rhs = matRef * pureRef;
656  qcomp dotProd = 0;
657  for (size_t i=0; i<rhs.size(); i++)
658  dotProd += conj(pureRef[i]) * rhs[i];
659 
660  REQUIRE( imag(dotProd) == Approx(0).margin(REAL_EPS) );
661  REQUIRE( calcFidelity(mat,pure) == Approx(real(dotProd)) );
662  }
663  SECTION( "unnormalised" ) {
664 
665  // test when both density matrix and pure state are unnormalised
666  QVector pureRef = getRandomQVector(1<<NUM_QUBITS);
667  QMatrix matRef = getRandomQMatrix(1<<NUM_QUBITS);
668  toQureg(pure, pureRef);
669  toQureg(mat, matRef);
670 
671  // real[ <pure|mat|pure> ] = real[ <pure| (Mat|pure>) ]
672  QVector rhs = matRef * pureRef;
673  qcomp dotProd = 0;
674  for (size_t i=0; i<rhs.size(); i++)
675  dotProd += conj(pureRef[i]) * rhs[i];
676 
677  REQUIRE( calcFidelity(mat,pure) == Approx(real(dotProd)) );
678  }
679  }
680  }
681  SECTION( "input validation" ) {
682 
683  SECTION( "dimensions" ) {
684 
685  // two state-vectors
686  Qureg vec2 = createQureg(vec.numQubitsRepresented + 1, QUEST_ENV);
687  REQUIRE_THROWS_WITH( calcFidelity(vec2,vec), Contains("Dimensions") && Contains("don't match") );
688  destroyQureg(vec2, QUEST_ENV);
689 
690  // density-matrix and state-vector
691  Qureg mat2 = createDensityQureg(vec.numQubitsRepresented + 1, QUEST_ENV);
692  REQUIRE_THROWS_WITH( calcFidelity(mat2,vec), Contains("Dimensions") && Contains("don't match") );
693  destroyQureg(mat2, QUEST_ENV);
694  }
695  SECTION( "density-matrices" ) {
696 
697  // two mixed statess
698  REQUIRE_THROWS_WITH( calcFidelity(mat,mat), Contains("Second argument must be a state-vector") );
699  }
700  }
701  destroyQureg(vec, QUEST_ENV);
702  destroyQureg(mat, QUEST_ENV);
703  destroyQureg(pure, QUEST_ENV);
704 }
705 
706 
707 
712 TEST_CASE( "calcHilbertSchmidtDistance", "[calculations]" ) {
713 
714  Qureg mat1 = createDensityQureg(NUM_QUBITS, QUEST_ENV);
715  Qureg mat2 = createDensityQureg(NUM_QUBITS, QUEST_ENV);
716 
717  SECTION( "correctness" ) {
718 
719  // perform these random tests 10 times
720  GENERATE( range(0,10) );
721 
722  SECTION( "density-matrix" ) {
723 
724  SECTION( "pure" ) {
725 
726  // create random |r1><r1| and |r2><r2| states
727  QVector r1 = getRandomStateVector(NUM_QUBITS);
728  QMatrix m1 = getKetBra(r1,r1);
729  toQureg(mat1, m1);
730  QVector r2 = getRandomStateVector(NUM_QUBITS);
731  QMatrix m2 = getKetBra(r2,r2);
732  toQureg(mat2, m2);
733 
734  // Tr{ (a-b)(a-b)^dagger } = sum_{ij} |a_{ij} - b_{ij}|^2
735  qreal tr = 0;
736  for (size_t i=0; i<m1.size(); i++)
737  for (size_t j=0; j<m1.size(); j++)
738  tr += pow(abs(m1[i][j] - m2[i][j]), 2);
739 
740  qreal res = calcHilbertSchmidtDistance(mat1, mat2);
741  REQUIRE( res == Approx(sqrt(tr)) );
742 
743  }
744  SECTION( "normalised" ) {
745 
746  QMatrix ref1 = getRandomDensityMatrix(NUM_QUBITS);
747  QMatrix ref2 = getRandomDensityMatrix(NUM_QUBITS);
748  toQureg(mat1, ref1);
749  toQureg(mat2, ref2);
750 
751  // Tr{ (a-b)(a-b)^dagger } = sum_{ij} |a_{ij} - b_{ij}|^2
752  qreal tr = 0;
753  for (size_t i=0; i<ref1.size(); i++)
754  for (size_t j=0; j<ref1.size(); j++)
755  tr += pow(abs(ref1[i][j] - ref2[i][j]), 2);
756 
757  qreal res = calcHilbertSchmidtDistance(mat1, mat2);
758  REQUIRE( res == Approx(sqrt(tr)) );
759  }
760  SECTION( "unnormalised" ) {
761 
762  // mat1 and mat2 are both random matrices
763  QMatrix ref1 = getRandomQMatrix(1<<NUM_QUBITS);
764  QMatrix ref2 = getRandomQMatrix(1<<NUM_QUBITS);
765  toQureg(mat1, ref1);
766  toQureg(mat2, ref2);
767 
768  // Tr{ (a-b)(a-b)^dagger } = sum_{ij} |a_{ij} - b_{ij}|^2
769  qreal tr = 0;
770  for (size_t i=0; i<ref1.size(); i++)
771  for (size_t j=0; j<ref1.size(); j++)
772  tr += pow(abs(ref1[i][j] - ref2[i][j]), 2);
773 
774  qreal res = calcHilbertSchmidtDistance(mat1, mat2);
775  REQUIRE( res == Approx(sqrt(tr)) );
776  }
777  }
778  }
779  SECTION( "input validation") {
780 
781  SECTION( "dimensions" ) {
782 
783  Qureg mat3 = createDensityQureg(NUM_QUBITS + 1, QUEST_ENV);
784  REQUIRE_THROWS_WITH( calcHilbertSchmidtDistance(mat1,mat3), Contains("Dimensions") && Contains("don't match") );
785  destroyQureg(mat3, QUEST_ENV);
786  }
787  SECTION( "state-vector" ) {
788 
789  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
790 
791  REQUIRE_THROWS_WITH( calcHilbertSchmidtDistance(vec,mat1), Contains("valid only for density matrices") );
792  REQUIRE_THROWS_WITH( calcHilbertSchmidtDistance(mat1,vec), Contains("valid only for density matrices") );
793  REQUIRE_THROWS_WITH( calcHilbertSchmidtDistance(vec,vec), Contains("valid only for density matrices") );
794 
795  destroyQureg(vec, QUEST_ENV);
796  }
797  }
798  destroyQureg(mat1, QUEST_ENV);
799  destroyQureg(mat2, QUEST_ENV);
800 }
801 
802 
803 
808 TEST_CASE( "calcInnerProduct", "[calculations]" ) {
809 
810  Qureg vec1 = createQureg(NUM_QUBITS, QUEST_ENV);
811  Qureg vec2 = createQureg(NUM_QUBITS, QUEST_ENV);
812 
813  SECTION( "correctness" ) {
814 
815  // perform these random tests 10 times
816  GENERATE( range(0,10) );
817 
818  SECTION( "state-vector" ) {
819 
820  SECTION( "normalised" ) {
821 
822  // <r1|r2> = sum_j conj(r1_j) * r2_j
823  QVector r1 = getRandomStateVector(NUM_QUBITS);
824  QVector r2 = getRandomStateVector(NUM_QUBITS);
825  qcomp prod = 0;
826  for (size_t i=0; i<r1.size(); i++)
827  prod += conj(r1[i]) * r2[i];
828 
829  toQureg(vec1, r1);
830  toQureg(vec2, r2);
831  Complex res = calcInnerProduct(vec1,vec2);
832 
833  REQUIRE( res.real == Approx(real(prod)) );
834  REQUIRE( res.imag == Approx(imag(prod)) );
835  }
836  SECTION( "unnormalised" ) {
837 
838  // <r1|r2> = sum_j conj(r1_j) * r2_j
839  QVector r1 = getRandomQVector(1<<NUM_QUBITS);
840  QVector r2 = getRandomQVector(1<<NUM_QUBITS);
841  qcomp prod = 0;
842  for (size_t i=0; i<r1.size(); i++)
843  prod += conj(r1[i]) * r2[i];
844 
845  toQureg(vec1, r1);
846  toQureg(vec2, r2);
847  Complex res = calcInnerProduct(vec1,vec2);
848 
849  REQUIRE( res.real == Approx(real(prod)) );
850  REQUIRE( res.imag == Approx(imag(prod)) );
851  }
852  }
853  }
854  SECTION( "input validation" ) {
855 
856  SECTION( "dimensions" ) {
857 
858  Qureg vec3 = createQureg(NUM_QUBITS + 1, QUEST_ENV);
859  REQUIRE_THROWS_WITH( calcInnerProduct(vec1,vec3), Contains("Dimensions") && Contains("don't match") );
860  destroyQureg(vec3, QUEST_ENV);
861  }
862  SECTION( "density-matrix" ) {
863 
864  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
865 
866  REQUIRE_THROWS_WITH( calcInnerProduct(vec1,mat), Contains("valid only for state-vectors") );
867  REQUIRE_THROWS_WITH( calcInnerProduct(mat,vec1), Contains("valid only for state-vectors") );
868  REQUIRE_THROWS_WITH( calcInnerProduct(mat,mat), Contains("valid only for state-vectors") );
869 
870  destroyQureg(mat, QUEST_ENV);
871  }
872  }
873  destroyQureg(vec1, QUEST_ENV);
874  destroyQureg(vec2, QUEST_ENV);
875 }
876 
877 
878 
883 TEST_CASE( "calcProbOfAllOutcomes", "[calculations]" ) {
884 
885  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
886  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
887 
888  SECTION( "correctness" ) {
889 
890  // generate all possible qubit arrangements
891  int numQubits = GENERATE_COPY( range(1,NUM_QUBITS+1) );
892  int* qubits = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numQubits) );
893 
894  int numOutcomes = 1<<numQubits;
895  qreal probs[numOutcomes];
896  QVector refProbs = QVector(numOutcomes);
897 
898  SECTION( "state-vector" ) {
899 
900  SECTION( "normalised" ) {
901 
902  QVector ref = getRandomStateVector(NUM_QUBITS);
903  toQureg(vec, ref);
904 
905  // prob is sum of |amp|^2 of basis states which encode outcome
906  for (size_t i=0; i<ref.size(); i++) {
907  int outcome = 0;
908  for (int q=0; q<numQubits; q++) {
909  int bit = (i >> qubits[q]) & 1;
910  outcome += bit * (1 << q);
911  }
912  refProbs[outcome] += pow(abs(ref[i]), 2);
913  }
914 
915  calcProbOfAllOutcomes(probs, vec, qubits, numQubits);
916  REQUIRE( areEqual(refProbs, probs) );
917  }
918  SECTION( "unnormalised" ) {
919 
920  QVector ref = getRandomQVector(1<<NUM_QUBITS);
921  toQureg(vec, ref);
922 
923  // prob is sum of |amp|^2 of basis states which encode outcome
924  for (size_t i=0; i<ref.size(); i++) {
925  int outcome = 0;
926  for (int q=0; q<numQubits; q++) {
927  int bit = (i >> qubits[q]) & 1;
928  outcome += bit * (1 << q);
929  }
930  refProbs[outcome] += pow(abs(ref[i]), 2);
931  }
932 
933  calcProbOfAllOutcomes(probs, vec, qubits, numQubits);
934  REQUIRE( areEqual(refProbs, probs) );
935  }
936  }
937  SECTION( "density-matrix" ) {
938 
939  SECTION( "normalised" ) {
940 
941  QMatrix ref = getRandomDensityMatrix(NUM_QUBITS);
942  toQureg(mat, ref);
943 
944  // prob is sum of diagonals which encode outcome
945  for (size_t i=0; i<ref.size(); i++) {
946  int outcome = 0;
947  for (int q=0; q<numQubits; q++) {
948  int bit = (i >> qubits[q]) & 1;
949  outcome += bit * (1 << q);
950  }
951  refProbs[outcome] += real(ref[i][i]);
952  }
953 
954  calcProbOfAllOutcomes(probs, mat, qubits, numQubits);
955  REQUIRE( areEqual(refProbs, probs) );
956  }
957  SECTION( "unnormalised" ) {
958 
959  QMatrix ref = getRandomQMatrix(1<<NUM_QUBITS);
960  toQureg(mat, ref);
961 
962  // prob is sum of diagonals which encode outcome
963  for (size_t i=0; i<ref.size(); i++) {
964  int outcome = 0;
965  for (int q=0; q<numQubits; q++) {
966  int bit = (i >> qubits[q]) & 1;
967  outcome += bit * (1 << q);
968  }
969  refProbs[outcome] += real(ref[i][i]);
970  }
971 
972  calcProbOfAllOutcomes(probs, mat, qubits, numQubits);
973  REQUIRE( areEqual(refProbs, probs) );
974  }
975  }
976  }
977  SECTION( "input validation" ) {
978 
979  int numQubits = 3;
980  int qubits[] = {0, 1, 2};
981  qreal probs[8];
982 
983  SECTION( "number of qubits" ) {
984 
985  numQubits = GENERATE( -1, 0, NUM_QUBITS+1 );
986  REQUIRE_THROWS_WITH( calcProbOfAllOutcomes(probs, mat, qubits, numQubits), Contains("Invalid number of target qubits") );
987  }
988  SECTION( "qubit indices" ) {
989 
990  qubits[GENERATE_COPY(range(0,numQubits))] = GENERATE( -1, NUM_QUBITS );
991  REQUIRE_THROWS_WITH( calcProbOfAllOutcomes(probs, mat, qubits, numQubits), Contains("Invalid target qubit") );
992  }
993  SECTION( "repetition of qubits" ) {
994 
995  qubits[GENERATE_COPY(1,2)] = qubits[0];
996  REQUIRE_THROWS_WITH( calcProbOfAllOutcomes(probs, mat, qubits, numQubits), Contains("qubits must be unique") );
997  }
998  }
999  destroyQureg(vec, QUEST_ENV);
1000  destroyQureg(mat, QUEST_ENV);
1001 }
1002 
1003 
1004 
1005 
1010 TEST_CASE( "calcProbOfOutcome", "[calculations]" ) {
1011 
1012  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
1013  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1014 
1015  SECTION( "correctness" ) {
1016 
1017  int target = GENERATE( range(0,NUM_QUBITS) );
1018  int outcome = GENERATE( 0, 1 );
1019 
1020  SECTION( "state-vector" ) {
1021 
1022  SECTION( "normalised" ) {
1023 
1024  QVector ref = getRandomStateVector(NUM_QUBITS);
1025  toQureg(vec, ref);
1026 
1027  // prob is sum of |amp|^2 of amplitudes where target bit is outcome
1028  qreal prob = 0;
1029  for (size_t ind=0; ind<ref.size(); ind++) {
1030  int bit = (ind >> target) & 1; // target-th bit
1031  if (bit == outcome)
1032  prob += pow(abs(ref[ind]), 2);
1033  }
1034 
1035  REQUIRE( calcProbOfOutcome(vec, target, outcome) == Approx(prob) );
1036  }
1037  SECTION( "unnormalised" ) {
1038 
1039  QVector ref = getRandomQVector(1<<NUM_QUBITS);
1040  toQureg(vec, ref);
1041 
1042  // prob is (sum of |amp|^2 of amplitudes where target bit is zero)
1043  // or 1 - (this) if outcome == 1
1044  qreal prob = 0;
1045  for (size_t ind=0; ind<ref.size(); ind++) {
1046  int bit = (ind >> target) & 1; // target-th bit
1047  if (bit == 0)
1048  prob += pow(abs(ref[ind]), 2);
1049  }
1050  if (outcome == 1)
1051  prob = 1 - prob;
1052 
1053  REQUIRE( calcProbOfOutcome(vec, target, outcome) == Approx(prob) );
1054  }
1055  }
1056  SECTION( "density-matrix" ) {
1057 
1058  SECTION( "pure" ) {
1059 
1060  // set mat to a random |r><r|
1061  QVector ref = getRandomStateVector(NUM_QUBITS);
1062  toQureg(mat, getKetBra(ref, ref));
1063 
1064  // calc prob of the state-vector
1065  qreal prob = 0;
1066  for (size_t ind=0; ind<ref.size(); ind++) {
1067  int bit = (ind >> target) & 1; // target-th bit
1068  if (bit == outcome)
1069  prob += pow(abs(ref[ind]), 2);
1070  }
1071 
1072  REQUIRE( calcProbOfOutcome(mat, target, outcome) == Approx(prob) );
1073  }
1074  SECTION( "mixed" ) {
1075 
1076  QMatrix ref = getRandomDensityMatrix(NUM_QUBITS);
1077  toQureg(mat, ref);
1078 
1079  // prob is sum of diagonal amps (should be real) where target bit is outcome
1080  qcomp tr = 0;
1081  for (size_t ind=0; ind<ref.size(); ind++) {
1082  int bit = (ind >> target) & 1; // target-th bit
1083  if (bit == outcome)
1084  tr += ref[ind][ind];
1085  }
1086  REQUIRE( imag(tr) == Approx(0).margin(REAL_EPS) );
1087 
1088  REQUIRE( calcProbOfOutcome(mat, target, outcome) == Approx(real(tr)) );
1089  }
1090  SECTION( "unnormalised" ) {
1091 
1092  QMatrix ref = getRandomQMatrix(1<<NUM_QUBITS);
1093  toQureg(mat, ref);
1094 
1095  // prob is (sum of real of diagonal amps where target bit is outcome)
1096  // or 1 - (this) if outcome == 1
1097  qreal tr = 0;
1098  for (size_t ind=0; ind<ref.size(); ind++) {
1099  int bit = (ind >> target) & 1; // target-th bit
1100  if (bit == 0)
1101  tr += real(ref[ind][ind]);
1102  }
1103  if (outcome == 1)
1104  tr = 1 - tr;
1105 
1106  REQUIRE( calcProbOfOutcome(mat, target, outcome) == Approx(tr) );
1107  }
1108  }
1109  }
1110  SECTION( "input validation" ) {
1111 
1112  SECTION( "qubit indices" ) {
1113 
1114  int target = GENERATE( -1, NUM_QUBITS );
1115  REQUIRE_THROWS_WITH( calcProbOfOutcome(vec, target, 0), Contains("Invalid target qubit") );
1116  }
1117  SECTION( "outcome value" ) {
1118 
1119  int outcome = GENERATE( -1, 2 );
1120  REQUIRE_THROWS_WITH( calcProbOfOutcome(vec, 0, outcome), Contains("Invalid measurement outcome") );
1121  }
1122  }
1123  destroyQureg(vec, QUEST_ENV);
1124  destroyQureg(mat, QUEST_ENV);
1125 }
1126 
1127 
1128 
1133 TEST_CASE( "calcPurity", "[calculations]" ) {
1134 
1135  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1136 
1137  SECTION( "correctness" ) {
1138 
1139  // perform the following random tests 10 times
1140  GENERATE( range(1,10) );
1141 
1142  SECTION( "density-matrix" ) {
1143 
1144  SECTION( "pure" ) {
1145 
1146  // pure states have unity purity
1147  initZeroState(mat);
1148  REQUIRE( calcPurity(mat) == 1 );
1149 
1150  // (try also a pure random L2-vector)
1151  QVector r1 = getRandomStateVector(NUM_QUBITS); // |r>
1152  QMatrix m1 = getKetBra(r1, r1); // |r><r|
1153  toQureg(mat, m1);
1154  REQUIRE( calcPurity(mat) == Approx(1) );
1155 
1156  }
1157  SECTION( "mixed" ) {
1158 
1159  // mixed states have 1/2^N < purity < 1
1160  QMatrix ref = getRandomDensityMatrix(NUM_QUBITS);
1161  toQureg(mat, ref);
1162  qreal purity = calcPurity(mat);
1163  REQUIRE( purity < 1 );
1164  REQUIRE( purity >= 1/pow(2.,NUM_QUBITS) );
1165 
1166  // compare to Tr(rho^2)
1167  QMatrix prod = ref*ref;
1168  qreal tr = 0;
1169  for (size_t i=0; i<prod.size(); i++)
1170  tr += real(prod[i][i]);
1171  REQUIRE( purity == Approx(tr) );
1172  }
1173  SECTION( "unnormalised" ) {
1174 
1175  // unphysical states give sum_{ij} |rho_ij|^2
1176  QMatrix ref = getRandomQMatrix(1<<NUM_QUBITS);
1177  qreal tot = 0;
1178  for (size_t i=0; i<ref.size(); i++)
1179  for (size_t j=0; j<ref.size(); j++)
1180  tot += pow(abs(ref[i][j]), 2);
1181 
1182  toQureg(mat, ref);
1183  REQUIRE( calcPurity(mat) == Approx(tot) );
1184  }
1185  }
1186  }
1187  SECTION( "input validation" ) {
1188 
1189  SECTION( "state-vector" ) {
1190 
1191  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
1192  REQUIRE_THROWS_WITH( calcPurity(vec), Contains("valid only for density matrices") );
1193  destroyQureg(vec, QUEST_ENV);
1194  }
1195  }
1196  destroyQureg(mat, QUEST_ENV);
1197 }
1198 
1199 
1200 
1205 TEST_CASE( "calcTotalProb", "[calculations]" ) {
1206 
1207  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
1208  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1209 
1210  SECTION( "correctness" ) {
1211 
1212  SECTION( "state-vector" ) {
1213 
1214  // normalised: prob(vec) = 1
1215  initPlusState(vec);
1216  REQUIRE( calcTotalProb(vec) == Approx(1) );
1217 
1218  // zero norm: prob(vec) = 0
1219  initBlankState(vec);
1220  REQUIRE( calcTotalProb(vec) == 0 );
1221 
1222  // random L2 state: prob(vec) = 1
1223  toQureg(vec, getRandomStateVector(NUM_QUBITS));
1224  REQUIRE( calcTotalProb(vec) == Approx(1) );
1225 
1226  // unnormalised: prob(vec) = sum_i |vec_i|^2
1227  initDebugState(vec);
1228  QVector ref = toQVector(vec);
1229  qreal refProb = 0;
1230  for (size_t i=0; i<ref.size(); i++)
1231  refProb += pow(abs(ref[i]), 2);
1232  REQUIRE( calcTotalProb(vec) == Approx(refProb) );
1233  }
1234  SECTION( "density-matrix" ) {
1235 
1236  // normalised: prob(mat) = 1
1237  initPlusState(mat);
1238  REQUIRE( calcTotalProb(mat) == Approx(1) );
1239 
1240  // zero norm: prob(mat) = 0
1241  initBlankState(mat);
1242  REQUIRE( calcTotalProb(mat) == 0 );
1243 
1244  // random density matrix: prob(mat) = 1
1245  toQureg(mat, getRandomDensityMatrix(NUM_QUBITS));
1246  REQUIRE( calcTotalProb(mat) == Approx(1) );
1247 
1248  // unnormalised: prob(mat) = sum_i real(mat_{ii})
1249  initDebugState(mat);
1250  QMatrix ref = toQMatrix(mat);
1251  qreal refProb = 0;
1252  for (size_t i=0; i<ref.size(); i++)
1253  refProb += real(ref[i][i]);
1254  REQUIRE( calcTotalProb(mat) == Approx(refProb) );
1255  }
1256  }
1257  SECTION( "input validation" ) {
1258 
1259  // no validation
1260  SUCCEED();
1261  }
1262  destroyQureg(vec, QUEST_ENV);
1263  destroyQureg(mat, QUEST_ENV);
1264 }
1265 
1266 
1267 
1272 TEST_CASE( "getAmp", "[calculations]" ) {
1273 
1274  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
1275 
1276  SECTION( "correctness" ) {
1277 
1278  SECTION( "state-vector" ) {
1279 
1280  initDebugState(vec);
1281  QVector ref = toQVector(vec);
1282 
1283  int ind = GENERATE( range(0,1<<NUM_QUBITS) );
1284  Complex amp = getAmp(vec,ind);
1285  REQUIRE( fromComplex(amp) == ref[ind] );
1286  }
1287  }
1288  SECTION( "input validation" ) {
1289 
1290  SECTION( "state index" ) {
1291 
1292  int ind = GENERATE( -1, 1<<NUM_QUBITS );
1293  REQUIRE_THROWS_WITH( getAmp(vec,ind), Contains("Invalid amplitude index") );
1294  }
1295  SECTION( "density-matrix" ) {
1296 
1297  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1298  REQUIRE_THROWS_WITH( getAmp(mat,0), Contains("valid only for state-vectors") );
1299  destroyQureg(mat, QUEST_ENV);
1300  }
1301  }
1302  destroyQureg(vec, QUEST_ENV);
1303 }
1304 
1305 
1306 
1311 TEST_CASE( "getDensityAmp", "[calculations]" ) {
1312 
1313  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1314 
1315  SECTION( "correctness" ) {
1316 
1317  SECTION( "density-matrix" ) {
1318 
1319  initDebugState(mat);
1320  QMatrix ref = toQMatrix(mat);
1321 
1322  int row = GENERATE( range(0,1<<NUM_QUBITS) );
1323  int col = GENERATE( range(0,1<<NUM_QUBITS) );
1324 
1325  Complex amp = getDensityAmp(mat,row,col);
1326  REQUIRE( fromComplex(amp) == ref[row][col] );
1327  }
1328  }
1329  SECTION( "input validation" ) {
1330 
1331  SECTION( "state index" ) {
1332 
1333  int ind = GENERATE( -1, 1<<NUM_QUBITS );
1334  REQUIRE_THROWS_WITH( getDensityAmp(mat,ind,0), Contains("Invalid amplitude index") );
1335  REQUIRE_THROWS_WITH( getDensityAmp(mat,0,ind), Contains("Invalid amplitude index") );
1336 
1337  }
1338  SECTION( "state-vector" ) {
1339 
1340  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
1341  REQUIRE_THROWS_WITH( getDensityAmp(vec,0,0), Contains("valid only for density matrices") );
1342  destroyQureg(vec, QUEST_ENV);
1343  }
1344  }
1345  destroyQureg(mat, QUEST_ENV);
1346 }
1347 
1348 
1349 
1354 TEST_CASE( "getImagAmp", "[calculations]" ) {
1355 
1356  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
1357 
1358  SECTION( "correctness" ) {
1359 
1360  SECTION( "state-vector" ) {
1361 
1362  initDebugState(vec);
1363  QVector ref = toQVector(vec);
1364 
1365  int ind = GENERATE( range(0,1<<NUM_QUBITS) );
1366  REQUIRE( getImagAmp(vec,ind) == imag(ref[ind]) );
1367  }
1368  }
1369  SECTION( "input validation" ) {
1370 
1371  SECTION( "state index" ) {
1372 
1373  int ind = GENERATE( -1, 1<<NUM_QUBITS );
1374  REQUIRE_THROWS_WITH( getImagAmp(vec,ind), Contains("Invalid amplitude index") );
1375  }
1376  SECTION( "density-matrix" ) {
1377 
1378  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1379  REQUIRE_THROWS_WITH( getImagAmp(mat,0), Contains("valid only for state-vectors") );
1380  destroyQureg(mat, QUEST_ENV);
1381  }
1382  }
1383  destroyQureg(vec, QUEST_ENV);
1384 }
1385 
1386 
1387 
1392 TEST_CASE( "getNumAmps", "[calculations]" ) {
1393 
1394  SECTION( "correctness" ) {
1395 
1396  // test >= NUM_QUBITS so as not to limit distribution size
1397  int numQb = GENERATE( range(NUM_QUBITS, NUM_QUBITS+10) );
1398 
1399  SECTION( "state-vector" ) {
1400 
1401  Qureg vec = createQureg(numQb, QUEST_ENV);
1402  REQUIRE( getNumAmps(vec) == (1<<numQb) );
1403  destroyQureg(vec, QUEST_ENV);
1404  }
1405  }
1406  SECTION( "input validation" ) {
1407 
1408  SECTION( "density-matrix" ) {
1409  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1410  REQUIRE_THROWS_WITH( getNumAmps(mat), Contains("valid only for state-vectors") );
1411  destroyQureg(mat, QUEST_ENV);
1412  }
1413  }
1414 }
1415 
1416 
1417 
1422 TEST_CASE( "getNumQubits", "[calculations]" ) {
1423 
1424  SECTION( "correctness" ) {
1425 
1426  // test >= NUM_QUBITS so as not to limit distribution size
1427  int numQb = GENERATE( range(NUM_QUBITS, NUM_QUBITS+10) );
1428 
1429  SECTION( "state-vector" ) {
1430 
1431  Qureg vec = createQureg(numQb, QUEST_ENV);
1432  REQUIRE( getNumQubits(vec) == numQb );
1433  destroyQureg(vec, QUEST_ENV);
1434  }
1435  SECTION( "density-matrix" ) {
1436 
1437  Qureg mat = createDensityQureg(numQb, QUEST_ENV);
1438  REQUIRE( getNumQubits(mat) == numQb );
1439  destroyQureg(mat, QUEST_ENV);
1440  }
1441  }
1442  SECTION( "input validation" ) {
1443 
1444  // no validation
1445  SUCCEED();
1446  }
1447 }
1448 
1449 
1450 
1455 TEST_CASE( "getProbAmp", "[calculations]" ) {
1456 
1457  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
1458 
1459  SECTION( "correctness" ) {
1460 
1461  SECTION( "state-vector" ) {
1462 
1463  initDebugState(vec);
1464  QVector ref = toQVector(vec);
1465 
1466  int ind = GENERATE( range(0,1<<NUM_QUBITS) );
1467  qreal refCalc = pow(abs(ref[ind]), 2);
1468  REQUIRE( getProbAmp(vec,ind) == Approx(refCalc) );
1469  }
1470  }
1471  SECTION( "input validation" ) {
1472 
1473  SECTION( "state index" ) {
1474 
1475  int ind = GENERATE( -1, 1<<NUM_QUBITS );
1476  REQUIRE_THROWS_WITH( getProbAmp(vec,ind), Contains("Invalid amplitude index") );
1477  }
1478  SECTION( "density-matrix" ) {
1479 
1480  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1481  REQUIRE_THROWS_WITH( getProbAmp(mat,0), Contains("valid only for state-vectors") );
1482  destroyQureg(mat, QUEST_ENV);
1483  }
1484  }
1485  destroyQureg(vec, QUEST_ENV);
1486 }
1487 
1488 
1489 
1494 TEST_CASE( "getRealAmp", "[calculations]" ) {
1495 
1496  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
1497 
1498  SECTION( "correctness" ) {
1499 
1500  SECTION( "state-vector" ) {
1501 
1502  initDebugState(vec);
1503  QVector ref = toQVector(vec);
1504 
1505  int ind = GENERATE( range(0,1<<NUM_QUBITS) );
1506  REQUIRE( getRealAmp(vec,ind) == real(ref[ind]) );
1507  }
1508  }
1509  SECTION( "input validation" ) {
1510 
1511  SECTION( "state index" ) {
1512 
1513  int ind = GENERATE( -1, 1<<NUM_QUBITS );
1514  REQUIRE_THROWS_WITH( getRealAmp(vec,ind), Contains("Invalid amplitude index") );
1515  }
1516  SECTION( "density-matrix" ) {
1517 
1518  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
1519  REQUIRE_THROWS_WITH( getRealAmp(mat,0), Contains("valid only for state-vectors") );
1520  destroyQureg(mat, QUEST_ENV);
1521  }
1522  }
1523  destroyQureg(vec, QUEST_ENV);
1524 }
1525 
1526 
getProbAmp
qreal getProbAmp(Qureg qureg, long long int index)
Get the probability of a state-vector at an index in the full state vector.
Definition: QuEST.c:932
initBlankState
void initBlankState(Qureg qureg)
Initialises a qureg to have all-zero-amplitudes.
Definition: QuEST.c:119
pauliOpType
pauliOpType
Codes for specifying Pauli operators.
Definition: QuEST.h:96
getFullOperatorMatrix
QMatrix getFullOperatorMatrix(int *ctrls, int numCtrls, int *targs, int numTargs, QMatrix op, int numQubits)
Takes a 2^numTargs-by-2^numTargs matrix op and a returns a 2^numQubits-by-2^numQubits matrix where op...
Definition: utilities.cpp:304
QUEST_ENV
QuESTEnv QUEST_ENV
The global QuESTEnv instance, to be created and destroyed once in this main(), so that the MPI enviro...
Definition: main.cpp:20
fromComplex
#define fromComplex(comp)
PAULI_Z
@ PAULI_Z
Definition: QuEST.h:96
getKroneckerProduct
QVector getKroneckerProduct(QVector b, QVector a)
Returns b (otimes) a.
Definition: utilities.cpp:143
calcTotalProb
qreal calcTotalProb(Qureg qureg)
A debugging function which calculates the probability of the qubits in qureg being in any state,...
Definition: QuEST.c:1143
calcProbOfAllOutcomes
void calcProbOfAllOutcomes(qreal *retProbs, Qureg qureg, int *qubits, int numQubits)
Populates outcomeProbs with the probabilities of every outcome of the sub-register contained in qubit...
Definition: QuEST.c:1176
PAULI_I
@ PAULI_I
Definition: QuEST.h:96
getDensityAmp
Complex getDensityAmp(Qureg qureg, long long int row, long long int col)
Get an amplitude from a density matrix at a given row and column.
Definition: QuEST.c:949
getRandomInt
int getRandomInt(int min, int max)
Returns a random integer between min (inclusive) and max (exclusive), from the uniform distribution.
Definition: utilities.cpp:526
getImagAmp
qreal getImagAmp(Qureg qureg, long long int index)
Get the imaginary component of the complex probability amplitude at an index in the state vector.
Definition: QuEST.c:925
syncDiagonalOp
void syncDiagonalOp(DiagonalOp op)
Update the GPU memory with the current values in op.real and op.imag.
Definition: QuEST.c:1531
destroyDiagonalOp
void destroyDiagonalOp(DiagonalOp op, QuESTEnv env)
Destroys a DiagonalOp created with createDiagonalOp(), freeing its memory.
Definition: QuEST.c:1524
NUM_QUBITS
#define NUM_QUBITS
The default number of qubits in the registers created for unit testing (both statevectors and density...
Definition: utilities.hpp:36
calcPurity
qreal calcPurity(Qureg qureg)
Calculates the purity of a density matrix, by the trace of the density matrix squared.
Definition: QuEST.c:1185
calcProbOfOutcome
qreal calcProbOfOutcome(Qureg qureg, int measureQubit, int outcome)
Gives the probability of a specified qubit being measured in the given outcome (0 or 1).
Definition: QuEST.c:1166
calcFidelity
qreal calcFidelity(Qureg qureg, Qureg pureState)
Calculates the fidelity of qureg (a state-vector or density matrix) against a reference pure state (n...
Definition: QuEST.c:1191
getRandomReal
qreal getRandomReal(qreal min, qreal max)
Returns a random real between min (inclusive) and max (exclusive), from the uniform distribution.
Definition: utilities.cpp:421
getKetBra
QMatrix getKetBra(QVector ket, QVector bra)
Returns the matrix |ket><bra|, with ith-jth element ket(i) conj(bra(j)), since |ket><bra| = sum_i a_i...
Definition: utilities.cpp:169
calcExpecDiagonalOp
Complex calcExpecDiagonalOp(Qureg qureg, DiagonalOp op)
Computes the expected value of the diagonal operator op for state qureg.
Definition: QuEST.c:1228
qreal
#define qreal
toQMatrix
QMatrix toQMatrix(ComplexMatrix2 src)
Returns a copy of the given 2-by-2 matrix.
Definition: utilities.cpp:1044
toQureg
void toQureg(Qureg qureg, QVector vec)
Initialises the state-vector qureg to have the same amplitudes as vec.
Definition: utilities.cpp:1201
PAULI_X
@ PAULI_X
Definition: QuEST.h:96
getAmp
Complex getAmp(Qureg qureg, long long int index)
Get the complex amplitude at a given index in the state vector.
Definition: QuEST.c:939
calcExpecPauliHamil
qreal calcExpecPauliHamil(Qureg qureg, PauliHamil hamil, Qureg workspace)
Computes the expected value of qureg under Hermitian operator hamil.
Definition: QuEST.c:1219
setRandomPauliSum
void setRandomPauliSum(qreal *coeffs, pauliOpType *codes, int numQubits, int numTerms)
Populates the coeffs array with random qreals in (-5, 5), and populates codes with random Pauli codes...
Definition: utilities.cpp:1229
calcHilbertSchmidtDistance
qreal calcHilbertSchmidtDistance(Qureg a, Qureg b)
Computes the Hilbert Schmidt distance between two density matrices a and b, defined as the Frobenius ...
Definition: QuEST.c:1237
QVector
std::vector< qcomp > QVector
A complex vector, which can be zero-initialised with QVector(numAmps).
Definition: utilities.hpp:60
DiagonalOp::imag
qreal * imag
The imaginary values of the 2^numQubits complex elements.
Definition: QuEST.h:310
toQVector
QVector toQVector(Qureg qureg)
Returns an equal-size copy of the given state-vector qureg.
Definition: utilities.cpp:1113
qcomp
#define qcomp
TEST_CASE
TEST_CASE("calcDensityInnerProduct", "[calculations]")
Definition: test_calculations.cpp:21
PauliHamil::pauliCodes
enum pauliOpType * pauliCodes
The Pauli operators acting on each qubit, flattened over every operator.
Definition: QuEST.h:281
getRealAmp
qreal getRealAmp(Qureg qureg, long long int index)
Get the real component of the complex probability amplitude at an index in the state vector.
Definition: QuEST.c:918
DiagonalOp
Represents a diagonal complex operator on the full Hilbert state of a Qureg.
Definition: QuEST.h:297
PAULI_Y
@ PAULI_Y
Definition: QuEST.h:96
PauliHamil
A Pauli Hamiltonian, expressed as a real-weighted sum of pauli products, and which can hence represen...
Definition: QuEST.h:277
getNumQubits
int getNumQubits(Qureg qureg)
Returns the number of qubits represented by qureg.
Definition: QuEST.c:908
calcInnerProduct
Complex calcInnerProduct(Qureg bra, Qureg ket)
Computes the inner product of two equal-size state vectors, given by.
Definition: QuEST.c:1150
destroyQureg
void destroyQureg(Qureg qureg, QuESTEnv env)
Deallocate a Qureg, freeing its memory.
Definition: QuEST.c:77
getRandomStateVector
QVector getRandomStateVector(int numQb)
Returns a random numQb-length L2-normalised state-vector from an undisclosed distribution.
Definition: utilities.cpp:468
initDebugState
void initDebugState(Qureg qureg)
Initialises qureg to be in the un-normalised, non-physical state with with -th complex amplitude give...
Definition: QuEST.c:1578
getRandomQMatrix
QMatrix getRandomQMatrix(int dim)
Returns a dim-by-dim complex matrix, where the real and imaginary value of each element are independe...
Definition: utilities.cpp:379
Qureg
Represents a system of qubits.
Definition: QuEST.h:322
getNumAmps
long long int getNumAmps(Qureg qureg)
Returns the number of complex amplitudes in a state-vector qureg.
Definition: QuEST.c:912
DiagonalOp::numElemsPerChunk
long long int numElemsPerChunk
The number of the 2^numQubits amplitudes stored on each distributed node.
Definition: QuEST.h:302
QMatrix
std::vector< std::vector< qcomp > > QMatrix
A complex square matrix.
Definition: utilities.hpp:49
calcExpecPauliSum
qreal calcExpecPauliSum(Qureg qureg, enum pauliOpType *allPauliCodes, qreal *termCoeffs, int numSumTerms, Qureg workspace)
Computes the expected value of a sum of products of Pauli operators.
Definition: QuEST.c:1210
sublists
Catch::Generators::GeneratorWrapper< int * > sublists(int *list, int len, int sublen)
Returns a Catch2 generator of every length-sublen sublist of length-len list, in increasing lexograph...
Definition: utilities.cpp:1488
Qureg::numQubitsRepresented
int numQubitsRepresented
The number of qubits represented in either the state-vector or density matrix.
Definition: QuEST.h:327
DiagonalOp::real
qreal * real
The real values of the 2^numQubits complex elements.
Definition: QuEST.h:308
Complex::real
qreal real
Definition: QuEST.h:105
createQureg
Qureg createQureg(int numQubits, QuESTEnv env)
Creates a state-vector Qureg object representing a set of qubits which will remain in a pure state.
Definition: QuEST.c:36
destroyPauliHamil
void destroyPauliHamil(PauliHamil h)
Destroy a PauliHamil instance, created with either createPauliHamil() or createPauliHamilFromFile().
Definition: QuEST.c:1414
Complex::imag
qreal imag
Definition: QuEST.h:106
getRandomDensityMatrix
QMatrix getRandomDensityMatrix(int numQb)
Returns a random numQb-by-numQb density matrix, from an undisclosed distribution, in a very mixed sta...
Definition: utilities.cpp:490
Complex
Represents one complex number.
Definition: QuEST.h:103
getRandomQVector
QVector getRandomQVector(int dim)
Returns a dim-length vector with random complex amplitudes in the square joining {-1-i,...
Definition: utilities.cpp:435
initZeroState
void initZeroState(Qureg qureg)
Initialise qureg into the zero state.
Definition: QuEST.c:113
calcDensityInnerProduct
qreal calcDensityInnerProduct(Qureg rho1, Qureg rho2)
Computes the Hilbert-Schmidt scalar product (which is equivalent to the Frobenius inner product of ma...
Definition: QuEST.c:1158
areEqual
bool areEqual(QVector a, QVector b)
Returns true if the absolute value of the difference between every amplitude in vectors a and b is le...
Definition: utilities.cpp:398
createDensityQureg
Qureg createDensityQureg(int numQubits, QuESTEnv env)
Creates a density matrix Qureg object representing a set of qubits which can enter noisy and mixed st...
Definition: QuEST.c:50
applyReferenceOp
void applyReferenceOp(QVector &state, int *ctrls, int numCtrls, int *targs, int numTargs, QMatrix op)
Modifies the state-vector state to be the result of applying the multi-target operator matrix op,...
Definition: utilities.cpp:728
createPauliHamil
PauliHamil createPauliHamil(int numQubits, int numSumTerms)
Dynamically allocates a Hamiltonian expressed as a real-weighted sum of products of Pauli operators.
Definition: QuEST.c:1398
initPlusState
void initPlusState(Qureg qureg)
Initialise qureg into the plus state.
Definition: QuEST.c:125
calcExpecPauliProd
qreal calcExpecPauliProd(Qureg qureg, int *targetQubits, enum pauliOpType *pauliCodes, int numTargets, Qureg workspace)
Computes the expected value of a product of Pauli operators.
Definition: QuEST.c:1201
createDiagonalOp
DiagonalOp createDiagonalOp(int numQubits, QuESTEnv env)
Creates a DiagonalOp representing a diagonal operator on the full Hilbert space of a Qureg.
Definition: QuEST.c:1518