% GENFIG_3 Reproduces figure 3.
%
% USAGE     
%
%       GENFIG_3(color, display, filename)
%
% INPUTS    
%
%       color     Color (true) or grayscale (false) (optional, default = true)
%       display   Displays (true) or does not display (false) the figure (optional, default = true)
%       filename  Saves the figure in eps format under the name 'filename' (optional, default = not saved)
%
% COPYRIGHT
%
%       This file is part of the Matlab code provided for the following
%       reproducible paper:
%
%       Olivier Roy and Martin Vetterli,
%       "Dimensionality Reduction for Distributed Estimation in the Infinite Dimensional Regime",
%       vol. 54, no. 4, pp. 1655-1669, April 2008.
%
%       This program is free software; you can redistribute it and/or modify it  
%       under the terms of the GNU General Public License as published by the     
%       Free Software Foundation; either version 2 of the License, or (at your    
%       option) any later version. This software is distributed in the hope that  
%       it will be useful, but without any warranty; without even the implied     
%       warranty of merchantability or fitness for a particular purpose.          
%       See the GNU General Public License for more details                       
%       (enclosed in the file GPL).   
% 
%       GNU General Public License,
%       Copyright (C) 2008,
%       Audiovisual Communications Laboratory (LCAV),
%       Ecole Polytechnique F閐閞ale de Lausanne (EPFL),
%       CH-1015 Lausanne.
%
% COMMENTS
%
%       Author: Olivier Roy
%       Latest modifications: April 2, 2008.
function genfig_3(color, display, filename)

% We set the default values for the input arguments
save = true;
if nargin < 3
    save = false;
end
if nargin < 2
    display = true;
end
if nargin < 1
    color = true;
end

% We set the simulation parameters
global R_S R_X R_SX
randn('state',0); % for reproducible results
options     = optimset('MaxFunEvals', 100000, 'MaxIter', 100000, 'TolFun', 1e-10);
rho_val     = [0:0.01:0.99];
sigma_val   = [0.00001:0.01:3 3];

nb_try      = 1;

% We compute the gap for different values of rho and sigma
dist_min_val = zeros(length(sigma_val),length(rho_val));
dist_algo1_val = zeros(length(sigma_val),length(rho_val));
gap_val = zeros(length(sigma_val),length(rho_val));
for i = 1:length(sigma_val)
    
    sigma = sigma_val(i);
    
    for j=1:length(rho_val)
    
        rho     = rho_val(j);
        R_S     = [1 rho; rho 1];
        R_SX0   = R_S;
        R_SX1   = R_S;
        R_SX    = [R_SX0 R_SX1];
        R_N     = sigma^2*eye(2);
        R_X0    = R_S + R_N;
        R_X1    = R_S + R_N;
        R_X0X1  = R_S;
        R_X     = [R_X0 R_X0X1; R_X0X1' R_X1];

        dist_min = inf;
        for (k=1:nb_try)
            [x_min, dist_min_try] = fminsearch(@dist, [randn,randn], options);
            dist_min = min(dist_min, dist_min_try);
        end

        dist_algo1 = algo1(R_S, R_X, R_SX, [2 2], [1 1]);
        
        dist_min_val(i,j) = dist_min;
        dist_algo1_val(i,j) = dist_algo1;
        gap_val(i,j) = dist_algo1 - dist_min;
        
    end
end

% We plot the figure
figure;
surf(rho_val, sigma_val, gap_val);
if (~color)
    map = gray(256);
    map = map(80:240,:);
    colormap(map);
else
    map = colormap('jet');
    map = map(8:end,:);
    colormap(map);
end
view(-133,36);
a = axis;
axis([0, 1, 0, sigma_val(end), 0, a(6)]); 
box off
grid on
shading flat

if save
    xlabel('xLabel');
    ylabel('yLabel');
    zlabel('zLabel');
    title('title');
    saveas(gcf, filename,'psc2');
end
if display
    xlabel('\rho');
    ylabel('\sigma');
    zlabel('\Delta MSE');
else
    close gcf
end


function dist_step = dist(x)
global R_S R_X R_SX

theta  = x(1);
vartheta = x(2);
Q_1 = [cos(theta) sin(theta)];
Q_2 = [cos(vartheta) sin(vartheta)];
C = [Q_1 zeros(1,2);zeros(1,2) Q_2];

dist_step = trace(R_S - R_SX*C'*inv(C*R_X*C')*C*R_SX');