key_rate_utils.cpp

namespace key_rate_utils{
	//general stuff
	double h_entropy(double x){
		//Classical Shannon entropy of a variable x
		double ret;
		if (fabs(x) < .000000001 || fabs(x) > 1-.0000001)
			ret = 0.0;
		else
			ret = -1*x*log2(x) - (1-x)*log2(1-x);
		return ret;
	}

	//thm 1 general:
	double s_entropy(vector<double> E1s, vector<double> E2s, vector<double> Ls, int CAD_level){
		//Computes the summation of terms in thm1 (given Es adjsted for cad)
		double ret = 0;
		for (int i = 0; i < E1s.size(); i++)
			ret+= calc_s_term(E1s[i], E2s[i], Ls[i], CAD_level);
		return ret;
	}

	double calc_s_term(E1s[i])(double E1, double E2, double Lx, int CAD_level){
		//Computes the terms of the summation in thm1
		double sumE, bin1, bin2;
		sumE = E1+E2;
		bin1 = h_entrpoy(E1/sumE);
		bin2 = h_entropy(calc_lambdax(E1,E2,big_Lx, CAD_level));
		return sumE*(bin1-bin2)/2; //TODO: paramaterize thihs 2 probably??
	}

	inline double calc_lambdax(double E1N, double E2N, double big_Lx, int CAD_level){
		//computes lambda_x according to thm1 given elements 
		//of E1s and E2s and Ls (adjusted for CAD)
    	return (1+sqrt(pow(E1N - E2N, 2) + 4*pow(big_Lx, 2*CAD_level))/(E1N+E2N))/2;
  	}

  	int H_AB(vector<double> Pij_key){
  		//Returns conditional entropy of Alice and Bobs systems given array of 
  		//prob A=0, B=0, in binary order adjusted for CAD
  		ret = 0;  
  		for(auto p : Pij_key)
  			ret += -1*p*log2(p);
  		
  		return ret - h_entropy(Pij_key[0]+Pij_key[2]);
  	}
  	  
  	void generateNoiseStatistics(bool SYM_FLAG double Q=.11){
	    /*
	     * Returns a vector of the form {NoiseA-B, Noise ARA, Noise ABA}. If sym is
	     * true it assumes it is a symmetric channel and creates statistis using the
	     * noise values Q (noise Z basis) and Qx (noise in Z basis), else uses the
	     * utils package
	     */
	    if(SYM_FLAG){
	      //do nothing currently. Can we get matrices that are less random?
	      double Qx = IND_FLAG ? 2*Q*(1-Q) : Q;
	      AtoB["P00"]  = 1- Q;
	      AtoB["P11"]  = 1- Q;
	      AtoB["P10"]  = Q;
	      AtoB["P01"]  = Q;
	      AtoB["PP0"]  = .5;
	      AtoB["PP1"]  = .5;
	      AtoB["P0P"]  = .5
	      AtoB["P1P"]  = .5

	      AtoAproj["P000"]  = 1 - Q;
	      AtoAproj["P111"]  = 1 - Q;
	      AtoAproj["P011"]  = 1 - Q;
	      AtoAproj["P100"]  = 1 - Q;
	      AtoAproj["P001"]  = Q;
	      AtoAproj["P110"]  = Q;
	      AtoAproj["P010"]  = Q;
	      AtoAproj["P101"]  = Q;
	      AtoAproj["PP0P"]  = .5; 
	      AtoAproj["P00P"]  = .5;
	      AtoAproj["P10P"]  = .5;
	      AtoAproj["P11P"]  = .5;
	      AtoAproj["PP1P"]  = .5;
	      AtoAproj["P01P"]  = .5;

	      //Bug Doc about the Qx here, doesn't make sense to me still
	      AtoAref["P0R0"]   = (1-Q)*(1-Q)+Q*Q; // 1-Qx;
	      AtoAref["P0R1"]   = 2*Q*(1-Q); //Qx;
	      AtoAref["P1R0"]   = 2*Q*(1-Q); //Qx;
	      AtoAref["P1R1"]   = (1-Q)*(1-Q)+Q*Q; // 1-Qx;-Qx;

	      AtoAref["P0RP"]  = .5;
	      AtoAref["P1RP"]  = .5;
	      AtoAref["PPRM"]  = Qx;

	      trueEs["E1"] = NAN;
	      trueEs["E2"] = NAN;
	      trueEs["E3"] = NAN;
	      trueEs["E4"] = NAN;


	    }else{

	      complex<double> one0(1.0,0.0); 
	      complex<double> oneoRoot2(1/sqrt(2),0.0);
	      complex<double> noneoRoot2(-1.0/sqrt(2),0.0);
	      complex<double> ioRoot2(0.0,1.0/sqrt(2));
	      complex<double> nioRoot2(0.0,-1.0/sqrt(2));

	      QMat zeroM(2,1);
	      zeroM.setElement(0,0,one0);
	      Qubit zero("0", zeroM);

	      QMat oneM(2,1);
	      oneM.setElement(1,0,one0);
	      Qubit one("1", oneM);

	      QMat plusM(2,1);
	      plusM.setElement(0,0,oneoRoot2);
	      plusM.setElement(1,0,oneoRoot2);
	      Qubit plus("P", plusM);

	      QMat minusM(2,1);
	      minusM.setElement(0,0,oneoRoot2);
	      minusM.setElement(1,0,noneoRoot2);
	      Qubit minus("M", minusM);

	      QMat zeroyM(2,1);
	      zeroyM.setElement(0,0,oneoRoot2);
	      zeroyM.setElement(1,0,ioRoot2);
	      Qubit zeroy("ZY", zeroyM);

	      QMat oneyM(2,1);
	      oneyM.setElement(0,0,oneoRoot2);
	      oneyM.setElement(1,0,nioRoot2);
	      Qubit oney("OY", oneyM);

	      vector<Qubit> iSet = {zero, one, plus};
	      vector<Qubit> jSet = {zero, one};
	      vector<Qubit> kSet = {zero, one, plus, minus};

	      QMat Uf = random_unitary_matrix(32);
	      QMat Ur = random_unitary_matrix(32);

	      map<string, double> AtoBStats = noiseStatsAtoB(Uf, iSet, jSet);
	      map<string, double> AtoBreftoAStats = noiseStatsAtoBreftoA(Uf, Ur, iSet, kSet);
	      map<string, double> AtoBprojtoAStats = noiseStatsAtoBprojtoA(Uf,Ur, iSet, jSet, kSet);

	      map<string, double> values = true_value(Uf, Ur, zero, one);

	      AtoB = AtoBStats;
	      AtoAref = AtoBreftoAStats;
	      AtoAproj = AtoBprojtoAStats;
	      trueEs = values;
	    } 
	}
}