% Script to perform spectral HMM estimation on the save files from EstimCluster.
%
% This routine runs the spectral algorithm with the specified parameters and displays the resulting estimators.
%%

%% Spectral parameters
forced_k = 3;               % 0 to let the algorithm choose k. Otherwise, takes k = forced_k.
m = 17;                     % 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).


%% Random number generator initialization
% Initialize with current time to avoid to get the same results for different instances.
rng('shuffle');

%% Load the simulation parameters from param_simulation.
[kmax, min_m, max_m, approx_type, annotation, stopeval] = param_simulation();
% Load the model parameters from param_HMM.
[Qstar, pistar, param, kast, n, lawtype, lawtypeannotation] = param_HMM();

%% Data initialization and chain generation    
% Chain generation
[X,Y] = func_generateHMM(Qstar, pistar, param, kast, n, lawtype);


%% Estimation
% True emission laws and their projection onto the approximation space of dimension m.
if strcmp(approx_type, 'trigo')
    [x, basePM, Ostar, fstar, ~] = func_emissionProjection_Trigo(m, Qstar, param, kast, lawtype);
elseif strcmp(approx_type, 'histo')
    [x, basePM, Ostar, fstar, ~] = func_emissionProjection_Histo(m, Qstar, param, kast, lawtype);
end

[k, ~, ~, ~, Ohat, ~, ~, ~, Nhat] = func_spectralEstimation(m, Y, n_rotations, SVthreshold, approx_type, forced_k);

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


%% Spectrum
Nstar = Ostar * diag(pistar) * Qstar * Ostar';
% Compute and save spectrum.
[~,Dtemp,~]=svd(Nhat);
Dhat = diag(Dtemp);
[~,Dtemp,~]=svd(Nstar);
Dstar = diag(Dtemp);

reg_range = ceil(m/2):m; 		% Regression on the smallest half of the singular values (instead of the Mreg smallest).
p = polyfit(reg_range',Dhat(reg_range),1);
noise_singular_values = p(2)+p(1)*(1:m);


%% ----------------------- PLOT SPECTRUM ---------------------
figure(1);
%-------------------------------------------------------------------------
start_index = 1;          % First singular value to be displayed.
end_index = 8;            % Last singular value to be displayed.
%-------------------------------------------------------------------------
range = start_index:end_index;
plot(range, Dhat(range));
hold on;
plot(range, Dstar(range), 'r-.o');
% plot(range, noise_singular_values(range), 'g-');
hold off;
title('Spectrum of N');
ylim([0 (1.2*Dstar(range(1)))]);

%% ---------------------- PLOT DENSITIES ---------------------
figure(2);
% Estimation of the densities
fhat = Ohat' * basePM;            % Spectral densities
fm = Ostar'*basePM;               % Projected densities

%-------------------------------------------------------------------------
permutSpect = [2 1 3 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);
for j=1:k
    SpectralRisk(j)=norm(Ostar(:,permutTheo(j))-Ohat(:,permutSpect(j)));
end

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

% Display densities.
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 (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],{'True density', 'L2 projection', 'Spectral'});

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