clear;
N = 100;
[X,Y] = meshgrid(linspace(0,1,N));
n = 8;
%Let's set the population density f(:) to be a pyramid distribution, for no
%particular reason; add a small term so that f(:) > 0
f = min( min(X(:), 1 - X(:)) , min(Y(:), 1 - Y(:)) ) + 0.01;

for i = 1:n
    mu = rand(1,2);
    S = rand(2,2); S(1,2) = 0;
    Sigma = S*S';
    %let's make sure the eigenvalues aren't too distorted, just for visual
    %purposes
    while(max(eigs(Sigma))/min(eigs(Sigma)) > 10)
        S = rand(2,2); S(1,2) = 0; Sigma = S*S';
    end
    u = mvnpdf([X(:) , Y(:)], mu, Sigma);
    U(:,i) = u(:);
    fU(:,i) = f.*u(:);
    %fU corresponds to f(x)*u_i(x)
end

disp('Solving the max-min problem...');
cvx_begin quiet;
variable Z(N^2,n);
variable t;
dual variable lambda1;
maximize t
subject to
lambda1: t <= diag( Z'*fU );
Z >= 0;
sum(Z,2) == 1;
cvx_end
[W,I] = max(Z,[],2);
I = reshape(I,N,N);
pcolor(X,Y,I);
axis equal; shading flat;
title('Min-max solution');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q = Z'*f;
%flogU corresponds to f(x)*log(u_i(x))
flogU = zeros(N^2,n);
for i = 1:n, flogU(:,i) = f.*log(U(:,i)); end;
disp('Solving the social problem with log utilities...');
cvx_begin quiet;
variable Z(N^2,n);
dual variable lambda2;
maximize sum(diag(Z'*flogU) )
subject to
lambda2: diag(Z'*repmat(f,1,n)) == q;
Z >= 0;
sum(Z,2) == 1;
cvx_end
figure;
[W,I] = max(Z,[],2);
I = reshape(I,N,N);
pcolor(X,Y,I);
axis equal; shading flat;
title('Social problem with log utilities');