
% Fist initializing the source covariance matrix. I assume that BST
% still does soemthing like this.

% ===== ADJUSTING SOURCE COVARIANCE MATRIX =====
% Adjusting Source Covariance matrix to make trace of L*C_J*L' equal to number of sensors.
display('Adjusting source covariance matrix.')
C_J = speye(sspl, sspl);
C_J = spdiags(w, 0, C_J);
trclcl = trace(L * C_J * L');
C_J = C_J * (rnkC_noise / trclcl);

% Below just some notes on how to do a time-dependent SNR approach, in case
% it's of interest.
% % % ===== WHITENED ERF/ERP DATA =====
% % % Compute whitened data.
% % Bw=W*B;
% % % Compute time-dependent SNR and lambda
% % SNR=sum(Bw.^2)./size(Bw,1); % whitened power SNR.
% % figure; plot(SNR)
% % SNR_mean=mean(SNR)
% % lambda2=SNR.^(-2);
% % num_times=length(lambda2);


% I'm going to assume that the regularization parameter is computed 
% the usual way based on SNR and that whitening is used the "standard" way.
% I see a lot of changes in function since I last edited it.
% So hopefully this is still compatible (in terms of SNR, regularization
% parameter, and whitening). Not sure who changed code or why.
if strcmpi(OPTIONS.InverseMethod, 'eloreta')
    %lambda2=mean(lambda2); % If using time-dependent SNR, use mean for
    %learning source covariance for computational economy...
    display('Computing eLORETA optimal source covariance matrix.')
    % Compute the a priori covariance matrix!
    numnonz=0;
    for k = 1:numL
        numnonz=numnonz+(spl(k)*numdipcomp(k));
    end
    sigmaj=ones(numnonz,1);
    sigmaj_i=sigmaj;
    sigmaj_j=sigmaj;
    changesigmajiter=inf;
    iter=0;
    costf=1;
    maxiter=100;
    cost=zeros(maxiter,1);
    while iter<maxiter && changesigmajiter>norm(sigmaj)*1e-10
        iter=iter+1;
        %%% Grammian approach
        G=(L*C_J*L' + eye(size(L,1)));
        [U,S,V] = svd(G,'econ');
        s = diag(S);
        ss = s ./ (s.^2 + lambda2);
        KernelR=L'*U*diag(ss)*U';
        %%%
        sigmaj_old=sigmaj;
        start=0;
        sigma_start=0;
        for k = 1:numL
            % Indices of dipole components.
            start = start + 1;
            endd = start + spl(k) - 1;
            if (numdipcomp(k) == 1)
                % Indices of all source covariance components .. !
                sigma_start=sigma_start+1;
                sigma_endd=sigma_start+spl(k)-1;
                R=sum(KernelR(start:endd,:).*L(:,start:endd)', 2);
                %%%
                sigmaj(sigma_start:sigma_endd)=R.^(-1/2);
                Rsqrtinv=1./sqrt(sigmaj(sigma_start:sigma_endd).*R);
                %%%
                if iter==1
                    sigmaj_i(sigma_start:sigma_endd)=start:endd;
                    sigmaj_j(sigma_start:sigma_endd)=start:endd;
                end
                sigma_start=sigma_endd;
                cost(iter)=cost(iter)+sum(abs(1-Rsqrtinv*sigmaj(sigma_start:sigma_endd)));
            elseif (numdipcomp(k)==3 || numdipcomp(k)==2)
                for spoint = start:numdipcomp(k):endd
                    % Indices of all source covariance components .. !
                    sigma_start=sigma_start+1;
                    sigma_endd=sigma_start+(numdipcomp(k)^2)-1;
                    indspt=spoint:spoint+numdipcomp(k)-1;
                    R=KernelR(indspt,:)*L(:,indspt);
                    if OPTIONS.gram
                        SIR=sqrtm(R);
                        SIR=pinv(SIR);
                        Rsqrtinv=pinv(sqrtm(full(C_J(indspt,indspt))*R));
                    end
                    sigmaj(sigma_start:sigma_endd)=SIR(:);
                    if iter==1
                        sigmaj_i_tmp=repmat(indspt(:),[1 numdipcomp(k)]);
                        sigmaj_i(sigma_start:sigma_endd)=sigmaj_i_tmp(:);
                        sigmaj_j_tmp=repmat(indspt(:)',[numdipcomp(k) 1]);
                        sigmaj_j(sigma_start:sigma_endd)=sigmaj_j_tmp(:);
                    end
                    sigma_start=sigma_endd;
                    cost(iter)=cost(iter)+sum(sum(abs(eye(numdipcomp(k))-Rsqrtinv)));
                end
            end
            start=endd;
        end  
        if iter==1  
            sigmaj_old=sigmaj;
            for sti=1:length(sigmaj_i)
                sigmaj_old(sti)=C_J(sigmaj_i(sti),sigmaj_j(sti));
            end
        end
        % display('Adjusting source covariance matrix.')
        C_J = sparse(sigmaj_i,sigmaj_j,sigmaj,sspl,sspl);
        trclcl = trace(L * C_J * L');
        C_J = C_J * (rnkC_noise / trclcl);
        sigmaj=sigmaj*(rnkC_noise / trclcl);
            
        changesigmajiter=norm(sigmaj-sigmaj_old);
        changesigmaj(iter)=changesigmajiter;
        figure(1),plot(changesigmaj(1:iter)); axis tight
        costf=cost(iter);
        %figure(1),plot(cost(1:iter)); axis tight
    end
    
    % Compute whitened eLORETA operator.
    if OPTIONS.gram
        G=(L*C_J*L' + eye(size(L,1)));
        [U,S,V] = svd(G,'econ');
        s = diag(S);
    end
    cost(iter+1:end)=[];
end
ss = s ./ (s.^2 + lambda2);
%ss = s ./ (s.^2 + lambda2(t)); % For time-dependent SNR, you could do this
%instead ...
Kernel_eLORETA=C_J*L'*U*diag(ss)*U';