-
Notifications
You must be signed in to change notification settings - Fork 22
/
UniversalQCompiler.m
4898 lines (4472 loc) · 248 KB
/
UniversalQCompiler.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
(* ::Package:: *)
(*
Copyright 2019 UniversalQCompiler (https://github.com/Q-Compiler/UniversalQCompiler)
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*)
BeginPackage["UniversalQCompiler`",{"QI`"}];
(*ToDos:
1. Add "do not simplify" option to methods, which allows to work with fixed circuit topologies.
2. Clear up unncessesary applications of Reverse[] (restructure the code to make it more redable)
3. Improve efficiency by removing For-loops and if statements with more efficient alternatives provided by Mathematica (and do not copy (long)lists)
4. Use ZXZ and XZX decompositions for single-qubit unitaries to simplify gate sequences
(for instance if one of the Zs in ZXZ is actually identity, the X might commute out of a neighbouring target gate)
5. For a gate sequence st, calling SimplifyGateList[XRGatesToCNOTRotations[CNOTRotationsToXXRGates[st]]] should give us back the same gate sequence st again.
6. Implementation of multi-controlled-Toffoli gates using ancillas to lower the C-NOT count.
7. [Raban] Extend Stateprepration scheme for small Schmidt rank
8. [Raban] Remove "\textnormal" for Latex outputs
9. In SimplifyGateList, reduce combinations of three rotation gates at the start of a circuit on ancilla qubits to two.
*)
(*Methods to handle and simplify gate sequence*)
(*CreateIsometryFromList::usage="CreateIsometryFromList[st,(n)] creates the operator corresponding to the list (optionally, the total number of qubits n can be determined)."
NCreateIsometryFromList::usage="CreateIsometryFromList[st,(n)] creates the operator (with numerical numbers) corresponding to the list (optionally, the total number of qubits n can be determined)."
CreateChannelFromList::usage="CreateChannelFromList[st, (n)] multiplies the gates in the list st on n qubits and outputs the channel represented."
NCreateChannelFromList::usage="NCreateChannelFromList[st, (n)] multiplies the gates in the list st on n qubits numericallyand outputs the channel represented."
CreateInstrumentFromList::usage="CreateInstrumentFromList[st, (n)] multiplies the gates in the list st on n qubits and outputs the instrument represented."
NCreateInstrumentFromList::usage="NCreateInstrumentFromList[st, (n)] multiplies the gates in the list st on n qubits numerically and outputs the instrument represented."
*)
CreatePOVMFromGateList::usage="CreatePOVMFromGateList[st, (n)] multiplies the gates in the list st on n qubits and outputs the POVM represented (optionally, the total number of qubits n can be determined)."
NCreatePOVMFromGateList::usage="NCreatePOVMFromGateList[st, (n)] multiplies the gates in the list st on n qubits numerically and outputs the POVM represented (optionally, the total number of qubits n can be determined)."
CreateOperationFromGateList::usage="CreateOperationFromGateList[st, (n)] generates the isometry/channel/instrument represented by st (optionally, the total number of qubits n can be determined)."
NCreateOperationFromGateList::usage="CreateOperationFromGateList[st, (n)] generates the isometry/channel/instrument represented by st numerically (optionally, the total number of qubits n can be determined)."
NGateList::usage="NGateList[st] makes the gate parameters numerical. [Helpful to speed up further processing of the gate list.]"
RelabelQubits::usage="RelabelQubits[st,numIn,numOut] relabels qubits in the gate sequence st (in list form) with qubit number in the list numIn with the qubit numbers given in the list numOut "
CNOTCount::usage="CNOTCount[st] returns the number of C-NOT gates in the gate sequence st."
InverseGateList::usage="InverseGateList[st] takes the inverse of a gate sequence st in list format."
AdjustAngle::usage="AdjustAngle[st] adjusts the angles in a gate sequence st to lie in [0,2 \[Pi]). Note that this may change the global phase by a factor -1."
ListFormToStr::usage="ListFormToStr[st] transforms a list in list format to a list containing the corresponding string representations of the gates."
SimplifyGateList::usage="SimplifyGateList[st] simplifies a gate sequence on in list format st by merging and commuting single-qubit and C-NOT gates (see the full documentation for more details)."
NSimplifyGateList::usage"NSimplifyGateList[st] simplifies a gate sequence (numerically) in list format st by merging and commuting single-qubit and C-NOT gates (see the full documentation for more details)."
NumberOfQubits::usage"NumberOfQubits[st] finds the largest qubit number featured in the list format st ."
(*Create gates in list form*)
CNOT::usage="CNOT[i,j] creates a C-NOT gate in list form with control qubit i and target qubit j."
CZ::usage="CZ[i,j] creates a controlled-Z gate in list form with control qubit i and target qubit j."
XX::usage="XX[phi,i,j] creates an XX gate in list form with parameter phi acting on qubits i,j."
RGate::usage="RGate[theta,phi,i] creates an R gate in list form with parameters theta,phi acting on qubit i."
Diag::usage="Diag[entr,act] creates a is a diagonal gate with diagonal entries given as a list entr and acting on the qubits listed in act."
Rx::usage="Rx[angle,act] creates an Rx gate in list form on the qubit act and with angle angle."
Ry::usage="Ry[diag,act] creates an Ry gate in list form on the qubit act and with angle angle."
Rz::usage="Rz[diag,act] creates an Rz gate in list form on the qubit act and with angle angle."
Mmt::usage="Mmt[act] creates a measurement in list form on the qubit act."
TrOut::usage="TrOut[act] creates a tracing out operation in list form on the qubit act."
Ancilla::usage="Ancilla[i_,act_] indicates that the qubit with number act is an ancilla starting in state |i> with i=0,1."
PostSelect::usage="PostSelect[i_,act_] indicates that the qubit with number act is postselected to state |i> with i=0,1."
(*Transform gates to matrices*)
MatrixFormOp::usage="MatrixFormOp[op] prints all the matrices in the list op in matrix form."
GateTypes::usage="GateTypes[] prints the convention for gate types used in this package."
ListFormToOp::usage="ListFormToOp[st,(n)] transforms a sequence of gates given in list from to a list containing the matrix repressentations of the gates (optionally, the total number of qubits n can be determined)."
RxM::usage="RxM[\[Alpha],i,n] creates a rotation gate around the x axis with angle \[Alpha] on qubit with number i, where the total number of qubits is n."
RyM::usage="RyM[\[Alpha],i,n] creates a rotation gate around the y axis with angle \[Alpha] on qubit with number i, where the total number of qubits is n."
RzM::usage="RzM[\[Alpha],i,n] creates a rotation gate around the z axis with angle \[Alpha] on qubit with number i, where the total number of qubits is n."
CNOTM::usage="CNOTM[i,j,n] creates a C-NOT gate on n qubits with control on the ith qubit and target on the jth qubit."
CZM::usage="CZM[i,j,n] creates a controlled-Z gate on n qubits with control on the ith qubit and target on the jth qubit."
XXM::usage="XXM[phi,{i,j},n] creates an XX gate with parameter phi acting on the i,j qubits with n qubits in total."
RGateM::usage="RGateM[theta,phi,i,n] creates an R-gate with parameters theta,phi acting on the ith qubit with n qubits in total."
DiagMat::usage="DiagMat[diag,act,n] creates a diagonal matrix with diagonal entries in the list diag representing a diagonal gate on n qubits acting on the qubits in the list act."
RzAngle::usage="RzAngle[Rz(\[Theta])] extracts the rotation angle \[Theta] from a rotation gate Rz (\[Theta])."
RyAngle::usage="RyAngle[Ry(\[Theta])] extracts the rotation angle \[Theta] from a rotation gate Ry (\[Theta])."
RxAngle::usage="RxAngle[Rx(\[Theta])] extracts the rotation angle \[Theta] from a rotation gate Rx (\[Theta])."
(*Apply gates efficiently*)
ApplyDiag::usage="TBA."
ApplyMCG::usage="TBA."
ApplyUCG::usage="TBA."
(*Visualization of ciruits*)
PrintCircuit::usage="PrintCircuit[st, (n)] prints a circuit st given in list form. Optionally, one can determine the total number of qubits n."
LatexQCircuit::usage="LatexQCircuit[st, (n)] provides an circuit in the form of the Latex package QCircuit. Optionally, one can determine the total number of qubits n."
(*Matrix decompositions*)
IsoToUnitary::usage="IsoToUnitary[v] expands an isometry v to an unitary by appending additional columns."
(*KronFactor::usage="KronFactor[m=a\[CircleTimes]b] for a 4\[Cross]4 matrix m (and two 2\[Cross]2 matrices a,b) returns c*a and 1/c*b for some constant c."*)
KronFactorUnitaryDim4::usage="KronFactorUnitaryDim4[u=a\[CircleTimes]b] for a 4\[Cross]4 unitary u (and two 2\[Cross]2 unitaries a,b) returns c*a and 1/c*b for some constant c."
KronFactorVector::usage="KronFactorVector[v=VecProd[a,b]] for a four-dimensional vector v, returns two two-dimensional vectors c*a and 1/c*b for some constant c."
VecProd::usage="VecProd[a,b] generates the tensor product of a,b."
XToYTransform::usage="TBA."
UnitaryEigenvalueDecomp::usage="TBA."
IsoToUnitarySpecial::usage="IsoToUnitarySpecial[v] extends the isometry v to a unitary such that the unitary has as many eigenvalues equal to 1 as possible."
(*Matrix decomposition*)
ZYZDecomposition::usage="ZYZDecomposition[u] for an single-qubit unitary u, returns {a,b,c,d} such that u=e^{i d}.R_z(c).R_y(b).R_z(a)."
ZXZDecomposition::usage="ZXZDecomposition[u] for an single-qubit unitary u, returns {a,b,c,d} such that u=e^{i d}.R_z(c).R_x(b).R_z(a)."
XYXDecomposition::usage="XYXDecomposition[u] for an single-qubit unitary u, returns {a,b,c,d} such that u=e^{i d}.R_x(c).R_y(b).R_x(a)."
XZXDecomposition::usage="XZXDecomposition[u] for an single-qubit unitary u, returns {a,b,c,d} such that u=e^{i d}.R_x(c).R_z(b).R_x(a)."
YXYDecomposition::usage="YXYDecomposition[u] for an single-qubit unitary u, returns {a,b,c,d} such that u=e^{i d}.R_y(c).R_x(b).R_y(a)."
YZYDecomposition::usage="YZYDecomposition[u] for an single-qubit unitary u, returns {a,b,c,d} such that u=e^{i d}.R_y(c).R_z(b).R_y(a)."
RQDecomposition::usage="RQDecomposition[m] for a square matrix m returns {r,q}, where r is an upper triangular matrix and q is an orthogonal/unitary matrix such that r.ConjugateTranspose[q] = m."
QLDecomposition::usage="QLDecomposition[m] for a square matrix m returns {q,l}, where l is a lower triangular matrix and q is an orthogonal/unitary matrix such that ConjugateTranspose[q].l = m."
CSD::usage="CSD[u] returns {m1,m2,m3}, the Cosine-Sine decomposition of u, i.e., u=m1.m2.m3."
(*Decompositions for single-qubit gates*)
ZYZDec::usage="ZYZDec[u,(action)] for an single-qubit unitary u acting on qubit action (default: action=1), returns the ZYZ decompostition in list form."
ZXZDec::usage="ZXZDec[u,(action)] for an single-qubit unitary u acting on qubit action (default: action=1), returns the ZXZ decompostition in list form."
XYXDec::usage="XYXDec[u,(action)] for an single-qubit unitary u acting on qubit action (default: action=1), returns the XYX decompostition in list form."
XZXDec::usage="XZXDec[u,(action)] for an single-qubit unitary u acting on qubit action (default: action=1), returns the XZX decompostition in list form."
YXYDec::usage="YXYDec[u,(action)] for an single-qubit unitary u acting on qubit action (default: action=1), returns the YXY decompostition in list form."
YZYDec::usage="YZYDec[u,(action)] for an single-qubit unitary u acting on qubit action (default: action=1), returns the YZY decompostition in list form."
(*Uniformly controlled rotations*)
DecUCY::usage="See the documentation."
DecUCZ::usage="See the documentation."
(*Quantum Shannon Decomposition*)
DemultiplexUC::usage="See the documentation."
QSD::usage="QSD[v,(action)] decomposes an isometry v from m qubits to n qubits into single-qubit gates and C-NOTs. Optionally, one can determine the qubits the isometry is acting on by providing their numers in teh list action (default: action=Range[n])."
(*Decompose diagonal gate*)
DecDiagGate::usage="DecDiagGate[{\!\(\*SubscriptBox[\(d\), \(1\)]\),...,\!\(\*SubscriptBox[\(d\), \(2^n\)]\)},action] decomposes a diagonal gate on n qubits listed in action (default:action=Range[n]) with diagonal entries \!\(\*SubscriptBox[\(d\), \(1\)]\),...,\!\(\*SubscriptBox[\(d\), \(2^n\)]\)."
(*Decompose two-qubit gates*)
DecUnitary2Qubits::usage="DecUnitary2Qubits[u,(action={n1,n2})] decomposes a two-qubit gate u acting on the two qubits n1 and n2 (default: action={1,2}) (see the full documentation for further options)."
(*Decompositions multi controlled Toffoli gates*)
DecToffoliMultiControlUpToDiagonal::usage="TBA."
DecToffoliUpToDiagonal::usage="TBA."
DecToffoli::usage="TBA."
DecToffoliMultiControl::usage="DecToffoliMultiControl[control,target,n] decomposes a controlled NOT gate on n>=3 qubits controlling on the qubits in the list control and acting on the target qubit target."
(*Decompositions multi controlled single-qubit gates*)
DecMCG::usage="DecToffoliMultiControl[u,control,target,n] decomposes a controlled single-unitary gate u on n qubits controlled on qubits in the list control and acting on the target qubit target."
DecMCSpecialUnitary::usage="TBA."
DecMCSpecialUnitaryUpToDiagonal::usage="TBA."
(*Methods to create and handle multi controlled gates*)
CreateMCToffoli::usage="CreateMCToffoli[control,target,n] creates a matrix corresponding to a multi controlled Toffoli gate on n qubits with controls given in the list control and target qubit nubmer target."
CreateMCG::usage="TBA."
CreateUCG::usage="TBA."
ExtractMCG::usage="TBA."
ExtractUCG::usage="TBA."
(*Methods to decompose uniformly controlled gates*)
DecUCGUpToDiagonal::usage="TBA."
(*Column-by-column decomposition*)
ColumnByColumnDec::usage="ColumnByColumnDec[v,(action)] decomposes the isometry v from m to n qubits on the qubits listed in action (default: action=Range[n]) into single-qubit rotations and C-NOT gates (see the full documentation for additional options)."
(*Decomposition for isometries on a small number of qubits*)
StatePrep1Qubit::usage="TBA."
StatePrep2Qubits::usage="TBA."
StatePrep3Qubits::usage="TBA."
DecIso12::usage="DecIso12[v,(action)] decomposes an isometry v from 1 to 2 qubits acting on the qubits whose numbers are given in the list action (default: action={1,2}) into a sequence of single-qubit and C-NOT gates using a bespoke optimization for this case."
(*Optimal decomposition of an isometry*)
DecIsometry::usage="DecIsometry[v,(action)] decomposes an isometry v from m to n qubits acting on the n qubits whose numbers are given in the list action (default: action=Range[n]) into a sequence of single-qubit and C-NOT gates using the decomposition scheme that achieves the lowest known C-NOT count."
DecIsometryGeneric::usage="DecIsometryGeneric[v,(action)] decomposes an isometry v from m to n qubits acting on the n qubits whose numbers are given in the list action (default: action=Range[n]) into a sequence of single-qubit and C-NOT gates using the decomposition scheme that achieves the lowest known C-NOT count for a generic isometry of the given dimensions."
(*Housholder decomposition*)
DenseHouseholderDec::usage="DenseHouseholderDec[iso] returns a circuit implementing the isometry iso using the dense Householder decomposition."
(*Decompositions for sparse isometries*)
SparseStatePreparation::usage="SparseStatePreparation[vec] returns a circuit implementing sparse state preparation for a sparse state vec."
SparseHouseholderDec::usage="SparseHouseholderDec[iso] returns a circuit implementing the isometry iso using the sparse Householder decomposition."
(*State preparation (Plesch and Brukner)*)
StatePreparation::usage="StatePreparation[v,(action)] prepares an n qubit state (i.e. a \!\(\*SuperscriptBox[\(2\), \(n\)]\)-dimensional normalized column vector) v on the qubits listed in action (default: action=Range[n]) into single-qubit and C-NOT gates (see the full documentation for more details)."
(*Knill's decomposition for isometries*)
KnillDec::usage="KnillDec[v, (action)] decomposes an isometry v from m to n qubits acting on the qubits whose numbers are listed in action (default: action=Range[n]) into a sequence of single-qubit and C-NOT gates based on Knill's decomposition (see the full documentation for more details)."
(*Methods for channels and POVMs*)
KrausToChoi::usage="KrausToChoi[chan] gives the Choi state corresponding to the channel chan."
ChoiToKraus::usage="ChoiToKraus[st] gives the channel corresponding to Choi state st."
MinimizeKrausRank::usage="MinimizeKrausRank[chan] gives another representation of the channel chan potentially with fewer Kraus operators."
StinespringQubit::usage="TBA."
POVMToIsometry::usage="TBA."
(*Decompose channels and POVMs in the quantum ciruit model*)
DecChannelInQCM::usage="TBA."
DecPOVMInQCM::usage="TBA."
IsoFromInstrument::usage="IsoFromInstrument[instr, (actionAndAncilla)] constructs an isometry from instrument instr, which corresponds to the action of the instrument on a wider system, before tracing out and measuring ancilla qubits.
Returns the isometry and the number of ancilla qubits used."
DecInstrumentInQCM::usage="DecInstrumentInQCM[instr, (actionAndAncilla)] decomposes instrument instr from m to n qubits into a sequence of gates using the decomposition scheme that achieves the lowest known C-NOT count.
The instrument uses q ancilla qubits whose numbers are given in the list actionAndAncilla[[1]] (default : range[q]), with q the maximum number of Kraus operators of a channel,
and it acts on the n qubits whose numbers are given in the list actionAndAncilla[[2]] (default: action=Range[q+1,q+n])."
(*Create random objects, such as circuits or instruments*)
PickRandomCircuitIsometry::usage="PickRandomCircuitIsometry[m,n,t] creates a random circuit from m qubits to n qubits with t CNOTs, where an arbitrary single qubit unitary is included at the start and after each CNOT. With Option TotGates->True, the value of t is the total number of gates and these are placed randomly."
PickRandomInstrument::usage="PickRandomInstrument[dim1, dim2, nChan, nKraus] chooses a random instrument from dimension dim1 to dim2 containing num1 channels, each containing num2 Kraus operators.
If nKraus = {nMin, nMax} is a list of two elements, each channel contains a random number n of Kraus operators such that nMin <= n <= nMax (uniform probability).
If nKraus = {...} is a list of three or more elements, each channel contains a random number n of Kraus operators such that n is in nKraus (uniform probability)."
RPickRandomInstrument::usage="RPickRandomInstrument[dim1, dim2, nChan, nKraus] behaves as PickRandomInstrument but gives real output."
FPickRandomInstrument::usage="FPickRandomInstrument[dim1, dim2, nChan, nKraus] behaves as PickRandomInstrument but gives rational output."
PickRandomSparseIsometry::usage="PickRandomSparseIsometry[dim1,dim2] generates a random sparse isometry from dimension dim1 to dim2."
RPickRandomSparseIsometry::usage="PickRandomSparseIsometry[dim1,dim2] generates a random sparse isometry from dimension dim1 to dim2."
PickRandomSparsePsi::usage="PickRandomSparsePsi[dim,s] generates a random pure state with dimension dim and with s non-zero elements."
RPickRandomSparsePsi::usage="RPickRandomSparsePsi[dim,s] generates a random real pure state with dimension dim and with s non-zero elements."
FPickRandomSparsePsi::usage="FPickRandomSparsePsi[dim,s,tol] generates a random analytic real pure state with dimension dim and with s non-zero elements."
(*Various helper methods*)
PrepareForQASM::usage="PrepareForQASM[gatelist] takes a list of gates and prepares it into a form suitable for use with the python script that converts to QASM."
RxRGateDecomp::usage="RxRGateDecomp[u] takes a single qubit unitary u and outputs (a,b,c,d) such that u is equal to Rx[a] followed by R[b,c] up to the phase E^(I*d)."
ReplaceCNOTWithXX::usage="ReplaceCNOTWithXX[st] takes a gate list and replaces all CNOT gates with XX gates and additional single qubit rotations."
ReplaceXXWithCNOT::usage="ReplaceXXWithCNOT[st] takes a gate list and replaces all XX gates with CNOTs and additional single qubit rotations."
CNOTRotationsToXXRGates::usage="CNOTRotationsToXXRGates[st] takes a gate list and replaces all CNOT and single-qubit rotations by XX and R gates."
XXRGatesToCNOTRotations::usage="XXRGatesToCNOTRotations[st] takes a gate list and replaces all XX and R gates with CNOTs and single-qubit rotations."
NearbyIsometry::usage="NearbyIsometry[iso] uses the singular value decomposition to generate an isometry near to iso."
PivotingDec::usage="PivotingDec[vec] returns a circuit which implements the pivoting algorithm for vec, along with the row qubits, column qubits, final column and implemented permutation of the entries of the state vector."
Begin["`Private`"];
(*Set convention for gate types *)
cnotType=0;xType=1;yType=2;zType=3;
measType=4;ancillaType=5;postselType=6;
czType=-1;diagType=-2;xxType=-10;rType=-11;
(* comment in the lines below to do a test of different labels *)
(*cnotType=12;xType=91;yType=92;zType=93;
measType=-5;ancillaType=-15;postselType=0;
czType=25;diagType=42;xxType=1;rType=2;*)
(*Set debug to True to run additional tests during running the methods.*)
debug=False;
(*Set analyzeAnalyticDecUnitary2Qubits to True to get out intermediate results for DecUnitary2Qubits.*)
analyzeAnalyticDecUnitary2Qubits=False;
(*Set analyzeAnalyticCCDec to True to get out intermediate results for ColumnByColumnDec.*)
analyzeAnalyticCCDec=False;
(*Set analyzeAnalyticCCDec to True to get out intermediate results for SimplifyGateList.*)
analyzeAnalyticSimplifyGateList=False;
(*Set analyzeAnalyticCCDec to True to get out intermediate results for QSD.*)
analyzeAnalyticQSD=False;
(*Switch off warnings that are generated by Null arguments (e.g., If[True,,a=5] would return a warning)*)
Off[Syntax::com]
(*----------- Methods to handle and simplify gate sequences in list representation (public)------------*)
(*If one adds support for additional gates to one of the following methods,
one should update all of them.*)
(*Displays each element in a list in matrix form*)
MatrixFormOp[m_]:=If[Length[Dimensions[m]]<=2,If[m=={},{},MatrixForm[m]],MatrixForm[#]&/@m]
Options[CreatePOVMFromGateList]={DropZero->"Last"};
CreatePOVMFromGateList[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:=Module[{mat,i},
IsListForm[st,"CreatePOVMFromGateList"];
If[OptionValue[DropZero]==="None",mat=CreateChannelFromList[DeleteCases[st,{measType,0,_}],n,{POVM->True,DropZero->False}];Map[CT[#] . #&,mat],
If[OptionValue[DropZero]==="All",mat=CreateChannelFromList[DeleteCases[st,{measType,0,_}],n,{POVM->True,DropZero->True}];Map[CT[#] . #&,mat],
If[OptionValue[DropZero]==="Last",mat=CreatePOVMFromGateList[st,n,DropZero->"None"];
For[i=Length[mat],i>=1,i--,If[Chop[mat[[i]]]==0*mat[[i]],mat=Drop[mat,-1],Break[]]];mat,Throw[StringForm["Unknown option for DropZero in CreatePOVMFromGateList. Valid options are None, All or Last"]]]]
]]
Options[NCreatePOVMFromGateList]={DropZero->"Last"};
NCreatePOVMFromGateList[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:=Module[{},
IsListForm[st,"NCreatePOVMFromGateList"];
CreatePOVMFromGateList[NGateList[st],n,{DropZero->OptionValue[DropZero]}]
]
(*Relabel qubits in list numIn in the list format input st with numbers in numOut. Note that numIn must be of the same length as numOut.*)
RelabelQubits[st_,numIn_,numOut_]:=Module[{pos,st2,act},
IsListForm[st,"RelabelQubits"];
st2=st;
Do[
Which[
st2[[i]][[1]]==diagType,
act=Flatten[If[MemberQ[numIn,#],numOut[[Position[numIn,#][[1]]]],#]&/@st2[[i]][[3]]];
st2=ReplacePart[st2,i->{diagType,st2[[i]][[2]],act}]
,
MemberQ[{czType,cnotType},st2[[i]][[1]]],st2=ReplacePart[st2,i->{st2[[i]][[1]],pos=Flatten[Position[numIn,st2[[i]][[2]]]];If[Length[pos]!=0,numOut[[pos]][[1]],st2[[i]][[2]]],pos=Flatten[Position[numIn,st2[[i]][[3]]]];If[Length[pos]!=0,numOut[[pos]][[1]],st2[[i]][[3]]]}],
MemberQ[{xType,yType,zType,measType,ancillaType,postselType},st2[[i]][[1]]],st2=ReplacePart[st2,i->{st2[[i]][[1]],st2[[i]][[2]],pos=Flatten[Position[numIn,st2[[i]][[3]]]];If[Length[pos]!=0,numOut[[pos]][[1]],st2[[i]][[3]]]}],
st2[[i]][[1]]==xxType,act=Flatten[If[MemberQ[numIn,#],numOut[[Position[numIn,#][[1]]]],#]&/@st2[[i]][[3]]];st2=ReplacePart[st2,i->{xxType,st2[[i]][[2]],act}],
st2[[i]][[1]]==rType,st2=ReplacePart[st2,i->{rType,st2[[i]][[2]],pos=Flatten[Position[numIn,st2[[i]][[3]]]];If[Length[pos]!=0,numOut[[pos]][[1]],st2[[i]][[3]]]}],
True,Throw[StringForm["Gate type is not supported by RelabelQubits."]]
],
{i,Length[st]}];
st2]
(*Counts the number of cnot gates in the list st*)
CNOTCount[st_]:=
Module[{i,count},(
IsListForm[st,"CNOTCount"];
count=0;
For[i=1,i<=Length[st],i++,
If[st[[i]][[1]]==cnotType,count++]
];
count
)]
(*
Reverses the elements in the list representation along with changing the signs of the angles in the rotation matrices.
This gives the list representation for the decomposition of a matrix that is the conjugate transpose of the one represented by the input
*)
InverseGateList[st_]:=
Module[{numstr,i,st2},(
IsListForm[st,"InverseGateList"];
st2=Reverse[st];
Do[
Which[
MemberQ[{diagType},st2[[i]][[1]]],st2=ReplacePart[st2,i->{st2[[i]][[1]],Conjugate[st2[[i]][[2]]],st2[[i]][[3]]}],
MemberQ[{czType,cnotType},st2[[i]][[1]]],,
MemberQ[{xType,yType,zType},st2[[i]][[1]]],st2=ReplacePart[st2,i->{st2[[i]][[1]],-st2[[i]][[2]],st2[[i]][[3]]}],
st2[[i]][[1]]==measType,Throw[StringForm["Measurements/Traces in InverseGateList[] are not reversable."]],
st2[[i]][[1]]==ancillaType||st2[[i]][[1]]==postselType,st2=ReplacePart[st2,i->{If[st2[[i]][[2]]==ancillaType,postselType,ancillaType],st2[[i]][[2]],st2[[i]][[3]]}],
st2[[i]][[1]]==xxType,st2=ReplacePart[st2,i->{xxType,-st2[[i]][[2]],st2[[i]][[3]]}],
st2[[i]][[1]]==rType,st2=ReplacePart[st2,i->{rType,{-st2[[i]][[2]][[1]],st2[[i]][[2]][[2]]},st2[[i]][[3]]}],
True,Throw[StringForm["Unknown gate `` in InverseGateList[] cannot be reversed.",st2[[i]]]]
],
{i,Length[st]}
];
st2
)]
(* Adjusts the angles of the gates in the list representation and adjust them to be between 0 and 2\[Pi] *)
AdjustAngle[st_]:=
Module[{numstr,i,st2},(
IsListForm[st,"AdjustAngle"];
If[Length[Dimensions[st]]==1,
If[st=={},
Return[{}],
Return[AdjustAngle[{st}][[1]]]
]
,
st2=st;
];
Do[
If[MemberQ[{xType,yType,zType},st2[[i]][[1]]],st2=ReplacePart[st2,i->{st2[[i]][[1]],AdjustAngleHelp[st2[[i]][[2]]],st2[[i]][[3]]}]],
{i,Length[st]}
];
st2
)]
GateTypes[]:=Module[{},Print["Gate types for UniversalQCompiler"];
Print["{",diagType,",diag,act}: diagonal gate with entries diag on the qubits listed in act"];
Print["{",czType,",n,m}: CZ with qubit n the control and m the target"];
Print["{",cnotType,",n,m}: CNOT with qubit n the control and m the target"];
Print["{",xType,",t,n}: x-rotation by angle t for qubit n"];
Print["{",yType,",t,n}: y-rotation by angle t for qubit n"];
Print["{",zType,",t,n}: z-rotation by angle t for qubit n"];
Print["{",measType,",0,n}: trace out qubit n"];
Print["{",measType,",1,n}: measure qubit n in the computational basis"];
Print["{",ancillaType,",0,n}: qubit n starts in state |0>"];
Print["{",ancillaType,",1,n}: qubit n starts in state |1>"];
Print["{",postselType,",0,n}: qubit n is postselected on |0>"];
Print["{",postselType,",1,n}: qubit n is postselected on |1>"];
Print["{",xxType,",t,{n,m}}: XX-gate with angle t on qubits n,m"];
Print["{",rType,",{t,p},n}: R-gate with angles t,p on qubit n"]]
(*Transforms a list in list format to a list containing the corresponding matrices*)
ListFormToOp[list_,n_:Null]:=Module[{numQubits=n},
IsListForm[list,"ListFormToOp"];
If[n==Null,numQubits=If[Dimensions[list]=={3},NumberOfQubits[{list}],NumberOfQubits[list]]];If[Dimensions[list]=={3},
Which[
MemberQ[{xType,yType,zType},list[[1]]],RotGateM[list[[2]],list[[1]],list[[3]],numQubits],
list[[1]]==cnotType,CNOTM[list[[2]],list[[3]],numQubits],
list[[1]]==czType,CZM[list[[2]],list[[3]],numQubits],
list[[1]]==xxType,XXM[list[[2]],list[[3]][[1]],list[[3]][[2]],numQubits],
list[[1]]==rType,RGateM[list[[2]][[1]],list[[2]][[2]],list[[3]],numQubits],
list[[1]]==diagType,If[Dimensions[list[[2]]]!= {},DiagMat[list[[2]],list[[3]],numQubits],Throw[StringForm["Unspecified diagonal gate cannot be represented as a matrix using ListFormToOp[]"]]],
list[[1]]==measType,Throw[StringForm["Measurements/Trace cannot be represented by matrices using ListFormToOp[]"]],
list[[1]]==ancillaType,Throw[StringForm["Starting in |0> or |1> cannot be represented by matrices using ListFormToOp[]"]],
list[[1]]==postselType,Throw[StringForm["Postselection cannot be represented by matrices using ListFormToOp[]"]]
],ListFormToOp[#,numQubits]&/@list]]
(*Transforms a list in list format to a list containing the corresponding string representations of the gates*)
ListFormToStr[list_]:=
Module[{},
IsListForm[list,"ListFormToStr"];
If[Length[Dimensions[list]]==1,If[list=={},{},ListFormToStrSingleGate[list]],ListFormToStrSingleGate/@list]
]
Options[SimplifyGateList]={FullSimp->True};
SimplifyGateList[st_,OptionsPattern[]]:=Module[{traceout,ancillain,ancillaout,out},
IsListForm[st,"SimplifyGateList"];If[Cases[st,{xxType,_,_}]==={}&&Cases[st,{rType,_,_}]==={},
traceout=Cases[st,{measType,_,_}];ancillain=Cases[st,{ancillaType,_,_}];ancillaout=Cases[st,{postselType,_,_}];
out=Join[ancillain,Reverse[SimplifyGateListReverseGateOrder[Reverse[DeleteCases[st,x_/;x[[1]]==measType||x[[1]]==ancillaType||x[[1]]==postselType]],OptionValue[FullSimp]]],ancillaout,traceout],
Throw[StringForm["Gatetypes ",xxType," and ",rType," cannot be used in SimplifyGateList[]"]]];out]
NSimplifyGateList[st_]:=Module[{},
IsListForm[st,"NSimplifyGateList"];
SimplifyGateList[NGateList[st]]
]
(*----------- Methods to handle and simplify gate sequences in list representation (private)------------*)
(*Create an n qubit isometry from list form. Multiplies the unitaries described in the list (in reversed order!) and outputs the first m columns*)
(* Use FullSimp\[Rule]False to avoid attempts to use FullSimplify *)
Options[CreateIsometryFromList]={FullSimp->True};
CreateIsometryFromList[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:=Module[{mat,mat2,i,k,ancillain,ancillainnums,ancillainvals,ancillaout,ancillaoutnums,ancillaoutvals,st2,id,rest,n1=n},
IsListForm[st];
If[n===Null,n1=NumberOfQubits[st]];ancillain=SortBy[Cases[st,{ancillaType,_,_}],Last];
If[ancillain==={},ancillainnums={},ancillainnums=Transpose[ancillain][[3]];ancillainvals=Transpose[ancillain][[2]]];ancillaout=SortBy[Cases[st,{postselType,_,_}],Last];
If[ancillaout==={},ancillaoutnums={},ancillaoutnums=Transpose[ancillaout][[3]];ancillaoutvals=Transpose[ancillaout][[2]]];
st2=DeleteCases[st,{x_/;x==ancillaType||x==postselType,_,_}];mat={{1}};k=0;
For[i=1,i<=n1,i++,mat=KroneckerProduct[mat,If[MemberQ[ancillainnums,i],k++;KetV[ancillainvals[[k]],2],IdentityMatrix[2]]]];
isAnalytic=True;
For[i=1,i<=Length[st2],i++,
If[isAnalyticGate[st2[[i]]],,isAnalytic=False];
If[isAnalytic,
If[OptionValue[FullSimp],mat=FullSimplifyNoRoots[ListFormToOp[st2[[i]],n1] . mat],mat=Simplify[ListFormToOp[st2[[i]],n1] . mat],Print["CreateIsometryFromList: Error"]],
mat=ListFormToOp[st2[[i]],n1] . mat
]
];
If[ancillaoutnums=={},mat,mat2={{1}};k=0;For[i=1,i<=n1,i++,(* usually ancillaoutvals will be all 1s, so create ket 0, sometimes (for instrument generation) we want to create ket 1 *)mat2=KroneckerProduct[mat2,If[MemberQ[ancillaoutnums,i],k++;KetV[ancillaoutvals[[k]],2],IdentityMatrix[2]]]];
CT[mat2] . mat]]
(*Create an n qubit isometry from list form. Multiplies the unitaries described in the list (numerically) and outputs the first m columns*)
NCreateIsometryFromList[st_,n_:Null]:=CreateIsometryFromList[NGateList[st],n]
Options[CreateChannelFromList]={POVM->False,DropZero->True,FullSimp->True};
CreateChannelFromList[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:=Module[{mat,i,traces,tracesnums,postsel,postselnums,posn,chanout,st2,dims,n1=n},If[Not[OptionValue[POVM]]&&MemberQ[st,{measType,1,_}],Print["CreateChannelFromList: measurement gate type found"]];
If[n===Null,n1=NumberOfQubits[st]];traces=Cases[st,{measType,_,_}];If[traces==={},tracesnums={},tracesnums=Transpose[traces][[3]]];postsel=Cases[st,{postselType,_,_}];If[postsel==={},postselnums={},postselnums=Transpose[postsel][[3]]];If[Dimensions[Intersection[tracesnums,postselnums]]=={0},,Print["CreateChannelFromList: Cannot postselect on zero and measure/trace on the same qubit."]];st2=DeleteCases[st,{measType,_,_}];
mat=CreateIsometryFromList[st2,n1,FullSimp->OptionValue[FullSimp]];
chanout={};posn={};dims={};For[i=1,i<=n1,i++,If[MemberQ[postselnums,i],,If[MemberQ[tracesnums,i],posn=Insert[posn,1,-1];dims=Insert[dims,{1,2},-1],posn=Insert[posn,2,-1];dims=Insert[dims,{2,2},-1]]]];For[i=0,i<=2^(Length[tracesnums])-1,i++,chanout=Insert[chanout,Tensor[BraV[i,2^(Length[tracesnums])],IdentityMatrix[2^(n1-Length[postselnums]-Length[tracesnums])],posn,dims] . mat,-1]];
If[OptionValue[DropZero],For[i=Length[chanout],i>=1,i--,If[Chop[chanout[[i]]]==0*chanout[[i]],chanout=Drop[chanout,{i}]]]];chanout]
Options[NCreateChannelFromList]={POVM->False,DropZero->True};
NCreateChannelFromList[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:=CreateIsometryFromList[NGateList[st],n,{POVM->OptionValue[POVM],DropZero->OptionValue[DropZero]}]
Options[CreateInstrumentFromList]={DropZero->True,FullSimp->True};(* using DropZero here prevents identification using Length[Dimensions[out]], where out is the output of CreateOperationFromGateList *)
CreateInstrumentFromList[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:=Module[{mat,i,j,traces,tracesnums,postsel,postselnums,posn,mmt,mmtnums,inst,chanout,st2,st3,digs,dims,n1=n},
If[n===Null,n1=NumberOfQubits[st]];traces=Cases[st,{measType,0,_}];If[traces==={},tracesnums={},tracesnums=Transpose[traces][[3]]];postsel=Cases[st,{postselType,_,_}];If[postsel==={},postselnums={},postselnums=Transpose[postsel][[3]]];mmt=Cases[st,{measType,1,_}];If[mmt==={},mmtnums={},mmtnums=Transpose[mmt][[3]]];If[Dimensions[Intersection[tracesnums,postselnums,mmtnums]]=={0},,Print["CreateInstrumentFromList: Cannot have combinations of postselect on zero/measure/trace on the same qubit."]];inst={};st2=DeleteCases[st,{measType,1,_}];
For[j=1,j<=2^(Length[mmtnums]),j++,digs=IntegerDigits[j-1,2,Length[mmtnums]];st3=st2;For[i=1,i<=Length[mmtnums],i++,st3=Insert[st3,{postselType,digs[[i]],mmtnums[[i]]},-1]];
inst=Insert[inst,CreateChannelFromList[st3,n1,{DropZero->OptionValue[DropZero],FullSimp->OptionValue[FullSimp]}],-1]];inst]
Options[NCreateInstrumentFromList]={DropZero->True};(* using DropZero here prevents identification using Length[Dimensions[out]], where out is the output of CreateOperationFromGateList *)
NCreateInstrumentFromList[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:=CreateInstrumentFromList[NGateList[st],n,{DropZero->OptionValue[DropZero]}]
(*ToDo: Improve efficiency by implementing the application of C-NOTs and single-qubit rotations efficiently*)
Options[CreateOperationFromGateList]={FullSimp->True};
CreateOperationFromGateList[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:=Module[{four0s,four1s},four1s=MemberQ[st,{measType,1,_}];
If[four1s,CreateInstrumentFromList[st,n,{DropZero->False,FullSimp->OptionValue[FullSimp]}],four0s=MemberQ[st,{measType,0,_}];
If[four0s,CreateChannelFromList[st,n,FullSimp->OptionValue[FullSimp]],CreateIsometryFromList[st,n,FullSimp->OptionValue[FullSimp]]]]]
NCreateOperationFromGateList[st_,n_: Null]:=Module[{four0s,four1s},four1s=MemberQ[st,{measType,1,_}];
If[four1s,NCreateInstrumentFromList[st,n,DropZero->False],four0s=MemberQ[st,{measType,0,_}];
If[four0s,NCreateChannelFromList[st,n],NCreateIsometryFromList[st,n]]]]
ListFormToStrSingleGate[list_]:=Switch[list[[1]],
diagType,ToString[StringForm["D(`1`)(`2`)",list[[2]],list[[3]] ]],
czType,CZNotationStr[list[[2]],list[[3]]],
cnotType,CNOTNotationStr[list[[2]],list[[3]]],
xType,ToString[StringForm["Rx(`2`)(`1`)",list[[3]],list[[2]] ]],
yType,ToString[StringForm["Ry(`2`)(`1`)",list[[3]],list[[2]] ]],
zType,ToString[StringForm["Rz(`2`)(`1`)",list[[3]],list[[2]] ]],
measType,If[list[[2]]==0,ToString[StringForm["Tr(`1`)",list[[3]] ]],ToString[StringForm["M(`1`)",list[[3]]]]],
ancillaType,ToString[StringForm["|`1`>(`1`)",list[[2]],list[[3]]]],
postselType,ToString[StringForm["<`1`|(`1`)",list[[2]],list[[3]]]],
xxType,ToString[StringForm["XX(`1`)(`2`,`3`))",list[[2]],list[[3,1]],list[[3,2]]]],
rType,ToString[StringForm["RGate(`1`,`2`)(`3`))",list[[2,1]],list[[2,2]],list[[3]]]],
_,Throw[StringForm["Unknown gate type found in ListFormToStr[]."]]
]
(*Input: st: gate sequence containing only single-qubit and C-NOT gates, n: total number of qubits.
Merges the single-qubit gates if this helps to reduce the gate count.
Commutes Rz rotation with the control and Rx rotations with the target of C-NOT gates. Moreover, two following C-NOT gates
are canceled out.
Output: Simplified gate sequence.*)
SimplifyGateListReverseGateOrder[st_,FullSimp_:True] := Module[{numQubits, maxNum,A,stNew,controlQubit,targetQubit,mergedContr,commutedContr,mergedTar,commutedTarg},(
numQubits=NumberOfQubits[st];
A=ConstantArray[{},numQubits]; (*List to save the collected single-qubit gates on the n qubits. The list of a qubit is updated after merging.*)
stNew={};
Do[
If[st[[i]][[1]]==cnotType,
(*case: C-NOT gate*)
controlQubit=st[[i]][[2]];
targetQubit=st[[i]][[3]];
If[analyzeAnalyticSimplifyGateList,Print["CNOT gate with control qubit ",controlQubit," and target qubit ", targetQubit," found."]];
If[analyzeAnalyticSimplifyGateList,Print["Call MergeAndCommuteSQG with input ",{A[[controlQubit]],False,True}]];
{mergedContr,commutedContr}=MergeAndCommuteSQG[A[[controlQubit]],False,True,FullSimp];
If[analyzeAnalyticSimplifyGateList,Print["Merged gates on control qubit ",controlQubit," are given by ",mergedContr]];
If[analyzeAnalyticSimplifyGateList,Print["Commuted gates on control qubit ",controlQubit," are given by ",commutedContr]];
If[analyzeAnalyticSimplifyGateList,Print["Call MergeAndCommuteSQG with input ",{A[[targetQubit]],True,False}]];
{mergedTar,commutedTarg}=MergeAndCommuteSQG[A[[targetQubit]],True,False,FullSimp];
If[analyzeAnalyticSimplifyGateList,Print["Merged gates on target qubit ",targetQubit," are given by ",mergedTar]];
If[analyzeAnalyticSimplifyGateList,Print["Commuted gates on target qubit ",targetQubit," are given by ",commutedTarg]];
A[[controlQubit]]={};
A[[targetQubit]]={};
(*Add merged gates:*)
stNew=Join[stNew,mergedContr];
stNew=Join[stNew,mergedTar];
(*Append C-NOT:*)
stNew=AppendCNOT[stNew,st[[i]]];
(*If two C-NOTs were cancelled out, we have to add the last single-qubit gates again to the array A*)
If[Length[stNew]>=1&&stNew[[-1]][[1]]!=cnotType,
If[analyzeAnalyticSimplifyGateList,Print["Call ExtractSQG with input ",{stNew,controlQubit}]];
{stNew,A[[controlQubit]]}=ExtractSQG[stNew,controlQubit];
If[analyzeAnalyticSimplifyGateList,Print["Call ExtractSQG with input ",{stNew,targetQubit}]];
{stNew,A[[targetQubit]]}=ExtractSQG[stNew,targetQubit];
];
(*Upade gate array:*)
A[[controlQubit]]=Join[A[[controlQubit]],commutedContr];
A[[targetQubit]]=Join[A[[targetQubit]],commutedTarg];
,
(*case: single-qubit gate acting on st[[i]][[3]]*)
AppendTo[A[[st[[i]][[3]]]],st[[i]]];
]
,{i,1,Length[st]}
];
(*Add the single-qubit gates at the beginning of the circuit to stNew*)
Do[
If[analyzeAnalyticSimplifyGateList,Print["Call MergeAndCommuteSQG with input ", {A[[i]],False,False}]];
stNew=Join[stNew,MergeAndCommuteSQG[A[[i]],False,False,FullSimp][[1]]]
,{i,1,numQubits}
];
If[analyzeAnalyticSimplifyGateList,Print["Call AdjustAngle with input ", stNew]];
AdjustAngle[stNew]
)]
(*Gets a sequence of single-qubit gates as an input and the information about the previous action of the two qubit gate: control or target of
the C-NOT gate. If "isTarget" and "isControl" are wrong, then we assume that we can not commute anything to the left in the circuit.
Outputs the merged single-qubit gates as well as the gate that can be commuted trough the previous two-qubit gate.*)
MergeAndCommuteSQG[stSQGsInp_,isTarget_,isControl_,FullSimp_:True] := Module[{angles,stSQGs,mergedGates,commutedGate,gate},(
If[isTarget&&isControl,Throw[StringForm["isTarget and isControl are set to True in the method MergeAndCommuteSQG[]. This setting is impossible and not allowed as an input."]];Return[{}]];
(*Merge gates if we have two single-qubit rotations of the same kind:*)
If[analyzeAnalyticSimplifyGateList,Print["Call MergeSameRot with input ",stSQGsInp]];
stSQGs=MergeSameRot[stSQGsInp,FullSimp];
mergedGates={};
commutedGate={};
If[Length[stSQGs]>= 1 && Length[stSQGs]<= 2,
(*Commute single-qubit gate if possible*)
If[isTarget&&stSQGs[[-1]][[1]]==xType,
commutedGate={stSQGs[[-1]]};
stSQGs=Delete[stSQGs,-1]
,
If[isControl&&stSQGs[[-1]][[1]]==zType,
commutedGate={stSQGs[[-1]]};
stSQGs=Delete[stSQGs,-1];
];
];
(*Define merged gates as the ones that were not commuted*)
mergedGates=stSQGs;
];
If[Length[stSQGs]>= 3,
If[analyzeAnalyticSimplifyGateList,Print["Call MultiplySQGate with input ",stSQGs]];
gate=MultiplySQGates[stSQGs,FullSimp];
If[isTarget,
If[analyzeAnalyticSimplifyGateList,Print["Call XYXDecomposition with input ",gate]];
angles=XYXDecomposition[gate];
If[analyzeAnalyticSimplifyGateList,Print["Call of XYXDecomposition finished. "]];
angles=Reverse[angles];
If[Chop@angles[[3]]==0,
mergedGates={};
commutedGate=MergeSameRot[{{xType,If[FullSimp,FullSimplifyNoRoots[Chop@angles[[4]]+Chop@angles[[2]]],Simplify[Chop@angles[[4]]+Chop@angles[[2]]]],stSQGs[[-1]][[3]]}},FullSimp];
,
mergedGates=MergeSameRot[{{xType,Chop@angles[[2]],stSQGs[[-1]][[3]]},{yType,Chop@angles[[3]],stSQGs[[-1]][[3]]}},FullSimp];
commutedGate=MergeSameRot[{{xType,Chop@angles[[4]],stSQGs[[-1]][[3]]}},FullSimp];
]
,
If[analyzeAnalyticSimplifyGateList,Print["Call ZYZDecomposition with input ",gate]];
angles=ZYZDecomposition[gate];
If[analyzeAnalyticSimplifyGateList,Print["Call of XYXDecomposition finished. "]];
angles=Reverse[angles];
If[isControl,
If[Chop@angles[[3]]==0,
mergedGates={};
commutedGate=MergeSameRot[{{zType,If[FullSimp,FullSimplifyNoRoots[Chop@angles[[4]]+Chop@angles[[2]]],Simplify[Chop@angles[[4]]+Chop@angles[[2]]]],stSQGs[[-1]][[3]]}},FullSimp]
,
mergedGates=MergeSameRot[{{zType,Chop@angles[[2]],stSQGs[[-1]][[3]]},{yType,Chop@angles[[3]],stSQGs[[-1]][[3]]}},FullSimp];
commutedGate=MergeSameRot[{{zType,Chop@angles[[4]],stSQGs[[-1]][[3]]}},FullSimp]
]
,
(*The following line is executed in the case where there is no control or target previously (e.g., at the start of the circuit)*)
mergedGates=MergeSameRot[{{zType,Chop@angles[[2]],stSQGs[[-1]][[3]]},{yType,Chop@angles[[3]],stSQGs[[-1]][[3]]},{zType,Chop@angles[[4]],stSQGs[[-1]][[3]]}},FullSimp];
];
];
];
{mergedGates,commutedGate}
)]
(*Merge the single-qubit rotations of the same kind in a list stSQGs of single-qubit rotations. Remove
rotations that rotate with angle zero.*)
MergeSameRot[stSQGs_,FullSimp_:True] := Module[{stSQGsSimplified,counter},(
If[Length[stSQGs]==0,
{}
,
If[Length[stSQGs]==1,
If[rotIsZero[stSQGs[[1]]],{},stSQGs]
,
stSQGsSimplified={};
Do[
If[Length[stSQGsSimplified]>=1,
If[stSQGsSimplified[[-1]][[1]]==stSQGs[[i]][[1]],
stSQGsSimplified[[-1]]={stSQGsSimplified[[-1]][[1]],If[FullSimp,FullSimplifyNoRoots[stSQGsSimplified[[-1]][[2]]+stSQGs[[i]][[2]]],Simplify[stSQGsSimplified[[-1]][[2]]+stSQGs[[i]][[2]]]],stSQGsSimplified[[-1]][[3]]};
If[rotIsZero[stSQGsSimplified[[-1]]],stSQGsSimplified=Delete[stSQGsSimplified,-1]]
,
If[rotIsZero[stSQGs[[i]]],,AppendTo[stSQGsSimplified,stSQGs[[i]]]]
]
,
If[rotIsZero[stSQGs[[i]]],,AppendTo[stSQGsSimplified,stSQGs[[i]]]]
]
,
{i,1,Length[stSQGs]}
];
stSQGsSimplified
]
]
)]
(*Append a C-NOT gate and cancel it if two are following each other*)
AppendCNOT[stIn_,cnot_]:=Module[{ctr,tar,st,sameContr,sameTarget},(
st=stIn;
ctr=cnot[[2]];
tar=cnot[[3]];
Do[
If[st[[-j]][[1]]!=cnotType&&(st[[-j]][[3]]==tar||st[[-j]][[3]]==ctr),AppendTo[st,cnot];Goto[End]];
If[st[[-j]][[1]]==cnotType&&(st[[-j]][[2]]==tar||st[[-j]][[2]]==ctr||st[[-j]][[3]]==tar||st[[-j]][[3]]==ctr),
(*If we have another C-NOT with the same target qubit, but another control qubit, or another C-NOT
with the same control qubit and anothertarget qubit they commute
(and we should search further for C-NOTs that may cancel.)*)
If[(st[[-j]][[2]]==ctr&&st[[-j]][[3]]!= tar)||(st[[-j]][[2]]!= ctr &&st[[-j]][[3]]== tar)
,
,
If[st[[-j]][[2]]==ctr&&st[[-j]][[3]]==tar,
st=Delete[st,{-j}];Goto[End],
(*Case wher the considered C-NOT cannot be cancelled and does not commute*)
AppendTo[st,cnot];
Goto[End]
]
]
]
,{j,1,Length[st]}];
AppendTo[st,cnot];
Label[End];
st
)];
(*Extracts the single-qubit gates at the qubit with number quNum until the next C-NOT gate*)
ExtractSQG[stIn_,quNum_]:=Module[{st,SQGs,toDelete},(
st=stIn;
SQGs={};
toDelete={};
Do[
If[(st[[-j]][[1]]==cnotType&&(st[[-j]][[2]]==quNum||st[[-j]][[3]]==quNum))
||(st[[-j]][[3]]==quNum)
,
If[st[[-j]][[1]]==cnotType,
Goto[End],
AppendTo[SQGs,st[[-j]]];
AppendTo[toDelete,{-j}]
]
]
,{j,1,Length[st]}];
Label[End];
st=Delete[st,toDelete];
{st,Reverse[SQGs]}
)];
(*---------------------------------Some methods to generate Matrices (public)---------------------------------*)
(*Rotation matrices*)
RxM[\[Alpha]_,i_:1,n_:1]:=RotGateM[\[Alpha],xType,i,If[n<=i,i,n]]
RyM[\[Alpha]_,i_:1,n_:1]:=RotGateM[\[Alpha],yType,i,If[n<=i,i,n]]
RzM[\[Alpha]_,i_:1,n_:1]:=RotGateM[\[Alpha],zType,i,If[n<=i,i,n]]
(*Generates a CNOT gate with control on i, target on j and total number of qubits n.
Note that for inputs i and j must be labelled according to the notation -least significant qubit is labelled n and most significant qubit is labelled 1. *)
CNOTM[i_,j_,n1_:Null]:=
Module[{n,iden,gate,sup,sdown},(
n=If[n1===Null,Max[i,j],n1];
sup={1,0};
sdown={0,1};
iden=Table[IdentityMatrix[2],{n}];
gate=KroneckerProduct@@ReplacePart[iden,i->Outer[Times,sup,sup]]+KroneckerProduct@@ReplacePart[iden,{i->Outer[Times,sdown,sdown],j->PauliMatrix[1]}];
gate
)]
(*Generates a controlled z gate with control on i and target on j, wher n is the total number of qubits (for controlled z, control and target don't matter, i.e, are interchangeable).
Note that for the input to this gate, i and j must be labelled according to the notation -least significant qubit is labelled n and most significant qubit is labelled 1. *)
CZM[i_,j_,n1_:Null]:=
Module[{n,iden,gate,sup,sdown},(
n=If[n1===Null,Max[i,j],n1];
sup={1,0};
sdown={0,1};
iden=Table[IdentityMatrix[2],{n}];
gate=KroneckerProduct@@ReplacePart[iden,i->Outer[Times,sup,sup]]+KroneckerProduct@@ReplacePart[iden,{i->Outer[Times,sdown,sdown],j->PauliMatrix[3]}];
gate
)]
XXM[phi_,i_:Null,j_:Null,n_:Null]:=Module[{out,sys,n1,i1=i,j1=j},If[i1===Null,If[j1===Null,i1=1;j1=2,If[j1==1,i1=2,i1=1]]];
If[j1===Null,If[i1==2,j1=1,j1=2]];
If[n===Null,n1=Max[i1,j1],n1=n];
out={{Cos[phi],0,0,-I*Sin[phi]},{0,Cos[phi],-I*Sin[phi],0},{0,-I*Sin[phi],Cos[phi],0},{-I*Sin[phi],0,0,Cos[phi]}};
sys=Insert[Insert[Table[1,n1-2],2,Min[i1,j1]],2,Max[i1,j1]];
Tensor[IdentityMatrix[2^n1/4],out,sys,Table[2,n1]]]
RGateM[th_,phi_,i_:1,n_:Null]:=Module[{out,sys,n1},
If[n===Null,n1=i,n1=n];
out={{Cos[th/2],-I*E^(-I*phi)*Sin[th/2]},{-I*E^(I*phi)*Sin[th/2],Cos[th/2]}};
sys=Insert[Table[1,n1-1],2,i];Tensor[IdentityMatrix[2^n1/2],out,sys,Table[2,n1]]]
(*Created the diagonal matrix with diagonal entries in diag and acting on the qubits listed in
the list act, where we consider n qubits in total.*)
DiagMat[diag_,act1_:Null,n1_:Null]:=
Module[{n,act,diagReordered,newPos,dimDiag,dimId,pos,i,actSorted},(
n=If[n1===Null,If[act1===Null,Log2[Length[diag]],Max[act1]],n1];
act=If[act1===Null,Range[n],act1];
actSorted=Sort[act];
If[Log2[Length[diag]]>1,
newPos=Table[Position[actSorted,i][[1,1]],{i,act}];
diagReordered=Flatten[ExchangeSystems[Transpose[{diag}],newPos,ConstantArray[2,Log2[Length[diag]]]]],
diagReordered=diag
];
dimDiag=2^(Length[actSorted]);
dimId=2^(n-Length[actSorted]);
pos={};
For[i=1,i<=n,i++,
If[MemberQ[actSorted,i],
AppendTo[pos,1],
AppendTo[pos,2]
]
];
Tensor[DiagonalMatrix[diagReordered], IdentityMatrix[dimId], pos, ConstantArray[2,n]]
)]
(*For an isometry from m to n qubits (with m<n), this method adds extra columns to the matrix to make it an unitary on n qubits.*)
IsoToUnitary[m_]:=
Module[{vec,mat,dim,norm},(
dim=Dimensions[m];
mat=m\[Transpose];
mat=Join[mat,Table[ConstantArray[0,dim[[1]]],{dim[[1]]-dim[[2]]}]];
vec=Orthogonalize[NullSpace[mat]];
(*norm=Norm/@vec;
vec=vec;*)
mat[[-Length[vec];;-1]]=ConjSimplify[vec];
mat\[Transpose]
)]
(*Take a rectangular matrix with more columns than rows and adds zero rows until it is square *)
FillZero[r_] :=
Module[{dr, dc}, {dr, dc} = Dimensions[r];
If[dc - dr > 0, Join[r, ConstantArray[0, {dc - dr, dc}]], r]]
(*Extract angles from rotation matrices*)
(*For a Rz(\[Theta]) matrix, returns the value of \[Theta]*)
RzAngle[u_]:=
Module[{th},(
th=2*Arg[u[[1,1]]];
th//Chop
)]
(*For a Ry(\[Theta]) matrix, returns the value of \[Theta]*)
RyAngle[u_]:=
Module[{th},(
th=2*ArcTan2[u[[1,1]],u[[1,2]]];
th//Chop
)]
(*For a Rx(\[Theta]) matrix, returns the value of \[Theta]*)
RxAngle[u_]:=
Module[{th},(
th=2*ArcTan2[u[[1,1]],Im@u[[1,2]]];
th//Chop
)]
(*----------------------------------------Apply gates efficiently (public)---------------------------------*)
(*Apply a diagonal gate diagGate={diagonalEntries,actionQubits,bits} efficiently to matrix.*)
(*ToDo: Improve efficiency by parallelizing the computation.*)
ApplyDiag[diagGate_,mat_]:=Module[{n,diag,i,rowNum,act,diagIndex},(
rowNum=Dimensions[mat][[1]];
act=diagGate[[2]];
n=diagGate[[3]];
(*Construct diagonal*)
diag={};
For[i=1,i<= rowNum,i++,
diagIndex=FromDigits[IntegerDigits[i-1,2,n][[act]],2]+1;
AppendTo[diag,diagGate[[1]][[diagIndex]]]
];
diag*mat
)]
(*Apply a multi controlled gate MCG={gate,oneControls,zeroControls,target,bits} efficiently to a matrix mat
(controls on zeros are also supported).*)
ApplyMCG[MCG_,mat_]:=Module[{inPos,j,con,gate,bits,target,control, ind1,ind2,mat2,index,diag,i,colNum,act,diagIndex,freeQubits,binaryInd,controlZero,controlOne},(
colNum=Dimensions[mat][[2]];
gate=MCG[[1]];
controlOne=MCG[[2]];
controlZero=MCG[[3]];
bits=MCG[[5]];
target=MCG[[4]];
control=Sort[Join[controlOne,controlZero]];
freeQubits=bits-Length[control]-1;
(*Construct a list of indexes where the controled qubits acts on (ignoring the bit of the target qubit)*)
binaryInd=Tuples[{1, 0}, freeQubits];
For[i=1,i<=Length[control],i++,
con=control[[i]];
If[con> target,inPos=con-1,inPos=con];
If[MemberQ[controlOne,con],
binaryInd=Insert[#,1,inPos]&/@binaryInd,
binaryInd=Insert[#,0,inPos]&/@binaryInd
];
];
mat2=mat;
Do[
ind1=FromDigits[Insert[index,0,target],2]+1;
ind2=FromDigits[Insert[index,1,target],2]+1;
For[j=1,j<=colNum,j++,
{{mat2[[ind1]][[j]]},{mat2[[ind2]][[j]]}}=gate . {{mat[[ind1]][[j]]},{mat[[ind2]][[j]]}};
];
,{index,binaryInd}];
mat2
)]
(*Efficient way to apply a uniformly controlled gate UCG={gates,control,target,bits} to a matrix mat.*)
(*ToDo: Improve efficiency further by parallelizing the loop in parallel (unfortunately, ParallelDo seems not to work properly).*)
ApplyUCG[UCG_,mat_]:=Module[{newPos,orderForSQGs,controlSorted,gatesOrdered,mat2,colNum,rowNum,bits,rep,gateIndex,gates,target, control,spacing,i},(
{gates,control,target,bits}=UCG;
controlSorted=Sort[control];
If[Length[control]>1,
newPos=Table[Position[controlSorted,i][[1,1]],{i,control}];
orderForSQGs=Flatten[ExchangeSystems[Transpose[{Range[2^Length[control]]}],newPos,ConstantArray[2,Length[control]]]];
gatesOrdered=gates[[Flatten[orderForSQGs]]],
gatesOrdered=gates
];
spacing=2^(bits-target);
rep=2^(bits-Length[control]-1);
mat2=ConstantArray[0,Dimensions[mat]];(*ToDo: Improve efficiency by not copying matrix*)
Do[
i=Floor[j/spacing]*spacing+j;
gateIndex=FromDigits[IntegerDigits[i,2,bits][[control]],2]+1;
{mat2[[i+1]],mat2[[i+1+spacing]]}=gates[[gateIndex]] . {mat[[i+1]],mat[[i+1+spacing]]};
,{j,0,2^(bits-1)-1}];
mat2
)]
(*----------------------------------------Visualization of cirucits (pucblic)---------------------------------*)
(*Export a gates sequence st (given in list form) to Latex using the package QCircuit.*)
Options[LatexQCircuit]={AnglePrecision->0};
LatexQCircuit[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:= Module[{draw=False},If[OptionValue[AnglePrecision]>0,draw=True];IsListForm[st,"LatexQCircuit"];qCircuitFromGridForm[computeGridForm[st, n],{DrawRotationAngles->draw,Digits->OptionValue[AnglePrecision]}]]
(*Print a gate sequence st (given in list form)*)
Options[PrintCircuit]={AnglePrecision->0};
PrintCircuit[st_,n:Except[_?OptionQ]:Null,OptionsPattern[]]:= Module[{draw=False},If[OptionValue[AnglePrecision]>0,draw=True];IsListForm[st,"PrintCircuit"];drawGraphFromGridForm[computeGridForm[st, n],{DrawRotationAngles->draw,Digits->OptionValue[AnglePrecision]}]]
(*----------------------------------------Visualization of cirucits (private)---------------------------------*)
computeGridForm[st_, n_:Null(*,maxGatesPerLine_:Infinity*)]:= Module[{postSelSetectingOps,initializingOps,gate,maxNum,numQubits,gatesOnEachWireSoFar,type, controlOrParameter, target,controlWireLength,
targetWireLength,controlGate,shortestWire,longestWireLength,blankMustBeLeftBlankGate,blankWhichCanBeOverwrittenGate,minIndex, maxIndex,index,leftMostPositionThatCanBeFilled,overWriteOrAppend,
positionForNewGate,positionsOfblackMustBeLeftBlankGateOnTarget,positionsOfblackMustBeLeftBlankGateOnControl,setActingQubits,positionsOfblackMustBeLeftBlankGate,wire},
(*set some relevant constants*)
initializingOps={};
postSelSetectingOps={};
controlGate = "control";
blankWhichCanBeOverwrittenGate = "blankWhichCanBeOverwritten";
blankMustBeLeftBlankGate = "mustBeLeftBlank";
leftMostPositionThatCanBeFilled[wire_] :=Module[{fp},fp=FirstPosition[wire,blankWhichCanBeOverwrittenGate,Length[wire] +1];(*We want to give back a number, but FirstPosition gives back a list in genearal, and the number one if the input is the empty set*)If[Length[fp]==0,,fp=fp[[1]]];fp];
(*Set the number of qubits if not already set as an input*)
maxNum=NumberOfQubits[st];
numQubits = Switch[n,
Null, maxNum,
_, n];
If[numQubits<maxNum,Print["Error in computeGridForm[]. The total number of qubits must be larger than the largest qubit number the input circuit is acting on. Null was returned."];Return[Null]];
(*Initialize grid*)
gatesOnEachWireSoFar = ConstantArray[{}, numQubits];
(*Add gates:*)
Do[
{type, controlOrParameter, target} = gate;
If[type == czType || type == cnotType||type == xxType || type == diagType,
(*Code for multiqubit gates*)
If[type==diagType ||type == xxType,
{minIndex, maxIndex}={Min[target],Max[target]}
,
{minIndex, maxIndex} = If[target < controlOrParameter, {target,controlOrParameter} , {controlOrParameter, target}];
];
positionForNewGate = Max[Map[leftMostPositionThatCanBeFilled, gatesOnEachWireSoFar[[minIndex;;maxIndex]]]];
(*Fill in the grid for the entries between the control and the target (or the first and last qubit the diagonal gate is acting on)*)
For[index = minIndex+1, index < maxIndex, index++,
gatesOnEachWireSoFar[[index]] = PadRightHelp[gatesOnEachWireSoFar[[index]] , positionForNewGate, blankWhichCanBeOverwrittenGate];
gatesOnEachWireSoFar[[index, positionForNewGate]] = blankMustBeLeftBlankGate;
];
(*Fill up all entries where the multi-qubit gate is acting on with "blankMustBeLeftBlankGate" until the current gate in two steps:*)
If[type==diagType ||type == xxType,
setActingQubits=target,
setActingQubits={controlOrParameter,target}
];
Do[
positionsOfblackMustBeLeftBlankGate=
(*Step 1: Replace "blankWhichCanBeOverwrittenGate" entries*)Flatten[Position[gatesOnEachWireSoFar[[actQubit]],blankWhichCanBeOverwrittenGate]];
gatesOnEachWireSoFar[[actQubit]][[positionsOfblackMustBeLeftBlankGate]]= blankMustBeLeftBlankGate;
(*Step 2: Fill up with " blankMustBeLeftBlankGate" entries*)
gatesOnEachWireSoFar[[actQubit]] = PadRightHelp[gatesOnEachWireSoFar[[actQubit]] , positionForNewGate, blankMustBeLeftBlankGate];
gatesOnEachWireSoFar[[actQubit,positionForNewGate]] = {type, controlOrParameter,target};
,
{actQubit,setActingQubits}];
,
(*Initializing and post selection operations*)
If[type==ancillaType,
AppendTo[initializingOps,gate]
,
If[type==postselType,
AppendTo[postSelSetectingOps,gate]
,
(*Code for single qubit gates*)
positionForNewGate = leftMostPositionThatCanBeFilled[gatesOnEachWireSoFar[[target]]];
If[positionForNewGate == Length[gatesOnEachWireSoFar[[target]]] +1,
AppendTo[gatesOnEachWireSoFar[[target]], {type, controlOrParameter,target}],
gatesOnEachWireSoFar[[target]][[positionForNewGate]]= {type, controlOrParameter,target}
]
]
]
]
,{gate, st}];
longestWireLength = Max[Map[Length, gatesOnEachWireSoFar]];
For[wire = 1, wire <= Length[gatesOnEachWireSoFar], wire ++,
gatesOnEachWireSoFar[[wire]] = PadRightHelp[gatesOnEachWireSoFar[[wire]],longestWireLength+2,"mustBeLeftBlank" ,1];
];
(*ToDo: Finish:
(*Restructure gridForm (use several rows for the representation)*)
gridLength=Length[gatesOnEachWireSoFar[[1]]];
numOfRows=Ceiling[gridLength/maxGatesPerLine];
Do[
If[i\[Equal]0,
grid=gridForm[[All,i*24+1;;(i+1)*24]],
grid=Join[grid,ConstantArray["mustBeLeftBlank",{1,24}]];
grid=Join[grid,gridForm[[All,i*24+1;;(i+1)*24]]]
],
{i,{0,1}}
]
*)
Return[{gatesOnEachWireSoFar,initializingOps,postSelSetectingOps}];
]
PadRightHelp[list_,position_,entry_,margin_:0]:=Module[{out},
If[position>=Length[list],
out=PadRight[list,position,entry,margin],
out=list;
];
out
]
Options[drawGraphFromGridForm]={DrawRotationAngles->False,Digits->2};
drawGraphFromGridForm[gridForm_,OptionsPattern[]]:= Module[{post,init,out,postSelSet,initSet,newTraceOuts,traceOutWires,positions,edgeRenderingFunction,actQubits,text, posTemp,numWires,longestWireLength,vertexCoordRules,edges,wire, gateIndex,sortTarget,box,vertexRenderingFunction,
vertexName,vertexCoordinates,vertices,vetexShapeFunction,edgeShapeFunction},
longestWireLength = Max[Map[Length, gridForm[[1]]]];
numWires = Length[gridForm[[1]]];
initSet=gridForm[[2]];
postSelSet=gridForm[[3]];
(*Create vertexCoordinates. As of mathematica 12 this information takes the form of a list rather than a function.
So we also explicitly create a list of vertices as well as one of their coordinates to ensure the two lists are ordered the same.*)
vertices = Flatten[Table[{w,g}, {w,1,numWires}, {g,1,longestWireLength}],1];
vertexCoordinates = Map[Function[{#[[2]],-#[[1]]}], vertices]; (*we get the y coordinate from the wire the gate is on lowest wire at the top*)
edges = {};
For[wire =1, wire <= Length[gridForm[[1]]], wire++,
(*edges between adjacent gates on the same wire*)
For[gateIndex = 1, gateIndex < Length[gridForm[[1]][[wire]]], gateIndex++,
AppendTo[edges, {wire, gateIndex}-> {wire, gateIndex+1}];
];
(*edges for controls*)
For[gateIndex = 1, gateIndex <= Length[gridForm[[1]][[wire]]], gateIndex++,
If[Head[gridForm[[1]][[wire, gateIndex]]]=== List &&
(gridForm[[1]][[wire, gateIndex]][[1]]==cnotType||gridForm[[1]][[wire, gateIndex]][[1]]==czType)&&
gridForm[[1]][[wire, gateIndex]][[3]]==wire,
AppendTo[edges, {wire, gateIndex}-> {gridForm[[1]][[wire, gateIndex]][[2]],gateIndex}];
];
];
];
vetexShapeFunction[{x_,y_},v_,{w_,h_}] := Module[{vsWireIndex=v[[1]], vsGateIndex=v[[2]],pos={x,y}, gate,type, controlOrParameter,target},
gate = gridForm[[1]][[vsWireIndex, vsGateIndex]];
(*If a qubits starts in a fixed state, draw \ket{0}*)
init=If[vsGateIndex==1&&MemberQ[initSet,{_,_,vsWireIndex}],
If[MemberQ[initSet,{_,0,vsWireIndex}],
{Black, Text["|0>", pos-{0.2,0}]},
{Black, Text["|1>", pos-{0.2,0}]}
],
{}
];
(*If a qubits are postSelSetected in a fixed state, draw \bra{0}*)
post= If[vsGateIndex==longestWireLength&&MemberQ[postSelSet,{_,_,vsWireIndex}],
If[MemberQ[postSelSet,{_,0,vsWireIndex}],
{Black, Text["<0|", pos+{0.2,0}]},
{Black, Text["<1|", pos+{0.2,0}]}
],
{}
];
out=Which[TrueQ[gate == "mustBeLeftBlank"],{},TrueQ[gate == "blankWhichCanBeOverwritten"],{},
{type, controlOrParameter,target} = gate;TrueQ[(type==cnotType ||type==czType)&& controlOrParameter==-pos[[2]]],
(*Control is recognized (from a gate type cnotType or czType)*)
{Black, Disk[pos, 0.1]},
type == xType,
{White, EdgeForm[Black], Rectangle[pos - {0.2,0.2}, pos + {0.2,0.2}], Black, Text[If[OptionValue[DrawRotationAngles],ToStringStandard[StringForm["Rx(`1`)",DrawCircuitNumberFormat[controlOrParameter,OptionValue[Digits]]]],"Rx"], pos]},
type == yType,
{White, EdgeForm[Black], Rectangle[pos - {0.2,0.2}, pos + {0.2,0.2}], Black, Text[If[OptionValue[DrawRotationAngles],ToStringStandard[StringForm["Ry(`1`)",DrawCircuitNumberFormat[controlOrParameter,OptionValue[Digits]]]],"Ry"], pos]},
type == zType,
{White, EdgeForm[Black], Rectangle[pos - {0.2,0.2}, pos + {0.2,0.2}], Black, Text[If[OptionValue[DrawRotationAngles],ToStringStandard[StringForm["Rz(`1`)",DrawCircuitNumberFormat[controlOrParameter,OptionValue[Digits]]]],"Rz"], pos]},
type == rType,
{White, EdgeForm[Black], Rectangle[pos - {0.2,0.2}, pos + {0.2,0.2}], Black, Text[If[OptionValue[DrawRotationAngles],ToStringStandard[StringForm["R(`1`,`2`)",DrawCircuitNumberFormat[controlOrParameter[[1]],OptionValue[Digits]],DrawCircuitNumberFormat[controlOrParameter[[2]],OptionValue[Digits]]]],"R"], pos]},
type == measType,
If[controlOrParameter==0,
{White, EdgeForm[Black], Rectangle[pos - {0.2,0.2}, pos + {0.2,0.2}], Black, Text["Tr", pos]},
{White, EdgeForm[Black], Rectangle[pos - {0.2,0.2}, pos + {0.2,0.2}], Black, Text["Mmt", pos]}
],
type == cnotType,
{White, EdgeForm[Black], Disk[pos , 0.2], Black, Line[{pos - {0,0.2}, pos + {0,0.2}}] , Line[{pos - {0.2,0}, pos + {0.2,0}}] },
type == czType,
{Black, Disk[pos, 0.05]},
type==diagType || type== xxType,
If[Min[target]==-pos[[2]],
(*Draw the whole diagonal or XX gate*)
If[Max[target]==Min[target],
{White, EdgeForm[Black], Rectangle[pos - {0.2,0.2}, pos + {0.2,0.2}], Black, If[type == diagType,Text["\[CapitalDelta]", pos],Text["XX", pos]]},
sortTarget=Sort[target];
(*Draw box for diagonal gate*)
box={White, EdgeForm[Black], Rectangle[{pos[[1]] -0.2,-sortTarget[[-1]]-0.2},{pos[[1]] +0.2,pos[[2]]+0.2}]};
(*Text for diag or XX gate*)
If[EvenQ[sortTarget[[-1]]-sortTarget[[1]]],