% This function calculates the adaptive threshold according to Cai and Liu
% Input:
% X - random matrix of size  NxT, where N is (most likely) smaller T
% SampleCov - sample covariance of the random matrix
% etaS - the shrinkage rule eta 
% delta - the hyperparameter to select the degree of shrinkage


function CovThresholded=adaptiveThresh(X,SampleCov,etaS,delta)

[N,T]=size(X);
muX = mean(X,2);
eta = etaS;
CovThresholded = zeros(N,N);
lambda = zeros(N,N);
summand = zeros(N,N,T);
diagVar = trace(SampleCov)/N;
diagVar = diag(diag(SampleCov));

for t=1:T
    demeanedX = X(:,t)-muX;
    summand(:,:,t)=(demeanedX*demeanedX'-SampleCov).^2; 
end

ThetaHat = mean(summand,3);

lambda = delta*sqrt(ThetaHat*log(N)/T);

CovThresholded = SampleCov.*max( (1-abs(lambda./SampleCov).^eta) , 0) + diagVar*diag(min( (abs(diag(lambda)./diag(SampleCov)).^eta) ,1));

end

%% Backup
%% different shrinkage if the diagonal goes below zero
%         if min(diag(CovThresholded))<=0
%             for i=1:N
%                 if  CovThresholded(i,i)<=0
%                     shrinkageRow=mean(CovThresholded(i,:)./SampleCov(i,:));
%                     shrinkageColumn=mean(CovThresholded(:,j)./SampleCov(:,j));
%                     shrinkage=sqrt(shrinkageRow*shrinkageColumn);
%                     CovThresholded(i,i)=shrinkage*SampleCov(i,i);
%                 end
%             end
%         end

%%  How to make the matrix PSD
% [~,flag] = chol(CovThresholded);
% if flag ~= 0
%     CovThresholded = 0.5*(CovThresholded'+CovThresholded);
%     [~,Sigma,V] = svd(CovThresholded);
%     H = V*Sigma*V';
%     CovThresholded = (CovThresholded + H)/2;
%     CovThresholded = 0.5*(CovThresholded'+CovThresholded);
%     p = 1;
%     k = 0;
%     while p~=0
%         k = k + 1;
%         if p~=0
%             mineig = min(eig(CovThresholded));
%             CovThresholded = CovThresholded + (-(k.^2)*mineig + 1e-12)*diagVar; %
%         end
%         [~,p] = chol(CovThresholded);
%     end
% end

