%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Insert title here" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: Marc Gauger / Alexander Knaub   
% date: 11.04.2014   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [h_x,r] = transinformation_fun_3d(T, svg_en, webdemo_en, filename)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, set up graphical output

% Create an invisible figure.
if (webdemo_en)
  fig = figure(1); set(fig, "visible", "off");
else
  fig = figure(1); set(fig, "visible", "on");
endif

landscape = "-S930,350";
portrait  = "-S640,480";

output_format = landscape; %portrait;     % default


set(0, "defaultlinelinewidth", 2);

if (webdemo_en)
  set(0, "defaultaxesfontsize", 8);
  set(0, "defaultaxesfontname", "Helvetica");
  set(0, "defaulttextfontsize", 8);
  set(0, "defaulttextfontname", "Helvetica");
  const_lw = 2;
elseif (svg_en)
  set(0, "defaultaxesfontsize", 6);
  set(0, "defaultaxesfontname", "Helvetica");
  set(0, "defaulttextfontsize", 6);
  set(0, "defaulttextfontname", "Helvetica");
  const_lw = 1;
else
  set(0, "defaultaxesfontsize", 10);
  set(0, "defaultaxesfontname", "Helvetica");
  set(0, "defaulttextfontsize", 8);
  set(0, "defaulttextfontname", "Helvetica");
  const_lw = 2;
endif

% END, set up graphical output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%sum of row have to be equal one.
T(:,4)=1-T(:,1)-T(:,2)-T(:,3);
T(T<1e-6)=0;

N=100;

x1 = linspace (0, 1 , N);
x2 = linspace (0, 1 , N);
x3 = linspace (0, 1 , N);

H_x_max= 0;
I_x_y_max = 0;

%for n1=1:N
%   for n2=1:N-n1
%       for n3=1:N-n1-n2
%           
%           x4 = 1-x1(n1)-x2(n2)-x3(n3);
%           if x1(n1)+x2(n2)+x3(n3)<=1
%               x=[x1(n1) x2(n2) x3(n3) x4];
%
%
%               y= x*T;
%
%               I = zeros(size(T,1),size(T,2),N);
%               Y_inv = zeros(size(T,1),size(T,2),N);
%
%               Y_inv = diag(1./y);
%               I=log(T*Y_inv)/log(2);
%               I_x_y=sum(x*(T.*I));
%               H_x=(log(1./x)/log(2))*x';
%               
%               if I_x_y > I_x_y_max
%                   I_x_y_max = I_x_y;
%                   I_x_y_x1 = x1(n1);
%                   I_x_y_x2 = x2(n2);
%                   I_x_y_x3 = x3(n3);
%                   I_x_y_x4 = x4;
%                   
%               end
%               if H_x > H_x_max
%                   H_x_max = H_x;
%                   H_x_x1 = x1(n1);
%                   H_x_x2 = x2(n2);
%                   H_x_x3 = x3(n3); 
%                   H_x_x4 = x4;                  
%               end
%               
%           end   
%       end
%   end 
%end

[C px] = BlahutArimoto(T');

%I_x_y_max
%I_x_y_x1
%I_x_y_x2
%I_x_y_x3
%I_x_y_x4


%H_x_max
%H_x_x1
%H_x_x2
%H_x_x3
%H_x_x4

axis off;
text(0.05,0.7, 'T= ');
text(0.1,0.7, [num2str(T)]);

%peak=max(max(I_x_y));
%[C,rowmaxarray]=max(I_x_y);
%[maxvalue,colmax]=max(C);
%rowmax=rowmaxarray(colmax);

x1_max=px(1);
x2_max=px(2);

%x1_max = x1(colmax);
%x2_max = x2(rowmax);
x3_max = 1-x1_max -x2_max;


%text(0.55,0.75, ['C= max(I(x,y)) = ' num2str(I_x_y_max) ' bit/symbol']);
text(0.55,0.75, ['C= max(I(X;Y)) = ' num2str(C) ' bit/symbol']);


text(0.55,0.5, 
['P[x_1] = ' num2str(px(1));...
'P[x_2] = ' num2str(px(2));...
'P[x_3] = ' num2str(px(3));...
'P[x_4] = ' num2str(px(4))]);


print(filename, "-dpng","-S640,480");

if (svg_en)
  print(["fig_",filename], "-dsvg",output_format);
endif

end

function [C,px] = BlahutArimoto(P)
% Blahut-Arimoto algorithm to compute channel capacities for DMCs
% 
% [C,px] = BlahutArimoto(P)
%
% Input is channel transition matrix P where:
%   Columns of P correspond to different input values x_i
%   Rows of P correspond to different output values y_j
% I.e., entry at row j and column i is P(y_j | x_i)
%
% Output is Capacity (in bits) and optimum input distribution
%
% Laurent Schmalen 2014	 


	 eps = 1e-9;
	 
	 % number of inputs
	 V = size(P,2);

	 % number of outputs
	 W = size(P,1);
	 
	 % initial guess for input distribution
	 rx = ones(V,1)/V;

	 abbruch = false;
	 more off;
	 while ~abbruch
	   temp = P*rx;
	   qxy = (repmat(rx(:)',W,1).*P)./repmat(temp(:),1,V);

	   temp = zeros(V,1);
	   for k=1:V
	       temp(k) = prod(qxy(:,k).^P(:,k));
	   end
	   
	   new_rx = temp./sum(temp);
	   if sum(abs(new_rx-rx)) < eps
	      abbruch = true;
	   end
	   rx = new_rx;
	 end

	 px = rx(:);

	 Pt = P.*log2(P);
	 Pt(isnan(Pt))=0;
	 
	 C = -(P*px).'*log2(P*px) + sum(Pt*px);
end

