% Script to perform spectral HMM estimation on the save files from EstimCluster.
%
% This routine takes a simulation file from the least squares estimation procedure
% "main_EstimCluster.m" and runs the spectral algorithm with the specified parameters.
% It then displays the spectral AND least squares estimators.
%%
clear all; close all;

%% Simulation parameters
simulation_file = % data file name %

forced_k = 0;               % 0 to let the algorithm choose k. Otherwise, takes k = forced_k.
m = 13;                     % Size of the projection space.
SVthreshold = 1.5;          % Threshold to state which singular values are significant.
n_rotations = 1;			% Number of generations of the rotation matrix (Step 4).

%% Load data.
load(simulation_file);
% Loaded variables: 'min_m', 'max_m', 'kmax', 'contrastvalue', 'initmethod', 'n', 'vartilde', 'backup_Mhat', 'backup_Phat', ...
%    'backup_Nhat', 'backup_lhat', 'backup_Otilde', 'backup_Qtilde', 'approx_type', 'X', 'Y', 'Qstar', 'param', ...
%    'kast', 'annotation', 'lawtype', 'lawtypeannotation', 'instanceNumber'

%% True emission laws and their projection onto Sm
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


%% Estimation
[k, pihat, Qhat, Mhat, Ohat, pihatZERO, QhatZERO, OhatZERO, Nhat] = func_spectralEstimation(m, Y, n_rotations, SVthreshold, approx_type, forced_k);

disp(['selected k: ' num2str(k)])


%% Plots
% Compute densities.
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));

% Estimation of the densities
fhat = Ohat' * basePM;            % Spectral densities
ftilde = Otilde' * basePM;        % CMA-ES densities
fm = Ostar'*basePM;               % Projected densities

%-------------------------------------------------------------------------
permutSquares = [2 3 1 4];        % Permutation of least squares densities in plots.
permutSpect = [3 2 1 4];          % Permutation of spectral densities in plots.
permutTheo = [1 2 3 4];           % Permutation of theoretical densities in plots.
%-------------------------------------------------------------------------

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

SpectralRisk = zeros(k,1);
ContrastRisk = zeros(k,1);
for j=1:k
    SpectralRisk(j)=norm(Ostar(:,permutTheo(j))-Ohat(:,permutSpect(j)));
    ContrastRisk(j)=norm(Ostar(:,permutTheo(j))-Otilde(:,permutSquares(j)));
end

strSpecRisk = num2str(SpectralRisk', '%d and ');
disp(['Spectral L2-risks: ' strSpecRisk(1:end-4)]);
strContRisk = num2str(ContrastRisk', '%d and ');
disp(['Emp. cont. L2-risks: ' strContRisk(1:end-4)]);

% Display densities.
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 (permutSquares(j) <= k)
        line4=stairs(x,ftilde(permutSquares(j),:),'gx-','MarkerSize',3);
    end
    if (permutSpect(j) <= k)
        line3=stairs(x,fhat(permutSpect(j),:),'r-.','LineWidth',2);
    end
    hold off
    title(['Emission law ' num2str(j)])
end

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

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