diff --git a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold index 4e4da1161..1fa7f5357 100644 --- a/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold +++ b/Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold @@ -1,7 +1,12 @@ #include +#include <../../Library/cRz/cRz.scaffold> #define M 4 //number of bottom register qubits -#define TROTTER 16 //number of bottom register qubits +#define TROTTER 8 + +#define _upper_sz 4 +#define _iteration 0 +#define _estimate 0 enum spin {alpha, beta}; const enum spin sigma[M] = { @@ -16,54 +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.000000; -// 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.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.0,0.0,0.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.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.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}} } @@ -77,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 @@ -139,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 @@ -158,19 +152,20 @@ int main () { qbit bottomregister[M]; - // 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); + cRz(topregister[0], bottomregister[t], eta*time_step); CNOT(bottomregister[c], bottomregister[t]); } @@ -234,24 +234,19 @@ 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 ); } } - // 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]); + } 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 -#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 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,69 @@ 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 +#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); + H (ctrl[1]); + + 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,12); // 0+8+4+0+0 // 01100 + + endian (width, b); + // assert_classical(b,5,6); // 00110 // + QFT (width, b); + // assert_superposition(b,5); + cADD ( 1, ctrl[0], ctrl[1], width, a, b ); + // assert_superposition(b,5); + iQFT (width, b); + // 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,12); // 0+8+4+0+0 // 01100 + assert_product ( ctrl, 2, b, 5 ); + + 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); + H (ctrl[0]); + + 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); + + 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,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 + +#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 + +#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= 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]); } diff --git a/Library/QFT/QFT_test.scaffassert b/Library/QFT/QFT_test.scaffassert new file mode 100644 index 000000000..448e9557f --- /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, 4, 5 ); + + QFT(width, reg); + + assert_superposition ( reg, 4 ); + + iQFT(width, reg); + + assert_classical ( reg, 4, 5 ); + + for ( int i=0; i 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/register_value_csv.py b/scripts/simulation/register_value_csv.py new file mode 100755 index 000000000..95654c0e3 --- /dev/null +++ b/scripts/simulation/register_value_csv.py @@ -0,0 +1,90 @@ +#!/usr/bin/env python +import sys, re, cmath, os, glob +from csv import DictWriter + +regex = re.compile("([a-zA-Z]+)([0-9]+)") +outcome = { 'probability':0, 'polar_r':0, 'polar_phi':0, 'rect_real':0, 'rect_imag':0 } +reg_map = [] + +# open QX source file and parse it for the register names +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 = ") +# print (reg_map) + +with open(sys.argv[2], 'a+') as csvfile: # open for appending + + writer = DictWriter(csvfile, fieldnames=list(outcome.keys())) + if not csvfile.read(): + writer.writeheader() + + # get the printouts from all the trials + filenames = glob.glob(os.path.splitext(sys.argv[1])[0]+'.trial_*.out') + print ("filenames = ") + print (filenames) + for filename in filenames: + with open(filename) as file: + + # Parse QX Simulator output format + state_lines = False # whether these lines are state amplitude lines + first_basis = True # whether this is the first state such that we need to find the global phase + phase_global = 0.0 + for line in file: + # print ("line = ") + # 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): + + 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'] + 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_map) = ") + # print (zip(strings[2][::-1],reg_map)) + for elem in zip(strings[2][::-1],reg_map): + outcome[elem[1][0]] += int(elem[0]) << elem[1][1] + + out_string = 'probability: {:f} polar: ({:f},{:f}) rect: ({:f},{:f}) |{:s}>'.format( + outcome['probability'], + outcome['polar_r'], + outcome['polar_phi'], + outcome['rect_real'], + outcome['rect_imag'], + strings[2], + ) + print("out_string = ") + 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 100644 index da1a3ede7..000000000 --- a/scripts/simulation/remove_global_phase.py +++ /dev/null @@ -1,48 +0,0 @@ -import sys -import re -from math import sqrt -from math import atan2 -from math import cos -from math import sin - -with open(sys.argv[1]) as infile, open(sys.argv[2], 'w') as outfile: - - state_lines = False - first_basis = True - phase_global = 0.0 - - # Parse QX Simulator output format - for line in infile: - if line.strip() == "--------------[quantum state]--------------": - print("--------------[quantum state]--------------") - state_lines = True - elif line.strip() == "-------------------------------------------": - print("-------------------------------------------") - state_lines = False - elif state_lines: - strings = re.findall(r"[-+]?\d*\.\d+|\d+", line) - - # grab real and imaginary parts of basis state amplitude - real = float(strings[0]) - imag = float(strings[1]) - magnitude = sqrt(real*real + imag*imag) - phase_local = atan2(imag,real) - - # if this is the first basis state then this is the global phase to factor out - if first_basis: - phase_global = phase_local - first_basis = False - - real_no_global = magnitude*cos(phase_local-phase_global) - imag_no_global = magnitude*sin(phase_local-phase_global) - out_string = 'mag {:f} ang {:f} ({:f},{:f}) |{:s}> +'.format( - magnitude, - phase_local-phase_global, - real_no_global, - imag_no_global, - strings[2] - ) - - print(out_string) - outfile.write(out_string) - outfile.write('\n') diff --git a/scripts/simulation/scaffassert.py b/scripts/simulation/scaffassert.py new file mode 100755 index 000000000..5844a9f73 --- /dev/null +++ b/scripts/simulation/scaffassert.py @@ -0,0 +1,100 @@ +#!/usr/bin/env python + +import sys, re, copy, os + +break_indx = 0 +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])) + + param_string = re.findall("^assert_classical\s*\((.*?)\)\s*;", line.strip()) + params = [param.strip() for param in param_string[0].split(',')] + + # 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 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( + breakpoint_name, + params[0], + params[1], + params[2])) + + 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])) + + param_string = re.findall("^assert_superposition\s*\((.*?)\)\s*;", line.strip()) + params = [param.strip() for param in param_string[0].split(',')] + + # 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 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( + 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: + breakpoints[break_indx].append(line) + +# 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 breakpoints[break_indx]: + outfile.write(line) + +# tell all the subsequent scripts how many breakpoints we found here +print (break_indx+1) diff --git a/scripts/simulation/simulate.sh b/scripts/simulation/simulate.sh deleted file mode 100755 index e5f006474..000000000 --- a/scripts/simulation/simulate.sh +++ /dev/null @@ -1,15 +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 - -./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 diff --git a/scripts/simulation/simulate_parallel.bash b/scripts/simulation/simulate_parallel.bash new file mode 100755 index 000000000..5dcec2f3f --- /dev/null +++ b/scripts/simulation/simulate_parallel.bash @@ -0,0 +1,110 @@ +#!/bin/bash + +export 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 *.breakpoint_*.out +rm *.breakpoint_*.err +rm *.breakpoint_*.csv +rm *.breakpoint_*.bash + +# split into breakpoints +BREAKPOINTS=$($SCAFFCC_PATH/scripts/simulation/scaffassert.py $1) +ENSEMBLE=1 + +# 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 + +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 + +# 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_serial.bash b/scripts/simulation/simulate_serial.bash new file mode 100755 index 000000000..847de23e6 --- /dev/null +++ b/scripts/simulation/simulate_serial.bash @@ -0,0 +1,62 @@ +#!/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=8 + +# 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 the ensemble of measurements into several runs of the simulator + 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 > \ + $TEST_NAME.breakpoint_${bpIndx}.trial_${ensIndx}.out + done + + # reduce the measurement results into a summary CSV file + $SCAFFCC_PATH/scripts/simulation/register_value_csv.py \ + $TEST_NAME.breakpoint_${bpIndx}.qc \ + $TEST_NAME.breakpoint_${bpIndx}.csv + + # do statistical tests on the summary CSV files + # check for the assertions as recorded in the bash file + 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 diff --git a/scripts/simulation/trace_distance.py b/scripts/simulation/trace_distance.py index b386f4657..d8446e4f3 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 @@ -27,4 +28,4 @@ imag_diff*imag_diff ) - print trace_distance + print (trace_distance)