%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .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(T, plottype, 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(:,3)=1-T(:,1)-T(:,2);
T(T<1e-6)=0;

N=101;

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

%H_x= zeros(N,N);
I_x_y = zeros(N,N);

for n1=1:N
   for n2=1:(N-n1)
       x3 = 1-x1(n1)-x2(n2);
       if x1(n1)+x2(n2)<=1
           x=[x1(n1) x2(n2) x3];
           
           
           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(T==0)=0;
           I_x_y(n1,n2)=sum(x*(T.*I));
           %H_x(n1,n2)=(log(1./x)/log(2))*x';
           
       end
   end 
end
subplot(2,1,1);
axis off;
text(0.05,0.5, 'T= ');
text(0.1,0.5, [num2str(T)]);

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

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

text(0.55,0.75, ['C= max(I(X;Y)) = ' num2str(peak) ' bit/symbol'
'P[x_1] = ' num2str(x1_max);...
'P[x_2] = ' num2str(x2_max);...
'P[x_3] = 1-P[x_1]-P[x_2] = ' num2str(x3_max)]);



subplot(2,1,2);

if (plottype=='heatmap')
      imagesc(x1, x2, I_x_y');
      xlabel({'P[x_1]'});
      ylabel({'P[x_2]'});
      set(gca,'ydir','normal');
      colorbar("EastOutside");
      title(colorbar,"I(X;Y) in bit/symbol");
      output_format = portrait;
elseif(plottype=='surface')
      surf(x1,x2,I_x_y');
      xlabel({'P[x_1]'});
      ylabel({'P[x_2]'});
      zlabel({'I(x,y)'});
      title({'I(X;Y) in bit/symbol'});
      view(-37.5,50);
      shading interp;
      colorbar("EastOutside");
      title(colorbar,"I(X;Y) in bit/symbol");
      output_format = portrait;
elseif(plottype=='meshz')
      meshz(x1,x2,I_x_y);
      xlabel({'P[x_1]'});
      ylabel({'P[x_2]'});
      zlabel({'I(x,y)'});
      title({'I(X;Y) in bit/symbol'});
      output_format = portrait;
end

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

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

end

