function res=rbreginvover(Amat, rhs, lambda, LTL, blocks, varargin)
%
% res=rbreginvover(Amat, rhs, lambda, LTL, blocks, varargin)
%
% Solve an overdetermined Gauss-Newton normal equation
%
% author: Qianqian Fang (q.fang <at> neu.edu)
% 
% input:
%     Amat: the left-hand-side matrices (a containers.Map object) at specified wavelengths 
%     rhs: the right-hand-side vectors for all sources (independent of wavelengths)
%     lambda: the Tikhonov regularization parameter
%     LTL: the regularization matrix in the form of the inverse of the
%         unknown-covariance (inv(C_mu)), which is L'*L
%     blocks: the dimensions of the 2D submatrices of Amat, needed if
%         size(LTL,1)==size(LTL,2) is not the same as size(Amat,2)
%
% output:
%     res: the least-square solution of the matrix equation
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details 
%
% -- this function is part of Redbird-m toolbox
%

% solve an overdetermined Gauss-Newton normal equation
%  delta_mu=inv(J'J + lambda*(L'L))*J'*(y-phi)

% if any node has no sensitivity, remove them from inversion
length0=size(Amat,2);
idx0=find(sum(Amat)~=0);
if(length(idx0)<length0)
    Amat=Amat(:,idx0);
    %TODO: need to shrink Lqr as well
end

emptydata=find(sum(Amat')~=0);
if(length(emptydata)<size(Amat,1))
    Amat=Amat(emptydata,:);
    rhs=rhs(emptydata);
end

if(nargin>=4 && ~isempty(LTL))
    nx=size(LTL,1);
    oldnx = length0/length(fieldnames(blocks));
    if (nx > length(idx0)/length(fieldnames(blocks)))
        Lidx = idx0(idx0<=nx);
        LTL = LTL(Lidx,Lidx);
        oldnx = nx;
        nx = size(LTL,1);
    end
end

rhs=Amat'*rhs(:);

Hess=Amat'*Amat; % Gauss-Hessian matrix, approximation to Hessian (2nd order)

% [Hess,Gdiag]=rbnormalizediag(Hess);

if(nargin<4 || isempty(LTL))
    Hess(1:1+size(Hess,1):end)=Hess(1:1+size(Hess,1):end)+lambda;
else
    if(size(Hess,1)==size(LTL,1))
        Hess=Hess+lambda*LTL;
    else
        nx=size(LTL,1);
        for i=1:nx:size(Hess,1)
            Hess(i:i+nx-1,i:i+nx-1)=Hess(i:i+nx-1,i:i+nx-1)+lambda*LTL;
        end
%         if(isempty(blocks)) % assume the Hess matrix size is multiples of LTL
%             for i=1:nx:size(Hess,1)
%                 Hess(i:i+nx-1,i:i+nx-1)=Hess(i:i+nx-1,i:i+nx-1)+lambda*LTL;
%             end
%         else
%             len=cumsum([1; structfun(@(x) x(2), blocks(1))]);
%             for i=1:length(blocks)
%                 if(nx==len(i+1)-len(i)+1)  % if the block size match LTL
%                     Hess(len(i):len(i+1)-1,len(i):len(i+1)-1)=Hess(len(i):len(i+1)-1,len(i):len(i+1)-1)+lambda*LTL;
%                 else % if size does not match, add lambda*I
%                     idx=sub2ind(size(Hess),len(i):len(i+1)-1,len(i):len(i+1)-1);
%                     Hess(idx)=Hess(idx)+lambda;
%                 end
%             end
%         end
     end
end
[Hess,Gdiag]=rbnormalizediag(Hess);
res=Gdiag(:).*rbfemsolve(Hess, Gdiag(:).*rhs, varargin{:});
% res = rbfemsolve(Hess,rhs,varargin{:});

if(length(idx0)<length0)
    res0=zeros(length0,size(res,2));
    res0(idx0,:)=res;
    res=res0;
end