% Here's my two cents ElisaM.
% Sorry this took so long. I'm just very busy, and I haven't been working
% on sLORETA for a while, even though I did write the sLORETA code for BST. 
% I now use eLORETA (not sLORETA) for simple distributed linear imaging instead, 
% and I'm working on COMPLETELY different methods for inverse imaging. 
% This is an actual program written in Matlab tongue, so you can copy paste
% and reproduce everything as many times as you want.

L=randn(10,100);  % Toy problem, "lead field" matrix is sampled from a Gaussian distribution, only 10 channels, 100 dipoles.
Ln=randn(10,13);  % Add a noise dictionary to get covariance matrix with non-diagonal structure (i.e., correlatted noise).
Jn=randn(13,1000); % The noise current density time series associated with the noise dictionary.
Bn=Ln*Jn;  % Simulated noise data.
Cn=cov(Bn')/100;  % Noise covariance matrix.
Cn=(Cn+Cn')/2; % Making sure it's perfectly symmetric.
W=inv(sqrtm(Cn)); % Whitening matrix.
LW=W*L; % Whitened lead field.
Sigma_J=eye(100,100)*(10/trace(LW*LW')); % A priori source covariance.


% Numerically more stable SVD approach. I think this is how I wrote for
% Brainstorm many years ago.
Rc=chol(Sigma_J,'lower'); % Cholesky used only in case we had nonzero off-diagonal elements, which we usually never do, at least in BST.
LWRc=LW*Rc;  % This is the matrix we need to SVD.
[U,S,V]=svd(LWRc,'econ');  % SVDing.
s=diag(S); % vector of singular values.
SNR=3; % Just a common safe default, could be time-dependent based on whitened data.
lambda2=SNR.^(-2); % Would be time dependent if SNR was time-dependent.
ss=s./(s.^2 + lambda2); % Same as above.
Kernel=Rc*V*diag(ss)*U';  % This is the MNE "imaging kernel".
R=Kernel*LW;  % Resolution matrix for MNE "imaging kernel"
sloretadiag=sqrt(sum(Kernel.*LW',2)); % computing the square root of diagonal of R above. Equivalently, sloretadiag=sqrt(diag(R));
KernelSLO=Kernel./repmat(sloretadiag(:),[1 10]); % sLORETA operator for the simpler case of fixed dipole orientations.


% Gram numerically noisy approach (i.e., inner product might accumulate error when there's a lot of dipoles, but it doesn't really matter in noisy reality)
T=Sigma_J*LW'*inv(LW*Sigma_J*LW' + lamda2*eye(10,10)); % Exactly equivalent to MNE Kernel via SVD,
%  which is the same as:
T_alt=(1/lambda2)*Sigma_J*LW'*inv(LW*(1/lambda2)*Sigma_J*LW' + eye(10,10));  % Exactly equivalent to MNE Kernel via SVD
corrcoef(Kernel,T)
norm(Kernel)
norm(T)
corrcoef(T,T_alt)
norm(T)
norm(T_alt)


T_oldschool=Sigma_J*L'*inv(L*Sigma_J*L' + lamda2*Cn); % Matrix is not identical to Kernel (as it should be) but it will produce the same 
%                                              exact current density vector when multiplied by data vector that you get if you multiply Kernel 
%                                              by whitened data vector.
%                                              Note that L is used not LW.
jK=Kernel*LW(:,1); % essentialy the point-spread function for first dipole
j_oldschool=T_oldschool*L(:,1);  % essentialy the point-spread function for first dipole
corrcoef(j,jK)  % perfect correlation
norm(jK) % same norms
norm(j) % same norms


RT=T*LW; % Matches R as it should.  Resolution matrix obviously has to since T and Kernel are identical.
R_oldschool=T_oldschool*L; % Also matches R as it should. Note that you have to multiply by L not LW.
corrcoef(R,R_oldschool)
norm(R)
norm(R_oldschool)
%  So all these different ways of computing the MNE resolution matrix are going
%  to give you the same exact thing. Therefore, they will give you the same
%  exact sLORETA 'inverse' operator. The important matrix here again, I
%  want to emphasize is the resolution matrxix (as far as I know).

% Now let's look at equation (15) of Pascual-Marqui 2002
S_j_PM=T_oldschool*(L*Sigma_J*L' + lambda2*Cn)*T_oldschool';
% Is it exactly equal to the resolution matrix?  Let's check.
corrcoef(R,S_j_PM)  % well you should get a perfect correlation of 1 (just a fast check as we hack along)
% But is there a scale difference?
norm(R)
norm(S_j_PM)  % Yes it seems there's a scale factor that is off.
% To get the exact same scale for the resolution matrix.
% You have to do this:  
S_j_alt=(1/Sigma_J(1,1))*T_oldschool*(L*Sigma_J*L' + lambda2*Cn)*T_oldschool';  % Note that the factor is 1/Sigma_J(1,1) not lambda2 as in your email, these are different.
% Let's check correlation
corrcoef(R,S_j_alt) 
norm(R)
norm(S_j_alt)  % Bingo they match now, which is what I think you had noticed.  You are very very veru observant.

% But what's important for sLORETA is the resolution matrix.
% Let's tale a look at the old school resolution matrix.
R_oldschool=T_oldschool*L;  % Note that you multiply by L not LW in this case.
corrcoef(R,R_oldschool)  % correlation of 1.
norm(R) % checking scales
norm(R_oldschool) % yes they match, and they are identical.

% So my two cents is that as far as I know what really matters for computing sLORETA
% is the MNE resolution matrix, which was done here in several different
% ways and gets the same exact result.  You were right that there's a
% scaling issue with equation (15) but that doesn't matter as far as I
% know. You have to use the resolution matrix (see Backus-Gilbert work for
% more info on resoution matrices and their optimization).

% AS A FINAL NOTE, I REALLY WANT TO EMPHASIZE THAT I PERSONALLY STRONGLY
% RECOMMEND ELORETA OVER ANY OTHER LINEAR DISTRIBURED METHOD INCLUDING
% SLORETA.  WITH ELORETA YOU GET ZERO LOCALIZATION BIAS AND THEY ARE
% CURRENT DENISTY SOLUTIONS, MEANING REAL FEASIBLE SOLUTIONS, NOT ACTIVITY
% MAPS WHICH IS WHAT SLORETA OR BEAMFORMERS GIVE YOU, YUCK.

% About question #2  I think you have to ask Francois but I think basically
% he probably multiplies the imaging Kernel by the whitener so that the
% data is used instead of the whitened data. If this is the case, it
% actually uses more memory and more floating point operations, but there
% might be a practical reason.  And as for the EEG case, well I have a
% feeling that the data is already rereferenced to the average 'electrode' so no need to
% keep doing that. The data is already probably stored properly referenced to the mean.
% Hope it makes sense.  Best regards.  And keep looking for
% mistakes as there could be some waiting to be corrected. 
% In any case, if there had been an error, it would have been a simple
% scaling and that wouldn't matter much for sLORETA as those activity maps
% are in arbitrary units, not units of current densisty.
% These are only my two cents. Anyone should feel more than free to look
% into this some more and possibly correct me if I'm wrong this time.
% It's very tricky business, so the more brains the merrier.  Brainstorm
% away. ?Rey Rene Ramirez
