% payoffs for the Shapley game
R = [0 0 1; 1 0 0; 0 1 0];
S = R';

% use celp to build constraints describing the region of correlated 
% equilibria
ntot = numel(R);
A = [celp(R); celp(S); eye(ntot)];
b = zeros(size(A,1),1);
E = ones(1,ntot);
d = -1;

% parameters for solving linear programs
par.verbose = 0;
par.epsilon = 1e-5;

% find some random corners of the set of correlated equilibria and
% display them.
for i = 1:5
    c = randn * R(:) + randn * S(:);
    [val, pr, du] = iqph(1e-4, c, A, b, E, d, par);
    P = reshape(pr, size(R,1), [])
    fprintf('P1 value %g\n', P(:)'*R(:));
    fprintf('P2 value %g\n\n', P(:)'*S(:));
end