function [belief,value]=fictitiousPlay(A,K)
% [belief,value]=fictitiousPlay(A,K)
%
% Performs fictitious play for a N-player mixed "matrix" game.
%
% Inputs:
% A - (N+1)-dimensional tensor that defines the game
%     outcomes. Specifically, 
%       A(i1,i2,...,iN,i) = outcome for player i, when 
%                 player 1 plays the pure policy i1
%                 player 2 plays the pure policy i2
%                   ...
%                 player N plays the pure policy iN
% K - length of the mixed path
%
% Outputs:
%
% belief - N-dimensional cell array with the player's beliefs at
%          the end of the path. Specifically, 
%                 belief{i} is the belief for player i, 
%          which is a probablity distribution over the simplex of with
%          dimension equal to size(A,i);
% value  - N-vector with the average outcomes of the game for the
%          beliefs at the end of the path. Specifically,
%                 value(i) is the average outcome for player i
    t0=clock();
    
    N=length(size(A))-1;
    %% compute random initialization for the mixed path
    belief=cell(N,1);
    for i=1:N
        if 1
            belief{i}=rand(size(A,i),1); % random belief;
        else
            belief{i}=ones(size(A,i),1); % uniform belief;
        end
        belief{i}=belief{i}/sum(belief{i}); % normalize to get distribution
    end
    %% iterate fictitious play
    for k=1:K
        if mod(k,1000)==0
            fprintf('%5d:\n',k);
        end
        %% compute probability-weighted outcomes 
        Ay=A;
        for i=1:N
            reps=size(A);
            reps(i)=1;
            shape=ones(1,N+1);
            shape(i)=length(belief{i});
            yi=reshape(belief{i},shape);
            yi=repmat(yi,reps);
            % for player i's outcomes do not multiply by its own
            % belief to eventually compute the best response
            str=['yi(',repmat(':,',1,N),num2str(i),')=1;'];
            eval(str);
            Ay=Ay.*yi;
        end % for i=1:N
        %% compute best response
        y=cell(N,1);
        for i=1:N 
            % average outcomes over actions of other players
            str=['Ay(',repmat(':,',1,N),num2str(i),')'];
            S=eval(str);
            for j=1:N
                if j~=i
                    S=sum(S,j);
                end
            end
            % compute best responses
            [~,j]=min(S,[],i); % pure best response
            y{i}=zeros(size(S,i),1);
            y{i}(j)=1;                % mixed best response
        end % for i=1:N
        %% update belief
        for i=1:N
            belief{i}=(k*belief{i}+y{i})/(k+1);
        end % for i=1:N
    end  % for k=1:K
    %% compute final value
    value=nan(N,1);
    % compute probability-weighted outcomes 
    Ay=A;
    for i=1:N
        reps=size(A);
        reps(i)=1;
        shape=ones(1,N+1);
        shape(i)=length(belief{i});
        yi=reshape(belief{i},shape);
        yi=repmat(yi,reps);
        Ay=Ay.*yi;
    end % for i=1:N
    for i=1:N 
        % average outcomes over actions of all players
        str=['Ay(',repmat(':,',1,N),num2str(i),')'];
        S=eval(str);
        for j=1:N
            S=sum(S,j);
        end
        value(i)=S;
    end % for i=1:N
    fprintf('   done (%.3f sec)\n',etime(clock(),t0));
end

%% 2-player
% outcome for player i
%    Ji=sum_j1 sum_j2 belief{1}(j1) belief{2}(j2) A(j1,j2,i)

% options for player 1
%    Ji(j1) = sum_j2 belief{2}(j2) A(j1,j2,1)
%           = sum_j2 Ay(j1,j2,1)
% where
%    Ay(j1,j2,1) = belief{2}(j2) A(j1,j2,1)

%% 3-player
% outcome for player i
%    Ji=sum_j1 sum_j2 sum_j3 belief{1}(j1) belief{2}(j2) belief{3}(j3) A(j1,j2,j3,i)

% options for player 1
%    Ji(j1) = sum_j2 sum_j3 belief{2}(j2) belief{3}(j3) A(j1,j2,j3,1)
%           = sum_j2 sum_j3 Ay(j1,j2,j3,1)
% where
%    Ay(j1,j2,j3,1) = belief{2}(j2) belief{3}(j3) A(j1,j2,j3,1)