% Script to perform Non-Parametric Hidden Markov Models Estimation 
%
% This routine allows to visualize the least squares estimators g_{k,m}
% computed and saved by "main_EstimCluster.m".
%%
clear all; close all; 

simulation_file = % Simulation data file name %

%% Model for which the estimators should be displayed.
k=3;
m=13;


%% Load data.
load(simulation_file);

%% Load data for subspace (k,m).
if strcmp(approx_type, 'trigo')
    [x, basePM, Ostar, fstar, Mth] = func_emissionProjection_Trigo(m, Qstar, param, kast, lawtype);
elseif strcmp(approx_type, 'histo')
    [x, basePM, Ostar, fstar, Mth] = func_emissionProjection_Histo(m, Qstar, param, kast, lawtype);
end

Otilde = squeeze(backup_Otilde(k, m-min_m+1, 1:m, 1:k));
Qtilde = squeeze(backup_Qtilde(k, m-min_m+1, 1:k, 1:k));

% Densities
ftilde = Otilde' * basePM;        % CMAES densities.
fm = Ostar'*basePM;               % Projected densities.

%-------------------------------------------------------------------------
permutEstim = [1 2 3 4 5];        % Permutation of the estimated densities.
permutTheo = [1 2 3 3 3];         % Permutation of the theoretical densities.
%-------------------------------------------------------------------------

kaff = max(k, kast);       % Number of densities to display.

ContrastRisk = zeros(k,1);
for j=1:k
    ContrastRisk(j)=norm(Ostar(:,permutTheo(j))-Otilde(:,permutEstim(j)));
end

strContRisk = num2str(ContrastRisk', '%d and ');
disp(['Emp. cont. L2-risks: ' strContRisk(1:end-4)]);
if (k==kast)
    Contrast_Qerror = func_dperm_Q(Qtilde, Qstar);
    disp(['Emp. cont. error on Q (in Frobenius norm) :' num2str(Contrast_Qerror)]);
end

% Graphs.
figure(4);
for j=1:kaff
    subplot(1,kaff,j)
    hold on
    if (permutTheo(j) <= kast)
        line1=plot(x,fstar(permutTheo(j),:),'b-');
        line2=stairs(x,fm(permutTheo(j),:),'k--');
    end
    if (permutEstim(j) <= k)
        line4=stairs(x,ftilde(permutEstim(j),:),'gx-','MarkerSize',3);
    end
    hold off
end

% Construct a Legend with the data from the sub-plots.
hL = legend([line1,line2,line4],{'True density', 'L2 projection', 'Least Squares'});

% Programatically move the Legend.
newPosition = [0.35 0.6 0.2 0.2];
newUnits = 'normalized';
set(hL,'Position', newPosition,'Units', newUnits);
