From 35350d23d015634d98fd26537c3c2e6c83ad8d0e Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Tue, 12 Nov 2019 14:16:34 -0800 Subject: [PATCH] QAOA and MAXCUT landscape plotting --- examples/qaoa.lisp | 206 ++++++++++++++++++++++++++++++++++++++++++ qvm-examples.asd | 4 +- src/wavefunction.lisp | 1 - 3 files changed, 209 insertions(+), 2 deletions(-) create mode 100644 examples/qaoa.lisp diff --git a/examples/qaoa.lisp b/examples/qaoa.lisp new file mode 100644 index 00000000..bc339a50 --- /dev/null +++ b/examples/qaoa.lisp @@ -0,0 +1,206 @@ +;;;; qaoa.lisp +;;;; +;;;; Author: Robert Smith + +(in-package #:qvm-examples) + +(defun qaoa (p cost-ham) + "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: + + 1. The Quil program (as a PARSED-PROGRAM) + + 2. A list of gamma parameters (as MEMORY-REFs) which parameterize the cost unitary. + + 3. A list of beta parameters (as MEMORY-REFs) which parameterize the driver unitary. + +COST-HAM should consist of commuting Pauli terms only. (This is *not* checked.)" + (check-type p (integer 1)) + ;; Every Pauli in COST-HAM should have the same number of qubits. + (let* ((n (cl-quil.clifford:num-qubits + (alexandria:extremum cost-ham #'> :key #'cl-quil.clifford:num-qubits))) + (quil (quil:parse-quil (format nil "DECLARE beta REAL[~D]; DECLARE gamma REAL[~D]" p p))) + ;; Parameters + (betas (loop :for i :below p :collect (quil:mref "beta" i))) + (gammas (loop :for i :below p :collect (quil:mref "gamma" i))) + ;; Hamiltonian + (driver-ham (loop :for q :below n + :collect (cl-quil.clifford:embed cl-quil.clifford:+X+ n (list q)))) + (isns (quil:with-inst () + ;; Initialize + (dotimes (q n) + (quil:inst "H" () q)) + (loop :for beta :in betas + :for gamma :in gammas + :do + (let ((beta (quil::make-delayed-expression nil nil beta)) + (gamma (quil::make-delayed-expression nil nil gamma))) + ;; Cost. All the terms are assumed to commute. + (dolist (pauli cost-ham) + (mapc #'quil:inst (cl-quil.clifford::exp-pauli pauli gamma))) + + ;; Driver + (dolist (pauli driver-ham) + (mapc #'quil:inst (cl-quil.clifford::exp-pauli pauli beta)))))))) + (setf (quil:parsed-program-executable-code quil) (coerce isns 'simple-vector)) + (values quil betas gammas))) + +(defun maxcut (graph) + "Produce a cost Hamiltonian for the graph GRAPH as a list of Pauli terms. + +GRAPH should be a list of edges, each represented as a pair (A B) of integer vertices." + (loop :with n := (1+ (loop :for (from to) :in graph + :maximize (max from to))) + :with ZZ := (cl-quil.clifford:pauli-from-string "ZZ") + :for vertex :in graph + :collect (cl-quil.clifford:embed ZZ n vertex))) + +(defun complete-graph (n) + "Build a complete graph of N vertices." + (loop :for i :from 1 :below n + :append (loop :for j :below i + :collect (list j i)))) + +(defun line-graph (n) + "Build a line graph of N vertices." + (loop :for i :below (1- n) + :collect (list i (1+ i)))) + +(defun cut-weight (graph cut) + "Given a graph GRAPH and a cut (a list of vertices), compute the cut weight." + (loop :for (from to) :in graph + :sum (logxor (ldb (byte 1 from) cut) + (ldb (byte 1 to) cut)))) + +;;; Test fixtures + +(defun qaoa-dump (graph file &key (shots 1000) (width 25)) + "Create a \"landscape plot\" for p=1 QAOA on the graph GRAPH. Write the output to the file FILE. + +The landscape will be a plot + + cut_cost : [0, pi) x [0, pi) -> R + gamma x beta |-> cut_cost(GRAPH) + +After FILE is created, one may plot the data with gnuplot. The following commands will plot (and live-update) the data: + + set view map + splot \"FILE\" u 1:2:3 w pm3d + pause 2 + reread" + (flet ((out (string) + (write-string string) + (terpri) + (finish-output))) + (out "writing program") + (multiple-value-bind (program betas gammas) (qaoa 1 (maxcut graph)) + (let* ((n (quil:qubits-needed program)) + (beta (first betas)) + (gamma (first gammas)) + (qvm (qvm:make-qvm n))) + (qvm:load-program qvm program :supersede-memory-subsystem t) + (out "compiling program") + (qvm::enable-all-qvm-optimizations) + (qvm::compile-loaded-program qvm) + (with-open-file (s file :direction ':output + :if-exists ':supersede + :if-does-not-exist ':create) + (out "run") + (let* ((start-range (qvm:flonum 0.0)) + (end-range (qvm:flonum pi)) + (step (/ (- end-range start-range) width))) + (loop :with time := (get-internal-real-time) + :with count := 0 + :for gamma-angle :from start-range :below end-range :by step :do + (loop :for beta-angle :from start-range :below end-range :by step :do + (setf (qvm::dereference-mref qvm beta) beta-angle + (qvm::dereference-mref qvm gamma) gamma-angle) + (qvm:run qvm) + (let ((samples (qvm::sample-wavefunction-multiple-times (qvm::amplitudes qvm) shots))) + (format s "~F ~F ~F~%" gamma-angle beta-angle + (loop :for bitstring :across samples + :sum (cut-weight graph bitstring) :into s + :finally (return (/ s shots))))) + (qvm::reset-quantum-state qvm) + (setf time (/ (- (get-internal-real-time) time) internal-time-units-per-second)) + (format t " ~D: ~F s (est ~F s left)~%" (incf count) time (* time (- (* width width) count))) + (setf time (get-internal-real-time))) + (out " line") + (terpri s) + (finish-output s)))))))) + +(defun produce-qvm-for-qaoa-problem (graph) + (multiple-value-bind (program betas gammas) (qaoa 1 (maxcut graph)) + (let* ((n (quil:qubits-needed program)) + (qvm (qvm:make-qvm n))) + (qvm:load-program qvm program :supersede-memory-subsystem t) + ;(qvm::enable-all-qvm-optimizations) + ;(qvm::compile-loaded-program qvm) + (values qvm (first betas) (first gammas))))) + +(defun qaoa-serial (graph width &key (shots 1000) + (start-range 0.0d0) + (end-range pi)) + (let ((results (make-array (list width width) :element-type 'qvm:flonum + :initial-element (qvm:flonum -100))) + (step (qvm:flonum (/ (- end-range start-range) width))) + (start-time (get-internal-real-time))) + (multiple-value-bind (qvm beta gamma) (produce-qvm-for-qaoa-problem graph) + (let ((time (get-internal-real-time))) + (dotimes (index (* width width)) + (multiple-value-bind (row col) (floor index width) + (let ((gamma-angle (* col step)) + (beta-angle (* row step))) + (setf (qvm::dereference-mref qvm beta) beta-angle + (qvm::dereference-mref qvm gamma) gamma-angle) + (qvm:run qvm) + (let ((samples (qvm::sample-wavefunction-multiple-times (qvm::amplitudes qvm) shots))) + (setf (aref results row col) (/ (loop :for bitstring :across samples + :sum (cut-weight graph bitstring)) + (qvm:flonum shots)))) + (qvm::reset-quantum-state qvm) + (setf time (/ (- (get-internal-real-time) time) internal-time-units-per-second)) + (setf time (get-internal-real-time))))))) + (values results (- (get-internal-real-time) start-time)))) + +(defun qaoa-parallel (graph width &key (shots 1000) + (start-range 0.0d0) + (end-range pi)) + (let ((results (make-array (list width width) :element-type 'qvm:flonum + :initial-element (qvm:flonum -100))) + (step (qvm:flonum (/ (- end-range start-range) width))) + (start-time (get-internal-real-time))) + (qvm::with-parallel-subdivisions (start end (* width width) + :num-divisions 12 + :force-parallel t) + (multiple-value-bind (qvm beta gamma) (produce-qvm-for-qaoa-problem graph) + (loop :for index :from start :below end :do + (multiple-value-bind (row col) (floor index width) + (let ((gamma-angle (* col step)) + (beta-angle (* row step))) + (setf (qvm::dereference-mref qvm beta) beta-angle + (qvm::dereference-mref qvm gamma) gamma-angle) + (qvm:run qvm) + (let ((samples (qvm::sample-wavefunction-multiple-times (qvm::amplitudes qvm) shots))) + (setf (aref results row col) (/ (loop :for bitstring :across samples + :sum (cut-weight graph bitstring)) + (qvm:flonum shots)))) + ;;(qvm::reset-quantum-state qvm) + (let ((amps (qvm::amplitudes qvm))) + (declare (optimize speed (safety 0) (debug 0) (space 0)) + (type qvm::quantum-state amps)) + (fill amps (qvm:cflonum 0)))))))) + (values results (- (get-internal-real-time) start-time)))) + +(defun compare-serial-vs-parallel () + (qvm::prepare-for-parallelization) + (sb-ext:gc :full t) + (loop :with width := 25 + :for n :from 2 :to 26 + :for time-serial := (nth-value 1 (qaoa-serial (line-graph n) width)) + :for time-parallel := (nth-value 1 (qaoa-parallel (line-graph n) width)) + :do (format t "~2Dq: ~10,3F ~20T| ~10,3F ~10T[~Dx]~%" n + (/ time-serial internal-time-units-per-second) + (/ time-parallel internal-time-units-per-second) + (round time-serial time-parallel)) + (sb-ext:gc :full t) + :collect (cons time-serial time-parallel))) diff --git a/qvm-examples.asd b/qvm-examples.asd index 430666e4..59efe6db 100644 --- a/qvm-examples.asd +++ b/qvm-examples.asd @@ -7,6 +7,7 @@ :author "Robert Smith , Juan M. Bello-Rivas " :license "Apache License 2.0 (See LICENSE.txt)" :depends-on ( + (:version #:cl-quil "1.13.1") ;; Nelder-Mead #:cl-grnm ;; Quantum Virtual Machine @@ -18,4 +19,5 @@ :serial t :components ((:file "package") (:file "qft") - (:file "vqe"))) + (:file "vqe") + (:file "qaoa"))) diff --git a/src/wavefunction.lisp b/src/wavefunction.lisp index f2882a3e..2c423817 100644 --- a/src/wavefunction.lisp +++ b/src/wavefunction.lisp @@ -312,7 +312,6 @@ up to amplitude ordering." (inline apply-operator)) (let* ((n (length wavefunction)) (k (nat-tuple-cardinality qubits)) - ;; KERNEL is playing double duty as a kernel and a boolean. (use-kernel? (<= n (expt 2 (qubit-limit-for-using-serial-kernels)))) (kernel (find-serial-kernel k))) (cond