final v1.0 files

[weak_simulation_stab_extent.git] / weaksim.c

diff --git a/weaksim.c b/weaksim.c
index 8ec622248094818ab474bac6557e69ab730d7dd4..b2c93f4284a8d26bd419472e66c95f770ccff168 100644 (file)
--- a/weaksim.c
+++ b/weaksim.c
@@ -11,6 +11,7 @@
 #include "measurepauli.h"
 #include "innerproduct.h"
 #include "sparsify.h"
+#include "omp.h"
 
 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
 
@@ -19,14 +20,18 @@ int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, i
 int main(int argc, char *argv[])
 {
 
-  if(argc != 4) {
-    printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=yes)\"\n");
+  if(argc != 5 && argc != 6) {
+    printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=t; 2=2t-1; 3=xi^3t/2)\" \"seed\" \"(optional) sample number prefactor\"\n");
     exit(0);
   }
 
   double additiveErrorDelta = atof(argv[1]);        // additive error delta
   double phi = PI*atof(argv[2]);//PI/4.0; // PI/4.0 is the T gate magic state
-  int coherentSampling = atoi(argv[3]);           // perform coherent sampling (true=1 or false=0)
+  int coherentSampling = atoi(argv[3]);           // perform coherent sampling (false=0; linear t=1; linear 2t-1=2; exponential xi^3t/2=3)
+  int seed = atoi(argv[4]); // random number seed
+  double samplePrefactor = 1.0;
+  if(argc == 6)
+    samplePrefactor = atof(argv[5]); // factor that multiplies k, the number of samples used
 
   int N;              // number of qubits
   scanf("%d", &N);
@@ -57,16 +62,17 @@ int main(int argc, char *argv[])
   long* stabStateIndices;
   int numStabStates;
 
-  srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+  //srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+  srand((unsigned)seed); // seeding the random number generator for sparsify()
 
-  if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling))
+  if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling, samplePrefactor))
     return 1;
 
-  //printf("checking: numStabStateIndices:\n");
-  //for(i=0; i<numStabStates; i++)
-  //  printf("%d ", stabStateIndices[i]);
-  //printf("\n");
-  //fflush(stdout);
+  /* printf("checking: numStabStateIndices:\n"); */
+  /* for(i=0; i<numStabStates; i++) */
+  /*   printf("%ld ", stabStateIndices[i]); */
+  /* printf("\n"); */
+  /* fflush(stdout); */
 
   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
   double complex Gamma[(int)numStabStates]; // prefactor in front of resultant state
@@ -248,6 +254,34 @@ int main(int argc, char *argv[])
 
   memcpy(origQ, Q, numStabStates*sizeof(int));
 
+  double complex amplitude = 0.0 + 0.0*I;
+  double complex lastamplitude = 0.0 + 0.0*I;
+
+  double complex probability = 1.0 + 0.0*I;
+
+  // the first measurement should be the identity
+  omega = 0;
+  for(i=0; i<N; i++) {
+    alpha[i] = 1; beta[i] = 0; gamma[i] = 0; delta[i] = 0;
+  }
+
+  for(j=0; j<numStabStates; j++) { // the kets
+    Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
+  }
+
+  double complex newamplitude;
+  
+  #pragma omp parallel
+  {
+  #pragma omp for schedule(static) collapse(2)
+  for(i=0; i<numStabStates; i++) { // the bras
+    for(j=0; j<numStabStates; j++) {
+      newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
+      amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+    }
+  }
+  }
+
   while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
   
     Paulicounter++;
@@ -273,14 +307,26 @@ int main(int argc, char *argv[])
 
     }
 
-  }
-
-  double complex amplitude = 0.0 + 0.0*I;
-  for(i=0; i<numStabStates; i++) { // the bras
-    for(j=0; j<numStabStates; j++) {
-      double complex newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
-      amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+    lastamplitude = amplitude;
+    amplitude = 0.0 + 0.0*I;
+    #pragma omp parallel
+    {
+    #pragma omp for schedule(static) collapse(2)
+    for(i=0; i<numStabStates; i++) { // the bras
+      for(j=0; j<numStabStates; j++) {
+       newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
+       amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+      }
+    }
     }
+    printf("lastamplitude: %lf %c %lf I\n", cabs(creal(lastamplitude)), cimag(lastamplitude+0.00000001)>0?'+':'-' , cabs(cimag(lastamplitude)));
+    printf("amplitude: %lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+
+    probability *= amplitude/lastamplitude;
+    //amplitude = amplitude/lastamplitude; // for NORMALIZE-SPARSIFY you normalize the algorithm after every Pauli measurement
+
+    printf("probability = %lf\n", cabs(probability));
+    printf("!\n");
   }
 
   //printf("numStabStates=%d\n", numStabStates);
@@ -295,6 +341,8 @@ int main(int argc, char *argv[])
   printf("\nabs(amplitude):\n");
   printf("%lf\n", cabs(amplitude));
 
+  printf("probability = %lf\n", cabs(probability));
+
   return 0;
 
 }