test_decoherence.cpp
Go to the documentation of this file.
1 
2 #include "catch.hpp"
3 #include "QuEST.h"
4 #include "utilities.hpp"
5 #include <random>
6 
9 #define PREPARE_TEST(qureg, ref) \
10  Qureg qureg = createDensityQureg(NUM_QUBITS, QUEST_ENV); \
11  initDebugState(qureg); \
12  QMatrix ref = toQMatrix(qureg);
13 
14 /* allows concise use of Contains in catch's REQUIRE_THROWS_WITH */
15 using Catch::Matchers::Contains;
16 
17 
18 
23 TEST_CASE( "mixDamping", "[decoherence]" ) {
24 
25  PREPARE_TEST(qureg, ref);
26 
27  SECTION( "correctness" ) {
28 
29  int target = GENERATE( range(0,NUM_QUBITS) );
30  qreal prob = getRandomReal(0, 1);
31  mixDamping(qureg, target, prob);
32 
33  // ref -> kraus0 ref kraus0^dagger + kraus1 ref kraus1^dagger
34  QMatrix kraus0{{1,0},{0,sqrt(1-prob)}};
35  QMatrix rho0 = ref;
36  applyReferenceOp(rho0, target, kraus0);
37  QMatrix kraus1{{0,sqrt(prob)},{0,0}};
38  QMatrix rho1 = ref;
39  applyReferenceOp(rho1, target, kraus1);
40  ref = rho0 + rho1;
41 
42  REQUIRE( areEqual(qureg, ref) );
43  }
44  SECTION( "validation ") {
45 
46  SECTION( "qubit index" ) {
47 
48  int target = GENERATE( -1, NUM_QUBITS );
49  REQUIRE_THROWS_WITH( mixDamping(qureg, target, 0), Contains("Invalid target") );
50 
51  }
52  SECTION( "probability" ) {
53 
54  REQUIRE_THROWS_WITH( mixDamping(qureg, 0, -.1), Contains("Probabilities") );
55  REQUIRE_THROWS_WITH( mixDamping(qureg, 0, 1.1), Contains("Probabilities") );
56  }
57  SECTION( "density-matrix" ) {
58 
59  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
60  REQUIRE_THROWS_WITH( mixDamping(vec, 0, 0), Contains("density matrices") );
61  destroyQureg(vec, QUEST_ENV);
62  }
63  }
64  destroyQureg(qureg, QUEST_ENV);
65 }
66 
67 
68 
73 TEST_CASE( "mixDensityMatrix", "[decoherence]" ) {
74 
75  Qureg qureg1 = createDensityQureg(NUM_QUBITS, QUEST_ENV);
76  Qureg qureg2 = createDensityQureg(NUM_QUBITS, QUEST_ENV);
77  initDebugState(qureg1);
78  initDebugState(qureg2);
79  QMatrix ref1 = toQMatrix(qureg1);
80  QMatrix ref2 = toQMatrix(qureg2);
81 
82  SECTION( "correctness" ) {
83 
84  // test p in {0, 1} and 10 random values in (0,1)
85  qreal prob = GENERATE( 0., 1., take(10, random(0.,1.)) );
86  mixDensityMatrix(qureg1, prob, qureg2);
87 
88  // ensure target qureg modified correctly
89  ref1 = (1-prob)*ref1 + (prob)*ref2;
90  REQUIRE( areEqual(qureg1, ref1) );
91 
92  // enure other qureg was not modified
93  REQUIRE( areEqual(qureg2, ref2) );
94  }
95  SECTION( "input validation" ) {
96 
97  SECTION( "probabilities") {
98 
99  qreal prob = GENERATE( -0.1, 1.1 );
100  REQUIRE_THROWS_WITH( mixDensityMatrix(qureg1, prob, qureg2), Contains("Probabilities") );
101  }
102  SECTION( "density matrices" ) {
103 
104  // one is statevec
105  Qureg state1 = createQureg(qureg1.numQubitsRepresented, QUEST_ENV);
106  REQUIRE_THROWS_WITH( mixDensityMatrix(qureg1, 0, state1), Contains("density matrices") );
107  REQUIRE_THROWS_WITH( mixDensityMatrix(state1, 0, qureg1), Contains("density matrices") );
108 
109  // both are statevec
110  Qureg state2 = createQureg(qureg1.numQubitsRepresented, QUEST_ENV);
111  REQUIRE_THROWS_WITH( mixDensityMatrix(state1, 0, state2), Contains("density matrices") );
112 
113  destroyQureg(state1, QUEST_ENV);
114  destroyQureg(state2, QUEST_ENV);
115  }
116  SECTION( "matching dimensions" ) {
117 
118  Qureg qureg3 = createDensityQureg(1 + qureg1.numQubitsRepresented, QUEST_ENV);
119  REQUIRE_THROWS_WITH( mixDensityMatrix(qureg1, 0, qureg3), Contains("Dimensions") );
120  REQUIRE_THROWS_WITH( mixDensityMatrix(qureg3, 0, qureg1), Contains("Dimensions") );
121  destroyQureg(qureg3, QUEST_ENV);
122  }
123  }
124  destroyQureg(qureg1, QUEST_ENV);
125  destroyQureg(qureg2, QUEST_ENV);
126 }
127 
128 
129 
134 TEST_CASE( "mixDephasing", "[decoherence]" ) {
135 
136  PREPARE_TEST(qureg, ref);
137 
138  SECTION( "correctness " ) {
139 
140  int target = GENERATE( range(0,NUM_QUBITS) );
141  qreal prob = getRandomReal(0, 1/2.);
142  mixDephasing(qureg, target, prob);
143 
144  // ref -> (1 - prob) ref + prob Z ref Z
145  QMatrix phaseRef = ref;
146  applyReferenceOp(phaseRef, target, QMatrix{{1,0},{0,-1}}); // Z ref Z
147  ref = ((1 - prob) * ref) + (prob * phaseRef);
148 
149  REQUIRE( areEqual(qureg, ref) );
150  }
151  SECTION( "validation ") {
152 
153  SECTION( "qubit index" ) {
154 
155  int target = GENERATE( -1, NUM_QUBITS );
156  REQUIRE_THROWS_WITH( mixDephasing(qureg, target, 0), Contains("Invalid target") );
157 
158  }
159  SECTION( "probability" ) {
160 
161  REQUIRE_THROWS_WITH( mixDephasing(qureg, 0, -.1), Contains("Probabilities") );
162  REQUIRE_THROWS_WITH( mixDephasing(qureg, 0, .6), Contains("probability") && Contains("cannot exceed 1/2") );
163  }
164  SECTION( "density-matrix" ) {
165 
166  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
167  REQUIRE_THROWS_WITH( mixDephasing(vec, 0, 0), Contains("density matrices") );
168  destroyQureg(vec, QUEST_ENV);
169  }
170  }
171  destroyQureg(qureg, QUEST_ENV);
172 }
173 
174 
175 
180 TEST_CASE( "mixDepolarising", "[decoherence]" ) {
181 
182  PREPARE_TEST(qureg, ref);
183 
184  SECTION( "correctness " ) {
185 
186  int target = GENERATE( range(0,NUM_QUBITS) );
187  qreal prob = getRandomReal(0, 3/4.);
188  mixDepolarising(qureg, target, prob);
189 
190  QMatrix xRef = ref;
191  applyReferenceOp(xRef, target, QMatrix{{0,1},{1,0}}); // X ref X
192  QMatrix yRef = ref;
193  applyReferenceOp(yRef, target, QMatrix{{0,-qcomp(0,1)},{qcomp(0,1),0}}); // Y ref Y
194  QMatrix zRef = ref;
195  applyReferenceOp(zRef, target, QMatrix{{1,0},{0,-1}}); // Z ref Z
196  ref = ((1 - prob) * ref) + ((prob/3.) * ( xRef + yRef + zRef));
197 
198  REQUIRE( areEqual(qureg, ref) );
199  }
200  SECTION( "validation ") {
201 
202  SECTION( "qubit index" ) {
203 
204  int target = GENERATE( -1, NUM_QUBITS );
205  REQUIRE_THROWS_WITH( mixDepolarising(qureg, target, 0), Contains("Invalid target") );
206 
207  }
208  SECTION( "probability" ) {
209 
210  REQUIRE_THROWS_WITH( mixDepolarising(qureg, 0, -.1), Contains("Probabilities") );
211  REQUIRE_THROWS_WITH( mixDepolarising(qureg, 0, .76), Contains("probability") && Contains("cannot exceed 3/4") );
212  }
213  SECTION( "density-matrix" ) {
214 
215  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
216  REQUIRE_THROWS_WITH( mixDepolarising(vec, 0, 0), Contains("density matrices") );
217  destroyQureg(vec, QUEST_ENV);
218  }
219  }
220  destroyQureg(qureg, QUEST_ENV);
221 }
222 
223 
224 
229 TEST_CASE( "mixMultiQubitKrausMap", "[decoherence]" ) {
230 
231  PREPARE_TEST(qureg, ref);
232 
233  // figure out max-num (inclusive) targs allowed by hardware backend
234  // (each node must contain as 2^(2*numTargs) amps)
235  int maxNumTargs = calcLog2(qureg.numAmpsPerChunk) / 2;
236 
237  SECTION( "correctness" ) {
238 
239  /* note that this function incurs a stack overhead when numTargs < 4,
240  * and a heap overhead when numTargs >= 4
241  */
242 
243  int numTargs = GENERATE_COPY( range(1,maxNumTargs+1) ); // inclusive upper bound
244 
245  // note this is very expensive to try every arrangement (2 min runtime for numTargs=5 alone)
246  int* targs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numTargs) );
247 
248  // try the min and max number of operators, and 2 random numbers
249  // (there are way too many to try all!)
250  int maxNumOps = (2*numTargs)*(2*numTargs);
251  int numOps = GENERATE_COPY( 1, maxNumOps, take(2,random(1,maxNumOps)) );
252 
253  // use a new random map
254  std::vector<QMatrix> matrs = getRandomKrausMap(numTargs, numOps);
255 
256  // create map in QuEST datatypes
257  ComplexMatrixN ops[numOps];
258  for (int i=0; i<numOps; i++) {
259  ops[i] = createComplexMatrixN(numTargs);
260  toComplexMatrixN(matrs[i], ops[i]);
261  }
262 
263  mixMultiQubitKrausMap(qureg, targs, numTargs, ops, numOps);
264 
265  // set ref -> K_i ref K_i^dagger
266  QMatrix matrRefs[numOps];
267  for (int i=0; i<numOps; i++) {
268  matrRefs[i] = ref;
269  applyReferenceOp(matrRefs[i], targs, numTargs, matrs[i]);
270  }
271  ref = getZeroMatrix(ref.size());
272  for (int i=0; i<numOps; i++)
273  ref += matrRefs[i];
274 
275  REQUIRE( areEqual(qureg, ref, 1E2*REAL_EPS) );
276 
277  // cleanup QuEST datatypes
278  for (int i=0; i<numOps; i++)
279  destroyComplexMatrixN(ops[i]);
280  }
281  SECTION( "input validation" ) {
282 
283  SECTION( "repetition of target" ) {
284 
285  // make valid targets
286  int targs[NUM_QUBITS];
287  for (int i=0; i<NUM_QUBITS; i++)
288  targs[i] = i;
289 
290  // duplicate one
291  int badInd = GENERATE( range(0,NUM_QUBITS) );
292  int copyInd = GENERATE_COPY( filter([=](int i){ return i!=badInd; }, range(0,NUM_QUBITS)) );
293  targs[badInd] = targs[copyInd];
294 
295  REQUIRE_THROWS_WITH( mixMultiQubitKrausMap(qureg, targs, NUM_QUBITS, NULL, 1), Contains("target qubits") && Contains("unique") );
296  }
297  SECTION( "qubit indices" ) {
298 
299  // make valid targets
300  int targs[NUM_QUBITS];
301  for (int i=0; i<NUM_QUBITS; i++)
302  targs[i] = i;
303 
304  // make one invalid
305  targs[GENERATE( range(0,NUM_QUBITS) )] = GENERATE( -1, NUM_QUBITS );
306 
307  REQUIRE_THROWS_WITH( mixMultiQubitKrausMap(qureg, targs, NUM_QUBITS, NULL, 1), Contains("Invalid target qubit") );
308  }
309  SECTION( "number of operators" ) {
310 
311  int numTargs = GENERATE_COPY( range(1,maxNumTargs+1) );
312  int maxNumOps = (2*numTargs)*(2*numTargs);
313  int numOps = GENERATE_REF( -1, 0, maxNumOps + 1 );
314 
315  // make valid targets to avoid triggering target validation
316  int targs[numTargs];
317  for (int i=0; i<numTargs; i++)
318  targs[i] = i;
319  REQUIRE_THROWS_WITH( mixMultiQubitKrausMap(qureg, targs, numTargs, NULL, numOps), Contains("operators may be specified") );
320  }
321  SECTION( "initialisation of operators" ) {
322 
323  /* compilers don't auto-initialise to NULL; the below circumstance
324  * only really occurs when 'malloc' returns NULL in createComplexMatrixN,
325  * which actually triggers its own validation. Hence this test is useless
326  * currently.
327  */
328 
329  int numTargs = NUM_QUBITS;
330  int numOps = (2*numTargs)*(2*numTargs);
331 
332  // no need to initialise ops, but set their attribs correct to avoid triggering other validation
333  ComplexMatrixN ops[numOps];
334  for (int i=0; i<numOps; i++)
335  ops[i].numQubits = numTargs;
336 
337  // make one of the max-ops explicitly NULL
338  ops[GENERATE_COPY( range(0,numTargs) )].real = NULL;
339 
340  // make valid targets to avoid triggering target validation
341  int targs[numTargs];
342  for (int i=0; i<numTargs; i++)
343  targs[i] = i;
344 
345  REQUIRE_THROWS_WITH( mixMultiQubitKrausMap(qureg, targs, numTargs, ops, numOps), Contains("ComplexMatrixN") && Contains("created") );
346  }
347  SECTION( "dimension of operators" ) {
348 
349  // make valid (dimension-wise) max-qubits Kraus map
350  int numTargs = NUM_QUBITS;
351  int numOps = (2*numTargs)*(2*numTargs);
352  ComplexMatrixN ops[numOps];
353  for (int i=0; i<numOps; i++)
354  ops[i] = createComplexMatrixN(numTargs);
355 
356  // make one have wrong-dimensions
357  int badInd = GENERATE_COPY( range(0,numTargs) );
358  destroyComplexMatrixN(ops[badInd]);
359  ops[badInd] = createComplexMatrixN(numTargs - 1);
360 
361  // make valid targets to avoid triggering target validation
362  int targs[numTargs];
363  for (int i=0; i<numTargs; i++)
364  targs[i] = i;
365 
366  REQUIRE_THROWS_WITH( mixMultiQubitKrausMap(qureg, targs, numTargs, ops, numOps), Contains("same number of qubits") );
367 
368  for (int i=0; i<numOps; i++)
369  destroyComplexMatrixN(ops[i]);
370  }
371  SECTION( "trace preserving" ) {
372 
373  int numTargs = GENERATE_COPY( range(1,maxNumTargs+1) );
374  int maxNumOps = (2*numTargs) * (2*numTargs);
375  int numOps = GENERATE_COPY( 1, 2, maxNumOps );
376 
377  // generate a valid map
378  std::vector<QMatrix> matrs = getRandomKrausMap(numTargs, numOps);
379  ComplexMatrixN ops[numOps];
380  for (int i=0; i<numOps; i++) {
381  ops[i] = createComplexMatrixN(numTargs);
382  toComplexMatrixN(matrs[i], ops[i]);
383  }
384 
385  // make only one invalid
386  ops[GENERATE_REF( range(0,numOps) )].real[0][0] = 0;
387 
388  // make valid targets to avoid triggering target validation
389  int targs[numTargs];
390  for (int i=0; i<numTargs; i++)
391  targs[i] = i;
392 
393  REQUIRE_THROWS_WITH( mixMultiQubitKrausMap(qureg, targs, numTargs, ops, numOps), Contains("trace preserving") );
394 
395  for (int i=0; i<numOps; i++)
396  destroyComplexMatrixN(ops[i]);
397  }
398  SECTION( "density-matrix" ) {
399 
400  Qureg statevec = createQureg(NUM_QUBITS, QUEST_ENV);
401 
402  // make valid targets to avoid triggering target validation
403  int targs[NUM_QUBITS];
404  for (int i=0; i<NUM_QUBITS; i++)
405  targs[i] = i;
406 
407  REQUIRE_THROWS_WITH( mixMultiQubitKrausMap(statevec, targs, NUM_QUBITS, NULL, 1), Contains("valid only for density matrices") );
408  destroyQureg(statevec, QUEST_ENV);
409 
410  }
411  SECTION( "operator fits in node" ) {
412 
413  // each node requires (2 numTargs)^2 amplitudes
414  int minAmps = (2*NUM_QUBITS) * (2*NUM_QUBITS);
415 
416  // make valid targets to avoid triggering target validation
417  int targs[NUM_QUBITS];
418  for (int i=0; i<NUM_QUBITS; i++)
419  targs[i] = i;
420 
421  // make a simple Identity map
422  ComplexMatrixN ops[] = {createComplexMatrixN(NUM_QUBITS)};
423  for (int i=0; i<(1<<NUM_QUBITS); i++)
424  ops[0].real[i][i] = 1;
425 
426  // fake a smaller qureg
427  qureg.numAmpsPerChunk = minAmps - 1;
428  REQUIRE_THROWS_WITH( mixMultiQubitKrausMap(qureg, targs, NUM_QUBITS, ops, 1), Contains("targets too many qubits") && Contains("cannot all fit") );
429 
430  destroyComplexMatrixN(ops[0]);
431  }
432  }
433  destroyQureg(qureg, QUEST_ENV);
434 }
435 
436 
437 
442 TEST_CASE( "mixPauli", "[decoherence]" ) {
443 
444  PREPARE_TEST(qureg, ref);
445 
446  SECTION( "correctness" ) {
447 
448  int target = GENERATE( range(0,NUM_QUBITS) );
449 
450  // randomly generate valid pauli-error probabilities
451  qreal probs[3];
452  qreal max0 = 1/2.; // satisfies p1 < 1 - py
453  probs[0] = getRandomReal(0, max0);
454  qreal max1 = (max0 - probs[0])/2.; // p2 can use half of p1's "unused space"
455  probs[1] = getRandomReal(0, max1);
456  qreal max2 = (max1 - probs[1])/2.; // p3 can use half of p2's "unused space"
457  probs[2] = getRandomReal(0, max2);
458 
459  // uniformly randomly assign probs (bound to target)
460  int inds[3] = {0,1,2};
461  std::shuffle(inds,inds+3, std::default_random_engine(1E5 * target));
462  qreal probX = probs[inds[0]]; // seed:target shows no variation
463  qreal probY = probs[inds[1]];
464  qreal probZ = probs[inds[2]];
465 
466  mixPauli(qureg, target, probX, probY, probZ);
467 
468  QMatrix xRef = ref;
469  applyReferenceOp(xRef, target, QMatrix{{0,1},{1,0}}); // X ref X
470  QMatrix yRef = ref;
471  applyReferenceOp(yRef, target, QMatrix{{0,-qcomp(0,1)},{qcomp(0,1),0}}); // Y ref Y
472  QMatrix zRef = ref;
473  applyReferenceOp(zRef, target, QMatrix{{1,0},{0,-1}}); // Z ref Z
474  ref = ((1 - probX - probY - probZ) * ref) +
475  (probX * xRef) + (probY * yRef) + (probZ * zRef);
476 
477  REQUIRE( areEqual(qureg, ref) );
478  }
479  SECTION( "input validation" ) {
480 
481  SECTION( "qubit index" ) {
482 
483  int target = GENERATE( -1, NUM_QUBITS );
484  REQUIRE_THROWS_WITH( mixPauli(qureg, target, 0, 0, 0), Contains("Invalid target") );
485 
486  }
487  SECTION( "probability" ) {
488 
489  int target = 0;
490 
491  // probs clearly must be in [0, 1]
492  REQUIRE_THROWS_WITH( mixPauli(qureg, target, -.1, 0, 0), Contains("Probabilities") );
493  REQUIRE_THROWS_WITH( mixPauli(qureg, target, 0, -.1, 0), Contains("Probabilities") );
494  REQUIRE_THROWS_WITH( mixPauli(qureg, target, 0, 0, -.1), Contains("Probabilities") );
495 
496  // max single-non-zero-prob is 0.5
497  REQUIRE_THROWS_WITH( mixPauli(qureg, target, .6, 0, 0), Contains("cannot exceed the probability") );
498  REQUIRE_THROWS_WITH( mixPauli(qureg, target, 0, .6, 0), Contains("cannot exceed the probability") );
499  REQUIRE_THROWS_WITH( mixPauli(qureg, target, 0, 0, .6), Contains("cannot exceed the probability") );
500 
501  // must satisfy px, py, pz < 1 - px - py - pz
502  REQUIRE_THROWS_WITH( mixPauli(qureg, target, .3, .3, .3), Contains("cannot exceed the probability") );
503  }
504  SECTION( "density-matrix" ) {
505 
506  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
507  REQUIRE_THROWS_WITH( mixPauli(vec, 0, 0, 0, 0), Contains("density matrices") );
508  destroyQureg(vec, QUEST_ENV);
509  }
510  }
511  destroyQureg(qureg, QUEST_ENV);
512 }
513 
514 
515 
520 TEST_CASE( "mixKrausMap", "[decoherence]" ) {
521 
522  PREPARE_TEST(qureg, ref);
523 
524  SECTION( "correctness" ) {
525 
526  int target = GENERATE( range(0,NUM_QUBITS) );
527  int numOps = GENERATE( range(1,5) ); // max 4 inclusive
528  std::vector<QMatrix> matrs = getRandomKrausMap(1, numOps);
529 
530  ComplexMatrix2 ops[numOps];
531  for (int i=0; i<numOps; i++)
532  ops[i] = toComplexMatrix2(matrs[i]);
533  mixKrausMap(qureg, target, ops, numOps);
534 
535  // set ref -> K_i ref K_i^dagger
536  QMatrix matrRefs[numOps];
537  for (int i=0; i<numOps; i++) {
538  matrRefs[i] = ref;
539  applyReferenceOp(matrRefs[i], target, matrs[i]);
540  }
541  ref = getZeroMatrix(ref.size());
542  for (int i=0; i<numOps; i++)
543  ref += matrRefs[i];
544 
545  REQUIRE( areEqual(qureg, ref, 10*REAL_EPS) );
546  }
547  SECTION( "input validation" ) {
548 
549  SECTION( "number of operators" ) {
550 
551  int numOps = GENERATE( 0, 5 );
552  REQUIRE_THROWS_WITH( mixKrausMap(qureg, 0, NULL, numOps), Contains("operators") );
553  }
554  SECTION( "trace preserving" ) {
555 
556  // valid Kraus map
557  int numOps = GENERATE( range(1,5) ); // max 4 inclusive
558  std::vector<QMatrix> matrs = getRandomKrausMap(1, numOps);
559  ComplexMatrix2 ops[numOps];
560  for (int i=0; i<numOps; i++)
561  ops[i] = toComplexMatrix2(matrs[i]);
562 
563  // make invalid
564  ops[GENERATE_REF( range(0,numOps) )].real[0][0] = 0;
565  REQUIRE_THROWS_WITH( mixKrausMap(qureg, 0, ops, numOps), Contains("trace preserving") );
566 
567  }
568  SECTION( "qubit index" ) {
569 
570  int target = GENERATE( -1, NUM_QUBITS );
571  REQUIRE_THROWS_WITH( mixKrausMap(qureg, target, NULL, 1), Contains("Invalid target qubit") );
572  }
573  SECTION( "density-matrix" ) {
574 
575  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
576  REQUIRE_THROWS_WITH( mixKrausMap(vec, 0, NULL, 1), Contains("density matrices") );
577  destroyQureg(vec, QUEST_ENV);
578  }
579  SECTION( "operators fit in node" ) {
580 
581  qureg.numAmpsPerChunk = 3; // min 4
582  REQUIRE_THROWS_WITH( mixKrausMap(qureg, 0, NULL, 1), Contains("targets too many qubits") );
583  }
584  }
585  destroyQureg(qureg, QUEST_ENV);
586 }
587 
588 
589 
594 TEST_CASE( "mixTwoQubitDephasing", "[decoherence]" ) {
595 
596  PREPARE_TEST(qureg, ref);
597 
598  SECTION( "correctness" ) {
599 
600  int targ1 = GENERATE( range(0,NUM_QUBITS) );
601  int targ2 = GENERATE_COPY( filter([=](int t){ return t!=targ1; }, range(0,NUM_QUBITS)) );
602  qreal prob = getRandomReal(0, 3/4.);
603 
604  mixTwoQubitDephasing(qureg, targ1, targ2, prob);
605 
606  // ref -> (1 - prob) ref + prob/3 (Z1 ref Z1 + Z2 ref Z2 + Z1 Z2 ref Z1 Z2)
607  QMatrix zMatr{{1,0},{0,-1}};
608  QMatrix z1Ref = ref;
609  applyReferenceOp(z1Ref, targ1, zMatr); // Z1 ref Z1
610  QMatrix z2Ref = ref;
611  applyReferenceOp(z2Ref, targ2, zMatr); // Z2 ref Z2
612  QMatrix z1z2Ref = ref;
613  applyReferenceOp(z1z2Ref, targ1, zMatr);
614  applyReferenceOp(z1z2Ref, targ2, zMatr); // Z1 Z2 ref Z1 Z2
615  ref = ((1 - prob) * ref) + (prob/3.) * (z1Ref + z2Ref + z1z2Ref);
616 
617  REQUIRE( areEqual(qureg, ref) );
618  }
619  SECTION( "input validation" ) {
620 
621  SECTION( "qubit indices" ) {
622 
623  int targ = GENERATE( -1, NUM_QUBITS );
624  REQUIRE_THROWS_WITH( mixTwoQubitDephasing(qureg, 0, targ, 0), Contains("Invalid target") );
625  REQUIRE_THROWS_WITH( mixTwoQubitDephasing(qureg, targ, 0, 0), Contains("Invalid target") );
626  }
627  SECTION( "target collision" ) {
628 
629  int targ = GENERATE( range(0,NUM_QUBITS) );
630  REQUIRE_THROWS_WITH( mixTwoQubitDephasing(qureg, targ, targ, 0), Contains("target") && Contains("unique") );
631  }
632  SECTION( "probability" ) {
633 
634  REQUIRE_THROWS_WITH( mixTwoQubitDephasing(qureg, 0, 1, -.1), Contains("Probabilities") );
635  REQUIRE_THROWS_WITH( mixTwoQubitDephasing(qureg, 0, 1, 3/4. + .01), Contains("probability") && Contains("cannot exceed 3/4") );
636  }
637  SECTION( "density-matrix" ) {
638 
639  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
640  REQUIRE_THROWS_WITH( mixTwoQubitDephasing(vec, 0, 1, 0), Contains("density matrices") );
641  destroyQureg(vec, QUEST_ENV);
642  }
643  }
644  destroyQureg(qureg, QUEST_ENV);
645 }
646 
647 
648 
653 TEST_CASE( "mixTwoQubitDepolarising", "[decoherence]" ) {
654 
655  PREPARE_TEST(qureg, ref);
656 
657  SECTION( "correctness" ) {
658 
659  int targ1 = GENERATE( range(0,NUM_QUBITS) );
660  int targ2 = GENERATE_COPY( filter([=](int t){ return t!=targ1; }, range(0,NUM_QUBITS)) );
661  qreal prob = getRandomReal(0, 15/16.);
662 
663  mixTwoQubitDepolarising(qureg, targ1, targ2, prob);
664 
665  QMatrix paulis[4] = {
666  QMatrix{{1,0},{0,1}}, // I
667  QMatrix{{0,1},{1,0}}, // X
668  QMatrix{{0,-qcomp(0,1)},{qcomp(0,1),0}}, // Y
669  QMatrix{{1,0},{0,-1}} // Z
670  };
671 
672  int targs[2] = {targ1, targ2};
673  QMatrix refInit = ref;
674  ref = (1 - (16/15.)*prob) * ref;
675  for (int i=0; i<4; i++) {
676  for (int j=0; j<4; j++) {
677  QMatrix term = refInit;
678  applyReferenceOp(term, targs, 2,
679  getKroneckerProduct(paulis[i], paulis[j]));
680  ref += (prob/15.) * term;
681  }
682  }
683 
684  REQUIRE( areEqual(qureg, ref, 1E4*REAL_EPS) );
685  }
686  SECTION( "input validation" ) {
687 
688  SECTION( "qubit indices" ) {
689 
690  int targ = GENERATE( -1, NUM_QUBITS );
691  REQUIRE_THROWS_WITH( mixTwoQubitDepolarising(qureg, 0, targ, 0), Contains("Invalid target") );
692  REQUIRE_THROWS_WITH( mixTwoQubitDepolarising(qureg, targ, 0, 0), Contains("Invalid target") );
693  }
694  SECTION( "target collision" ) {
695 
696  int targ = GENERATE( range(0,NUM_QUBITS) );
697  REQUIRE_THROWS_WITH( mixTwoQubitDepolarising(qureg, targ, targ, 0), Contains("target") && Contains("unique") );
698  }
699  SECTION( "probability" ) {
700 
701  REQUIRE_THROWS_WITH( mixTwoQubitDepolarising(qureg, 0, 1, -.1), Contains("Probabilities") );
702  REQUIRE_THROWS_WITH( mixTwoQubitDepolarising(qureg, 0, 1, 15/16. + .01), Contains("probability") && Contains("cannot exceed 15/16") );
703  }
704  SECTION( "density-matrix" ) {
705 
706  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
707  REQUIRE_THROWS_WITH( mixTwoQubitDepolarising(vec, 0, 1, 0), Contains("density matrices") );
708  destroyQureg(vec, QUEST_ENV);
709  }
710  }
711  destroyQureg(qureg, QUEST_ENV);
712 }
713 
714 
715 
720 TEST_CASE( "mixTwoQubitKrausMap", "[decoherence]" ) {
721 
722  PREPARE_TEST(qureg, ref);
723 
724  SECTION( "correctness" ) {
725 
726  int targ1 = GENERATE( range(0,NUM_QUBITS) );
727  int targ2 = GENERATE_COPY( filter([=](int t){ return t!=targ1; }, range(0,NUM_QUBITS)) );
728  int numOps = GENERATE( range(1,17) ); // max 16 inclusive
729  std::vector<QMatrix> matrs = getRandomKrausMap(2, numOps);
730 
731  ComplexMatrix4 ops[numOps];
732  for (int i=0; i<numOps; i++)
733  ops[i] = toComplexMatrix4(matrs[i]);
734  mixTwoQubitKrausMap(qureg, targ1, targ2, ops, numOps);
735 
736  // set ref -> K_i ref K_i^dagger
737  int targs[2] = {targ1, targ2};
738  QMatrix matrRefs[numOps];
739  for (int i=0; i<numOps; i++) {
740  matrRefs[i] = ref;
741  applyReferenceOp(matrRefs[i], targs, 2, matrs[i]);
742  }
743  ref = getZeroMatrix(ref.size());
744  for (int i=0; i<numOps; i++)
745  ref += matrRefs[i];
746 
747  REQUIRE( areEqual(qureg, ref, 10*REAL_EPS) );
748  }
749  SECTION( "input validation" ) {
750 
751  SECTION( "number of operators" ) {
752 
753  int numOps = GENERATE( 0, 17 );
754  REQUIRE_THROWS_WITH( mixTwoQubitKrausMap(qureg, 0,1, NULL, numOps), Contains("operators") );
755  }
756  SECTION( "trace preserving" ) {
757 
758  // valid Kraus map
759  int numOps = GENERATE( range(1,16) );
760  std::vector<QMatrix> matrs = getRandomKrausMap(2, numOps);
761  ComplexMatrix4 ops[numOps];
762  for (int i=0; i<numOps; i++)
763  ops[i] = toComplexMatrix4(matrs[i]);
764 
765  // make only one of the ops at a time invalid
766  ops[GENERATE_REF( range(0,numOps) )].real[0][0] = 0;
767  REQUIRE_THROWS_WITH( mixTwoQubitKrausMap(qureg, 0,1, ops, numOps), Contains("trace preserving") );
768  }
769  SECTION( "target collision" ) {
770 
771  int target = GENERATE( range(0,NUM_QUBITS) );
772  REQUIRE_THROWS_WITH( mixTwoQubitKrausMap(qureg, target, target, NULL, 1), Contains("target qubits") && Contains("unique") );
773  }
774  SECTION( "qubit index" ) {
775 
776  int target = GENERATE( -1, NUM_QUBITS );
777  REQUIRE_THROWS_WITH( mixTwoQubitKrausMap(qureg, 0,target, NULL, 1), Contains("Invalid target qubit") );
778  REQUIRE_THROWS_WITH( mixTwoQubitKrausMap(qureg, target,0, NULL, 1), Contains("Invalid target qubit") );
779  }
780  SECTION( "density-matrix" ) {
781 
782  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
783  REQUIRE_THROWS_WITH( mixTwoQubitKrausMap(vec, 0,1, NULL, 1), Contains("density matrices") );
784  destroyQureg(vec, QUEST_ENV);
785  }
786  SECTION( "operators fit in node" ) {
787 
788  qureg.numAmpsPerChunk = 15; // min 16
789  REQUIRE_THROWS_WITH( mixTwoQubitKrausMap(qureg, 0,1, NULL, 1), Contains("targets too many qubits") );
790  }
791  }
792  destroyQureg(qureg, QUEST_ENV);
793 }
794 
mixDephasing
void mixDephasing(Qureg qureg, int targetQubit, qreal prob)
Mixes a density matrix qureg to induce single-qubit dephasing noise.
Definition: QuEST.c:1250
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
ComplexMatrix4::real
qreal real[4][4]
Definition: QuEST.h:177
getKroneckerProduct
QVector getKroneckerProduct(QVector b, QVector a)
Returns b (otimes) a.
Definition: utilities.cpp:143
mixDepolarising
void mixDepolarising(Qureg qureg, int targetQubit, qreal prob)
Mixes a density matrix qureg to induce single-qubit homogeneous depolarising noise.
Definition: QuEST.c:1272
mixDamping
void mixDamping(Qureg qureg, int targetQubit, qreal prob)
Mixes a density matrix qureg to induce single-qubit amplitude damping (decay to 0 state).
Definition: QuEST.c:1283
destroyComplexMatrixN
void destroyComplexMatrixN(ComplexMatrixN m)
Destroy a ComplexMatrixN instance created with createComplexMatrixN()
Definition: QuEST.c:1369
getRandomKrausMap
std::vector< QMatrix > getRandomKrausMap(int numQb, int numOps)
Returns a random Kraus map of #numOps 2^numQb-by-2^numQb operators, from an undisclosed distribution.
Definition: utilities.cpp:578
createComplexMatrixN
ComplexMatrixN createComplexMatrixN(int numQubits)
Allocate dynamic memory for a square complex matrix of any size, which can be passed to functions lik...
Definition: QuEST.c:1348
mixKrausMap
void mixKrausMap(Qureg qureg, int target, ComplexMatrix2 *ops, int numOps)
Apply a general single-qubit Kraus map to a density matrix, as specified by at most four Kraus operat...
Definition: QuEST.c:1314
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
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
mixTwoQubitDepolarising
void mixTwoQubitDepolarising(Qureg qureg, int qubit1, int qubit2, qreal prob)
Mixes a density matrix qureg to induce two-qubit homogeneous depolarising noise.
Definition: QuEST.c:1291
ComplexMatrix4
Represents a 4x4 matrix of complex numbers.
Definition: QuEST.h:175
ComplexMatrixN
Represents a general 2^N by 2^N matrix of complex numbers.
Definition: QuEST.h:186
qreal
#define qreal
toQMatrix
QMatrix toQMatrix(ComplexMatrix2 src)
Returns a copy of the given 2-by-2 matrix.
Definition: utilities.cpp:1044
calcLog2
unsigned int calcLog2(long unsigned int num)
returns log2 of numbers which must be gauranteed to be 2^n
Definition: QuEST_validation.c:328
toComplexMatrixN
void toComplexMatrixN(QMatrix qm, ComplexMatrixN cm)
Initialises cm with the values of qm.
Definition: utilities.cpp:1033
qcomp
#define qcomp
toComplexMatrix4
ComplexMatrix4 toComplexMatrix4(QMatrix qm)
Returns a ComplexMatrix4 copy of QMatix qm.
Definition: utilities.cpp:1027
destroyQureg
void destroyQureg(Qureg qureg, QuESTEnv env)
Deallocate a Qureg, freeing its memory.
Definition: QuEST.c:77
mixTwoQubitKrausMap
void mixTwoQubitKrausMap(Qureg qureg, int target1, int target2, ComplexMatrix4 *ops, int numOps)
Apply a general two-qubit Kraus map to a density matrix, as specified by at most sixteen Kraus operat...
Definition: QuEST.c:1324
ComplexMatrixN::real
qreal ** real
Definition: QuEST.h:189
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
Qureg
Represents a system of qubits.
Definition: QuEST.h:322
mixTwoQubitDephasing
void mixTwoQubitDephasing(Qureg qureg, int qubit1, int qubit2, qreal prob)
Mixes a density matrix qureg to induce two-qubit dephasing noise.
Definition: QuEST.c:1260
ComplexMatrix2::real
qreal real[2][2]
Definition: QuEST.h:139
PREPARE_TEST
#define PREPARE_TEST(qureg, ref)
Prepares a density matrix in the debug state, and the reference QMatrix.
Definition: test_decoherence.cpp:9
QMatrix
std::vector< std::vector< qcomp > > QMatrix
A complex square matrix.
Definition: utilities.hpp:49
toComplexMatrix2
ComplexMatrix2 toComplexMatrix2(QMatrix qm)
Returns a ComplexMatrix2 copy of QMatix qm.
Definition: utilities.cpp:1021
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
mixDensityMatrix
void mixDensityMatrix(Qureg combineQureg, qreal otherProb, Qureg otherQureg)
Modifies combineQureg to become (1-prob)combineProb + prob otherQureg.
Definition: QuEST.c:1012
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
mixPauli
void mixPauli(Qureg qureg, int qubit, qreal probX, qreal probY, qreal probZ)
Mixes a density matrix qureg to induce general single-qubit Pauli noise.
Definition: QuEST.c:1303
getZeroMatrix
QMatrix getZeroMatrix(size_t dim)
Returns a dim-by-dim square complex matrix, initialised to all zeroes.
Definition: utilities.cpp:153
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
TEST_CASE
TEST_CASE("mixDamping", "[decoherence]")
Definition: test_decoherence.cpp:23
mixMultiQubitKrausMap
void mixMultiQubitKrausMap(Qureg qureg, int *targets, int numTargets, ComplexMatrixN *ops, int numOps)
Apply a general N-qubit Kraus map to a density matrix, as specified by at most (2N)^2 Kraus operators...
Definition: QuEST.c:1334
ComplexMatrix2
Represents a 2x2 matrix of complex numbers.
Definition: QuEST.h:137