key_rate_utilsb92.cpp

#include "key_rate_utilsb92.hpp"

//namespace key_rate_utils{
	//general stuff
	using namespace Algebra;
	double safe_log2(double x){
		if (fabs(x) < .00000000001)
			return 0;
		else 
			return log2(x);
	}
	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*safe_log2(x) - (1-x)*safe_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+= (E1s[i] + E1s[i]) * calc_s_term(E1s[i], E2s[i], Ls[i], CAD_level);
		double N =0;
		for (auto& n : E1s)
    		N += n;
    	for (auto& n : E2s)
    		N += n;
		return ret/N;
	}

	double calc_s_term(double E1, double E2, double Lx, int CAD_level){
		//Computes the terms of the summation in thm1
		if (E1 <= 0 or E2 <= 0){
			return 0;
		}
		double sumE, bin1, bin2;
		sumE = E1+E2;
		bin1 = h_entropy(E1/sumE);
		bin2 = h_entropy(calc_lambdax(E1,E2,Lx, CAD_level));
		double ret = bin1-bin2; 
		return ret;
	}

	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)
		//cout << "E1N" << E1N << endl;
		//cout << "E2N" << E2N<< endl;
		//cout << "Lx" << big_Lx<< endl;
		//cout << "CAD_level" << CAD_level<< endl;
		return .5 + sqrt(
						pow(E1N - E2N,2) +
						4*pow(big_Lx, 2*CAD_level)
					)/(2*E1N+2*E2N);

  	}

  	double 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
  		double ret = 0;  
  		for(auto p : Pij_key)
  			ret += -1*p*safe_log2(p);
  			//cout << p <<"*safe_log2" << p <<"="<<p*safe_log2(p) << endl;
  			
  		//cout << ret << endl;
  		//cout << "ret " << ret << endl;
  		//cout << "B " << Pij_key[0]+Pij_key[2] << endl;
  		//cout << ret - h_entropy(Pij_key[0]+Pij_key[2]) << endl;
  		return ret - h_entropy(Pij_key[0]+Pij_key[2]);
  	}
  	  
  	void generateNoiseStatistics(bool SYM_FLAG, double alph, double Q){
	    /*
	     * 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
	      AtoB["P10"] = Q;
	      AtoB["P01"] = Q;
	      AtoB["P11"] = 1 - AtoB["P10"];
	      AtoB["P00"] = 1 - AtoB["P01"];
	      // fix these they should be different
	      AtoB["PA0"]  = Q + (1-2*Q)*alph*alph;
	      AtoB["P0A"]  = Q + (1-2*Q)*alph*alph;
	      AtoB["P1A"]  = Q + (1-2*Q)*(1-alph*alph);
	      AtoB["PAB"]  = Q;

	      //key rounds
	      AtoB["P0B"] = 1 - AtoB["P0A"] ;
	      AtoB["PA1"] = 1 - AtoB["PA0"];


	    }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;
	    } 
	}
//}