% ERROR CORRECTION FOR QUANTUM ANNEALING ARCHITECTURE
% requires: nothing
% author: Fernando Pastawski
% license: GPL3
% If this code is usefull to you, please cite the accompanying article

% This script benchmarks the perfomance of decoding information from the 
% quantum annealing architecture proposed by Lechner et al. assuming i.i.d. 
% measurement errors.
% The decoding function is included inline and may be easily extracted and 
% used independently for the actual decoding of information from an experimets.


% nruns = number of runs to average over to obtain perfomace estimate
nruns=50; % set to 5000 for the corresponding article which took 8 hours on a desktop
% Different number of logical variables
NNvalues = [2,4,6,8,10,12,14,16,18,20, 30, 40];
%NNvalues = [200];
% Different values for the physical bit flip probability
pvalues = [0.05, 0.07, 0.1, 0.15, 0.2,0.3, 0.4];
%pvalues = [0.4];

fail =zeros(length(NNvalues),length(pvalues));

% Sweep over the different values NN chosen for the number of logical variables
for x =1:length(NNvalues);
    NN = NNvalues(x);
    
    % Sweep over the different values p chosen for the measurement error probability
    for w = 1:length(pvalues);
        p = pvalues(w);
        
        % Number of iterations made to accumulate statistics
        for l=1:nruns;
            
            % mu0 = ones(NN) is the correct parity data
            % Due to the symmetry and linearity of the error model it is
            % sufficient to analyze the decoder for this input.
            % Generate a random NN by NN array with all ones except for randomly chosen
            % positions which contain -1 with probability p
            mu = 2*(rand(NN) >= p)-1; % mu contains our current guess -1s are errors but we don't know that
            per = p*(ones(NN)-eye(NN)); % per contains the likelihood for each of the NN*NN entries to be incorrect 
            
            % Make distribution symmetric and acertain that self-parity is a certainty;
            for i=1:NN;
                mu(i,i)=1; % Diagonal elements correspond to no error
                for j=i+1:NN
                    mu(j,i)=mu(i,j); % Only half the indices are used (i<j)
                end;
            end;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % HERE BEGGINS BELIEF PROPAGATION DECODING ALGORITH
            % functions mu = BPdecoder(N, p, nits, mu)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            % nits = Number of belief propagation iterations
            nits=5;
            for it=1:nits %Number of iterations;
                Pm1  = (per-1/2).* mu + 1/2; % Likelihood of being -1  ( per or 1-per )
                Pp1 = (mu+1)/2 - per.*(mu);  % Likelihood of being +1  ( per or 1-per )
                for i=1:NN;
                    for j=i+1:NN % j > i
                        for k = [1:i-1,i+1:j-1,j+1:NN] % k ~= i and k ~= j
                            % Calculate the likelihood (unnormalized) of (i,j) being pm1
                            % given neighborhood and associated trust value
                            Pm1(i,j) = Pm1(i,j) * ... % Previous likelihood of -1 at the variable
                                ( (mu(i,k)*mu(k,j)==-1)*((1-per(i,k))*(1-per(k,j))+per(i,k)*per(k,j)) ... % Probability of a consistent check (0 or 2 errors)
                                +  (mu(i,k)*mu(k,j)== 1)*((1-per(i,k))*per(k,j)+per(i,k)*(1-per(k,j))) ); % Probability of an inconsistent check (1 error)
                            Pp1(i,j) = Pp1(i,j) * ...% Previous likelihood of +1 at the variable 
                                ( (mu(i,k)*mu(k,j)==1)*((1-per(i,k))*(1-per(k,j))+per(i,k)*per(k,j)) ... % Probability of a consistent check ( 0 or 2 errors)
                                +  (mu(i,k)*mu(k,j)== -1)*((1-per(i,k))*per(k,j)+per(i,k)*(1-per(k,j))));% Probability of an inconsistent check (1 error)
                        end
                    end
                end;
                % Update our current best guess in mu and 
                % calculate normalized marginal distributions
                for i=1:NN;
                    for j=i+1:NN;
                        mu(i,j) =  2*( Pp1(i,j) >= Pm1(i,j) )-1; % Recalculate our favorite value
                        mu(j,i) = mu(i,j); 
                        % Recalculate relative probabilities use eps to
                        % regulate overconfidence
                        per(i,j) = max(min( Pm1(i,j), Pp1(i,j))/(Pp1(i,j)+Pm1(i,j)),eps); 
                        per(j,i) = per(i,j);
                    end;
                end;
            end;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%
            % end  
            % Here ends the function estimating the most likely mu
            % through belief propagations.
            % mu(NN,j) contains all the parities with respect to b(NN)
            %%%%%%%%%%%%%%%%%%%%%%%%
            
            % Without danger of impressision we have assumed for
            % benchmarking to start with the all +1 state
            % Hence it is a failure when this is not recovered.
            if (sum(mu(NN,:))~= NN)
                fail(x,w) = fail(x,w) + 1;
            end;
        end;
    end;
end;

fail=fail/nruns;
toc

% Generate a plot with one curve per value of N used
P = pvalues % Abreviation, just for brevity in next line
figure
plot(P,fail(1,:),'r',P,fail(2,:),'b--o',P,fail(3,:),'r--', ...
     P,fail(4,:),'b',P,fail(5,:),'r--o',P,fail(6,:),'b--', ...
     P,fail(7,:),'r',P,fail(8,:),'b--o',P,fail(9,:),'r--', ...
     P,fail(10,:),'b', P,fail(10,:),'r--o',P,fail(12,:),'b--')

% Generate a plot with one curve per measurement error probability analyzed
NN = NNvalues % Abreviation overwriting a variable with a different use, just for brevity
figure
plot(NN,fail(:,1),'r',NN,fail(:,2),'b--o', NN,fail(:,3),'r--', ...
     NN,fail(:,4),'b',NN,fail(:,5),'r--o',NN,fail(:,6),'b--',  NN,fail(:,7),'r')