function sums3d

% Code to make 3d plots of mutual information and conditional mutual 
% information for various detector configurations from the experiment. 
%
% Experiment set up:
%
%                      D3      D2
%                      |       |
%   Laser - > - - - - -/- - - -/- - - D1
%                     (s1)    (s2)
%
%
% This works for both coherent (Laser) and Thermal (SLED) settings.

%These spinners are actually the wrong way round. As labelled in the code. 
%No they are the right way around, I just got confused about which order I
%was putting data into the array. ****Double check this with the net*******
%Which just means that the x and y axis on my graphs are the wrong way
%around. This would make sense.


%Spinner Data
s1 = 31; %start Spinner 1 (to D3)
s2 = 11; %start Spinner 2 (to D1/2)
r1 = 49; %range Spinner 1
r2 = 49; %range Spinner 2
stp = 1; %step size
f1 = s1 + r1*stp; %finish position of Spinner 1
f2 = s2 + r2*stp; %finish position of Spinner 2

%File names
floc = 'C:\Users\MyUser\Documents\PhD\UniDimensionalCentralBroadcast\PythonControlled\';
fna = '8-11-16-set1-50x50-step1'; %This is the THERMAL data
% fna = '10-11-16-140mA-50x50-step1'; %This is the COHERENT data

%Calculation method
Imeth = 'variances'; %Use this one if no modulation
% Imeth = 'slicing'; %Use ths one if there was modulation.

%Save?
ISsave = 1; %Save everything
% ISsave = 0; %Save nothing.

fprintf('Initialising...\n')
% Start initialising
n = 0; %Spinner 1 counter

% Big 3d datasets
I12 = zeros(r1*stp+1,r2*stp+1); %I(D1;D2)
I13 = zeros(r1*stp+1,r2*stp+1); %I(D1;D3)
I23 = zeros(r1*stp+1,r2*stp+1); %I(D2;D3)
% I12g3 = zeros(r1*stp+1,r2*stp+1); %I(D1;D2|D3)
% I13g2 = zeros(r1*stp+1,r2*stp+1); %I(D1;D3|D2)
% I23g1 = zeros(r1*stp+1,r2*stp+1); %I(D2;D3|D1)

%g^{(2)}(0)
g12 = zeros(r1*stp+1,r2*stp+1);
g13 = zeros(r1*stp+1,r2*stp+1);
g23 = zeros(r1*stp+1,r2*stp+1);

% full data from each detector
D1_3d = zeros(r1*stp+1,r2*stp+1);
D2_3d = zeros(r1*stp+1,r2*stp+1);
D3_3d = zeros(r1*stp+1,r2*stp+1);

% Average intensities
x = zeros(1,r1*stp+1); %D3
y = zeros(1,r2*stp+1); %D2
z = zeros(1,r2*stp+1); %D1

fprintf('Starting loop...\n')
for SPIN1 = s1:stp:f1; % Run through spinner 1 data
    %Indexing
    n = n+1; %Spinner 1 counter
    m = 0; %Spinner 2 counter
    fprintf('Spinner 1: %d\n', SPIN1)
    
    for SPIN2 = s2:stp:f2; % Run through spinner 2 data
        fprintf('.')
        
        %index
        m = m+1; %Spinner 2 counter
        
        %read file
        fname = strcat(floc, fna, '/S', num2str(SPIN1), 'S', num2str(SPIN2), '.csv');
        Stuff = dlmread(fname);
        
        D1 = (Stuff(1,:)); %+ Stuff(4,:) + Stuff(7,:))/3;
        D2 = (Stuff(2,:)); %+ Stuff(5,:) + Stuff(8,:))/3;
        D3 = (Stuff(3,:)); %+ Stuff(6,:) + Stuff(9,:))/3;
        % Three sets of data were taken for each experimental setting
        % Here we are(/were) averaging across them. 

 
        %Store intensity data
        x(n) = (x(n) + mean(D3)); 
        D3_3d(n,m) = mean(D3);
        y(m) = (y(m) + mean(D2)); 
        D2_3d(n,m) = mean(D2);
        z(m) = (z(m) + mean(D1)); 
        D1_3d(n,m) = mean(D1);
        
        if strcmp(Imeth, 'slicing') == 1
        
            %Rough Slice
            if mean(D3)>1 || mean(D2)>1 || mean(D1)>1
                fprintf('You need to find a better method of rough slicing\n')
            end
            D3(D3>=mean(D3)) = 1;
            D3(D3<mean(D3)) = 0;
            D2(D2>=mean(D2)) = 1;
            D2(D2<mean(D2)) = 0;
            D1(D1>=mean(D1)) = 1;
            D1(D1<mean(D1)) = 0;
            
            %Find Mutual & Cond. Mutual Info.
            I12(n,m) = MutualInformation(D1, D2);
            I13(n,m) = MutualInformation(D1, D3);
            I23(n,m) = MutualInformation(D2, D3);
%             I12g3(n,m) = ConditionalMutualInformation(D1, D2, D3);
%             I13g2(n,m) = ConditionalMutualInformation(D1, D3, D2);
%             I23g1(n,m) = ConditionalMutualInformation(D2, D3, D1);
        
        elseif strcmp(Imeth, 'variances') == 1
            % Plot scattergraphs of raw data
%             figure(1)
%             clf
%             hold on
%             scatter(D1,D2, '.r')
%             xlabel('Detector 1 (V)')
%             ylabel('Detector 2 (V)')
%             title('Raw detector values')
%             
%             figure(2)
%             clf
%             hold on
%             scatter(D1,D3, '.b')
%             xlabel('Detector 1 (V)')
%             ylabel('Detector 3 (V)')
%             title('Raw detector values')
%             
%             figure(3)
%             clf
%             hold on
%             scatter(D2,D3, '.g')
%             xlabel('Detector 2 (V)')
%             ylabel('Detector 3 (V)')
%             title('Raw detector values')
            
            %Find the variances
            V2 = var(D2);
            V3 = var(D3);
            
            %Find the conditional variances:
            % 1) Find the mean of the x-axis
            M1 = mean(D1);
            M2 = mean(D2);
            
            % 2) Find the  line of points (scope precision) closest to mean
            [~, idx1] = min(abs(D1-M1));
            [~, idx2] = min(abs(D2-M2));
            MM1 = D1(idx1);
            MM2 = D2(idx2);
            
            % 3) Find the values of these points (y axis)
            P2M1 = D2(D1 == MM1);
            P3M1 = D3(D1 == MM1);
            P3M2 = D3(D2 == MM2);
%             figure(1)
%             scatter(D1(D1==MM1), P2M1, 'b.')
%             if ISsave == 1
%                 fname = strcat('sums3d_2/Scatter1_', num2str(n), '_', num2str(m));
%                 saveas(gcf,fname,'png')
%                 saveas(gcf,fname,'fig')
%             end
%             figure(2)
%             scatter(D1(D1==MM1), P3M1, 'r.')
%             if ISsave == 1
%                fname = strcat('sums3d_2/Scatter2_', num2str(n), '_', num2str(m));
%                 saveas(gcf,fname,'png')
%                 saveas(gcf,fname,'fig') 
%             end
%             figure(3)
%             scatter(D2(D2==MM2), P3M2, 'm.')
%             if ISsave == 1
%                fname = strcat('sums3d_2/Scatter3_', num2str(n), '_', num2str(m));
%                 saveas(gcf,fname,'png')
%                 saveas(gcf,fname,'fig') 
%             end
            
            
            
            % 4) Find the variance of these. (V_{B|A})
            V2g1 = var(P2M1);
            V3g1 = var(P3M1);
            V3g2 = var(P3M2);
            
            %Calculate the mutual information
            %Sometimes the noise will mean it goes slightly negative, that
            %is perfectly normal, due to the way we are calculating this,
            %but some people go crazy about negative mutual informations, 
            %so these checks happen. Take them out if you like, they never
            %go very negative. 
            i12 = log2(V2/V2g1);
%             if i12 < 0
%                 i12 = 0;
%             end
            i13 = log2(V3/V3g1);
%             if i13 < 0
%                 i13 = 0;
%             end
            i23 = log2(V3/V3g2);
%             if i23 < 0
%                 i23 = 0;
%             end
            I12(n,m) = i12;
            I13(n,m) = i13;
            I23(n,m) = i23;
            
        else
            fprintf('Youve probably had a typo, please check Imeth.\n')
        end
        
        %Calculate the g(2)(0)
        g12(n,m) = 1 + mean((D1-mean(D1)).*(D2-mean(D2)))...
            /(mean(D1)*mean(D2));
        g13(n,m) = 1 + mean((D1-mean(D1)).*(D3-mean(D3)))...
            /(mean(D1)*mean(D3));
        g23(n,m) = 1 + mean((D2-mean(D2)).*(D3-mean(D3)))...
            /(mean(D2)*mean(D3));
      
    end
    fprintf('\n')
    
    % Degrees for some graphs
    s2deg = s2:stp:f2; 
    
    % plot data as spinner 1 varies
    figure(4)
    clf
    hold on
    plot(s2deg, I12(n,:), 'r')
    plot(s2deg, I13(n,:), 'b')
    plot(s2deg, I23(n,:), 'g')
%     plot(s2deg, I12g3(n,:), 'm')
%     plot(s2deg, I13g2(n,:), 'c')
%     plot(s2deg, I23g1(n,:), 'y')
    legend('I(D1;D2)', 'I(D1;D3)', 'I(D2;D3)',...
        'I(D1;D2|D3)', 'I(D1;D3|D2)', 'I(D2;D3|D1)')
    xlabel('Spinner 2 position (degrees)')
    ylabel('Mutual Information')
    title(['m = ', num2str(n)])
    if ISsave == 1
        fname = strcat('sums3d_2/I_', num2str(n));
        saveas(gcf,fname,'png')
        saveas(gcf,fname,'fig')
    end

end

fprintf('Starting Graphs...\n')

% Normailse Intensities
x = x - min(x);
y = y - min(y);

x = x./max(x);
y = y./max(y);


% Degrees for some graphs
s1deg = s1:stp:f1;

% Graphs of Intensity for each player vs spinner degrees
figure(1)
clf
hold on
surf(s1deg, s2deg, D1_3d)
xlabel('Spinner 1 (deg)')
ylabel('Spinner 2 (deg)')
zlabel('Detector intensity (V)')
title('Detector 1')
if ISsave == 1
    fname = strcat('sums3d_2/D1_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(2)
clf
hold on
surf(s1deg, s2deg, D2_3d)
xlabel('Spinner 1 (deg)')
ylabel('Spinner 2 (deg)')
zlabel('Detector intensity (V)')
title('Detector 2')
if ISsave == 1
    fname = strcat('sums3d_2/D2_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(3)
clf
hold on
surf(s1deg, s2deg, D3_3d)
xlabel('Spinner 1 (deg)')
ylabel('Spinner 2 (deg)')
zlabel('Detector intensity (V)')
title('Detector 3')
if ISsave == 1
    fname = strcat('sums3d_2/D3_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

% Graphs of the MI and CMI vs spinner degrees
figure(4)
clf
surf(s2deg, s1deg, I12)
xlabel('Spinner 2 position (deg)')
ylabel('Spinner 1 position (deg)')
zlabel('I(D1;D2)')
if ISsave == 1
    fname = strcat('sums3d_2/I12_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(5)
clf
surf(s2deg, s1deg, I13)
xlabel('Spinner 2 position (deg)')
ylabel('Spinenr 1 position (deg)')
zlabel('I(D1;D3)')
if ISsave == 1
    fname = strcat('sums3d_2/I13_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(6)
clf
surf(s2deg, s1deg, I23)
xlabel('Spinner 2 position (deg)')
ylabel('Spinenr 1 position (deg)')
zlabel('I(D2;D3)')
if ISsave == 1
    fname = strcat('sums3d_2/I23_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

% figure(7)
% clf
% surf(s2deg, s1deg, I12g3)
% xlabel('Spinner 2 position (deg)')
% ylabel('Spinenr 1 position (deg)')
% zlabel('I(D1;D2|D3)')
% 
% figure(8)
% clf
% surf(s2deg, s1deg, I13g2)
% xlabel('Spinner 2 position (deg)')
% ylabel('Spinenr 1 position (deg)')
% zlabel('I(D1;D3|D2)')
% 
% figure(9)
% clf
% surf(s2deg, s1deg, I23g1)
% xlabel('Spinner 2 position (deg)')
% ylabel('Spinenr 1 position (deg)')
% zlabel('I(D2;D3|D1)')

% Graphs of MI and CMI vs intensity
figure(10)
clf
surf(y, x, I12)
xlabel('D3 Intensity')
ylabel('D1 Intensity')
zlabel('I(D1;D2)')
if ISsave == 1
    fname = strcat('sums3d_2/I12_int');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(11)
clf
surf(y, x, I13)
xlabel('D3 Intensity')
ylabel('D1 Intensity')
zlabel('I(D1;D3)')
if ISsave == 1
    fname = strcat('sums3d_2/I13_int');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(12)
clf
surf(y, x, I23)
xlabel('D3 Intensity')
ylabel('D1 Intensity')
zlabel('I(D2;D3)')
if ISsave == 1
    fname = strcat('sums3d_2/I23_int');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

% figure(13)
% clf
% surf(y, x, I12g3)
% xlabel('D3 Intensity')
% ylabel('D1 Intensity')
% zlabel('I(D1;D2|D3)')
% 
% figure(14)
% clf
% surf(y, x, I13g2)
% xlabel('D3 Intensity')
% ylabel('D1 Intensity')
% zlabel('I(D1;D3|D2)')
% 
% figure(15)
% clf
% surf(y, x, I23g1)
% xlabel('D3 Intensity')
% ylabel('D1 Intensity')
% zlabel('I(D2;D3|D1)')

% Graphs of Intensity for each player vs intensity (debugging only)

figure(16)
clf
hold on
surf(y, x, D1_3d)
xlabel('D3 Intensity')
ylabel('D1 Intensity')
title('D1')
if ISsave == 1
    fname = strcat('sums3d_2/D1_int');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(17)
clf
hold on
surf(y, x, D2_3d)
xlabel('D3 Intensity')
ylabel('D1 Intensity')
title('D2')
if ISsave == 1
    fname = strcat('sums3d_2/D2_int');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(18)
clf
hold on
surf(y, x, D3_3d)
xlabel('D3 Intensity')
ylabel('D1 Intensity')
title('D3')
if ISsave == 1
    fname = strcat('sums3d_2/D3_int');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

% Graphs of Intensity vs spinner degrees (debugging only)

figure(19)
clf
hold on
plot(x, s1deg, 'r')
xlabel('D1 Intensity')
ylabel('Spinner 1 (deg)')
if ISsave == 1
    fname = strcat('sums3d_2/deg1_int1');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(20)
clf
hold on
plot(y, s2deg, 'b')
xlabel('D3 Intensity')
ylabel('Spinner 2 (deg)')
if ISsave == 1
    fname = strcat('sums3d_2/deg2_int3');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(21)
clf
hold on
surf(s1deg, s2deg, g12)
xlabel('Spinner 1 (deg)')
ylabel('Spinner 2 (deg)')
zlabel('g2')
title('g12')
if ISsave == 1
    fname = strcat('sums3d_2/g12_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(22)
clf
hold on
surf(s1deg, s2deg, g13)
xlabel('Spinner 1 (deg)')
ylabel('Spinner 2 (deg)')
zlabel('g2')
title('g13')
if ISsave == 1
    fname = strcat('sums3d_2/g13_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(24)
clf
hold on
surf(s1deg, s2deg, g23)
xlabel('Spinner 1 (deg)')
ylabel('Spinner 2 (deg)')
zlabel('g2')
title('g23')
if ISsave == 1
    fname = strcat('sums3d_2/g23_deg');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(25)
clf
hold on
surf(y, x, g12)
xlabel('y')
ylabel('x')
zlabel('g2')
title('g12')
if ISsave == 1
    fname = strcat('sums3d_2/g12_i');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(26)
clf
hold on
surf(y, x, g13)
xlabel('y')
ylabel('x')
zlabel('g2')
title('g13')
if ISsave == 1
    fname = strcat('sums3d_2/g13_i');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

figure(27)
clf
hold on
surf(y, x, g23)
xlabel('y')
ylabel('x')
zlabel('g2')
title('g23')
if ISsave == 1
    fname = strcat('sums3d_2/g23_i');
    saveas(gcf,fname,'png')
    saveas(gcf,fname,'fig')
end

if ISsave == 1
fprintf('saving...\n')
fname2 = 'sums3d_2/data3d_noav.mat';
save(fname2, 'x', 'y', 'z', 'I12', 'I13', 'I23',...
    'D1_3d', 'D2_3d', 'D3_3d', 'g12', 'g13', 'g23');
end

fprintf('...Finished.\n')
return

function [Iab] = MutualInformation(A, B)
%Finds mutual information from binary.

if length(A) ~= length(B)
    fprintf('Please ensure data is of the same length.\n')
end

L = length(A);

% Prob of one thing happening
pa1 = length(find(A == 0))/L;
pa0 = 1 - pa1;
pb1 = length(find(B == 0))/L;
pb0 = 1 - pb1;

% Prob of two things happeing
pab00 = length(find(A == 0 & B == 0))/L;
pab01 = length(find(A == 0 & B == 1))/L;
pab10 = length(find(A == 1 & B == 0))/L;
pab11 = length(find(A == 1 & B == 1))/L;

% Finding Iab by summing each contribution.
Iab = 0;
if pab00~=0
    Iab =  Iab + pab00*log2(pab00/pa0/pb0);
end
if pab01~=0
    Iab = Iab + pab01*log2(pab01/pa0/pb1);
end
if pab10~=0
    Iab = Iab + pab10*log2(pab10/pa1/pb0);
end
if pab11~=0
    Iab = Iab + pab11*log2(pab11/pa1/pb1);
end

return


function [Iabc] = ConditionalMutualInformation(A, B, C)

if length(A) ~= length(B) || length(A) ~= length(C) || length(B) ~= length(C)
    fprintf('Please ensure data is of the same length.\n')
end
L = length(C);

% Prob of one thing happening
pc1 = length(find(C == 0))/L;
pc0 = 1 - pc1;


% Prob of two things happeing
pac00 = length(find(A == 0 & C == 0))/L;
pac01 = length(find(A == 0 & C == 1))/L;
pac10 = length(find(A == 1 & C == 0))/L;
pac11 = length(find(A == 1 & C == 1))/L;
pbc00 = length(find(B == 0 & C == 0))/L;
pbc01 = length(find(B == 0 & C == 1))/L;
pbc10 = length(find(B == 1 & C == 0))/L;
pbc11 = length(find(B == 1 & C == 1))/L;


% Prob of three things happening
pabc000 = length(find(A == 0 & B == 0 & C == 0))/L;
pabc001 = length(find(A == 0 & B == 0 & C == 1))/L;
pabc010 = length(find(A == 0 & B == 1 & C == 0))/L;
pabc011 = length(find(A == 0 & B == 1 & C == 1))/L;
pabc100 = length(find(A == 1 & B == 0 & C == 0))/L;
pabc101 = length(find(A == 1 & B == 0 & C == 1))/L;
pabc110 = length(find(A == 1 & B == 1 & C == 0))/L;
pabc111 = length(find(A == 1 & B == 1 & C == 1))/L;

Iabc = 0;
if pabc000~=0
   Iabc = Iabc + pabc000*log2(pc0*pabc000/pac00/pbc00); 
end
if pabc001~=0
   Iabc = Iabc + pabc001*log2(pc1*pabc001/pac01/pbc01); 
end
if pabc010~=0
   Iabc = Iabc + pabc010*log2(pc0*pabc010/pac00/pbc10); 
end
if pabc011~=0
   Iabc = Iabc + pabc011*log2(pc1*pabc011/pac01/pbc11);
end
if pabc100~=0
   Iabc = Iabc + pabc100*log2(pc0*pabc100/pac10/pbc00); 
end
if pabc101~=0
   Iabc = Iabc + pabc101*log2(pc1*pabc101/pac11/pbc01); 
end
if pabc110~=0
   Iabc = Iabc + pabc110*log2(pc0*pabc110/pac10/pbc10); 
end
if pabc111~=0
   Iabc = Iabc + pabc111*log2(pc1*pabc111/pac11/pbc11); 
end
return