function  blahutArimoto
% accuracy
givenEpsilon = 0.000001;
% beta = 10;

% init points
numXs =100;
Xs = rand(numXs,2) * 0.8 + 0.1 ;
% Xs(1,1) = 0;
% Xs(1,2) = 0.2;
% Xs(2,1) = 0;
% Xs(2,2) = 0.7;


% init representatives
numTs = 2;
%Ts = rand(numTs,2)* 0.8 + 0.1;
Ts(1,1) = 0;
Ts(1,2) = 0.5;
Ts(2,1) = 1;
Ts(2,2) = 1;


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

% probabilities for t (initial guess)
probTs = ones(numTs,1);
probTs = probTs./sum(probTs);

% show initial problem
f=figure; hold on;
colormapSize = 1024;
colormap(gray(colormapSize));
hold on;
title(['start']);
c = colormap;
for x = 1:numXs
    plot(Xs(x,1), Xs(x,2), 'o','MarkerEdgeColor','g','MarkerFaceColor','w','MarkerSize',probXs(x)* 15*numXs+ 1);
end

myColormap = [0 0 0; 1 1 1];
for t = 1:numTs
    plot(Ts(t,1), Ts(t,2), 's','MarkerEdgeColor','r','MarkerFaceColor',myColormap(t,:),'MarkerSize',probTs(t)*15* numTs + 1);
end
saveas(f,['blahutArimoto_start.tif']);
close(f);   

for beta = [1:100]        
	% probabilities for t (initial guess)
	probTs = ones(numTs,1);
	probTs = probTs./sum(probTs);
	
	% conditional probabilities 
	probTsGivenXs = zeros(numTs,numXs);
	prevProbTsGivenXs = zeros(numTs,numXs);
	iterateFlag = 1;
	iterationsCount = 0;
	frameCounter = 1;
	
	% iterate
	while (iterateFlag)
        % calc new conditional 
        prevProbTsGivenXs = probTsGivenXs;
        
        for x = 1:numXs
            for t = 1:numTs
                probTsGivenXs(t,x) = probTs(t) * exp(-beta * calcDist(Xs(x,:),Ts(t,:)));
            end
            probTsGivenXs(:,x) = probTsGivenXs(:,x)./sum(probTsGivenXs(:,x));
        end
        
	
        % calc new prob on t
        for t = 1:numTs
            probTs(t) =  probTsGivenXs(t,:) * probXs; 
        end
	
	
        % converegence?
        if (max(probTsGivenXs - prevProbTsGivenXs)  < givenEpsilon)
            iterateFlag = 0;
        end
	
        if (beta == 1 | beta == 10 | beta == 100 | iterateFlag == 0)
            % show current approximation    
            f=figure; hold on;
            colormapSize = 1024;
            colormap(gray(colormapSize));
            hold on;
            title(['iteration: ' num2str(iterationsCount) ', beta: ' num2str(beta)]);
            c = colormap;
            for x = 1:numXs
                plot(Xs(x,1), Xs(x,2), 'o','MarkerEdgeColor','g','MarkerFaceColor',c(floor((1-probTsGivenXs(1,x)) * (colormapSize-1) + 1),:),'MarkerSize',probXs(x)* 15*numXs+ 1);
            end
            
            myColormap = [0 0 0; 1 1 1];
            for t = 1:numTs
                plot(Ts(t,1), Ts(t,2), 's','MarkerEdgeColor','r','MarkerFaceColor',myColormap(t,:),'MarkerSize',probTs(t)*15* numTs + 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 = calcDist(a,b)
    dist = sqrt(sum((a-b).^2));
return