From d8b00af22be39a14a14939de10cd79fdbc9dc7f5 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Mon, 24 Sep 2018 18:52:48 -0400 Subject: [PATCH 01/24] H2 simulation code --- .../ground_state_estimation.h2.scaffold | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold index 4e4da1161..355ec37d1 100644 --- a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold +++ b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold @@ -1,7 +1,7 @@ #include #define M 4 //number of bottom register qubits -#define TROTTER 16 //number of bottom register qubits +#define TROTTER 64 enum spin {alpha, beta}; const enum spin sigma[M] = { @@ -162,7 +162,7 @@ int main () // PrepZ ( bottomregister[i], 1 ); // } PrepZ ( bottomregister[3], 1 ); - PrepZ ( bottomregister[2], 1 ); + PrepZ ( bottomregister[2], 0 ); PrepZ ( bottomregister[1], 1 ); PrepZ ( bottomregister[0], 0 ); @@ -217,6 +217,7 @@ int main () } + // TODO: something in here causing Trotter to not converge for (int t=3; t>=0; t--) { for (int c=t-1; c>=0; c--) { CNOT(bottomregister[c], bottomregister[t]); @@ -230,6 +231,7 @@ int main () } // two electron operators: excitation-exitation operators + // TODO: something in here causing Trotter to not converge DoubleExcitationOperator ( topregister[0], bottomregister, From 464bd7a3778c8f188872f9159473f4a5c754a72f Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Thu, 4 Oct 2018 15:15:14 -0400 Subject: [PATCH 02/24] Script for iterative phase estimation --- iterative_phase_estimation_qchem.py | 39 +++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100755 iterative_phase_estimation_qchem.py diff --git a/iterative_phase_estimation_qchem.py b/iterative_phase_estimation_qchem.py new file mode 100755 index 000000000..b0651ee26 --- /dev/null +++ b/iterative_phase_estimation_qchem.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python +import sys +import re +from subprocess import check_output + +upper_sz = 9 +estimate = 0 + +for iteration in range(0, upper_sz): + + with open(sys.argv[1]) as infile, open(sys.argv[1]+'.ipe', 'w') as outfile: + + for line in infile: + if "#define _upper_sz" in line: + line = re.sub('\d+', str(upper_sz), line) + if "#define _iteration" in line: + line = re.sub('\d+', str(iteration), line) + if "#define _estimate" in line: + line = re.sub('\d+', str(estimate), line) + outfile.write(line) + + sim_out = check_output(['./scripts/simulation/simulate.sh', './Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold.ipe']) + + state_lines = False + for line in sim_out.split('\n'): + if line.strip() == "--------------[quantum state]--------------": + print("--------------[quantum state]--------------") + state_lines = True + elif line.strip() == "-------------------------------------------": + print("-------------------------------------------") + state_lines = False + elif state_lines: + print(line) + strings = re.findall(r"[-+]?\d*\.\d+|\d+", line) + if (float(strings[0]) > 0.5): + meas_val = int(strings[4][1]) + print(meas_val) + estimate += meas_val * pow(2,iteration) + print(estimate) From 8aa9fa7ab424984fc74fe702953115791d613cd8 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Thu, 4 Oct 2018 15:15:38 -0400 Subject: [PATCH 03/24] Ground state estimation --- .../ground_state_estimation.h2.scaffold | 60 ++++--- .../ground_state_estimation.m10.scaffold | 164 +++++++++--------- 2 files changed, 122 insertions(+), 102 deletions(-) diff --git a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold index 355ec37d1..b957c2094 100644 --- a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold +++ b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold @@ -1,7 +1,11 @@ #include #define M 4 //number of bottom register qubits -#define TROTTER 64 +#define TROTTER 4 + +#define _upper_sz 4 +#define _iteration 0 +#define _estimate 0 enum spin {alpha, beta}; const enum spin sigma[M] = { @@ -31,9 +35,12 @@ const enum spin sigma[M] = { // h_double[2][3][3][2]=h_double[3][2][2][3]=0.697398010; // // h_double[0][2][0][2]=h_double[1][3][1][3]=0.181287518; -// h_double[2][0][2][0]=h_double[3][1][3][1]=0.000000; +// h_double[2][0][2][0]=h_double[3][1][3][1]=0.181287518; // h_double[0][3][1][2]=h_double[0][1][3][2]=0.181287518; +// according to 1.4768229 +// h_double[2][1][3][0]=h_double[2][3][1][0]=0.181287518; + const double h_single[M][M] = { {-1.252477495, 0.0, 0.0, 0.0}, {0.0, -1.252477495, 0.0, 0.0}, @@ -56,14 +63,14 @@ const double h_double[M][M][M][M] = {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.181287518},{0.0,0.0,0.0,0.0},{0.0,0.663472101,0.0,0.0}} }, { - {{0.0,0.0,0.663472101,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, + {{0.0,0.0,0.663472101,0.0},{0.0,0.0,0.0,0.0},{0.181287518,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.663472101,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.697398010,0.0}} }, { {{0.0,0.0,0.0,0.663472101},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, - {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.663472101},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, + {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.663472101},{0.0,0.0,0.0,0.0},{0.0,0.181287518,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.697398010},{0.0,0.0,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}} } @@ -161,28 +168,38 @@ int main () // for ( int i=0; i=0; t--) { for (int c=t-1; c>=0; c--) { CNOT(bottomregister[c], bottomregister[t]); double eta = 0.25 * h_double[c][t][t][c]; if (sigma[c]==sigma[t]) eta -= 0.25*h_double[c][t][c][t]; - ControlledRotation(topregister[0], bottomregister[t], 2*eta*time_step); + ControlledRotation(topregister[0], bottomregister[t], eta*time_step); CNOT(bottomregister[c], bottomregister[t]); } @@ -236,10 +251,10 @@ int main () topregister[0], bottomregister, 0, 1, 2, 3, - -time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 4.0, - -time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 4.0, - time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 4.0, - time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 4.0 + -time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 8.0, + -time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 8.0, + time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 8.0, + time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 8.0 ); } } @@ -250,6 +265,11 @@ int main () // Rz(topregister[0],-M_PI/8);// // Rz(topregister[0],-M_PI/4);// // Rz(topregister[0],-M_PI/2);// + for (int c=0; c<_iteration; c++) { + if ((_estimate>>c)&1) { + Rz ( topregister[0], -M_PI/pow(2,_iteration-c) ); + } + } H(topregister[0]); // MeasZ(bottomregister[0]); // MeasZ(bottomregister[1]); diff --git a/Algorithms/Ground_State_Estimation/ground_state_estimation.m10.scaffold b/Algorithms/Ground_State_Estimation/ground_state_estimation.m10.scaffold index eacb0cf1b..7144a2226 100644 --- a/Algorithms/Ground_State_Estimation/ground_state_estimation.m10.scaffold +++ b/Algorithms/Ground_State_Estimation/ground_state_estimation.m10.scaffold @@ -4,14 +4,14 @@ * * Oana Theogarajan, UCSB * -* +* * Implements the QCS GFI for the Ground State Estimation Algorithm * * December 2011 * *****************************************************************************************/ // The ground state is stored in M qubits coresponding to spin-orbital states (the quantum register qbit bottomregister[M]) -// There are M/2 spatial orbitals, each with spin up and down +// There are M/2 spatial orbitals, each with spin up and down // The qubits are indexed by p=0 to M-1 corresponding to i,j i=0 to M/2-1, j=+,- (i=spatial quantum number, j=spin quantum number) // So qubits 0,1,2,3,4,5,.... correspond to states with quantum numbers 0+,0-,1+,1-,2+,2-,.... @@ -34,10 +34,10 @@ #define Emax 2.2 //dummy value for now - we don't know it, not in the GFIs #define DeltaE 0.0184 */ -#define tau 341.47746 //ajavadia +#define tau 341.47746 //ajavadia // file containing one electron integrals as a one dimensional array - one element per line -//const char *integralsFile1="hmol1.txt"; +//const char *integralsFile1="hmol1.txt"; // file containing two electron integrals as a one dimensional array - one element per line //const char *integralsFile2="hmols2.txt"; @@ -45,7 +45,7 @@ #define SZ2 20 // hmol1.txt -double const hmol1[SZ1][SZ1] = +double const hmol1[SZ1][SZ1] = { {-350.29076436161165,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}, {-1.3305426082250455e-11,-350.29076949405697,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}, @@ -9053,7 +9053,7 @@ module Swap(qbit q1,qbit q2) // Implements a controlled rotation // inputs are the target qubit, control qubit and the rotation angle -module ControlledRotation(qbit target, qbit control, double theta) +module ControlledRotation(qbit target, qbit control, double theta) { Rz(target, theta / 2.0); CNOT(control, target); @@ -9062,9 +9062,9 @@ module ControlledRotation(qbit target, qbit control, double theta) } // Implement the Y gate defined in the GFIs Y=Rx(-PI/2), we'll name it Y_Rx -module Y_Rx(qbit target) +module Y_Rx(qbit target) { -// G(target,-PI/4); //indistinguishable global phase + // G(target,-PI/4); //indistinguishable global phase Z(target); H(target); S(target); @@ -9084,7 +9084,7 @@ module Y_Rx_Dagger(qbit target) // Controlled Phase gate // controlled-[e^(i*theta) 0;0 e^(i*theta)] is equivalent to the phase shift gate [1 0; 0 e^(i*theta)] on the control qubit // can be written as a product of G(theta/2) and Rz(theta); G is physically unobservable so it can be neglected -module ControlledPhase(qbit target, qbit control, double theta) +module ControlledPhase(qbit target, qbit control, double theta) { Rz(control, theta); } @@ -9102,10 +9102,10 @@ module ControlledRotationPiMinus(qbit target, qbit control,int j) // Cascade CNOT between qubits p and q (with p>q) and controlled z rotation of angle theta on qubit q // The control qubit is passed as a parameter and it will be one of the top register qubits -module CascadeControlledRotation(qbit reg[M], int p, int q,double theta,qbit control) -{ +module CascadeControlledRotation(qbit reg[M], int p, int q,double theta,qbit control) +{ int i; - + for (i=p;i>q;i--) CNOT(reg[i],reg[i-1]); ControlledRotation(reg[q],control,theta); @@ -9117,9 +9117,9 @@ module CascadeControlledRotation(qbit reg[M], int p, int q,double theta,qbit con // p>q>r>s // Cascade CNOT between p and q, CNOT between q and r and again cascade CNOT between r and s; controlled rotation of angle theta on qubit s // The control qubit is passed as a parameter and it will be one of the top register qubits -module CascadeControlledRotation_pqrs(qbit reg[M], int p, int q,int r, int s, double theta,qbit control) +module CascadeControlledRotation_pqrs(qbit reg[M], int p, int q,int r, int s, double theta,qbit control) { - int i; + int i; for (i=p;i>q;i--) CNOT(reg[i],reg[i-1]); CNOT(reg[p],reg[r]); @@ -9137,9 +9137,9 @@ module CascadeControlledRotation_pqrs(qbit reg[M], int p, int q,int r, int s, do // p>r>q // Cascade CNOT between p and r, CNOT between r and q and controlled rotation of angle theta on qubit q // The control qubit is passed as a parameter and it will be one of the top register qubits -module CascadeControlledRotation_prq(qbit reg[M], int p, int q,int r, double theta,qbit control) +module CascadeControlledRotation_prq(qbit reg[M], int p, int q,int r, double theta,qbit control) { - int i; + int i; for (i=p;i>r;i--) CNOT(reg[i],reg[i-1]); CNOT(reg[r],reg[q]); @@ -9151,10 +9151,10 @@ module CascadeControlledRotation_prq(qbit reg[M], int p, int q,int r, double the // q>p>r // CNOT between q and p, cascade CNOT between p and r and controlled rotation of angle theta on qubit r -// The control qubit is passed as a parameter and it will be one of the top register qubits -module CascadeControlledRotation_qpr(qbit reg[M], int p, int q,int r, double theta,qbit control) -{ - int i; +// The control qubit is passed as a parameter and it will be one of the top register qubits +module CascadeControlledRotation_qpr(qbit reg[M], int p, int q,int r, double theta,qbit control) +{ + int i; CNOT(reg[q],reg[p]); for (i=p;i>r;i--) CNOT(reg[i],reg[i-1]); @@ -9165,10 +9165,10 @@ module CascadeControlledRotation_qpr(qbit reg[M], int p, int q,int r, double the } -// Controlled version of the number operator "circuit primitive" on qubit p of the bottom register -// The control qubit is passed as a parameter and it will be one of the top register qubits -module NumberOperator(qbit target,double theta, qbit control) -{ +// Controlled version of the number operator "circuit primitive" on qubit p of the bottom register +// The control qubit is passed as a parameter and it will be one of the top register qubits +module NumberOperator(qbit target,double theta, qbit control) +{ ControlledPhase(target,control,theta/2.0); ControlledRotation(target,control,theta); } @@ -9176,7 +9176,7 @@ module NumberOperator(qbit target,double theta, qbit control) //p>q // Controlled version of the Excitation Operator Template circuit primitive -// The control qubit is passed as a parameter and it will be one of the top register qubits +// The control qubit is passed as a parameter and it will be one of the top register qubits module ExcitationOperatorTemplate(qbit reg[M], int p, int q,double theta, qbit control) { H(reg[q]); @@ -9188,7 +9188,7 @@ module ExcitationOperatorTemplate(qbit reg[M], int p, int q,double theta, qbit c Y_Rx(reg[p]); CascadeControlledRotation(reg,p,q,theta,control); Y_Rx_Dagger(reg[q]); - Y_Rx_Dagger(reg[p]); + Y_Rx_Dagger(reg[p]); } @@ -9199,9 +9199,9 @@ module ExcitationOperator(qbit reg[M], int p, int q, double theta,qbit control) } //controlled version of the coulomb Exchange Operator circuit primitive on qubits p and q -// Assumes p>q - will be called only for qubits with p>q -module CoulombExchgOperator(qbit reg[M],int p, int q, double theta, qbit control) -{ +// Assumes p>q - will be called only for qubits with p>q +module CoulombExchgOperator(qbit reg[M],int p, int q, double theta, qbit control) +{ ControlledPhase(reg[p],control,theta/2); ControlledRotation(reg[q],control,theta); ControlledRotation(reg[p],control,theta); @@ -9211,7 +9211,7 @@ module CoulombExchgOperator(qbit reg[M],int p, int q, double theta, qbit control } //Implements controlled version of Number Excitation Operator circuit primitive for the case p>q>r -module NrExcitationOperator_PQR(qbit reg[M],int p, int q,int r, double theta, qbit control) +module NrExcitationOperator_PQR(qbit reg[M],int p, int q,int r, double theta, qbit control) { ExcitationOperatorTemplate(reg,p,r,theta,control); @@ -9225,12 +9225,12 @@ module NrExcitationOperator_PQR(qbit reg[M],int p, int q,int r, double theta, q Y_Rx(reg[p]); CascadeControlledRotation_pqrs(reg,p,q+1,q-1,r,-theta,control); Y_Rx_Dagger(reg[r]); - Y_Rx_Dagger(reg[p]); + Y_Rx_Dagger(reg[p]); } //Implements controlled version of Number Excitation Operator circuit primitive for the case p>r>q -module NrExcitationOperator_PRQ(qbit reg[M],int p, int q,int r, double theta, qbit control) +module NrExcitationOperator_PRQ(qbit reg[M],int p, int q,int r, double theta, qbit control) { ExcitationOperatorTemplate(reg,p,r,theta,control); @@ -9244,12 +9244,12 @@ module NrExcitationOperator_PRQ(qbit reg[M],int p, int q,int r, double theta, q Y_Rx(reg[p]); CascadeControlledRotation_prq(reg,p,q,r,-theta,control); Y_Rx_Dagger(reg[r]); - Y_Rx_Dagger(reg[p]); + Y_Rx_Dagger(reg[p]); } //Implements controlled version of Number Excitation Operator circuit primitive for the case q>p>r -module NrExcitationOperator_QPR(qbit reg[M],int p, int q,int r, double theta, qbit control) +module NrExcitationOperator_QPR(qbit reg[M],int p, int q,int r, double theta, qbit control) { ExcitationOperatorTemplate(reg,p,r,theta,control); @@ -9263,13 +9263,13 @@ module NrExcitationOperator_QPR(qbit reg[M],int p, int q,int r, double theta, q Y_Rx(reg[p]); CascadeControlledRotation_qpr(reg,p,q,r,-theta,control); Y_Rx_Dagger(reg[r]); - Y_Rx_Dagger(reg[p]); + Y_Rx_Dagger(reg[p]); } -// controlled version of Number Excitation Operator circuit primitive +// controlled version of Number Excitation Operator circuit primitive // Assumes p>r and handles separately the 3 different circuits for the 3 cases p>q>r, p>r>q, q>p>r -module NrExcitationOperator(qbit reg[M],int p, int q,int r, double theta, qbit control) +module NrExcitationOperator(qbit reg[M],int p, int q,int r, double theta, qbit control) { if (qq>r>s -module DoubleExcitationOperator(qbit reg[M],int p, int q,int r, int s, double theta1, double theta2, double theta3, double theta4, qbit control) +module DoubleExcitationOperator(qbit reg[M],int p, int q,int r, int s, double theta1, double theta2, double theta3, double theta4, qbit control) { //1 H(reg[s]); @@ -9374,9 +9374,9 @@ module DoubleExcitationOperator(qbit reg[M],int p, int q,int r, int s, double th Y_Rx_Dagger(reg[p]); } -// The controlled version of PauliH1 module form the GFI pseudocode +// The controlled version of PauliH1 module form the GFI pseudocode // Includes the Hamiltonian terms involving one electron integrals -module PauliH1(qbit bottomregister[M],qbit control) +module PauliH1(qbit bottomregister[M],qbit control) { int p,q,pp,qq; double theta; @@ -9404,7 +9404,7 @@ module PauliH1(qbit bottomregister[M],qbit control) //get the spatial indices for the integrals pp=p/2; qq=q/2; - + //theta=hpq[p,q]*tau/100; //read h[p,q] from file using the corresponding spatial indices pp,qq // theta=readline(integralsFile1,Index2to1Dimension(pp,qq)); @@ -9415,13 +9415,13 @@ module PauliH1(qbit bottomregister[M],qbit control) theta=theta*((p-q+1)%2); if(theta!=0) - ExcitationOperator(bottomregister,p,q,theta,control); + ExcitationOperator(bottomregister,p,q,theta,control); } - } + } } -// The controlled version of PauliH2 module form the GFI pseudocode +// The controlled version of PauliH2 module form the GFI pseudocode // Includes the Hamiltonian terms involving two electron integrals module PauliH2(qbit bottomregister[M],qbit control) { @@ -9436,22 +9436,22 @@ module PauliH2(qbit bottomregister[M],qbit control) pp=p/2; qq=q/2; // theta=(1/2)*(h[p,q,p,q]-h[p,q,q,p])*tau/100; - + //read h[p,q,p,q] from file using the corresponding spatial indices pp,qq,rr // tmp1=readline(integralsFile2,Index4to1Dimension(pp,qq,pp,qq)); tmp1=hmol2[pp][qq][pp][qq]; - + //read h[p,q,q,p] from file using the corresponding spatial indices pp,qq,rr // tmp2=readline(integralsFile2,Index4to1Dimension(pp,qq,qq,pp)); tmp2=hmol2[pp][qq][qq][pp]; - + //only apply module if theta not zero, make theta=0 if spins are not equal for il and jk in h_ijkl //for h[p,q,q,p] the spins are always equal since il=pp and jk=qq, so it is always non-zero - theta1=(1/2.0)*(tmp1*((p-q+1)%2)-tmp2)*tau/100; + theta1=(1/2.0)*(tmp1*((p-q+1)%2)-tmp2)*tau/100; if(theta1!=0) { - CoulombExchgOperator(bottomregister,p,q,theta1,control); - } + CoulombExchgOperator(bottomregister,p,q,theta1,control); + } } } for (r=0;rp>r, p>q>r, p>r>q) @@ -9488,11 +9488,11 @@ module PauliH2(qbit bottomregister[M],qbit control) } } } - } + } } for (s=0;s Date: Wed, 10 Oct 2018 10:23:02 -0400 Subject: [PATCH 04/24] Simplified version of Grover's for simulation. --- .../Square_Root/square_root.n4.scaffold | 231 ++++++++++++++++++ 1 file changed, 231 insertions(+) create mode 100644 Algorithms/Square_Root/square_root.n4.scaffold diff --git a/Algorithms/Square_Root/square_root.n4.scaffold b/Algorithms/Square_Root/square_root.n4.scaffold new file mode 100644 index 000000000..97346bec3 --- /dev/null +++ b/Algorithms/Square_Root/square_root.n4.scaffold @@ -0,0 +1,231 @@ +/***************************************************************************** + * + * Grover's Search Algorithm + * + * Implements Grover's Search algorithm in the Scaffold programming language + * + * February 2013 + * + *****************************************************************************/ + +#include + +#define n 4 // problem size (log of database size) + +#define pi 3.141592653589793238462643383279502884197 + +// Module prototypes +void diffuse(qbit q[n], qbit x[n-1]); +void Sqr(qbit a[n], qbit b[n]); +void EQxMark(qbit b[n], qbit x[n-1]); + +module cRz ( qbit ctrl, qbit target, const double angle ) { + /* cRz identity matrix: + [ [ 1 0 0 0 ] + [ 0 1 0 0 ] + [ 0 0 e^(-i*angle/2) 0 ] + [ 0 0 0 e^(i*angle/2) ] ] + */ + + S(ctrl); + CNOT(ctrl, target); + Rz(target,-angle/2); + CNOT(ctrl, target); + Rz(target,angle/2); +} + +/*************************** + Diffusion Module +***************************/ +void diffuse(qbit q[n], qbit x[n-1]) { + int j; + // local qbit x[n-1]; // scratch register + // No local registers yet + // qbit x[n-1]; // scratch register + + // Hadamard applied to q + // No forall yet + // forall(j=0; j 0; j--) + Toffoli(x[j-1], q[j+1], x[j]); + Toffoli(q[1], q[0], x[0]); + + // Restore q + for(j = 0; j < n; j++) + X(q[j]); + + // Complete the diffusion + // No forall yet + // forall(j=0; j 0; j--) + Toffoli(x[j-1], b[j+1], x[j]); + Toffoli(b[1], b[0], x[0]); + + // Restore b + for(j = 0; j < n; j++) + if(j!=1) X(b[j]); +} + +/*********************************************** + Squaring a(x) = a0 + a1*x + ... + a_(n-1)*x^(n-1) + over the ring GF(2)[x] mod (x^n + x + 1). + + Result placed in the n-qubit register b +************************************************/ +void Sqr(qbit a[n], qbit b[n]) { + + int i; + int k; + + // Using forall indicates the CNOT's are independent + // So these instructions can be done in parallel + // No forall yet + // forall(i=0; i<=(n-1)/2; i++) { + for(i=0; i<=(n-1)/2; i++) { + k = 2*i; + CNOT(a[i], b[k]); + } + + // Would it make a difference to split this into two loops? + // Maybe since deleaving would be two forall loops! + // Or perhaps it is okay to just replace the for with a forall. + for(i=(n+1)/2; i Date: Wed, 10 Oct 2018 15:20:57 -0400 Subject: [PATCH 05/24] Synching Shor's implementation --- Algorithms/Shors/Shors_15/shors.n3.scaffold | 73 +++++++++++++-------- 1 file changed, 47 insertions(+), 26 deletions(-) diff --git a/Algorithms/Shors/Shors_15/shors.n3.scaffold b/Algorithms/Shors/Shors_15/shors.n3.scaffold index 1c6868e33..22445f48a 100644 --- a/Algorithms/Shors/Shors_15/shors.n3.scaffold +++ b/Algorithms/Shors/Shors_15/shors.n3.scaffold @@ -1,4 +1,4 @@ -//--- Scaffold Code for Shor’s Algorithm ---// +// --- Scaffold Code for Shor’s Algorithm ---// //Input x -> H -> f(x) -> iQFT //f(x) = a^x mod N @@ -9,9 +9,13 @@ #include "../cMODMUL/cMODMUL.scaffold" #include "../../../Library/cSWAP/cSWAP.scaffold" -#define _width 5 // one extra than number of bits in N -#define _N 15 // number to be factorized -// #define _a 7 // randomly chosen s.t. gcd(a,N) = 1 + +#define _a_val 7 // randomly chosen s.t. gcd(a,N) = 1 +#define _a_inv 13 // randomly chosen s.t. gcd(a,N) = 1 + +#define _upper_sz 3 // total iterations of phase estimation +#define _iteration 2 // current iteration index +#define _estimate 0 // tally of less significant bits module cUa ( qbit ctrl, @@ -32,50 +36,67 @@ module cUa ( // b[0] is least significant bit int main ( ) { - const unsigned int upper_sz = 5; - qbit upper[upper_sz]; - for (int i=0; i>i)&1 ); } + // prepare ancilla b to zero const unsigned int b_val = 0; - qbit b[_width]; - for ( int i=0; i<_width; i++ ) { + qbit b[width]; + for ( int i=0; i>i)&1 ); } qbit ancilla[1]; PrepZ ( ancilla[0], 0 ); - endian (_width, lower); - endian (_width, b); - cUa ( upper[0], _width, 7, 13, lower, b, _N, ancilla[0] ); - cUa ( upper[1], _width, 4, 4, lower, b, _N, ancilla[0] ); - for (int i=2; i>c)&1) { + // Rz ( upper[0], -M_PI/pow(2,_iteration-c) ); + // } + // } + // H(upper[0]); - for (int i=0; i Date: Wed, 10 Oct 2018 15:32:43 -0400 Subject: [PATCH 06/24] Synching simulation scripts --- .../iterative_phase_estimation_qchem.py | 39 ++++++++++++ .../iterative_phase_estimation_shors.py | 62 +++++++++++++++++++ scripts/simulation/remove_global_phase.py | 5 +- scripts/simulation/trace_distance.py | 1 + 4 files changed, 106 insertions(+), 1 deletion(-) create mode 100755 scripts/simulation/iterative_phase_estimation_qchem.py create mode 100755 scripts/simulation/iterative_phase_estimation_shors.py diff --git a/scripts/simulation/iterative_phase_estimation_qchem.py b/scripts/simulation/iterative_phase_estimation_qchem.py new file mode 100755 index 000000000..e34920d54 --- /dev/null +++ b/scripts/simulation/iterative_phase_estimation_qchem.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python +import sys +import re +from subprocess import check_output + +upper_sz = 10 +estimate = 0 + +for iteration in range(0, upper_sz): + + with open(sys.argv[1]) as infile, open(sys.argv[1]+'.ipe', 'w') as outfile: + + for line in infile: + if "#define _upper_sz" in line: + line = re.sub('\d+', str(upper_sz), line) + if "#define _iteration" in line: + line = re.sub('\d+', str(iteration), line) + if "#define _estimate" in line: + line = re.sub('\d+', str(estimate), line) + outfile.write(line) + + sim_out = check_output(['./scripts/simulation/simulate.sh', './Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold.ipe']) + + state_lines = False + for line in sim_out.split('\n'): + if line.strip() == "--------------[quantum state]--------------": + print("--------------[quantum state]--------------") + state_lines = True + elif line.strip() == "-------------------------------------------": + print("-------------------------------------------") + state_lines = False + elif state_lines: + print(line) + strings = re.findall(r"[-+]?\d*\.\d+|\d+", line) + if (float(strings[0]) > 0.5): + meas_val = int(strings[4][1]) + print(meas_val) + estimate += meas_val * pow(2,iteration) + print(estimate) diff --git a/scripts/simulation/iterative_phase_estimation_shors.py b/scripts/simulation/iterative_phase_estimation_shors.py new file mode 100755 index 000000000..a60b47668 --- /dev/null +++ b/scripts/simulation/iterative_phase_estimation_shors.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python +import sys +import re +from subprocess import check_output + +def egcd(a, b): + if a == 0: + return (b, 0, 1) + else: + g, y, x = egcd(b % a, a) + return (g, x - (b // a) * y, y) + +def modinv(a, m): + g, x, y = egcd(a, m) + if g != 1: + raise Exception('modular inverse does not exist') + else: + return x % m + +a_val = [1,1,4,7] +a_inv = [1,1,4,13] +upper_sz = 4 +estimate = 0 + +for iteration in range(0, upper_sz): + + with open(sys.argv[1]) as infile, open(sys.argv[1]+'.ipe', 'w') as outfile: + + for line in infile: + if "#define _a_val" in line: + line = re.sub('\d+', str(a_val[iteration]), line) + if "#define _a_inv" in line: + line = re.sub('\d+', str(a_inv[iteration]), line) + if "#define _upper_sz" in line: + line = re.sub('\d+', str(1), line) + if "#define _iteration" in line: + line = re.sub('\d+', str(iteration), line) + if "#define _estimate" in line: + line = re.sub('\d+', str(estimate), line) + outfile.write(line) + + sim_out = check_output(['./scripts/simulation/simulate.sh', './Algorithms/Shors/Shors_15/shors.n3.scaffold.ipe']) + + state_lines = False + for line in sim_out.split('\n'): + if line.strip() == "--------------[quantum state]--------------": + print("--------------[quantum state]--------------") + state_lines = True + elif line.strip() == "-------------------------------------------": + print("-------------------------------------------") + state_lines = False + elif state_lines: + print(line) + strings = re.findall(r"[-+]?\d*\.\d+|\d+", line) + if (float(strings[0]) > 0.5): + meas_val = int(strings[4][-1]) + print(meas_val) + estimate += meas_val * pow(2,iteration) + print(estimate) + + # a_val = (a_val*a_val) % 15 + # a_inv = modinv(a_val, 15) diff --git a/scripts/simulation/remove_global_phase.py b/scripts/simulation/remove_global_phase.py index da1a3ede7..234ff9489 100644 --- a/scripts/simulation/remove_global_phase.py +++ b/scripts/simulation/remove_global_phase.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python + import sys import re from math import sqrt @@ -43,6 +45,7 @@ strings[2] ) - print(out_string) + if (magnitude>0.00390625): + print(out_string) outfile.write(out_string) outfile.write('\n') diff --git a/scripts/simulation/trace_distance.py b/scripts/simulation/trace_distance.py index b386f4657..c9b9c8f77 100644 --- a/scripts/simulation/trace_distance.py +++ b/scripts/simulation/trace_distance.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python import sys import re from math import sqrt From 5cb64c7c7a7ba224943ab7bd585b20b9446ef175 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Wed, 10 Oct 2018 15:34:03 -0400 Subject: [PATCH 07/24] Moving scripts to simulation folder --- iterative_phase_estimation_qchem.py | 39 ----------------------------- 1 file changed, 39 deletions(-) delete mode 100755 iterative_phase_estimation_qchem.py diff --git a/iterative_phase_estimation_qchem.py b/iterative_phase_estimation_qchem.py deleted file mode 100755 index b0651ee26..000000000 --- a/iterative_phase_estimation_qchem.py +++ /dev/null @@ -1,39 +0,0 @@ -#!/usr/bin/env python -import sys -import re -from subprocess import check_output - -upper_sz = 9 -estimate = 0 - -for iteration in range(0, upper_sz): - - with open(sys.argv[1]) as infile, open(sys.argv[1]+'.ipe', 'w') as outfile: - - for line in infile: - if "#define _upper_sz" in line: - line = re.sub('\d+', str(upper_sz), line) - if "#define _iteration" in line: - line = re.sub('\d+', str(iteration), line) - if "#define _estimate" in line: - line = re.sub('\d+', str(estimate), line) - outfile.write(line) - - sim_out = check_output(['./scripts/simulation/simulate.sh', './Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold.ipe']) - - state_lines = False - for line in sim_out.split('\n'): - if line.strip() == "--------------[quantum state]--------------": - print("--------------[quantum state]--------------") - state_lines = True - elif line.strip() == "-------------------------------------------": - print("-------------------------------------------") - state_lines = False - elif state_lines: - print(line) - strings = re.findall(r"[-+]?\d*\.\d+|\d+", line) - if (float(strings[0]) > 0.5): - meas_val = int(strings[4][1]) - print(meas_val) - estimate += meas_val * pow(2,iteration) - print(estimate) From 8aa9e8dfeb4d45c155f0282fe03c9624574d24ca Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Fri, 12 Oct 2018 17:20:37 -0400 Subject: [PATCH 08/24] Synching scripts --- scripts/simulation/iterative_phase_estimation_qchem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/simulation/iterative_phase_estimation_qchem.py b/scripts/simulation/iterative_phase_estimation_qchem.py index e34920d54..b007781e0 100755 --- a/scripts/simulation/iterative_phase_estimation_qchem.py +++ b/scripts/simulation/iterative_phase_estimation_qchem.py @@ -3,7 +3,7 @@ import re from subprocess import check_output -upper_sz = 10 +upper_sz = 13 estimate = 0 for iteration in range(0, upper_sz): From a014f220055ea2c3a46d9ab9601084040dc87f02 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Fri, 12 Oct 2018 17:25:45 -0400 Subject: [PATCH 09/24] Cleanup --- Library/QFT/QFT.scaffold | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Library/QFT/QFT.scaffold b/Library/QFT/QFT.scaffold index b579c112d..470f26a9f 100644 --- a/Library/QFT/QFT.scaffold +++ b/Library/QFT/QFT.scaffold @@ -15,7 +15,7 @@ module QFT ( const unsigned int width, qbit qbits[] ) { for ( int t = 0; t < width; t++ ) { H(qbits[t]); for ( int c = t+1; c < width; c++ ) { - cRz ( qbits[c], qbits[t], M_PI/(pow(2,c-t)) ); + cRz ( qbits[c], qbits[t], M_PI/pow(2,c-t) ); } } @@ -27,7 +27,7 @@ module iQFT ( const unsigned int width, qbit qbits[] ) { for ( int t = width-1; t >= 0; t-- ) { for ( int c = width-1; c > t; c-- ) { - cRz ( qbits[c], qbits[t], -M_PI/(pow(2,c-t)) ); + cRz ( qbits[c], qbits[t], -M_PI/pow(2,c-t) ); } H(qbits[t]); } From bf2636f42a7e039b463d6872ff6c9c57bfe4a33b Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Fri, 12 Oct 2018 18:03:06 -0400 Subject: [PATCH 10/24] Correct version of H2 ground state estimation, with results matching the Lanyon & Whitfield 2010 paper --- .../ground_state_estimation.h2.scaffold | 139 +++++++----------- 1 file changed, 56 insertions(+), 83 deletions(-) diff --git a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold index b957c2094..bbc6ed721 100644 --- a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold +++ b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold @@ -1,4 +1,5 @@ #include +#include <../../Library/cRz/cRz.scaffold> #define M 4 //number of bottom register qubits #define TROTTER 4 @@ -20,57 +21,57 @@ const enum spin sigma[M] = { // info=0.00,0.7414 // nuc=0.713776188 // Ehf=-1.116685636 -// 0,0=-1.252477495 -// 1,1=-0.475934275 -// 0,0,0,0=0.674493166 -// 0,1,0,1=0.181287518 -// 0,1,1,0=0.663472101 +// 0,0=-1.252477495 -1.2563390710389655 +// 1,1=-0.475934275 -0.4718960093502949 +// 0,0,0,0=0.674493166 0.6757101541858549 +// 0,1,0,1=0.181287518 0.18093119996471013 +// 0,1,1,0=0.663472101 0.664581729691277 // 1,1,1,1=0.697398010 -// h_double[0][1][1][0]=h_double[1][0][0][1]=0.674493166; -// h_double[0][2][2][0]=h_double[2][0][0][2]=0.663472101; -// h_double[0][3][3][0]=h_double[3][0][0][3]=0.663472101; -// h_double[1][2][2][1]=h_double[2][1][1][2]=0.663472101; -// h_double[1][3][3][1]=h_double[3][1][1][3]=0.663472101; +// h_double[0][1][1][0]=h_double[1][0][0][1]=0.6757101541858549; +// h_double[0][2][2][0]=h_double[2][0][0][2]=0.664581729691277; +// h_double[0][3][3][0]=h_double[3][0][0][3]=0.664581729691277; +// h_double[1][2][2][1]=h_double[2][1][1][2]=0.664581729691277; +// h_double[1][3][3][1]=h_double[3][1][1][3]=0.664581729691277; // h_double[2][3][3][2]=h_double[3][2][2][3]=0.697398010; // -// h_double[0][2][0][2]=h_double[1][3][1][3]=0.181287518; -// h_double[2][0][2][0]=h_double[3][1][3][1]=0.181287518; -// h_double[0][3][1][2]=h_double[0][1][3][2]=0.181287518; +// h_double[0][2][0][2]=h_double[1][3][1][3]=0.18093119996471013; +// h_double[2][0][2][0]=h_double[3][1][3][1]=0.18093119996471013; +// h_double[0][3][1][2]=h_double[0][1][3][2]=0.18093119996471013; // according to 1.4768229 -// h_double[2][1][3][0]=h_double[2][3][1][0]=0.181287518; +// h_double[2][1][3][0]=h_double[2][3][1][0]=0.18093119996471013; const double h_single[M][M] = { - {-1.252477495, 0.0, 0.0, 0.0}, - {0.0, -1.252477495, 0.0, 0.0}, - {0.0, 0.0, -0.475934275, 0.0}, - {0.0, 0.0, 0.0, -0.475934275} + {-1.2563390710389655, 0.0, 0.0, 0.0}, + {0.0, -1.2563390710389655, 0.0, 0.0}, + {0.0, 0.0, -0.4718960093502949, 0.0}, + {0.0, 0.0, 0.0, -0.4718960093502949} }; const double h_double[M][M][M][M] = { { {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, - {{0.0,0.0,0.0,0.0},{0.674493166,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.181287518,0.0}}, - {{0.0,0.0,0.181287518,0.0},{0.0,0.0,0.0,0.0},{0.663472101,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, - {{0.0,0.0,0.0,0.0},{0.0,0.0,0.181287518,0.0},{0.0,0.0,0.0,0.0},{0.663472101,0.0,0.0,0.0}} + {{0.0,0.0,0.0,0.0},{0.6757101541858549,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.18093119996471013,0.0}}, + {{0.0,0.0,0.18093119996471013,0.0},{0.0,0.0,0.0,0.0},{0.664581729691277,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, + {{0.0,0.0,0.0,0.0},{0.0,0.0,0.18093119996471013,0.0},{0.0,0.0,0.0,0.0},{0.664581729691277,0.0,0.0,0.0}} }, { - {{0.0,0.674493166,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, + {{0.0,0.6757101541858549,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, - {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.663472101,0.0,0.0},{0.0,0.0,0.0,0.0}}, - {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.181287518},{0.0,0.0,0.0,0.0},{0.0,0.663472101,0.0,0.0}} + {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.664581729691277,0.0,0.0},{0.0,0.0,0.0,0.0}}, + {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.18093119996471013},{0.0,0.0,0.0,0.0},{0.0,0.664581729691277,0.0,0.0}} }, { - {{0.0,0.0,0.663472101,0.0},{0.0,0.0,0.0,0.0},{0.181287518,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, - {{0.0,0.0,0.0,0.0},{0.0,0.0,0.663472101,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, + {{0.0,0.0,0.664581729691277,0.0},{0.0,0.0,0.0,0.0},{0.18093119996471013,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, + {{0.0,0.0,0.0,0.0},{0.0,0.0,0.664581729691277,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.697398010,0.0}} }, { - {{0.0,0.0,0.0,0.663472101},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, - {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.663472101},{0.0,0.0,0.0,0.0},{0.0,0.181287518,0.0,0.0}}, + {{0.0,0.0,0.0,0.664581729691277},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}, + {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.664581729691277},{0.0,0.0,0.0,0.0},{0.0,0.18093119996471013,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.697398010},{0.0,0.0,0.0,0.0}}, {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}} } @@ -84,21 +85,7 @@ module ControlledPhase ( qbit target, double theta ) { - Rz(control, theta); -} - -// Implements a controlled rotation -// inputs are the target qubit, control qubit and the rotation angle -module ControlledRotation ( - qbit control, - qbit target, - double theta -) { - // Rz (control, theta / 2.0); - CNOT (control, target); - Rz (target, theta / 2.0); - CNOT (control, target); - Rz (target, -1.0 * theta / 2.0); + // Rz(control, theta); } // Implement the Y gate defined in the GFIs Y=Rx(-PI/2), we'll name it Y_Rx @@ -146,10 +133,10 @@ module DoubleExcitationOperator ( } CNOT(reg[p], reg[q]); CNOT(reg[q], reg[r]); CNOT(reg[r], reg[s]); switch (i) { // M operator (as defined in A1) - case 0: ControlledRotation(control, reg[s], theta0); break; - case 1: ControlledRotation(control, reg[s], theta1); break; - case 2: ControlledRotation(control, reg[s], theta2); break; - case 3: ControlledRotation(control, reg[s], theta3); break; + case 0: cRz(control, reg[s], theta0); break; + case 1: cRz(control, reg[s], theta1); break; + case 2: cRz(control, reg[s], theta2); break; + case 3: cRz(control, reg[s], theta3); break; } CNOT(reg[r], reg[s]); CNOT(reg[q], reg[r]); CNOT(reg[p], reg[q]); switch (i) { // M^dagger operator is identical since matrices are Hermitian @@ -165,15 +152,14 @@ int main () { qbit bottomregister[M]; - // for ( int i=0; i=0; t--) { for (int c=t-1; c>=0; c--) { @@ -239,14 +223,13 @@ int main () double eta = 0.25 * h_double[c][t][t][c]; if (sigma[c]==sigma[t]) eta -= 0.25*h_double[c][t][c][t]; - ControlledRotation(topregister[0], bottomregister[t], eta*time_step); + cRz(topregister[0], bottomregister[t], eta*time_step); CNOT(bottomregister[c], bottomregister[t]); } } // two electron operators: excitation-exitation operators - // TODO: something in here causing Trotter to not converge DoubleExcitationOperator ( topregister[0], bottomregister, @@ -259,21 +242,11 @@ int main () } } - // Rz(topregister[0],-M_PI/64);// - // Rz(topregister[0],-M_PI/32);// - // Rz(topregister[0],-M_PI/16);// - // Rz(topregister[0],-M_PI/8);// - // Rz(topregister[0],-M_PI/4);// - // Rz(topregister[0],-M_PI/2);// for (int c=0; c<_iteration; c++) { if ((_estimate>>c)&1) { Rz ( topregister[0], -M_PI/pow(2,_iteration-c) ); } } H(topregister[0]); - // MeasZ(bottomregister[0]); - // MeasZ(bottomregister[1]); - // MeasZ(bottomregister[2]); - // MeasZ(bottomregister[3]); - // MeasZ(topregister[0]); + } From 41c8905882eaa13c8c106b47b31fc2e52c9e248a Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Thu, 15 Nov 2018 20:22:59 -0500 Subject: [PATCH 11/24] Moving to cluster simulation --- .../Ground_State_Estimation/ground_state_estimation.h2.scaffold | 2 +- Algorithms/Shors/Shors_15/shors.n3.scaffold | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold index bbc6ed721..1fa7f5357 100644 --- a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold +++ b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold @@ -2,7 +2,7 @@ #include <../../Library/cRz/cRz.scaffold> #define M 4 //number of bottom register qubits -#define TROTTER 4 +#define TROTTER 8 #define _upper_sz 4 #define _iteration 0 diff --git a/Algorithms/Shors/Shors_15/shors.n3.scaffold b/Algorithms/Shors/Shors_15/shors.n3.scaffold index 22445f48a..dd9452d73 100644 --- a/Algorithms/Shors/Shors_15/shors.n3.scaffold +++ b/Algorithms/Shors/Shors_15/shors.n3.scaffold @@ -80,6 +80,8 @@ int main ( endian (_upper_sz, upper); iQFT(_upper_sz, upper); + endian (_upper_sz, upper); + // or you can do iterative phase estimation instead of iQFT: // for (int c=0; c<_iteration; c++) { // if ((_estimate>>c)&1) { From 5c82e54885e7e888265d4eeafc64b16a491319c7 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Thu, 15 Nov 2018 20:26:22 -0500 Subject: [PATCH 12/24] Moving to cluster simulation --- Library/QFT/QFT_test.scaffold | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Library/QFT/QFT_test.scaffold b/Library/QFT/QFT_test.scaffold index 03bc800f0..3d57f3403 100644 --- a/Library/QFT/QFT_test.scaffold +++ b/Library/QFT/QFT_test.scaffold @@ -1,16 +1,21 @@ #include "QFT.scaffold" -#define width 8 +#define width 4 int main () { qbit reg[width]; for ( int i = 0; i < width; i++ ) { + // PrepZ(reg[i], 0); PrepZ(reg[i], (i+1)%2); } QFT(width, reg); - // iQFT(width, reg); + + for ( int i=0; i Date: Thu, 15 Nov 2018 20:29:52 -0500 Subject: [PATCH 13/24] moving to cluster simulation --- scripts/simulation/assert_integer.py | 28 ++++++ scripts/simulation/assert_product.py | 40 +++++++++ scripts/simulation/assert_uniform.py | 22 +++++ scripts/simulation/chisquare.py | 65 ++++++++++++++ scripts/simulation/fidelity.py | 44 ++++++++++ .../iterative_phase_estimation_qchem.py | 2 +- scripts/simulation/remove_global_phase.py | 88 ++++++++++++------- scripts/simulation/scaffassert.py | 47 ++++++++++ scripts/simulation/simulate.sh | 23 +++-- 9 files changed, 319 insertions(+), 40 deletions(-) create mode 100755 scripts/simulation/assert_integer.py create mode 100755 scripts/simulation/assert_product.py create mode 100755 scripts/simulation/assert_uniform.py create mode 100644 scripts/simulation/chisquare.py create mode 100644 scripts/simulation/fidelity.py mode change 100644 => 100755 scripts/simulation/remove_global_phase.py create mode 100755 scripts/simulation/scaffassert.py diff --git a/scripts/simulation/assert_integer.py b/scripts/simulation/assert_integer.py new file mode 100755 index 000000000..3876912a7 --- /dev/null +++ b/scripts/simulation/assert_integer.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python + +import sys +import csv + +import warnings +warnings.filterwarnings("ignore", message="numpy.dtype size changed") +from scipy import stats + +tallies_model = {} +tallies_dut = {} + +qubit_count = int(sys.argv[2]) +integer_val = int(sys.argv[3]) + +for val in range(1<0.00390625): + (outcome['polar_r'],outcome['polar_phi']) = cmath.polar(complex(float(strings[0]),float(strings[1]))) + + # if (outcome['polar_r']>0.0): + if (outcome['polar_r']>1.0/256.0): + + # if this is the first basis state then this is the global phase to factor out + if first_basis: + phase_global = outcome['polar_phi'] + first_basis = False + + outcome['polar_phi'] -= phase_global + rect = cmath.rect(outcome['polar_r'],outcome['polar_phi']) + outcome['rect_real'] = rect.real + outcome['rect_imag'] = rect.imag + + # print zip(strings[2][::-1],reg_mapping) + for elem in zip(strings[2][::-1],reg_map): + outcome[elem[1][0]] += int(elem[0]) << elem[1][1] + + out_string = 'polar: ({:f},{:f}) rect: ({:f},{:f}) |{:s}>'.format( + outcome['polar_r'], + outcome['polar_phi'], + outcome['rect_real'], + outcome['rect_imag'], + strings[2], + ) print(out_string) - outfile.write(out_string) - outfile.write('\n') + writer.writerow(outcome) diff --git a/scripts/simulation/scaffassert.py b/scripts/simulation/scaffassert.py new file mode 100755 index 000000000..a1670f7a4 --- /dev/null +++ b/scripts/simulation/scaffassert.py @@ -0,0 +1,47 @@ +#!/usr/bin/env python + +import sys, re, copy, os + +check_indx = 0 +checkpoints = [[]] + +with open(sys.argv[1]) as source: + for line in source: + + if "assert_integer" in line: + # start the next checkpoint + checkpoints.append(copy.copy(checkpoints[check_indx])) + + param_string = re.findall("assert_integer\s*\((.*?)\)\s*;", line.strip()) + params = [param.strip() for param in param_string[0].split(',')] + + # close out this checkpoint + checkpoints[check_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) + checkpoints[check_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) + checkpoints[check_indx].append('}\n') + + check_indx+=1 + + if "assert_uniform" in line: + # start the next checkpoint + checkpoints.append(copy.copy(checkpoints[check_indx])) + + param_string = re.findall("assert_uniform\s*\((.*?)\)\s*;", line.strip()) + params = [param.strip() for param in param_string[0].split(',')] + + # close out this checkpoint + checkpoints[check_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) + checkpoints[check_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) + checkpoints[check_indx].append('}\n') + + check_indx+=1 + + else: + checkpoints[check_indx].append(line) + +# generate checkpoints +for check_indx in range(len(checkpoints)): + out_name = os.path.basename(sys.argv[1]) + str(check_indx) + '.scaffold' + with open(out_name, 'w') as outfile: + for line in checkpoints[check_indx]: + outfile.write(line) diff --git a/scripts/simulation/simulate.sh b/scripts/simulation/simulate.sh index e5f006474..046f0d1c2 100755 --- a/scripts/simulation/simulate.sh +++ b/scripts/simulation/simulate.sh @@ -7,9 +7,22 @@ then tar -xvf ../qx_simulator_linux_x86_64.tar.gz -C ../ fi -./scaffold.sh -c $1 -./scaffold.sh -s $1 source=${1##*/} -# echo ${source%.scaffold} -../qx_simulator_linux_x86_64/qx_simulator_1.0.beta_linux_x86_64 ${source%.scaffold}.qc > ${source%.scaffold}.global_phase.output -python scripts/simulation/remove_global_phase.py ${source%.scaffold}.global_phase.output ${source%.scaffold}.no_global_phase.output + +./scaffold.sh -c $1 + +scripts/simulation/scaffassert.py $1 +# ./scaffold.sh -s $1 +# # echo ${source%.scaffold} +# +# # map +# rm ${source%.scaffold}.csv +# for i in `seq 1 1`; +# do +# ../qx_simulator_linux_x86_64/qx_simulator_1.0.beta_linux_x86_64 ${source%.scaffold}.qc | scripts/simulation/remove_global_phase.py ${source%.scaffold}.qc ${source%.scaffold}.csv +# done +# +# # reduce +# # scripts/simulation/assert_uniform.py ${source%.scaffold}.csv 4 +# # scripts/simulation/assert_integer.py ${source%.scaffold}.csv 5 30 +# scripts/simulation/assert_product.py ${source%.scaffold}.csv From 7e171ad022f2c310bf711a743920168e4d67b854 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Thu, 15 Nov 2018 20:31:07 -0500 Subject: [PATCH 14/24] Moving to cluster simulation --- Library/QFT/QFT_test.scaffassert | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 Library/QFT/QFT_test.scaffassert diff --git a/Library/QFT/QFT_test.scaffassert b/Library/QFT/QFT_test.scaffassert new file mode 100644 index 000000000..361c6bf0b --- /dev/null +++ b/Library/QFT/QFT_test.scaffassert @@ -0,0 +1,30 @@ +#include "QFT.scaffold" +#define width 4 + +int main () { + qbit reg[width]; + + for ( int i = 0; i < width; i++ ) { + // PrepZ(reg[i], 0); + PrepZ(reg[i], (i+1)%2); + } + + assert_classical ( reg, width, 10 ); + + QFT(width, reg); + + assert_superposition ( reg, width ); + + iQFT(width, reg); + + assert_classical ( reg, width, 10 ); + + for ( int i=0; i Date: Wed, 28 Nov 2018 14:58:45 -0500 Subject: [PATCH 15/24] Simulation scripts for assertions testing --- scripts/simulation/assert_integer.py | 2 + scripts/simulation/assert_product.py | 18 ++++ scripts/simulation/assert_uniform.py | 2 + scripts/simulation/chisquare.py | 65 ------------- scripts/simulation/register_value_csv.py | 80 ++++++++++++++++ scripts/simulation/remove_global_phase.py | 71 -------------- scripts/simulation/scaffassert.py | 58 ++++++------ scripts/simulation/simulate.bash | 107 ++++++++++++++++++++++ scripts/simulation/simulate.sh | 28 ------ 9 files changed, 240 insertions(+), 191 deletions(-) delete mode 100644 scripts/simulation/chisquare.py create mode 100755 scripts/simulation/register_value_csv.py delete mode 100755 scripts/simulation/remove_global_phase.py create mode 100755 scripts/simulation/simulate.bash delete mode 100755 scripts/simulation/simulate.sh diff --git a/scripts/simulation/assert_integer.py b/scripts/simulation/assert_integer.py index 3876912a7..d3c9aa6e8 100755 --- a/scripts/simulation/assert_integer.py +++ b/scripts/simulation/assert_integer.py @@ -26,3 +26,5 @@ print tallies_model print tallies_dut print stats.chisquare(tallies_dut.values(), f_exp=tallies_model.values()) + +# chisquare([16, 18, 16, 14, 12, 12], f_exp=[16, 16, 16, 16, 16, 8]) diff --git a/scripts/simulation/assert_product.py b/scripts/simulation/assert_product.py index 02a739e1c..5b1495112 100755 --- a/scripts/simulation/assert_product.py +++ b/scripts/simulation/assert_product.py @@ -38,3 +38,21 @@ print cross_tab print stats.chi2_contingency(cross_tab) + + +# chi2_stat, p_val, dof, ex = stats.chi2_contingency(shor) +# +# print("===Chi2 Stat===") +# print(chi2_stat) +# print("\n") +# +# print("===Degrees of Freedom===") +# print(dof) +# print("\n") +# +# print("===P-Value===") +# print(p_val) +# print("\n") +# +# print("===Contingency Table===") +# print(ex) diff --git a/scripts/simulation/assert_uniform.py b/scripts/simulation/assert_uniform.py index e6f4ecf82..bc4ee36fb 100755 --- a/scripts/simulation/assert_uniform.py +++ b/scripts/simulation/assert_uniform.py @@ -20,3 +20,5 @@ print tallies print stats.chisquare(tallies.values()) + +# chisquare([16, 18, 16, 14, 12, 12], f_exp=[16, 16, 16, 16, 16, 8]) diff --git a/scripts/simulation/chisquare.py b/scripts/simulation/chisquare.py deleted file mode 100644 index b850943cf..000000000 --- a/scripts/simulation/chisquare.py +++ /dev/null @@ -1,65 +0,0 @@ -import warnings -warnings.filterwarnings("ignore", message="numpy.dtype size changed") -# warnings.filterwarnings("ignore", message="numpy.ufunc size changed") -# -# from scipy.stats import chisquare -from scipy import stats - -import numpy as np - -# a1 = [1/8, 0.001, 1/8, 0.001, 1/8, 0.001, 1/8, 0.001] -# a2 = [1/64, 1/64, 1/64, 1/64, 1/64, 1/64, 1/64, 1/64] -# a3 = [1/64, 1/64, 1/64, 1/64, 1/64, 1/64, 1/64, 1/64] -# a4 = [1/64, 1/64, 1/64, 1/64, 1/64, 1/64, 1/64, 1/64] -# a5 = [1/64, 1/64, 1/64, 1/64, 1/64, 1/64, 1/64, 1/64] - -# a1 = [8, 8, 8, 8] -# a2 = [8, 8, 8, 8] -# a3 = [8, 8, 8, 8] -# a4 = [8, 8, 8, 8] -# a5 = [8, 8, 8, 8] -# a1 = [8, 0, 8, 0, 8, 0, 8, 0] -# a2 = [1, 1, 1, 1, 1, 1, 1, 1] -# a3 = [1, 1, 1, 1, 1, 1, 1, 1] -# a4 = [1, 1, 1, 1, 1, 1, 1, 1] -# a5 = [1, 1, 1, 1, 1, 1, 1, 1] - -# a1 = [8,1,1,1,1] -# a2 = [0,1,1,1,1] -# a3 = [8,1,1,1,1] -# a4 = [0,1,1,1,1] -# a5 = [8,1,1,1,1] -# a6 = [0,1,1,1,1] -# a7 = [8,1,1,1,1] -# a8 = [0,1,1,1,1] - - -# shor = 32*np.array([a2, a3, a4, a5]) -# shor = 32*np.array([a1, a2, a3, a4, a5]) -# shor = np.array([a1, a2, a3, a4, a5, a6, a7, a8]) -# shor = np.transpose(np.array([a1, a2, a3, a4, a5])) - - -# -print stats.chisquare([1,6]) -# chisquare([16, 18, 16, 14, 12, 12], f_exp=[16, 16, 16, 16, 16, 8]) -# -# chi2_stat, p_val, dof, ex = stats.chi2_contingency(shor) -# -# print("===Chi2 Stat===") -# print(chi2_stat) -# print("\n") -# -# print("===Degrees of Freedom===") -# print(dof) -# print("\n") -# -# print("===P-Value===") -# print(p_val) -# print("\n") -# -# print("===Contingency Table===") -# print(ex) -# scipy.stats.entropy(pk, qk=None, base=None)[source] - -print "done" diff --git a/scripts/simulation/register_value_csv.py b/scripts/simulation/register_value_csv.py new file mode 100755 index 000000000..e79f42789 --- /dev/null +++ b/scripts/simulation/register_value_csv.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python +import sys, re, cmath, csv, os, glob + +regex = re.compile("([a-zA-Z]+)([0-9]+)") +outcome = { 'polar_r':0, 'polar_phi':0, 'rect_real':0, 'rect_imag':0 } +reg_map = [] + +# open QX source file +with open(sys.argv[1], 'r') as qcfile: + for line in qcfile: + if 'map q' in line: + print line + words = line.split() + reg_index_pair = regex.match(words[2]) + register = reg_index_pair.group(1) + index = int(reg_index_pair.group(2)) + + reg_map.append((register,index)) + outcome[register] = 0 + +# print reg_map + +with open(sys.argv[2], 'a+') as csvfile: + + writer = csv.DictWriter(csvfile, fieldnames=list(outcome.keys())) + if not csvfile.read(): + writer.writeheader() + + filenames = glob.glob(os.path.splitext(sys.argv[1])[0]+'.trial_*.out') + print filenames + for filename in filenames: + with open(filename) as file: + + # Parse QX Simulator output format + state_lines = False + first_basis = True + phase_global = 0.0 + for line in file: + print line + if line.strip() == "--------------[quantum state]--------------": + print("--------------[quantum state]--------------") + state_lines = True + elif state_lines and line.strip() == "-------------------------------------------": + print("-------------------------------------------") + state_lines = False + elif state_lines: + strings = re.findall(r"[-+]?\d*\.\d+|\d+", line) + + # zeros an outcome structure + outcome = dict.fromkeys(outcome, 0) + + # grab real and imaginary parts of basis state amplitude + (outcome['polar_r'],outcome['polar_phi']) = cmath.polar(complex(float(strings[0]),float(strings[1]))) + + # if (outcome['polar_r']>0.0): + if (outcome['polar_r']>1.0/256.0): + + # if this is the first basis state then this is the global phase to factor out + if first_basis: + phase_global = outcome['polar_phi'] + first_basis = False + + outcome['polar_phi'] -= phase_global + rect = cmath.rect(outcome['polar_r'],outcome['polar_phi']) + outcome['rect_real'] = rect.real + outcome['rect_imag'] = rect.imag + + # print zip(strings[2][::-1],reg_mapping) + for elem in zip(strings[2][::-1],reg_map): + outcome[elem[1][0]] += int(elem[0]) << elem[1][1] + + out_string = 'polar: ({:f},{:f}) rect: ({:f},{:f}) |{:s}>'.format( + outcome['polar_r'], + outcome['polar_phi'], + outcome['rect_real'], + outcome['rect_imag'], + strings[2], + ) + print(out_string) + writer.writerow(outcome) diff --git a/scripts/simulation/remove_global_phase.py b/scripts/simulation/remove_global_phase.py deleted file mode 100755 index b00647ef7..000000000 --- a/scripts/simulation/remove_global_phase.py +++ /dev/null @@ -1,71 +0,0 @@ -#!/usr/bin/env python -import sys, re, cmath, csv - -regex = re.compile("([a-zA-Z]+)([0-9]+)") -outcome = { 'polar_r':0, 'polar_phi':0, 'rect_real':0, 'rect_imag':0 } -reg_map = [] - -# open QX source file -with open(sys.argv[1], 'r') as qcfile: - for line in qcfile: - words = line.split() - if words[0]=='map': - reg_index_pair = regex.match(words[2]) - register = reg_index_pair.group(1) - index = int(reg_index_pair.group(2)) - - outcome[register] = 0 - reg_map.append((register,index)) - # print reg_map - -with open(sys.argv[2], 'a+') as csvfile: - - writer = csv.DictWriter(csvfile, fieldnames=list(outcome.keys())) - if not csvfile.read(): - writer.writeheader() - - # Parse QX Simulator output format - state_lines = False - first_basis = True - phase_global = 0.0 - for line in sys.stdin: - if line.strip() == "--------------[quantum state]--------------": - print("--------------[quantum state]--------------") - state_lines = True - elif state_lines and line.strip() == "-------------------------------------------": - print("-------------------------------------------") - state_lines = False - elif state_lines: - strings = re.findall(r"[-+]?\d*\.\d+|\d+", line) - - outcome = dict.fromkeys(outcome, 0) - - # grab real and imaginary parts of basis state amplitude - (outcome['polar_r'],outcome['polar_phi']) = cmath.polar(complex(float(strings[0]),float(strings[1]))) - - # if (outcome['polar_r']>0.0): - if (outcome['polar_r']>1.0/256.0): - - # if this is the first basis state then this is the global phase to factor out - if first_basis: - phase_global = outcome['polar_phi'] - first_basis = False - - outcome['polar_phi'] -= phase_global - rect = cmath.rect(outcome['polar_r'],outcome['polar_phi']) - outcome['rect_real'] = rect.real - outcome['rect_imag'] = rect.imag - - # print zip(strings[2][::-1],reg_mapping) - for elem in zip(strings[2][::-1],reg_map): - outcome[elem[1][0]] += int(elem[0]) << elem[1][1] - - out_string = 'polar: ({:f},{:f}) rect: ({:f},{:f}) |{:s}>'.format( - outcome['polar_r'], - outcome['polar_phi'], - outcome['rect_real'], - outcome['rect_imag'], - strings[2], - ) - print(out_string) - writer.writerow(outcome) diff --git a/scripts/simulation/scaffassert.py b/scripts/simulation/scaffassert.py index a1670f7a4..92194f3cd 100755 --- a/scripts/simulation/scaffassert.py +++ b/scripts/simulation/scaffassert.py @@ -2,46 +2,50 @@ import sys, re, copy, os -check_indx = 0 -checkpoints = [[]] +break_indx = 0 +breakpoints = [[]] with open(sys.argv[1]) as source: for line in source: - if "assert_integer" in line: - # start the next checkpoint - checkpoints.append(copy.copy(checkpoints[check_indx])) + if re.findall("^assert_", line.strip()): - param_string = re.findall("assert_integer\s*\((.*?)\)\s*;", line.strip()) - params = [param.strip() for param in param_string[0].split(',')] + if "assert_classical" in line: + # start the next breakpoint + breakpoints.append(copy.copy(breakpoints[break_indx])) - # close out this checkpoint - checkpoints[check_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) - checkpoints[check_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) - checkpoints[check_indx].append('}\n') + param_string = re.findall("^assert_classical\s*\((.*?)\)\s*;", line.strip()) + params = [param.strip() for param in param_string[0].split(',')] - check_indx+=1 + # close out this breakpoint + breakpoints[break_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) + breakpoints[break_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) + breakpoints[break_indx].append('}\n') - if "assert_uniform" in line: - # start the next checkpoint - checkpoints.append(copy.copy(checkpoints[check_indx])) + break_indx+=1 - param_string = re.findall("assert_uniform\s*\((.*?)\)\s*;", line.strip()) - params = [param.strip() for param in param_string[0].split(',')] + elif "assert_superposition" in line: + # start the next breakpoint + breakpoints.append(copy.copy(breakpoints[break_indx])) - # close out this checkpoint - checkpoints[check_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) - checkpoints[check_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) - checkpoints[check_indx].append('}\n') + param_string = re.findall("^assert_superposition\s*\((.*?)\)\s*;", line.strip()) + params = [param.strip() for param in param_string[0].split(',')] - check_indx+=1 + # close out this breakpoint + breakpoints[break_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) + breakpoints[break_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) + breakpoints[break_indx].append('}\n') + + break_indx+=1 else: - checkpoints[check_indx].append(line) + breakpoints[break_indx].append(line) -# generate checkpoints -for check_indx in range(len(checkpoints)): - out_name = os.path.basename(sys.argv[1]) + str(check_indx) + '.scaffold' +# generate breakpoints +for break_indx in range(len(breakpoints)): + out_name = os.path.splitext(os.path.basename(sys.argv[1]))[0] + '.breakpoint_' + str(break_indx+1) + '.scaffold' with open(out_name, 'w') as outfile: - for line in checkpoints[check_indx]: + for line in breakpoints[break_indx]: outfile.write(line) + +print (break_indx+1) diff --git a/scripts/simulation/simulate.bash b/scripts/simulation/simulate.bash new file mode 100755 index 000000000..836c45665 --- /dev/null +++ b/scripts/simulation/simulate.bash @@ -0,0 +1,107 @@ +#!/bin/bash + +SCAFFCC_PATH=/n/fs/qdb/ScaffCC +TEST_NAME=${1%.scaffassert} + +# CLEANUP +$SCAFFCC_PATH/scaffold.sh -c $1 +rm *.breakpoint.*.scaffold +rm sbatch.bash +rm *.breakpoint.*.ll +rm *.breakpoint.*.tmp +rm *.breakpoint.*.resources +rm *.breakpoint.*.qasmh +rm *.breakpoint.*.qasmf +rm *.breakpoint.*.qc +rm *.out +rm *.err +rm *.csv + +# split into breakpoints +BREAKPOINTS=$($SCAFFCC_PATH/scripts/simulation/scaffassert.py $1) +ENSEMBLE=3 + +# compilation without slurm: +# for ((i=1;i<=$BREAKPOINTS;i++)) +# do +# $SCAFFCC_PATH/scaffold.sh -s $TEST_NAME.breakpoint_$i.scaffold +# done + +# get QX Simulator +if [ ! -e $SCAFFCC_PATH/../qx_simulator_linux_x86_64/qx_simulator_1.0.beta_linux_x86_64 ] +then + wget -nc -P $SCAFFCC_PATH/../ http://quantum-studio.net/qx_simulator_linux_x86_64.tar.gz + tar -xvf $SCAFFCC_PATH/../qx_simulator_linux_x86_64.tar.gz -C $SCAFFCC_PATH/../ +fi + +SBATCH_SCRIPT="sbatch.bash" +/bin/cat <$SBATCH_SCRIPT +#!/bin/bash + +#SBATCH --job-name=$TEST_NAME.simulate +#SBATCH --time=00:01:00 + +#SBATCH --mail-type=ALL +#SBATCH --mail-user=yipengh@cs.princeton.edu + +echo "HEAD =" $HOSTNAME +echo "LOCAL =" \$HOSTNAME +echo "SLURM_ARRAY_JOB_ID =" \$SLURM_ARRAY_JOB_ID +echo "SLURM_ARRAY_TASK_ID =" \$SLURM_ARRAY_TASK_ID +echo "SLURM_CLUSTER_NAME =" \$SLURM_CLUSTER_NAME +echo "SLURM_CPUS_PER_TASK =" \$SLURM_CPUS_PER_TASK +echo "SLURM_JOB_ACCOUNT =" \$SLURM_JOB_ACCOUNT +echo "SLURM_JOB_ID =" \$SLURM_JOB_ID +echo "SLURM_JOB_NAME =" \$SLURM_JOB_NAME +echo "SLURM_JOB_NODELIST =" \$SLURM_JOB_NODELIST +echo "SLURM_JOB_NUM_NODES =" \$SLURM_JOB_NUM_NODES +echo "SLURM_JOB_PARTITION =" \$SLURM_JOB_PARTITION +echo "SLURM_JOB_UID =" \$SLURM_JOB_UID +echo "SLURM_JOB_USER =" \$SLURM_JOB_USER +echo "SLURM_STEP_ID =" \$SLURM_STEP_ID + +###################### +# Begin work section # +###################### + +SCRATCH_PATH=/scratch/$USER/$TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID +mkdir -p \$SCRATCH_PATH +cd \$SCRATCH_PATH + +# COMPILE +# compile Scaffold to QASM to QX input +$SCAFFCC_PATH/scaffold.sh -s \$SLURM_SUBMIT_DIR/$TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID.scaffold +cp $TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID.qc \$SLURM_SUBMIT_DIR +cd \$SLURM_SUBMIT_DIR + +# noise model +# error_model depolarizing_channel, 0.001 + +# SIMULATE +# map +srun \ +--ntasks=$ENSEMBLE \ +--output=$TEST_NAME.breakpoint_%a.trial_%n.%A.%N.out \ +--error=$TEST_NAME.breakpoint_%a.trial_%n.%A.%N.err \ +$SCAFFCC_PATH/../qx_simulator_linux_x86_64/qx_simulator_1.0.beta_linux_x86_64 $TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID.qc + +# reduce +$SCAFFCC_PATH/scripts/simulation/register_value_csv.py \ +$TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID.qc \ +$TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID.csv + +# scripts/simulation/assert_uniform.py ${source%.scaffold}.csv 4 +# scripts/simulation/assert_integer.py ${source%.scaffold}.csv 5 30 +# scripts/simulation/assert_product.py ${source%.scaffold}.csv + +# users are expected to clean up their files at the end of each job +rm -r \$SCRATCH_PATH + +EOM + +sbatch \ +--array=1-$BREAKPOINTS \ +-N $((ENSEMBLE)) \ +--output=$TEST_NAME.breakpoint_%a.compile.%A.%N.out \ +--error=$TEST_NAME.breakpoint_%a.compile.%A.%N.err \ +$SBATCH_SCRIPT diff --git a/scripts/simulation/simulate.sh b/scripts/simulation/simulate.sh deleted file mode 100755 index 046f0d1c2..000000000 --- a/scripts/simulation/simulate.sh +++ /dev/null @@ -1,28 +0,0 @@ -#!/bin/sh - -# get QX Simulator -if [ ! -e ../qx_simulator_linux_x86_64/qx_simulator_1.0.beta_linux_x86_64 ] -then - wget -nc -P ../ http://quantum-studio.net/qx_simulator_linux_x86_64.tar.gz - tar -xvf ../qx_simulator_linux_x86_64.tar.gz -C ../ -fi - -source=${1##*/} - -./scaffold.sh -c $1 - -scripts/simulation/scaffassert.py $1 -# ./scaffold.sh -s $1 -# # echo ${source%.scaffold} -# -# # map -# rm ${source%.scaffold}.csv -# for i in `seq 1 1`; -# do -# ../qx_simulator_linux_x86_64/qx_simulator_1.0.beta_linux_x86_64 ${source%.scaffold}.qc | scripts/simulation/remove_global_phase.py ${source%.scaffold}.qc ${source%.scaffold}.csv -# done -# -# # reduce -# # scripts/simulation/assert_uniform.py ${source%.scaffold}.csv 4 -# # scripts/simulation/assert_integer.py ${source%.scaffold}.csv 5 30 -# scripts/simulation/assert_product.py ${source%.scaffold}.csv From faaf30d5a6d980c3e4fc62919a69878c52e4be20 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Thu, 29 Nov 2018 11:52:34 -0500 Subject: [PATCH 16/24] Updating simulation scripts for assertions testing --- scripts/simulation/assert_classical.py | 31 ++++++++++ scripts/simulation/assert_integer.py | 30 ---------- scripts/simulation/assert_product.py | 70 ++++++++++++---------- scripts/simulation/assert_superposition.py | 20 +++++++ scripts/simulation/assert_uniform.py | 24 -------- scripts/simulation/register_value_csv.py | 11 ++-- scripts/simulation/scaffassert.py | 43 +++++++++++++ scripts/simulation/simulate.bash | 27 +++++---- 8 files changed, 156 insertions(+), 100 deletions(-) create mode 100755 scripts/simulation/assert_classical.py delete mode 100755 scripts/simulation/assert_integer.py create mode 100755 scripts/simulation/assert_superposition.py delete mode 100755 scripts/simulation/assert_uniform.py diff --git a/scripts/simulation/assert_classical.py b/scripts/simulation/assert_classical.py new file mode 100755 index 000000000..6b77f293b --- /dev/null +++ b/scripts/simulation/assert_classical.py @@ -0,0 +1,31 @@ +#!/usr/bin/env python + +import sys, csv, warnings +import numpy as np +warnings.filterwarnings("ignore", message="numpy.dtype size changed") +from scipy import stats + +qubit_count = int(sys.argv[3]) +integer_val = int(sys.argv[4]) + +# generate golden model +tallies_model = {} +tallies_dut = {} +for val in range(1<0.0): if (outcome['polar_r']>1.0/256.0): + outcome['probability'] = outcome['polar_r'] * outcome['polar_r'] + # if this is the first basis state then this is the global phase to factor out if first_basis: phase_global = outcome['polar_phi'] @@ -69,7 +71,8 @@ for elem in zip(strings[2][::-1],reg_map): outcome[elem[1][0]] += int(elem[0]) << elem[1][1] - out_string = 'polar: ({:f},{:f}) rect: ({:f},{:f}) |{:s}>'.format( + out_string = 'probability: {:f} polar: ({:f},{:f}) rect: ({:f},{:f}) |{:s}>'.format( + outcome['probability'], outcome['polar_r'], outcome['polar_phi'], outcome['rect_real'], diff --git a/scripts/simulation/scaffassert.py b/scripts/simulation/scaffassert.py index 92194f3cd..859f49ec0 100755 --- a/scripts/simulation/scaffassert.py +++ b/scripts/simulation/scaffassert.py @@ -22,6 +22,15 @@ breakpoints[break_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) breakpoints[break_indx].append('}\n') + # generate script for checking the breakpoint + breakpoint_name = os.path.splitext(os.path.basename(sys.argv[1]))[0] + '.breakpoint_' + str(break_indx+1) + with open(breakpoint_name + '.bash', 'w') as outfile: + outfile.write('$SCAFFCC_PATH/scripts/simulation/assert_classical.py {:s}.csv {:s} {:s} {:s}'.format( + breakpoint_name, + params[0], + params[1], + params[2])) + break_indx+=1 elif "assert_superposition" in line: @@ -36,6 +45,40 @@ breakpoints[break_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) breakpoints[break_indx].append('}\n') + # generate script for checking the breakpoint + breakpoint_name = os.path.splitext(os.path.basename(sys.argv[1]))[0] + '.breakpoint_' + str(break_indx+1) + with open(breakpoint_name + '.bash', 'w') as outfile: + outfile.write('$SCAFFCC_PATH/scripts/simulation/assert_superposition.py {:s}.csv {:s} {:s}'.format( + breakpoint_name, + params[0], + params[1])) + + break_indx+=1 + + elif "assert_product" in line: + # start the next breakpoint + breakpoints.append(copy.copy(breakpoints[break_indx])) + + param_string = re.findall("^assert_product\s*\((.*?)\)\s*;", line.strip()) + params = [param.strip() for param in param_string[0].split(',')] + + # close out this breakpoint + breakpoints[break_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) + breakpoints[break_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) + breakpoints[break_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[3])) + breakpoints[break_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[2])) + breakpoints[break_indx].append('}\n') + + # generate script for checking the breakpoint + breakpoint_name = os.path.splitext(os.path.basename(sys.argv[1]))[0] + '.breakpoint_' + str(break_indx+1) + with open(breakpoint_name + '.bash', 'w') as outfile: + outfile.write('$SCAFFCC_PATH/scripts/simulation/assert_product.py {:s}.csv {:s} {:s} {:s} {:s}'.format( + breakpoint_name, + params[0], + params[1], + params[2], + params[3])) + break_indx+=1 else: diff --git a/scripts/simulation/simulate.bash b/scripts/simulation/simulate.bash index 836c45665..1395a78de 100755 --- a/scripts/simulation/simulate.bash +++ b/scripts/simulation/simulate.bash @@ -1,25 +1,26 @@ #!/bin/bash -SCAFFCC_PATH=/n/fs/qdb/ScaffCC +export SCAFFCC_PATH=/n/fs/qdb/ScaffCC TEST_NAME=${1%.scaffassert} # CLEANUP $SCAFFCC_PATH/scaffold.sh -c $1 -rm *.breakpoint.*.scaffold +rm *.breakpoint_*.scaffold rm sbatch.bash -rm *.breakpoint.*.ll -rm *.breakpoint.*.tmp -rm *.breakpoint.*.resources -rm *.breakpoint.*.qasmh -rm *.breakpoint.*.qasmf -rm *.breakpoint.*.qc -rm *.out -rm *.err -rm *.csv +rm *.breakpoint_*.ll +rm *.breakpoint_*.tmp +rm *.breakpoint_*.resources +rm *.breakpoint_*.qasmh +rm *.breakpoint_*.qasmf +rm *.breakpoint_*.qc +rm *.breakpoint_*.out +rm *.breakpoint_*.err +rm *.breakpoint_*.csv +rm *.breakpoint_*.bash # split into breakpoints BREAKPOINTS=$($SCAFFCC_PATH/scripts/simulation/scaffassert.py $1) -ENSEMBLE=3 +ENSEMBLE=8 # compilation without slurm: # for ((i=1;i<=$BREAKPOINTS;i++)) @@ -90,6 +91,8 @@ $SCAFFCC_PATH/scripts/simulation/register_value_csv.py \ $TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID.qc \ $TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID.csv +bash $TEST_NAME.breakpoint_\$SLURM_ARRAY_TASK_ID.bash + # scripts/simulation/assert_uniform.py ${source%.scaffold}.csv 4 # scripts/simulation/assert_integer.py ${source%.scaffold}.csv 5 30 # scripts/simulation/assert_product.py ${source%.scaffold}.csv From de5751a6e5ecc22d1f22920d0ad7de8e4f895753 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Thu, 29 Nov 2018 11:53:52 -0500 Subject: [PATCH 17/24] Assertions testing for QFT --- Library/QFT/QFT_test.scaffassert | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Library/QFT/QFT_test.scaffassert b/Library/QFT/QFT_test.scaffassert index 361c6bf0b..448e9557f 100644 --- a/Library/QFT/QFT_test.scaffassert +++ b/Library/QFT/QFT_test.scaffassert @@ -9,15 +9,15 @@ int main () { PrepZ(reg[i], (i+1)%2); } - assert_classical ( reg, width, 10 ); + assert_classical ( reg, 4, 5 ); QFT(width, reg); - assert_superposition ( reg, width ); + assert_superposition ( reg, 4 ); iQFT(width, reg); - assert_classical ( reg, width, 10 ); + assert_classical ( reg, 4, 5 ); for ( int i=0; i Date: Thu, 29 Nov 2018 11:54:35 -0500 Subject: [PATCH 18/24] Assertions testing scripts --- Algorithms/QFT/qft.scaffold | 42 --- .../Shors/Shors_15/shors.n3.scaffassert | 116 +++++++++ Algorithms/Shors/cADD/cADD_test.scaffassert | 48 ++++ .../Shors/cMODADD/cMODADD_test.scaffassert | 62 +++++ .../Shors/cMODMUL/cMODMUL_test.scaffassert | 50 ++++ .../Square_Root/square_root.n4.scaffassert | 239 ++++++++++++++++++ 6 files changed, 515 insertions(+), 42 deletions(-) delete mode 100644 Algorithms/QFT/qft.scaffold create mode 100644 Algorithms/Shors/Shors_15/shors.n3.scaffassert create mode 100644 Algorithms/Shors/cADD/cADD_test.scaffassert create mode 100644 Algorithms/Shors/cMODADD/cMODADD_test.scaffassert create mode 100644 Algorithms/Shors/cMODMUL/cMODMUL_test.scaffassert create mode 100644 Algorithms/Square_Root/square_root.n4.scaffassert diff --git a/Algorithms/QFT/qft.scaffold b/Algorithms/QFT/qft.scaffold deleted file mode 100644 index b81428949..000000000 --- a/Algorithms/QFT/qft.scaffold +++ /dev/null @@ -1,42 +0,0 @@ -#include -#define pi 3.141592653589793238462643383279502884197 - -#define n 5 - -module cRz ( qbit ctrl, qbit target, const double angle ) { - /* cRz identity matrix: - [ [ 1 0 0 0 ] - [ 0 1 0 0 ] - [ 0 0 e^(-i*angle/2) 0 ] - [ 0 0 0 e^(i*angle/2) ] ] - */ - - Rz(ctrl, angle/2); - CNOT(ctrl, target); - Rz(target, -angle/2); - CNOT(ctrl, target); - Rz(target, angle/2); -} - -module qft ( qbit bit[n] ) { - H(bit[0]); - for ( int i = 1; i < n; i++ ) { - for ( int j = 0; j < i; j++ ) { - cRz ( bit[j], bit[i], pi/(pow(2,i-j)) ); - } - H(bit[i]); - } -} - -int main () { - qbit reg[n]; - - for ( int i = 0; i < n; i++ ) { - PrepZ(reg[i], 0); - } - qft(reg); - for ( int i = 0; i < n; i++ ) { - H(reg[i]); - } - return 0; -} diff --git a/Algorithms/Shors/Shors_15/shors.n3.scaffassert b/Algorithms/Shors/Shors_15/shors.n3.scaffassert new file mode 100644 index 000000000..10f9202e7 --- /dev/null +++ b/Algorithms/Shors/Shors_15/shors.n3.scaffassert @@ -0,0 +1,116 @@ +// --- Scaffold Code for Shor’s Algorithm ---// +//Input x -> H -> f(x) -> iQFT +//f(x) = a^x mod N + +// minimal qubit implementation as described in +// Circuit for Shor’s algorithm using 2n+3 qubits +// Stephane Beauregard +// https://arxiv.org/abs/quant-ph/0205095v3 + +#include "../cMODMUL/cMODMUL.scaffold" +#include "../../../Library/cSWAP/cSWAP.scaffold" + +#define _a_val 7 // randomly chosen s.t. gcd(a,N) = 1 +#define _a_inv 13 // randomly chosen s.t. gcd(a,N) = 1 + +#define _upper_sz 3 // total iterations of phase estimation +#define _iteration 2 // current iteration index +#define _estimate 0 // tally of less significant bits + +module cUa ( + qbit ctrl, + const unsigned int width, + const unsigned int a, + const unsigned int a_inv, + qbit x[], + qbit b[], + unsigned int N, + qbit ancilla +) { + cMODMUL ( ctrl, width, a, x, b, N, ancilla ); + cSWAPs ( ctrl, width, x, b ); + ciMODMUL ( ctrl, width, a_inv, x, b, N, ancilla ); +} + +// b[width-1] is most significant bit +// b[0] is least significant bit +int main ( +) { + const unsigned int width = 5; // one extra than number of bits in N + const unsigned int N = 15; // number to be factorized + + // prepare upper register + qbit upper[_upper_sz]; + for (int i=0; i<_upper_sz; i++) { + PrepZ ( upper[i], 0 ); + H (upper[i]); + } + // the above is just a degenerate case of QFT + // QFT(_upper_size, upper); + // assert_superposition(upper,3); + + // prepare lower register to any value + const unsigned int lower_val = 1; + qbit lower[width]; + for ( int i=0; i>i)&1 ); + } + // assert_classical(lower,5,1); + + // prepare ancilla b to zero + const unsigned int b_val = 0; + qbit b[width]; + for ( int i=0; i>i)&1 ); + } + // assert_classical(b,5,0); + + qbit ancilla[1]; + PrepZ ( ancilla[0], 0 ); + + endian (width, lower); + endian (width, b); + + // cUa ( upper[0], width, _a_val, _a_inv, lower, b, N, ancilla[0] ); + cUa ( upper[0], width, 7, 13, lower, b, N, ancilla[0] ); + // cUa ( upper[0], width, 7, 12, lower, b, N, ancilla[0] ); + if (1<_upper_sz) cUa ( upper[1], width, 4, 4, lower, b, N, ancilla[0] ); + for (int i=2; i<_upper_sz; i++) { + cUa ( upper[i], width, 1, 1, lower, b, N, ancilla[0] ); + } + + endian (width, b); + // assert_classical(b,5,0); + + endian (width, lower); + // assert_classical(lower,5,13); // 1, 4, 7, 13 + + endian (_upper_sz, upper); + iQFT(_upper_sz, upper); + endian (_upper_sz, upper); + + // assert_classical(upper,3,6); // 0, 2, 4, 6 + assert_product(upper,3,b,5); + assert_product(upper,3,lower,5); + assert_product(b,5,lower,5); + + // or you can do iterative phase estimation instead of iQFT: + // for (int c=0; c<_iteration; c++) { + // if ((_estimate>>c)&1) { + // Rz ( upper[0], -M_PI/pow(2,_iteration-c) ); + // } + // } + // H(upper[0]); + + for (int i=0; i<_upper_sz; i++) { + // MeasZ(upper[i]); + } + + for ( int i=0; i + +#define width 5 // one extra than number of bits in N + +// b[width-1] is most significant bit +// b[0] is least significant bit +int main () { + + qbit ctrl[2]; + PrepZ (ctrl[0], 1); + PrepZ (ctrl[1], 1); + assert_classical(ctrl,2,3); + + const unsigned int a = 15; + + const unsigned int b_val = 15; + qbit b[width]; + for ( int i=0; i>i)&1 ); + } + assert_classical(b,5,15); + + endian (width, b); + assert_classical(b,5,30); + + QFT (width, b); + assert_superposition(b,5); + + cADD ( 2, ctrl[0], ctrl[1], width, a, b ); + assert_superposition(b,5); + + iQFT (width, b); + assert_classical(b,5,15); + + endian (width, b); + assert_classical(b,5,30); + + for ( int i=0; i + +#define width 5 // one extra than number of bits in N +#define N 37 // one extra than number of bits in N + +// b[_n-1] is most significant bit +// b[0] is least significant bit +int main () { + + qbit ctrl[2]; + PrepZ (ctrl[0], 1); + PrepZ (ctrl[1], 1); + assert_classical(ctrl,2,3); + + const unsigned int a = 22; + + const unsigned int b_val = 29; + qbit b[width]; + for ( int i=0; i>i)&1 ); + } + // 29 to binary + // 16+8+4+1 + // 11101 + assert_classical(b,5,29); + + qbit ancilla[1]; + PrepZ ( ancilla[0], 0 ); + + endian (width, b); + // 10111 to decimal + // 16+4+2+1 + // 23 + assert_classical(b,5,23); + + QFT (width, b); + assert_superposition(b,5); + + cMODADD ( ctrl[0], ctrl[1], width, a, b, N, ancilla[0] ); + assert_superposition(b,5); + + iQFT (width, b); + // 14 to binary + // 8+4+2 + // 01110 + assert_classical(b,5,14); + + endian (width, b); + assert_classical(b,5,14); + + for ( int i=0; i +#define width 5 // one extra than number of bits in N +#define N 15 + +// b[width-1] is most significant bit +// b[0] is least significant bit +int main () { + + qbit ctrl[1]; + PrepZ (ctrl[0], 1); + + const unsigned int a = 6; + + const unsigned int x_val = 6; + qbit x[width]; + for ( int i=0; i>i)&1 ); + } + assert_classical(x,5,6); + + const unsigned int b_val = 7; + qbit b[width]; + for ( int i=0; i>i)&1 ); + } + assert_classical(b,5,7); + + qbit ancilla[1]; + PrepZ ( ancilla[0], 0 ); + + endian (width, x); + endian (width, b); + cMODMUL ( ctrl[0], width, a, x, b, N, ancilla[0] ); + endian (width, b); + endian (width, x); + + assert_classical(b,5,13); + + for ( int i=0; i + +#define n 4 // problem size (log of database size) + +#define pi 3.141592653589793238462643383279502884197 + +// Module prototypes +void diffuse(qbit q[n], qbit x[n-1]); +void Sqr(qbit a[n], qbit b[n]); +void EQxMark(qbit b[n], qbit x[n-1]); + +module cRz ( qbit ctrl, qbit target, const double angle ) { + /* cRz identity matrix: + [ [ 1 0 0 0 ] + [ 0 1 0 0 ] + [ 0 0 e^(-i*angle/2) 0 ] + [ 0 0 0 e^(i*angle/2) ] ] + */ + + S(ctrl); + CNOT(ctrl, target); + Rz(target,-angle/2); + CNOT(ctrl, target); + Rz(target,angle/2); +} + +/*************************** + Diffusion Module +***************************/ +void diffuse(qbit q[n], qbit x[n-1]) { + int j; + // local qbit x[n-1]; // scratch register + // No local registers yet + // qbit x[n-1]; // scratch register + + // Hadamard applied to q + // No forall yet + // forall(j=0; j 0; j--) + Toffoli(x[j-1], q[j+1], x[j]); + Toffoli(q[1], q[0], x[0]); + + // Restore q + for(j = 0; j < n; j++) + X(q[j]); + + // Complete the diffusion + // No forall yet + // forall(j=0; j 0; j--) + Toffoli(x[j-1], b[j+1], x[j]); + Toffoli(b[1], b[0], x[0]); + + // Restore b + for(j = 0; j < n; j++) + if(j!=1) X(b[j]); +} + +/*********************************************** + Squaring a(x) = a0 + a1*x + ... + a_(n-1)*x^(n-1) + over the ring GF(2)[x] mod (x^n + x + 1). + + Result placed in the n-qubit register b +************************************************/ +void Sqr(qbit a[n], qbit b[n]) { + + int i; + int k; + + // Using forall indicates the CNOT's are independent + // So these instructions can be done in parallel + // No forall yet + // forall(i=0; i<=(n-1)/2; i++) { + for(i=0; i<=(n-1)/2; i++) { + k = 2*i; + CNOT(a[i], b[k]); + } + + // Would it make a difference to split this into two loops? + // Maybe since deleaving would be two forall loops! + // Or perhaps it is okay to just replace the for with a forall. + for(i=(n+1)/2; i Date: Tue, 18 Dec 2018 18:59:09 -0500 Subject: [PATCH 19/24] Ensemble size change --- scripts/simulation/simulate.bash | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/simulation/simulate.bash b/scripts/simulation/simulate.bash index 1395a78de..f86964650 100755 --- a/scripts/simulation/simulate.bash +++ b/scripts/simulation/simulate.bash @@ -20,7 +20,7 @@ rm *.breakpoint_*.bash # split into breakpoints BREAKPOINTS=$($SCAFFCC_PATH/scripts/simulation/scaffassert.py $1) -ENSEMBLE=8 +ENSEMBLE=16 # compilation without slurm: # for ((i=1;i<=$BREAKPOINTS;i++)) From f0efa01dd15c8cc8e18b3d83ccad8dc32e861eb5 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Tue, 18 Dec 2018 19:08:23 -0500 Subject: [PATCH 20/24] Changing testing script --- Algorithms/Shors/cADD/cADD_test.scaffassert | 40 ++++++++++++--------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/Algorithms/Shors/cADD/cADD_test.scaffassert b/Algorithms/Shors/cADD/cADD_test.scaffassert index b9c52548c..c2d22836e 100644 --- a/Algorithms/Shors/cADD/cADD_test.scaffassert +++ b/Algorithms/Shors/cADD/cADD_test.scaffassert @@ -4,41 +4,49 @@ // https://arxiv.org/abs/quant-ph/0205095v3 #include - #define width 5 // one extra than number of bits in N - // b[width-1] is most significant bit // b[0] is least significant bit int main () { qbit ctrl[2]; PrepZ (ctrl[0], 1); + H (ctrl[0]); PrepZ (ctrl[1], 1); - assert_classical(ctrl,2,3); - - const unsigned int a = 15; + H (ctrl[1]); - const unsigned int b_val = 15; + const unsigned int a = 13; + const unsigned int b_val = 12; qbit b[width]; for ( int i=0; i>i)&1 ); } - assert_classical(b,5,15); + // assert_classical(b,5,12); // 0+8+4+0+0 // 01100 endian (width, b); - assert_classical(b,5,30); - + // assert_classical(b,5,6); // 00110 // QFT (width, b); - assert_superposition(b,5); - - cADD ( 2, ctrl[0], ctrl[1], width, a, b ); - assert_superposition(b,5); - + // assert_superposition(b,5); + cADD ( 1, ctrl[0], ctrl[1], width, a, b ); + // assert_superposition(b,5); iQFT (width, b); - assert_classical(b,5,15); + // assert_classical(b,5,19); // 10011 // 16+2+1 + endian (width, b); + // assert_classical(b,5,25); // 25 = 16+8+0+0+1 + assert_product ( ctrl, 2, b, 5 ); + // now reverse compute it + endian (width, b); + // assert_classical(b,5,19); // 10011 // 16+2+1 + QFT (width, b); + // assert_superposition(b,5); + ciADD ( 1, ctrl[0], ctrl[1], width, a, b ); + // assert_superposition(b,5); + iQFT (width, b); + // assert_classical(b,5,6); // 00110 // endian (width, b); - assert_classical(b,5,30); + // assert_classical(b,5,12); // 0+8+4+0+0 // 01100 + assert_product ( ctrl, 2, b, 5 ); for ( int i=0; i Date: Tue, 18 Dec 2018 19:08:49 -0500 Subject: [PATCH 21/24] Changing testing script --- .../Shors/cMODMUL/cMODMUL_test.scaffassert | 20 +++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/Algorithms/Shors/cMODMUL/cMODMUL_test.scaffassert b/Algorithms/Shors/cMODMUL/cMODMUL_test.scaffassert index d770fd847..b9c086b42 100644 --- a/Algorithms/Shors/cMODMUL/cMODMUL_test.scaffassert +++ b/Algorithms/Shors/cMODMUL/cMODMUL_test.scaffassert @@ -13,22 +13,24 @@ int main () { qbit ctrl[1]; PrepZ (ctrl[0], 1); + H (ctrl[0]); - const unsigned int a = 6; + const unsigned int a = 7; + const unsigned int a_inv = 13; const unsigned int x_val = 6; qbit x[width]; for ( int i=0; i>i)&1 ); } - assert_classical(x,5,6); + // assert_classical(x,5,6); const unsigned int b_val = 7; qbit b[width]; for ( int i=0; i>i)&1 ); } - assert_classical(b,5,7); + // assert_classical(b,5,7); qbit ancilla[1]; PrepZ ( ancilla[0], 0 ); @@ -39,7 +41,17 @@ int main () { endian (width, b); endian (width, x); - assert_classical(b,5,13); + // assert_classical(b,5,4); + assert_product ( ctrl, 1, b, 5 ); + + endian (width, x); + endian (width, b); + cMODMUL ( ctrl[0], width, a_inv, x, b, N, ancilla[0] ); + endian (width, b); + endian (width, x); + + // assert_classical(b,5,7); + assert_product ( ctrl, 1, b, 5 ); for ( int i=0; i Date: Wed, 12 Jun 2019 16:09:40 -0400 Subject: [PATCH 22/24] Adding more comments on how this script splits the Scaffassert program into several Scaffold programs --- scripts/simulation/scaffassert.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/scripts/simulation/scaffassert.py b/scripts/simulation/scaffassert.py index 859f49ec0..5844a9f73 100755 --- a/scripts/simulation/scaffassert.py +++ b/scripts/simulation/scaffassert.py @@ -3,13 +3,17 @@ import sys, re, copy, os break_indx = 0 -breakpoints = [[]] +breakpoints = [[]] # the Scaffold code to be generated for each breakpoint +# open the Scaffassert file as input with open(sys.argv[1]) as source: + # for all the lines of code in the Scaffassert program... for line in source: + # decide if the line is any of the three types of assertions at all if re.findall("^assert_", line.strip()): + # syntax: assert_classical ( quantum_register, register_width, assertion_value ) if "assert_classical" in line: # start the next breakpoint breakpoints.append(copy.copy(breakpoints[break_indx])) @@ -17,12 +21,12 @@ param_string = re.findall("^assert_classical\s*\((.*?)\)\s*;", line.strip()) params = [param.strip() for param in param_string[0].split(',')] - # close out this breakpoint + # close out this breakpoint by doing measurement breakpoints[break_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) breakpoints[break_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) breakpoints[break_indx].append('}\n') - # generate script for checking the breakpoint + # generate script for checking this breakpoint after simulation and measurement breakpoint_name = os.path.splitext(os.path.basename(sys.argv[1]))[0] + '.breakpoint_' + str(break_indx+1) with open(breakpoint_name + '.bash', 'w') as outfile: outfile.write('$SCAFFCC_PATH/scripts/simulation/assert_classical.py {:s}.csv {:s} {:s} {:s}'.format( @@ -33,6 +37,7 @@ break_indx+=1 + # syntax: assert_superposition ( quantum_register, register_width ) elif "assert_superposition" in line: # start the next breakpoint breakpoints.append(copy.copy(breakpoints[break_indx])) @@ -40,12 +45,12 @@ param_string = re.findall("^assert_superposition\s*\((.*?)\)\s*;", line.strip()) params = [param.strip() for param in param_string[0].split(',')] - # close out this breakpoint + # close out this breakpoint by doing measurement breakpoints[break_indx].append('for ( int i=0 ; i<{:s} ; i++ )\n'.format(params[1])) breakpoints[break_indx].append(' MeasZ ( {:s}[i] );\n'.format(params[0])) breakpoints[break_indx].append('}\n') - # generate script for checking the breakpoint + # generate script for checking this breakpoint after simulation and measurement breakpoint_name = os.path.splitext(os.path.basename(sys.argv[1]))[0] + '.breakpoint_' + str(break_indx+1) with open(breakpoint_name + '.bash', 'w') as outfile: outfile.write('$SCAFFCC_PATH/scripts/simulation/assert_superposition.py {:s}.csv {:s} {:s}'.format( @@ -91,4 +96,5 @@ for line in breakpoints[break_indx]: outfile.write(line) +# tell all the subsequent scripts how many breakpoints we found here print (break_indx+1) From 90efb093186563640cf1f3da043c828702dce443 Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Wed, 12 Jun 2019 16:33:58 -0400 Subject: [PATCH 23/24] Assertion simulation script for simulations on one machine, with no parallelization. --- scripts/simulation/simulate_serial.bash | 59 +++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100755 scripts/simulation/simulate_serial.bash diff --git a/scripts/simulation/simulate_serial.bash b/scripts/simulation/simulate_serial.bash new file mode 100755 index 000000000..15dafd8df --- /dev/null +++ b/scripts/simulation/simulate_serial.bash @@ -0,0 +1,59 @@ +#!/bin/bash + +export SCAFFCC_PATH=/n/fs/qdb/ScaffCC +TEST_NAME=${1%.scaffassert} + +# CLEANUP +$SCAFFCC_PATH/scaffold.sh -c $1 +rm *.breakpoint_*.scaffold +rm *.breakpoint_*.ll +rm *.breakpoint_*.tmp +rm *.breakpoint_*.resources +rm *.breakpoint_*.qasmh +rm *.breakpoint_*.qasmf +rm *.breakpoint_*.qc +rm *.breakpoint_*.out +rm *.breakpoint_*.err +rm *.breakpoint_*.csv +rm *.breakpoint_*.bash + +# split into breakpoints +BREAKPOINTS=$($SCAFFCC_PATH/scripts/simulation/scaffassert.py $1) +# number of measurements we want for each breakpoint +ENSEMBLE=1 + +# get QX Simulator +if [ ! -e $SCAFFCC_PATH/../qx_simulator_linux_x86_64/qx_simulator_1.0.beta_linux_x86_64 ] +then + wget -nc -P $SCAFFCC_PATH/../ http://quantum-studio.net/qx_simulator_linux_x86_64.tar.gz + tar -xvf $SCAFFCC_PATH/../qx_simulator_linux_x86_64.tar.gz -C $SCAFFCC_PATH/../ +fi + +###################### +# Begin work section # +###################### + +for ((bpIndx=1;bpIndx<=$BREAKPOINTS;bpIndx++)) +do + # compile into QASM code + $SCAFFCC_PATH/scaffold.sh -s $TEST_NAME.breakpoint_${bpIndx}.scaffold + + # SIMULATE the whole ensemble + # map + for ((ensIndx=1;ensIndx<=$ENSEMBLE;ensIndx++)) + do + $SCAFFCC_PATH/../qx_simulator_linux_x86_64/qx_simulator_1.0.beta_linux_x86_64 $TEST_NAME.breakpoint_${bpIndx}.qc + done + + # reduce + $SCAFFCC_PATH/scripts/simulation/register_value_csv.py \ + $TEST_NAME.breakpoint_${bpIndx}.qc \ + $TEST_NAME.breakpoint_${bpIndx}.csv + + bash $TEST_NAME.breakpoint_${bpIndx}.bash + + # scripts/simulation/assert_uniform.py ${source%.scaffold}.csv 4 + # scripts/simulation/assert_integer.py ${source%.scaffold}.csv 5 30 + # scripts/simulation/assert_product.py ${source%.scaffold}.csv + +done From 0a98828fd2e7bf290f5fbd0c02e4b88bf367af1a Mon Sep 17 00:00:00 2001 From: Yipeng Huang Date: Thu, 13 Jun 2019 13:42:57 -0400 Subject: [PATCH 24/24] Improved comments on simulation scripts. --- scripts/simulation/assert_classical.py | 10 +++---- scripts/simulation/assert_product.py | 22 +++++++------- scripts/simulation/assert_superposition.py | 5 +++- scripts/simulation/fidelity.py | 2 +- scripts/simulation/register_value_csv.py | 29 ++++++++++++------- .../{simulate.bash => simulate_parallel.bash} | 2 +- scripts/simulation/simulate_serial.bash | 13 +++++---- scripts/simulation/trace_distance.py | 2 +- 8 files changed, 49 insertions(+), 36 deletions(-) rename scripts/simulation/{simulate.bash => simulate_parallel.bash} (99%) diff --git a/scripts/simulation/assert_classical.py b/scripts/simulation/assert_classical.py index 6b77f293b..dae97cba0 100755 --- a/scripts/simulation/assert_classical.py +++ b/scripts/simulation/assert_classical.py @@ -14,9 +14,9 @@ for val in range(1<