function [k, pihat, Qhat, Mhat, Ohat, pihatZERO, QhatZERO, OhatZERO, Nhat] = func_spectralEstimation(m,Y,n_rotations, SVthreshold, approx_type, forced_k)
% Spectral algorithm (with optional selection of k) for the estimation of HMM parameters.

temp = size(Y);
n = temp(2);

%% Part I: SPECTRAL ALGORITHM
disp(['Spectral estimation... (m=' num2str(m) ')']);
% Step 1: Empirical estimators
if strcmp(approx_type, 'trigo')
    [Mhat, Phat, Nhat, lhat] = func_empiricalEstimators_Trigo(m, Y);
elseif strcmp(approx_type, 'histo')
    [Mhat, Phat, Nhat, lhat] = func_empiricalEstimators_Histo(m, Y);
end

%% Estimation of k
[~,Dhat,U1hat]=svd(Nhat);

if(forced_k)
    % Force choice of k
    k = forced_k;
else
    D = diag(Dhat);

    % Regression
    D_range = ceil(m/2):m; 				% Regression on the smallest half of the singular values (instead of the Mreg smallest).
    p = polyfit(D_range',D(D_range),1);

    noise_singular_values = p(2)+p(1)*(1:m);    % Estimated pure noise singular values.

    ratio = D ./ noise_singular_values';
    i=1;
    while ratio(i) > SVthreshold
       i = i+1;
    end
    k = i-1;
end

%% Step 2 : SVD
Uhat=U1hat(:,1:k);              % Keep the k largest sing. val.

%% Step 3 and 4 - multiple Theta
ecartMin = 0;
for auxVar = 1:n_rotations
    Theta = func_uniformUnitary(k);
    eta=Uhat*Theta;

    Bhat=zeros(k,k,m);
    A = zeros(k,k,k,m);
    for b=1:m
        Bhat(:,:,b)=(Uhat'*Phat*Uhat)\(Uhat'*squeeze(Mhat(:,b,:))*Uhat);
        for j=1:k
            A(j,:,:,b)=eta(b,j).*Bhat(:,:,b);
        end
    end
    Chat = sum(A,4);

    % Step 5 and 6
    [Rhat,~]=eig(squeeze(Chat(1,:,:)));

    Lhat=zeros(k,k);
    for b=1:k
        for j=1:k
            Lhat(j,:)=diag(Rhat\(squeeze(Chat(j,:,:)))*Rhat);
        end
    end

    minimum = abs(Lhat(1,1) - Lhat(1,2));
    for k1 = 1:k
        for k2=1:k
            for k3=1:k
                if (k2 ~= k3)
                    minimum = min( [ minimum (abs(Lhat(k1,k2) - Lhat(k1,k3))) ] );
                end
            end
        end
    end

    if (auxVar == 1)
        OhatZERO = Uhat*Theta*Lhat;
        pihatZERO=(Uhat'*OhatZERO)\Uhat'*lhat;
        pihatZERO=func_simplexProjection(pihatZERO);         % proj onto Delta_k
        QhatZERO=(Uhat'*OhatZERO*diag(pihatZERO))\Uhat'*Nhat*Uhat/(OhatZERO'*Uhat);
        pihatZERO = func_stationaryLaw(QhatZERO);
    end
    if (minimum < ecartMin || auxVar == 1)
        ecartMin = minimum;
        LhatMin = Lhat;
        ThetaMin = Theta;
        disp('Changement de Theta')
    end
end
Ohat=Uhat*ThetaMin*LhatMin;

%% Step 7 - in plots

%% Step 8 and 9
pihat=(Uhat'*Ohat)\Uhat'*lhat;
pihat=func_simplexProjection(pihat);
Qhat=(Uhat'*Ohat*diag(pihat))\Uhat'*Nhat*Uhat/(Ohat'*Uhat);

for j=1:k
    Qhat(j,:) = func_simplexProjection(Qhat(j,:)')';
end
pihat = func_stationaryLaw(Qhat);

disp(['...done']);
