#include "matrix.h"
#include "supplement.h"
+#include "supplement2.h"
#include "sparsify.h"
/****************************************************************************
// Note: Make sure you seed the random number generator before calling this!
// i.e. srand((unsigned)time(NULL));
-int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling) {
+int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling, double samplePrefactor) {
int i, j, k, l;
double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
//printf("L1Norm squared (= stabilizer extent) = %lf\n", pow(L1Norm,2.0));
+ // fraction of optimal T supplemental states
+ //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/((double)T);
+ //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/(4.0+pow(delta,2.0))/((double)T);
+ double alpha;
+
+ // fraction of T supplemental states that we actually need for optimal scaling
+ // if first term is kept from Seddon et al.
+ //double fractionSupplement = alpha*0.25*delta*pow(L1Norm,2.0)/((double)T);
+ // if first term isn't kept
+ double fractionSupplement;
+
printf("delta = %lf\n", delta);
- //printf("L1Norm = %lf\n", L1Norm);
+ printf("L1Norm = %lf\n", L1Norm);
double gamma;
int increment;
- if(coherentSampling) {
+ if(coherentSampling!=0) {
+ if(coherentSampling==2) {
+ alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(2*T-1));
+ fractionSupplement = alpha*1.0;
+ *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(2*T-1))+1)))*(floor(fractionSupplement*(2*T-1))+1);
+ increment = floor(fractionSupplement*(2*T-1))+1;
+ } else if(coherentSampling==1) {
+ alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(T));
+ fractionSupplement = alpha*1.0;
+ *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(T))+1)))*(floor(fractionSupplement*(T))+1);
+ increment = floor(fractionSupplement*(T))+1;
+ }
+
printf("Performing coherent sampling!\n");
- gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
- // since we will be taking a multiple of (T+1) samples, we round numStabStates up to the nearest multiple of (T+1)
- *numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
- increment = T+1;
- } else {
+ printf("Using an optimal fraction of T supplemental states!\n");
+ printf("Optimal fraction = %lf\n", fractionSupplement);
+ // since we will be taking a multiple of (fractionSupplement*T+1) samples, we round numStabStates up to the nearest multiple of (fractionSupplement*T+1)
+ // if first term is kept from Seddon et al.
+ //*numStabStates = ceil((8.0/(8.0*delta-sqrt(2.0)*sqrt(16.0-16.0*alpha+pow(alpha,3.0)*pow(delta,4.0))))*pow(L1Norm,2.0)/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1);
+ // if first term isn't kept
+ //*numStabStates = ceil(0.0005*(sqrt(2.0)*pow(L1Norm,2.0)/delta-0.5*fractionSupplement*T/pow(delta,2.0)+pow(fractionSupplement*T/L1Norm,2.0)/8.0/sqrt(2.0)/pow(delta,3.0))/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1);
+ //*numStabStates = ceil(sqrt(2.0)*(2.0+pow(delta,2.0))*pow(L1Norm,2.0)/delta/(4.0+pow(delta,2.0))/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1);
+
+ printf("increment= %d\n", increment);
+ //printf("prebinned numStabStates=%d\n", (int)ceil((8.0/(8.0*delta-sqrt(2.0)*sqrt(16.0-16.0*alpha+pow(alpha,3.0)*pow(delta,4.0))))*pow(L1Norm,2.0)));
+ //printf("prebinned numStabStates=%d\n", (int)ceil(0.0005*(sqrt(2.0)*pow(L1Norm,2.0)/delta-0.5*fractionSupplement*T/pow(delta,2.0)+pow(fractionSupplement*T/L1Norm,2.0)/8.0/sqrt(2.0)/pow(delta,3.0))));
+ //printf("prebinned numStabStates=%d\n", (int)ceil(sqrt(2.0)*(2.0+pow(delta,2.0))*pow(L1Norm,2.0)/delta/(4.0+pow(delta,2.0))));
+ printf("prebinned numStabStates=%d\n", (int)ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta));
+ fflush(stdout);
+ if(fractionSupplement > 1.0) {
+ printf("Too many supplemental states requested!\n");
+ printf("Note: the correlated sampling is currently only implemented for a maximum of (2t-1) supplemental correlated states!\n");
+ return 1;
+ }
+ } else if(coherentSampling==0) {
//gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
+ // SPARSIFY
//*numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
- gamma = 1.0;
- *numStabStates = ceil((pow(L1Norm,2.0)-gamma)/pow(delta,2.0));
+ // if first term is kept from Seddon et al.
+ //*numStabStates = ceil((2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0));
+ // if first term isn't kept
+ *numStabStates = ceil(samplePrefactor*(1.0*2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0));
+ // multiply the factor of 2.0 by 1.0 to get the prior state-of-the-art
+ // multiply the factor of 2.0 by 0.0 to get the tighter bound on the prior state-of-the-art
increment = 1;
+ } else {
+ printf("Error: Unknown value for \"coherent sampling\" flag");
+ return 1;
}
}
(*stabStateIndices)[i] = index;
//printf("first index = %d\n", index);
- if(coherentSampling)
- for(k=1; k<=T; k++)
+ if(coherentSampling==2)
+ //for(k=1; k<=T; k++)
+ //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index
+ supplement2(&(*stabStateIndices)[i], T, fractionSupplement);
+ else if(coherentSampling==1)
+ //for(k=1; k<=T; k++)
//(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index
- supplement(&(*stabStateIndices)[i], T);
+ supplement(&(*stabStateIndices)[i], T, fractionSupplement);
break;
}
//printf("%lf\n", sum);
projects / weak_simulation_stab_extent.git / blobdiff