function [relValues, lambda, policy] = FindPolicyPI_AVE(contribArray, transArray, ...
    feasibleActions, initialPolicy, optDir, verboseFlag)
numStates = length(initialPolicy);
oldPolicy = initialPolicy(:);
newPolicy = zeros(numStates,1);
stopit = false;
count = 1;
while (~stopit)
    % build vector and matrix
    [rhs,matrix] = makeArrays(oldPolicy);
    % evaluate current policy
    auxSol = matrix\rhs;
    oldRelVal = [auxSol(1:(numStates-1));0];
    lambda = auxSol(numStates);
    % improve policy
    if strcmpi(optDir, 'min')
        for j = 1:numStates
            newPolicy(j) = findMin(j);
        end
    else
        for j = 1:numStates
            newPolicy(j) = findMax(j);
        end
    end
    if verboseFlag
        printPolicy(oldPolicy);
    end
    if all(oldPolicy == newPolicy)
        stopit = true;
    else
        count = count + 1;
        oldPolicy = newPolicy;
    end
end
relValues = oldRelVal;
policy = oldPolicy;

% nested functions 
    % build system of linear equations
    function [outVet, outMatrix] = makeArrays(policy)
        % we assume that the relative value of the last state is zero
        % the last variable is the average contribution 
        outVet = zeros(numStates,1);
        auxMatrix = zeros(numStates,numStates);
        for k=1:numStates
            outVet(k) = contribArray(k,policy(k));
            auxMatrix(k,:) = transArray(k,:,policy(k));
        end  % for
        auxMatrix = eye(numStates) - auxMatrix;
        outMatrix = [auxMatrix(:,1:(numStates-1)), ones(numStates,1)];
    end
    % print policy
    function printPolicy(vet)
       fprintf(1,'Policy at step %d:  ',count);
       aux = [sprintf('%d,', vet(1:(end-1))), sprintf('%d', vet(end))];
       fprintf(1,'%s\n',aux);
    end
    % find min cost action
    function best = findMin(k)
        best = NaN;
        valBest = inf;
        act = feasibleActions{k};
        for a = 1:length(act)
            aux = contribArray(k,act(a)) + dot(transArray(k,:,act(a)), oldRelVal);
            if aux < valBest 
                valBest = aux;
                best = act(a);
            end
        end
    end
    % find max reward action
    function best = findMax(k)
        best = NaN;
        valBest = -inf;
        act = feasibleActions{k};
        for a = 1:length(act)
            aux = contribArray(k,act(a)) + dot(transArray(k,:,act(a)), oldRelVal);
            if aux > valBest 
                valBest = aux;
                best = act(a);
            end
        end
    end

end
