clear all

load ../solns/1063_acoustic_brain-5mm/SolutionStats.mat
s5 = Soln;
load ../solns/1063_acoustic_brain-3mm/SolutionStats.mat
s3 = Soln;
load ../solns/1063_acoustic_brain-5mm/AngleAssingments.mat
P5 = Prob;
load ../solns/1063_acoustic_brain-3mm/AngleAssingments.mat
P3 = Prob;

clear Soln Prob;

PDF = struct( ...
   't1', ['BalancedAvg'], ...
   't2', ['BalancedMax'], ...
   't3', ['PrimalAvg'], ...
   't4', ['PrimalMax'], ...
   't5', ['DualAvg'], ...
   't6', ['DualMax'], ...
   't7', ['InteriorAvg'], ...
   't8', ['InteriorMax'], ...
   't9', ['SetCover1'], ...
   't10', ['SetCover2'], ...
   't11', ['Scoring'], ...
   't12', ['EntropyAvgBalanced']);

for i = 1:12
   % Vectors of Solution Times for the Angle Assignments
   PDFtime5(i) = ...
     eval(['s5.AngleAssignment.',eval([sprintf('PDF.t%d',i)]),'.Time']);
   PDFtime3(i) = ...
     eval(['s3.AngleAssignment.',eval([sprintf('PDF.t%d',i)]),'.Time']);

   % Each method gets its own matrix, cols are in order of Angle
   % value stated above, with each col containing the selected angles.
   VQAngles5(1:9,i) = ...
      eval(['s5.BeamSelector.VQ.',eval([sprintf('PDF.t%d',i)]),'.B9.Angles']);
   VQValue5(i) = ...
      eval(['s5.BeamSelector.VQ.',eval([sprintf('PDF.t%d',i)]),'.B9.JudgementValue']);
   VQAngles3(1:9,i) = ...
      eval(['s3.BeamSelector.VQ.',eval([sprintf('PDF.t%d',i)]),'.B9.Angles']);
   VQValue3(1:9,i) = ...
      eval(['s3.BeamSelector.VQ.',eval([sprintf('PDF.t%d',i)]),'.B9.JudgementValue']);
   SC1Angles5(1:9,i) = ...
      eval(['s5.BeamSelector.SC1.',eval([sprintf('PDF.t%d',i)]),'.B9.Angles']);
   SC1Value5(1:9,i) = ...
      eval(['s5.BeamSelector.SC1.',eval([sprintf('PDF.t%d',i)]),'.B9.JudgementValue']);
   SC1Angles3(1:9,i) = ...
      eval(['s3.BeamSelector.SC1.',eval([sprintf('PDF.t%d',i)]),'.B9.Angles']);
   SC1Value3(1:9,i) = ...
      eval(['s3.BeamSelector.SC1.',eval([sprintf('PDF.t%d',i)]),'.B9.JudgementValue']);
   SAngles5(1:9,i) = ...
      eval(['s5.BeamSelector.S.',eval([sprintf('PDF.t%d',i)]),'.B9.Angles']);
   SValue5(1:9,i) = ...
      eval(['s5.BeamSelector.S.',eval([sprintf('PDF.t%d',i)]),'.B9.JudgementValue']);
   SAngles3(1:9,i) = ...
      eval(['s3.BeamSelector.S.',eval([sprintf('PDF.t%d',i)]),'.B9.Angles']);
   SValue3(1:9,i) = ...
      eval(['s3.BeamSelector.S.',eval([sprintf('PDF.t%d',i)]),'.B9.JudgementValue']);
end
%5*SAngles5(:,11)

%
% Lets Build a multilayer Angles Matrix
%     layer 1 - VQ
%     layer 2 - SC1
%     layer 3 - S
%
Angles5(:,:,1) = VQAngles5;
Angles5(:,:,2) = SC1Angles5;
Angles5(:,:,3) = SAngles5;
Angles3(:,:,1) = VQAngles3;
Angles3(:,:,2) = SC1Angles3;
Angles3(:,:,3) = SAngles3;

% Calculate the differences between the angles 
for k = 1:3
   A5(:,:,k) = sort(Angles5(:,:,k));
   A3(:,:,k) = sort(Angles3(:,:,k));
   for i = 1:8
      D5(i,:,k) = A5(i+1,:,k) - A5(i,:,k);
      D3(i,:,k) = A3(i+1,:,k) - A3(i,:,k);
   end
   D5(9,:,k) = mod(A5(1,:,k) - A5(9,:,k), 72);
   D3(9,:,k) = mod(A3(1,:,k) - A3(9,:,k), 72);
end
D5(:,7,1)'
D3(:,7,1)'

% Calculate nearest neighbor differences between 3 and 5mm
Q1(:,:,:) = abs(D5(:,:,:) - D3(:,:,:));
for k = 1:3
   for j = 1:12
      if A5(1,j,k) <= A3(1,j,k)
         Q2(:,j,k) = abs(D5(:,j,k) - [D3(9,j,k) ; D3(1:8,j,k)]);
      else
         Q2(:,j,k) = abs(D5(:,j,k) - [D3(2:9,j,k) ; D3(1,j,k)]);
      end
   end
end
Q1(:,7,1)'
Q2(:,7,1)'
sum(Q1(:,7,1))
sum(Q2(:,7,1))

% Q contains the smallest differences between 3 and 5mm
for k = 1:3
   for j = 1:12
      if sum(Q1(:,j,k)) <= sum(Q2(:,j,k))
         Q(:,j,k) = Q1(:,j,k);
      else
         Q(:,j,k) = Q2(:,j,k);
      end
   end
end

% calculate the Average shift and the vaiance of the shifts
% NOTE: this is in degrees.
for k = 1:3
   for j = 1:12
      Avg(j,k) = mean(5*Q(:,j,k));
      Vrc(j,k) = var(5*Q(:,j,k));
   end
end

% Calculate the Probabilities
for j = 1:12
   AngleValues5(:,j) = eval(['P5.', eval(sprintf('PDF.t%d',j)) ] );
   AngleValues3(:,j) = eval(['P3.', eval(sprintf('PDF.t%d',j)) ] );
end
for k=1:3
   for j=1:12
      PRB3(k,j) = sum(AngleValues3(Angles3(:,j,k),j)) / sum(AngleValues3(:,j));
      PRB5(k,j) = sum(AngleValues5(Angles5(:,j,k),j)) / sum(AngleValues5(:,j));
   end
end

% Calculate the difference in function values
for j = 1:12
   PDFtime5(1,j) = ...
      eval(['s5.BeamSelector.VQ.',eval([sprintf('PDF.t%d',j)]),'.B9.JudgementValue']);
   PDFtime5(2,j) = ...
      eval(['s5.BeamSelector.SC1.',eval([sprintf('PDF.t%d',j)]),'.B9.JudgementValue']);
   PDFtime5(3,j) = ...
      eval(['s5.BeamSelector.S.',eval([sprintf('PDF.t%d',j)]),'.B9.JudgementValue']);
   PDFtime3(1,j) = ...
      eval(['s3.BeamSelector.VQ.',eval([sprintf('PDF.t%d',j)]),'.B9.JudgementValue']);
   PDFtime3(2,j) = ...
      eval(['s3.BeamSelector.SC1.',eval([sprintf('PDF.t%d',j)]),'.B9.JudgementValue']);
   PDFtime3(3,j) = ...
      eval(['s3.BeamSelector.S.',eval([sprintf('PDF.t%d',j)]),'.B9.JudgementValue']);
end
