function probTsGivenXs = informationBottleneck(beta)
% accuracy
givenEpsilon = 10^(-8);

% read the image witht he conditional probabilites
probYsGivenXsPrivate = imread('face1.tif');

% create probabilities from picture
maxIntensity = double(max(max(probYsGivenXsPrivate)));
probYsGivenXsPrivateProb = double(probYsGivenXsPrivate)./maxIntensity;
% figure; imagesc(probYsGivenXsPrivate);

picSize = size(probYsGivenXsPrivate,1);

% inits relevant variable
numYs = 2;


% init points
numXs =300;
Xs = floor(rand(numXs,2) * (picSize-1)) + 1 ;
% Xs(1,1) = 0;
% Xs(1,2) = 0.2;
% Xs(2,1) = 0;
% Xs(2,2) = 0.7;

% fill conditional probabilities y given x
probYsGivenXs = zeros(numYs,numXs);
for x = 1:numXs
    probYsGivenXs(1,x) = probYsGivenXsPrivateProb(Xs(x,2),Xs(x,1));
    probYsGivenXs(2,x) = 1- probYsGivenXs(1,x);
end

% init representatives
numTs = 2;

% probabilities for x
 probXs = rand(numXs,1);
 probXs = probXs./sum(probXs);

% joint probaiblities x, y
probYsXs = probYsGivenXs.* repmat(probXs',numYs,1); 

% for beta = [0.1 2 5 5.5 6.0 6.2 6.7 7 9 13 16 20] % use 2.8 for face

    % conditional probabilities 
	probTsGivenXs = rand(numTs,numXs);
    probTsGivenXs = probTsGivenXs./repmat(sum(probTsGivenXs),numTs,1);
	prevProbTsGivenXs = zeros(numTs,numXs);

    
	iterateFlag = 1;
	iterationsCount = 0;
	frameCounter = 1;
	
	% iterate
	while (iterateFlag)
        % calc new prob on t
        probTs =  probTsGivenXs * probXs; 
        
        zeroProbIdx = find(probTs == 0);
        % calc new conditional prob y given t
        probYsGivenTs =  ((probTsGivenXs * probYsXs')./ repmat(probTs,1,numYs))'; 
        probYsGivenTs(:,zeroProbIdx) = 0;
        
        % calc new conditional t given x
        prevProbTsGivenXs = probTsGivenXs;
        
        for x = 1:numXs
            for t = 1:numTs
                probTsGivenXs(t,x) = probTs(t) * exp(-beta * calcKullbackDist(probYsGivenXs(:,x),probYsGivenTs(:,t)));
            end
            probTsGivenXs(:,x) = probTsGivenXs(:,x)./sum(probTsGivenXs(:,x));
        end
	
%         % calc I(xt)
%         Ixt = 0;
%         for x = 1:numXs
%             for t = 1:numTs
%                 probXT = probXs(x) * probTsGivenXs(t,x);
%                 if ((probXs(x) * probTs(t) == 0) & (probXT ~= 0))
%                     Ixt = 1000;
%                 elseif (probXT ~= 0)
%                     Ixt = Ixt +  probXT * log(probXT/ (probXs(x) * probTs(t)));
%                 end
%             end
%         end
%         
%         % calc I(ty)
%         Ity = 0;
%         probYs = probYsGivenTs * probTs;
%         for y = 1:numYs
%             for t = 1:numTs
%                 probTY = probTs(t) * probYsGivenTs(y,t);
%                 if ((probYs(y) * probTs(t) == 0) & (probTY ~= 0))
%                     Ity = 1000;
%                 elseif (probTY ~= 0)
%                     Ity = Ity +  probTY * log(probTY/ (probYs(y) * probTs(t)));
%                 end
%             end
%         end
%         F = Ixt - beta * Ity
        
        
        
        
        % converegence?
        iterateFlag = 0;
        maxJs = 0;
        for x = 1:numXs
            pTag = (probTsGivenXs(:,x) + prevProbTsGivenXs(:,x))./2;
            Js = 0.5 * (calcKullbackDist(probTsGivenXs(:,x),pTag) + calcKullbackDist(prevProbTsGivenXs(:,x),pTag));
            if (Js > givenEpsilon)
                iterateFlag = 1;
            end
            if (Js > maxJs)
                maxJs = Js;
            end
        end % check Js for x
        maxJs;
        if (iterationsCount > 201)
            iterateFlag = 0;
            disp(['not convergent']);
        end
	
%         if (beta == 1 | beta == 10 | beta == 100 | iterateFlag == 0)
         if ( iterateFlag == 0)
            % show current approximation    
            f=figure; 
            colormapSize = 256;
            %colormap(autumn(colormapSize));
%             colormap(gray(colormapSize));
            imagesc(probYsGivenXsPrivateProb); 
            hold on;
            title(['iteration: ' num2str(iterationsCount) ', beta: ' num2str(beta)]);
            c = colormap(gray(colormapSize));
            colormap(bone(colormapSize));
            colorbar;
            for x = 1:numXs
                plot(Xs(x,1), Xs(x,2), 'o','MarkerEdgeColor','g','MarkerFaceColor',c(floor(probTsGivenXs(1,x) * (colormapSize-1)) + 1,:),'MarkerSize',probXs(x)* 15*numXs+ 1);
            end
            
%             frames(frameCounter) = getframe;
%             saveas(f,['blahutArimoto_'  num2str(beta) '_frame_' num2str(frameCounter) '.tif']);
            frameCounter = frameCounter + 1;
%             close(f); 
        end
        iterationsCount = iterationsCount + 1;
	end% while: iteration
% end% for: beta
      
return

function dist = calcKullbackDist(p,q)
if (min(q) == 0)
    dist =1000;
else
    dist = 0;
    for i = 1:length(p)
        if (p(i) ~= 0)
            dist = dist + p(i).*log(p(i)./q(i));    
        end
    end
end
return