mdp_policy_iteration(P, R, discount, policy0, max_iter, eval_type)

function [V, policy, iter, cpu_time] = mdp_policy_iteration(P, R, discount, policy0, max_iter, eval_type)


% mdp_policy_iteration  Resolution of discounted MDP 
%                       with policy iteration algorithm
% Arguments ---------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA) = transition matrix 
%              P could be an array with 3 dimensions or 
%              a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or 
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse  
%   discount = discount rate, in ]0, 1[
%   policy0(S) = starting policy, optional 
%   max_iter = maximum number of iteration to be done, upper than 0, 
%              optional (default 1000)
%   eval_type = type of function used to evaluate policy: 
%              0 for mdp_eval_policy_matrix, else mdp_eval_policy_iterative
%              optional (default 0)
% Evaluation --------------------------------------------------------------
%   V(S)   = value function 
%   policy(S) = optimal policy
%   iter     = number of done iterations
%   cpu_time = used CPU time
%--------------------------------------------------------------------------
% In verbose mode, at each iteration, displays the number 
% of differents actions between policy n-1 and n

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification, 
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice, 
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice, 
%      this list of conditions and the following disclaimer in the documentation 
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors 
%      may be used to endorse or promote products derived from this software 
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE 
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.


cpu_time = cputime;

global mdp_VERBOSE;
if ~exist('mdp_VERBOSE'); mdp_VERBOSE=0; end;

% check of arguments
if discount <= 0 | discount >= 1
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0; 1[')
    disp('--------------------------------------------------------')
elseif ~iscell(P) & nargin > 3 & size(policy0) ~= size(P,1)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: policy0 must have the same dimension as P')
    disp('--------------------------------------------------------')
elseif iscell(P) & nargin > 3 & size(policy0) ~= size(P{1},1)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: policy0 must have the same dimension as P')
    disp('--------------------------------------------------------')
elseif nargin > 4 & max_iter <= 0
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: The maximum number of iteration must be upper than 0')
    disp('--------------------------------------------------------')
else
    
    if iscell(P)
        S = size(P{1},1);
        A = length(P);
    else
        S = size(P,1);   
        A = size(P,3);
    end
    
    PR = mdp_computePR(P,R);

    % initialization of optional arguments
    if nargin < 6; eval_type = 0; end;
    if nargin < 5; max_iter = 1000; end;
    if nargin < 4;
        % initialization of policy: 
        % the one wich maximizes the expected immediate reward
        [Vunused, policy0] = mdp_bellman_operator(P,PR,discount,zeros(S,1));
    end;
    
    if mdp_VERBOSE; disp('  Iteration  Number_of_different_actions'); end;
    
    iter = 0;
    policy = policy0;
    is_done = false;
    while ~is_done
        iter = iter + 1;
        if  (eval_type==0)   
            V = mdp_eval_policy_matrix(P,PR,discount,policy);         
        else
            V = mdp_eval_policy_iterative(P,PR,discount,policy);
        end;
        [Vnext, policy_next] = mdp_bellman_operator(P,PR,discount,V);
        
	n_different = sum(policy_next ~= policy);
        if mdp_VERBOSE; disp(['       ' num2str(iter) '                 '  num2str(n_different)]); end;
        
        if policy_next==policy | iter == max_iter
            is_done = true; 
        else
            policy = policy_next;
        end;
    end;
    
end; 

cpu_time = cputime - cpu_time;