Skip to content

Commit 35350d2

Browse files
Robert Smithstylewarning
Robert Smith
authored and
stylewarning
committed
QAOA and MAXCUT landscape plotting
1 parent c7d1b15 commit 35350d2

File tree

3 files changed

+209
-2
lines changed

3 files changed

+209
-2
lines changed

‎examples/qaoa.lisp

Lines changed: 206 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,206 @@
1+
;;;; qaoa.lisp
2+
;;;;
3+
;;;; Author: Robert Smith
4+
5+
(in-package #:qvm-examples)
6+
7+
(defun qaoa (p cost-ham)
8+
"Produce a QAOA program for the list of Pauli (CL-QUIL.CLIFFORD:PAULI) terms COST-HAM representing a cost Hamiltonian to optimize. Return three values:
9+
10+
1. The Quil program (as a PARSED-PROGRAM)
11+
12+
2. A list of gamma parameters (as MEMORY-REFs) which parameterize the cost unitary.
13+
14+
3. A list of beta parameters (as MEMORY-REFs) which parameterize the driver unitary.
15+
16+
COST-HAM should consist of commuting Pauli terms only. (This is *not* checked.)"
17+
(check-type p (integer 1))
18+
;; Every Pauli in COST-HAM should have the same number of qubits.
19+
(let* ((n (cl-quil.clifford:num-qubits
20+
(alexandria:extremum cost-ham #'> :key #'cl-quil.clifford:num-qubits)))
21+
(quil (quil:parse-quil (format nil "DECLARE beta REAL[~D]; DECLARE gamma REAL[~D]" p p)))
22+
;; Parameters
23+
(betas (loop :for i :below p :collect (quil:mref "beta" i)))
24+
(gammas (loop :for i :below p :collect (quil:mref "gamma" i)))
25+
;; Hamiltonian
26+
(driver-ham (loop :for q :below n
27+
:collect (cl-quil.clifford:embed cl-quil.clifford:+X+ n (list q))))
28+
(isns (quil:with-inst ()
29+
;; Initialize
30+
(dotimes (q n)
31+
(quil:inst "H" () q))
32+
(loop :for beta :in betas
33+
:for gamma :in gammas
34+
:do
35+
(let ((beta (quil::make-delayed-expression nil nil beta))
36+
(gamma (quil::make-delayed-expression nil nil gamma)))
37+
;; Cost. All the terms are assumed to commute.
38+
(dolist (pauli cost-ham)
39+
(mapc #'quil:inst (cl-quil.clifford::exp-pauli pauli gamma)))
40+
41+
;; Driver
42+
(dolist (pauli driver-ham)
43+
(mapc #'quil:inst (cl-quil.clifford::exp-pauli pauli beta))))))))
44+
(setf (quil:parsed-program-executable-code quil) (coerce isns 'simple-vector))
45+
(values quil betas gammas)))
46+
47+
(defun maxcut (graph)
48+
"Produce a cost Hamiltonian for the graph GRAPH as a list of Pauli terms.
49+
50+
GRAPH should be a list of edges, each represented as a pair (A B) of integer vertices."
51+
(loop :with n := (1+ (loop :for (from to) :in graph
52+
:maximize (max from to)))
53+
:with ZZ := (cl-quil.clifford:pauli-from-string "ZZ")
54+
:for vertex :in graph
55+
:collect (cl-quil.clifford:embed ZZ n vertex)))
56+
57+
(defun complete-graph (n)
58+
"Build a complete graph of N vertices."
59+
(loop :for i :from 1 :below n
60+
:append (loop :for j :below i
61+
:collect (list j i))))
62+
63+
(defun line-graph (n)
64+
"Build a line graph of N vertices."
65+
(loop :for i :below (1- n)
66+
:collect (list i (1+ i))))
67+
68+
(defun cut-weight (graph cut)
69+
"Given a graph GRAPH and a cut (a list of vertices), compute the cut weight."
70+
(loop :for (from to) :in graph
71+
:sum (logxor (ldb (byte 1 from) cut)
72+
(ldb (byte 1 to) cut))))
73+
74+
;;; Test fixtures
75+
76+
(defun qaoa-dump (graph file &key (shots 1000) (width 25))
77+
"Create a \"landscape plot\" for p=1 QAOA on the graph GRAPH. Write the output to the file FILE.
78+
79+
The landscape will be a plot
80+
81+
cut_cost : [0, pi) x [0, pi) -> R
82+
gamma x beta |-> cut_cost(GRAPH)
83+
84+
After FILE is created, one may plot the data with gnuplot. The following commands will plot (and live-update) the data:
85+
86+
set view map
87+
splot \"FILE\" u 1:2:3 w pm3d
88+
pause 2
89+
reread"
90+
(flet ((out (string)
91+
(write-string string)
92+
(terpri)
93+
(finish-output)))
94+
(out "writing program")
95+
(multiple-value-bind (program betas gammas) (qaoa 1 (maxcut graph))
96+
(let* ((n (quil:qubits-needed program))
97+
(beta (first betas))
98+
(gamma (first gammas))
99+
(qvm (qvm:make-qvm n)))
100+
(qvm:load-program qvm program :supersede-memory-subsystem t)
101+
(out "compiling program")
102+
(qvm::enable-all-qvm-optimizations)
103+
(qvm::compile-loaded-program qvm)
104+
(with-open-file (s file :direction ':output
105+
:if-exists ':supersede
106+
:if-does-not-exist ':create)
107+
(out "run")
108+
(let* ((start-range (qvm:flonum 0.0))
109+
(end-range (qvm:flonum pi))
110+
(step (/ (- end-range start-range) width)))
111+
(loop :with time := (get-internal-real-time)
112+
:with count := 0
113+
:for gamma-angle :from start-range :below end-range :by step :do
114+
(loop :for beta-angle :from start-range :below end-range :by step :do
115+
(setf (qvm::dereference-mref qvm beta) beta-angle
116+
(qvm::dereference-mref qvm gamma) gamma-angle)
117+
(qvm:run qvm)
118+
(let ((samples (qvm::sample-wavefunction-multiple-times (qvm::amplitudes qvm) shots)))
119+
(format s "~F ~F ~F~%" gamma-angle beta-angle
120+
(loop :for bitstring :across samples
121+
:sum (cut-weight graph bitstring) :into s
122+
:finally (return (/ s shots)))))
123+
(qvm::reset-quantum-state qvm)
124+
(setf time (/ (- (get-internal-real-time) time) internal-time-units-per-second))
125+
(format t " ~D: ~F s (est ~F s left)~%" (incf count) time (* time (- (* width width) count)))
126+
(setf time (get-internal-real-time)))
127+
(out " line")
128+
(terpri s)
129+
(finish-output s))))))))
130+
131+
(defun produce-qvm-for-qaoa-problem (graph)
132+
(multiple-value-bind (program betas gammas) (qaoa 1 (maxcut graph))
133+
(let* ((n (quil:qubits-needed program))
134+
(qvm (qvm:make-qvm n)))
135+
(qvm:load-program qvm program :supersede-memory-subsystem t)
136+
;(qvm::enable-all-qvm-optimizations)
137+
;(qvm::compile-loaded-program qvm)
138+
(values qvm (first betas) (first gammas)))))
139+
140+
(defun qaoa-serial (graph width &key (shots 1000)
141+
(start-range 0.0d0)
142+
(end-range pi))
143+
(let ((results (make-array (list width width) :element-type 'qvm:flonum
144+
:initial-element (qvm:flonum -100)))
145+
(step (qvm:flonum (/ (- end-range start-range) width)))
146+
(start-time (get-internal-real-time)))
147+
(multiple-value-bind (qvm beta gamma) (produce-qvm-for-qaoa-problem graph)
148+
(let ((time (get-internal-real-time)))
149+
(dotimes (index (* width width))
150+
(multiple-value-bind (row col) (floor index width)
151+
(let ((gamma-angle (* col step))
152+
(beta-angle (* row step)))
153+
(setf (qvm::dereference-mref qvm beta) beta-angle
154+
(qvm::dereference-mref qvm gamma) gamma-angle)
155+
(qvm:run qvm)
156+
(let ((samples (qvm::sample-wavefunction-multiple-times (qvm::amplitudes qvm) shots)))
157+
(setf (aref results row col) (/ (loop :for bitstring :across samples
158+
:sum (cut-weight graph bitstring))
159+
(qvm:flonum shots))))
160+
(qvm::reset-quantum-state qvm)
161+
(setf time (/ (- (get-internal-real-time) time) internal-time-units-per-second))
162+
(setf time (get-internal-real-time)))))))
163+
(values results (- (get-internal-real-time) start-time))))
164+
165+
(defun qaoa-parallel (graph width &key (shots 1000)
166+
(start-range 0.0d0)
167+
(end-range pi))
168+
(let ((results (make-array (list width width) :element-type 'qvm:flonum
169+
:initial-element (qvm:flonum -100)))
170+
(step (qvm:flonum (/ (- end-range start-range) width)))
171+
(start-time (get-internal-real-time)))
172+
(qvm::with-parallel-subdivisions (start end (* width width)
173+
:num-divisions 12
174+
:force-parallel t)
175+
(multiple-value-bind (qvm beta gamma) (produce-qvm-for-qaoa-problem graph)
176+
(loop :for index :from start :below end :do
177+
(multiple-value-bind (row col) (floor index width)
178+
(let ((gamma-angle (* col step))
179+
(beta-angle (* row step)))
180+
(setf (qvm::dereference-mref qvm beta) beta-angle
181+
(qvm::dereference-mref qvm gamma) gamma-angle)
182+
(qvm:run qvm)
183+
(let ((samples (qvm::sample-wavefunction-multiple-times (qvm::amplitudes qvm) shots)))
184+
(setf (aref results row col) (/ (loop :for bitstring :across samples
185+
:sum (cut-weight graph bitstring))
186+
(qvm:flonum shots))))
187+
;;(qvm::reset-quantum-state qvm)
188+
(let ((amps (qvm::amplitudes qvm)))
189+
(declare (optimize speed (safety 0) (debug 0) (space 0))
190+
(type qvm::quantum-state amps))
191+
(fill amps (qvm:cflonum 0))))))))
192+
(values results (- (get-internal-real-time) start-time))))
193+
194+
(defun compare-serial-vs-parallel ()
195+
(qvm::prepare-for-parallelization)
196+
(sb-ext:gc :full t)
197+
(loop :with width := 25
198+
:for n :from 2 :to 26
199+
:for time-serial := (nth-value 1 (qaoa-serial (line-graph n) width))
200+
:for time-parallel := (nth-value 1 (qaoa-parallel (line-graph n) width))
201+
:do (format t "~2Dq: ~10,3F ~20T| ~10,3F ~10T[~Dx]~%" n
202+
(/ time-serial internal-time-units-per-second)
203+
(/ time-parallel internal-time-units-per-second)
204+
(round time-serial time-parallel))
205+
(sb-ext:gc :full t)
206+
:collect (cons time-serial time-parallel)))

‎qvm-examples.asd

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
:author "Robert Smith <[email protected]>, Juan M. Bello-Rivas <[email protected]>"
88
:license "Apache License 2.0 (See LICENSE.txt)"
99
:depends-on (
10+
(:version #:cl-quil "1.13.1")
1011
;; Nelder-Mead
1112
#:cl-grnm
1213
;; Quantum Virtual Machine
@@ -18,4 +19,5 @@
1819
:serial t
1920
:components ((:file "package")
2021
(:file "qft")
21-
(:file "vqe")))
22+
(:file "vqe")
23+
(:file "qaoa")))

‎src/wavefunction.lisp

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -312,7 +312,6 @@ up to amplitude ordering."
312312
(inline apply-operator))
313313
(let* ((n (length wavefunction))
314314
(k (nat-tuple-cardinality qubits))
315-
;; KERNEL is playing double duty as a kernel and a boolean.
316315
(use-kernel? (<= n (expt 2 (qubit-limit-for-using-serial-kernels))))
317316
(kernel (find-serial-kernel k)))
318317
(cond

0 commit comments

Comments
 (0)